CN104573269A - 一种基于强耦合整体技术的索膜结构抗风设计方法 - Google Patents

Abstract

本发明涉及一种基于强耦合整体技术的索膜结构抗风设计方法，建立索膜结构的初始形态；建立流固耦合体系控制方程和流固交界面处的耦合条件，对其进行时空离散，得到流固耦合体系的强耦合整体方程；采用SSTK-ω湍流模型模拟湍流，得到考虑湍流模型的流固耦合体系的强耦合整体方程，并采用Newton方法将其进行线性化；构造线性化后的强耦合整体方程矩阵的预处理器，得到带有预处理器的流固耦合体系的强耦合整体方程，采用Krylov子空间投影法对其进行求解，得到强耦合整体方程流体压力、流体速度、索膜结构位移和线弹性有限元模型的位移；从而进行索膜结构的抗风设计，得到抗风设计参数。

Description

一种基于强耦合整体技术的索膜结构抗风设计方法

技术领域

本发明属于建筑行业索膜结构抗风分析与设计领域，具体涉及一种基于强耦合整体技术的索膜结构抗风设计方法。

背景技术

目前国内外索膜结构抗风设计理论尚不成熟。风与索膜结构的流固耦合作用是索膜结构抗风设计中的难点问题，流固耦合作用研究的滞后已成为制约索膜结构风振理论完善和发展的主要因素。沈世钊院士也曾指出，对于索膜结构的抗风设计，重点是解决风与索膜结构的流固耦合问题。因此采用合理方法研究风与索膜结构的流固耦合问题对于完善索膜结构的风振设计理论、指导工程实践具有重要的理论意义和实用价值。

目前工程界已经可以通过使用相对成熟的有限元软件对索膜结构在风荷载下的应力、变形效应进行精确计算。同时，也能够通过使用计算流体动力学软件对索膜结构所受的风荷载进行准确估计。然而，对于流固耦合计算的数值方法和应用却少之又少。流固耦合作用分析的数值模拟技术包括弱耦合分区技术、强耦合分区技术和强耦合整体技术。弱耦合分区技术是在每一时间步内先对流体控制方程和结构控制方程分别独立求解，然后将作用在流体域模型上的气动力荷载通过流体和结构的交界面传递给结构域模型，用来预测结构的位移，然后再将结构的位移作为新的荷载传递给流体域，此过程如此反复，直到结果收敛于指定值。强耦合分区技术就是在流体域和结构域各自求解器的基础上再增加一次迭代循环，在每一个时间步求解非线性方程组，计算出全场变量的值。强耦合整体技术是指在每一时间步内对流体控制方程和结构控制方程同时联立求解，即在每一时间步内流体域和结构域互相传递气动力和结构位移的过程是同时进行的；狭义的说，就是对整个流固耦合问题用一个方程组表示，并进行求解。

目前对于索膜结构抗风设计中的流固耦合作用的分析主要采用弱耦合分区和强耦合分区技术，原因是它们可以利用现有的流体或结构计算软件，但这两种方法由于流体网格和结构网格的划分不一致，结构的反应要滞后于流体域，从而导致交界面处的计算能量耗散，引起数值计算不稳定，最终影响结果的准确性。并且不论流体和结构单个求解器的准确性有多高，这种方法也只能达到一阶准确性。而且对于索膜结构这类流体和结构的耦合作用很强的情况，弱耦合分区和强耦合分区技术非常容易发散，计算准确性会受到较大影响。

尽管采用强耦合整体技术计算流固耦合问题的研究还非常有限，但强耦合整体技术在稳定性和准确性方面却显示出其优势。目前国内外对于索膜结构抗风设计中流固耦合问题研究的主要发展趋势是力争采用更准确高效的计算方法和手段解决复杂的流固耦合计算，为索膜结构的抗风设计提供更可靠的理论或方法，因此强耦合整体技术可以成为研究风与索膜结构流固耦合作用的有利工具。

采用强耦合整体技术对索膜结构进行抗风设计，首先需要重点解决的问题之一就是处理空气流体域和结构域交界面处的数据信息传递，即流体域计算域和结构域计算域的计算结果的相互交换，完成下一时间步计算的网格更新过程。需要注意的是，对交界面处信息的处理要适应强耦合整体技术的实施。并在处理流固交界面处的数据信息传递之后，采用强耦合整体技术需要采用一个非线性方程组表示空气流体和索膜结构形成的流固耦合体系，即需要推导出耦合体系的非线性强耦合整体方程，

强耦合整体技术所采用的非线性方程组通常需要采用Newton方法线性化后求解。在求解该方程过程中，会耗费在对雅可比矩阵的反复集成和相应线性体系解的Newton修正中，这往往会导致该问题的计算量是不可接受的。因此在使用强耦合整体技术时对方程的高效求解也是采用强耦合整体技术需要解决的重点问题。构造预处理器是求解线性化后的非线性方程的一种有效方法，因此在强耦合整体方程的基础上构造合适的预处理器，使得方程的求解过程简化，并节约大量计算机时也是非常关键的步骤。目前的流固耦合计算的预处理器都是分别单独的针对流体域和结构域的，在使用过程中分别对流体域和结构域使用这些预处理器，这样在很大程度上就会降低计算效率。因此需要针对所要采用的强耦合整体技术的特点，构造相应的预处理器，以实现整体隐式方程的求解的高效性和稳定性。

索膜结构属于低矮的大跨度结构，位于大气边界层中湍流强度较高的近地区域，索膜结构周围的流场会使其发生变形的缘故，索膜结构常常会遇到包括分离、再附、漩涡脱落大尺度流瞬时流动过程。因此需要采用湍流模型以便更准确的进行索膜结构的抗风设计。这里选用合适的湍流模型并将其融入提出的整体隐式方程中，也是非常重要的环节，也是索膜结构抗风计算的主要技术难点之一。

发明内容

针对现有技术的不足，目前索膜结构抗风设计中的流固耦合分析的分区数值模拟技术存在易发散、准确性差的缺点，为索膜结构抗风设计带来安全隐患，本发明提出一种基于强耦合整体技术的索膜结构抗风设计方法。将线弹性有限元模型应用于空气流体与索膜结构形成的耦合系统中，采用线弹性有限元模型处理流体和结构交界面处的数据信息交换传递，并基于线弹性有限元模型方程提出一种适应于强耦合整体技术的网格更新方法；采用流体方程、结构方程和线弹性有限元模型方程的变分形式推导出强耦合整体方程；采用Newton方法线性化强耦合整体方程，构造强耦合整体方程的预处理器，采用两方程的SST k-ω湍流模型模拟湍流，并将其融入整体隐式方程中，给出了计算膜结构风振中流固耦合计算的实现过程和步骤，为索膜结构抗风设计提供更准确可靠的理论方法和造作简便的计算平台。

一种基于强耦合整体技术的索膜结构抗风设计方法，包括以下步骤：

步骤1：根据实际膜结构材料的物理性质和实际膜结构的初始预应力，通过找形方法建立索膜结构的初始形态；

所述的通过找形方法建立索膜结构有限元模型具体方法为：采用小杨氏模量逐步刚度逼近的方法建立索膜结构的初始形态；

步骤2：建立流固耦合体系控制方程和流固交界面处的耦合条件，流固耦合体系控制方程包括：流体域控制方程、结构域控制方程和线弹性有限元模型控制方程；流固交界面处的耦合条件包括流固交界面线弹性有限元模型和结构的位移耦合条件、流固交界面处流体和结构的应力耦合条件、流固交界面处流体和结构的速度耦合条件；

所述的建立流固耦合体系控制方程的具体方法为：采用纳维叶-斯托克斯方程建立流体域控制方程，采用拉格朗日方程建立结构域控制方程，采用动力平衡方程的半离散有限元方程的准静态形式建立线弹性有限元模型控制方程；

所述的采用线弹性有限元模型处理流体域变形的具体方法为：将流体网格视为线弹性固体，利用流体单元刚化矩阵中的刚化函数得到线弹性有限元模型单元的单元刚化矩阵，得到流固耦合体系的网格更新公式，对流固耦合体系网格进行实时更新，采用控制网格单元的纵横比控制更新后网格质量；

步骤3：采用线弹性有限元模型处理流体域变形，利用线弹性有限元模型控制方程对流固耦合体系网格进行实时更新；

步骤4：对流固耦合体系的控制方程和流固交界面处的耦合条件进行时空离散，得到流固耦合体系的各个控制方程的变分弱形式和流固交界面处的各耦合条件的变分弱形式，合并各个控制方程和耦合条件的变分弱形式，得到流固耦合体系的强耦合整体方程；

所述的对流固耦合体系的控制方程和流固交界面处的耦合条件进行时空离散的具体方法为：采用伽辽金有限元法对流固耦合体系控制方程、流固交界面处线弹性有限元模型和结构的位移耦合条件、流固交界面处流体和结构的应力耦合条件进行空间离散，采用隐式有限差分方法对流固交界面处流体和结构的速度耦合条件进行时间离散；

步骤5：采用SST K-ω湍流模型模拟湍流，建立SST K-ω湍流模型控制方程，对SSTK-ω湍流模型控制方程进行空间离散，得到SST K-ω湍流模型控制方程的变分弱形式，用SST K-ω湍流模型的变分弱形式替换流固耦合体系的强耦合整体方程中的流体域控制方程的变分弱形式，得到考虑湍流模型的流固耦合体系的强耦合整体方程；

所述的对SST K-ω湍流模型控制方程进行空间离散的具体方法为：采用伽辽金有限元法对SST K-ω湍流模型的控制方程进行空间离散，即将滑移边界条件作为湍流的边界，将结构壁面的法线方向速度作为垂直于壁面的流体速度，得到离散后的SST K-ω湍流模型的控制方程；

步骤6：采用Newton方法将考虑湍流模型的流固耦合体系的强耦合整体方程进行线性化，得到线性化的流固耦合体系的强耦合整体方程；

步骤7：构造线性化后的强耦合整体方程矩阵的预处理器，得到带有预处理器的流固耦合体系的强耦合整体方程；

所述的构造线性化后的强耦合整体方程矩阵的预处理器的具体方法为：构造线性化后的强耦合整体方程矩阵的分块表达，在雅可比矩阵中使用块状矩阵结构，使得预处理矩阵模块化，得到线性化后的强耦合整体方程矩阵的预处理器；

步骤8：采用Krylov子空间投影法求解带有预处理器的流固耦合体系的强耦合整体方程，得到强耦合整体方程流体压力、流体速度、索膜结构位移和线弹性有限元模型的位移；

步骤9：利用流体压力、流体速度和索膜结构位移进行索膜结构的抗风设计，得到抗风设计参数：索膜结构风压系数、索膜结构风振响应和索膜结构气动特性。

本发明的有益效果是：

本发明引入线弹性有限元模型处理流体域变形，根据线弹性有限元模型控制方程提出了流固耦合体系网格的更新方法，解决了流体域和结构域的数据交换传递问题，为实现强耦合整体技术奠定了基础；采用流体、结构以及线弹性有限元模型控制方程三者的变分弱形式实现了流固耦合体系的整体耦合，推导出了强耦合整体方程，实现了流体域和结构域的整体强耦合；在强耦合整体方程中引入了湍流模型，采用两方程的SST k-ω湍流模型，并采用伽辽金有限元法进行离散，给出了SST k-ω湍流模型的变分弱形式，进一步完整了强耦合整体技术对索膜结构进行抗风设计的计算框架；推导了强耦合整体方程的预处理器，在雅可比矩阵中使用特殊的块结构，使得预处理矩阵模块化，大大提高了强耦合整体方程的求解效率和计算稳定性；通过计算得到索膜结构的风压系数、风振响应、膜结构气动特性重要抗风设计参数，为索膜结构抗风设计提供了简便的计算工具。解决了索膜结构抗风设计的传统分区法计算中的收敛性和稳定性问题，提高了计算准确性，并将湍流的影响考虑在内。使得计算准确性相比以往分区方法有了显著的提高，而且计算速度较快，使这项技术从科学研究走向实际应用成为可能。

附图说明

图1为本发明具体实施方式中基于强耦合整体技术的索膜结构抗风设计方法的流程图；

图2为本发明具体实施方式中采用小杨氏模量逐步刚度逼近的方法建立索膜结构的初始形态示意图；

图3为本发明具体实施方式中流固耦合体系网格进行实时更新方法的流程图；

图4为本发明具体实施方式中鞍型膜结构表面网格划分图；

图5为本发明具体实施方式中鞍型膜结构的几何尺寸图和鞍型膜结构分区图；

其中，(a)为本发明具体实施方式中鞍型膜结构的几何尺寸图；(b)为本发明具体实施方式中鞍型膜结构分区图；

图6为本发明具体实施方式中不同风向角下索膜结构中点的平均风压系数时程图；

其中，(a)为0°风向角下中点的平均风压系数时程图；(b)为90°风向角下中点的平均风压系数时程图；

图7为本发明具体实施方式中无量纲附加空气质量的变化规律。

具体实施方式

下面结合附图对本发明具体实施方式加以详细的说明。

一种基于强耦合整体技术的索膜结构抗风设计方法，如图1所示，包括以下步骤：

步骤1：根据实际膜结构材料的物理性质和实际膜结构的初始预应力，采用小杨氏模量逐步刚度逼近的方法建立索膜结构的初始形态。

本实施方式中，确定索膜结构的初始形态属于给定预应力分布的形状确定问题，即已知实际索膜结构的初始预应力，通过找形得到膜的初始几何状态。利用杨氏模量逐步刚度逼近方法的找形可以通过ANSYS软件完成：输入索膜结构的材料性能模量，此时将索膜结构的杨氏模量(即弹性模量)输入的非常小(一股比实际索膜结构的杨氏模量小3个数量级)，然后从初始几何形状开始迭代，得到膜结构的找形结果时，再将索膜结构的杨氏模量设定为真实值，在此基础上进行自平衡迭代求解，获得索膜结构的初始形状。本发明采用小杨氏模量逐步刚度逼近的方法建立索膜结构的初始形态如图2所示。

步骤2：建立流固耦合体系控制方程和流固交界面处的耦合条件，流固耦合体系控制方程包括：流体域控制方程、结构域控制方程和线弹性有限元模型控制方程；流固交界面处的耦合条件包括流固交界面线弹性有限元模型和结构的位移耦合条件、流固交界面处流体和结构的应力耦合条件、流固交界面处流体和结构的速度耦合条件。

建立流固耦合体系控制方程的具体方法为：采用纳维叶-斯托克斯方程建立流体域控制方程，采用拉格朗日方程建立结构域控制方程，采用动力平衡方程的半离散有限元方程的准静态形式建立线弹性有限元模型控制方程。

采用纳维叶-斯托克斯方程建立的流体控制方程如式(1)、(2)、(3)所示：

$\▿\·v_f=0---\left(1\right)$

${\mathrm{\ρ}}_f(\frac{{\∂v}_f}{\∂t}+v_f\·{\▿v}_f)=\▿\·{\mathrm{\σ}}_f+f_f^B---\left(2\right)$

${\mathrm{\σ}}_f=\mathrm{\μ}[{\▿v}_f+{\left({\▿v}_f\right)}^T]-\mathrm{pI}---\left(3\right)$

其中，v_f为流体流动速度，为空间梯度，ρ_f为流体密度，σ_f为流体的全应力张量，即流体的压力和流体的粘性力，μ为流体粘度，p是流体压力，为流体体积力，I为单位矩阵，t为时间，为对流体流动速度空间梯度的转置。

采用拉格朗日方程建立的结构域的控制方程如式(4)、(5)、(6)、(7)、(8)所示：

$\▿\·{\mathrm{\σ}}_l+f_s=0---\left(4\right)$

σ_l＝Fσ_p (5)

$E=\frac{1}{2}(F^{T}F-I)---\left(6\right)$

$F=\▿x---\left(7\right)$

σ_p＝λ₁tr(E)+2λ₂E (8)

f_s是索膜结构受的体积力，σ_l为第二Piola-Lagrange应力张量，σ_p为第二Piola-Kirchhoff应力张量，E为Green-Lagrange应变张量，F为未变形或初始位形上的结构位移梯度张量，x＝ζe_x+ηe_y是位移，ζ、η为检验函数，e_x和e_y分别是x方向的单位向量和y方向的单位向量，λ₁和λ₂是常数，υ为索膜结构泊松比，tr(E)为单位均值的迹。

采用动力平衡方程的半离散有限元方程的准静态形式建立的线弹性有限元模型控制方程如式(9)、(10)所示：

$M_1{\stackrel{\·\·}{d}}_1+D_1{\stackrel{\·}{d}}_1+K_1{d}_1=0---\left(9\right)$

d₁＝d₀(在流固交界面Γ_FSI上) (10)

其中，M₁为流体的虚拟质量，D₁为流体的阻尼，K₁为流体的刚度矩阵，K₁通常是基于整个网格都是均一材料模型计算得到的，d₁是待确定的流体网格位移；d₀表示流固交界面Γ_FSI上的由结构网格的实际运动所引起的边界位移。

流固交界面处的耦合条件包括：流固交界面处线弹性有限元模型和结构的位移耦合条件、流固交界面处流体和结构的应力耦合条件、流固交界面处流体和结构的速度耦合条件。

流固交界面处线弹性有限元模型和结构的位移耦合条件如公式(11)所示：

其中，为流固交界面处线弹性有限元模型的位移，为流固交界面处结构的位移，Γ_FSI为流固交界面。

流固交界面处流体和结构的应力耦合条件如公式(12)所示：

其中，为流固交界面处的结构柯西(Cauchy)应力，为流固交界面处的流体柯西(Cauchy)应力，n_s为结构边界面的外法线方向单位向量，n_f为流体边界面的外法线方向单位向量，n_s＝-n_f。

流固交界面处流体和结构的速度耦合条件如式(13)所示：

其中，为交界面处流体的速度，为交界面处结构的速度。

步骤3：采用线弹性有限元模型处理流体域变形，利用线弹性有限元模型控制方程对流固耦合体系网格进行实时更新。

所述的采用线弹性有限元模型处理流体域变形的具体方法为：将流体网格视为线弹性固体，利用流体单元刚化矩阵中的刚化函数得到线弹性有限元模型单元的单元刚化矩阵，得到流固耦合体系的网格更新公式，对流固耦合体系网格进行实时更新，采用控制网格单元的纵横比控制更新后网格质量。

采用线弹性有限元模型处理流体域变形，利用线弹性有限元模型控制方程得到流固耦合体系网格更新公式，对流固耦合体系网格进行实时更新的方法，如图3所示。

步骤3.1：计算t_n时刻的流固耦合体系网格位移d₁(t_n)，如式(14)所示：

d₁(t_n)＝x(t_n)-x(t_n-1) (14)

其中，x(t_n)是流固耦合体系网格在t_n时刻的坐标，x(t_n-1)是流固耦合体系网格在t_n-1时刻的坐标。

步骤3.2：采用线性弹簧生成线弹性有限元模型控制方程的准静态方程，得到线弹性有限元模型的基本有限元方程如式(15)所示，得到线弹性有限元模型单元的单元刚化矩阵如式(16)所示：

其中，线弹性有限元模型单元的刚化矩阵，B为形状函数的导数矩阵，D为线弹性有限元模型单元的本构矩阵，为线弹性有限元模型单元雅可比行列式，为线弹性有限元模型网格移动的单元刚化函数，Ω为计算边界。

得到根据线弹性有限元模型的流体耦合体系网格更新的单元刚化函数如式(17)所示：

其中，ε₁为线弹性有限元模型网格移动的单元的最大主应变，ε₃线弹性有限元模型网格移动的单元的最小主应变；a是待确定的常数。

步骤3.3：根据公式(16)和公式(17)计算t_n+1时刻的线弹性有限元模型单元的刚化矩阵

步骤3.4：判断流固耦合体系网格单元内切球的半径r和外接球的半径R之比r/R是否小于10^-3，若是，则执行步骤3.5，否则返回步骤3.2。

步骤3.5：将流固耦合体系在t_n时刻得到的流体压力作为t_n+1时刻的单元外力F₂(d₂(t_n+1))，求解方程得到t_n+1时刻的流固耦合体系网格位置d₂(t_n+1)。

对流固耦合体系网格进行更新后，应该控制产生过大纵横比的单元数量，本实施方式中，流固耦合体系网格单元采用3-D四面体单元，采用r/R的比值作为单元的纵横比，r和R分别是流固耦合体系某网格单元的内切球和外接球的半径。这里认为当网格单元内切球的半径r和外接球的半径R之比r/R小于10^-3时，计算结果收敛。否则重新进行计算，直网格单元内切球的半径r和外接球的半径R之比r/R小于10^-3时为止。

步骤4：对流固耦合体系的控制方程和流固交界面处的耦合条件进行时空离散，得到流固耦合体系的各个控制方程的变分弱形式和流固交界面处的各耦合条件的变分弱形式，合并各个控制方程和耦合条件的变分弱形式，得到流固耦合体系的强耦合整体方程。

对流固耦合体系的控制方程和流固交界面处的耦合条件进行时空离散的具体方法为：采用伽辽金有限元法对流固耦合体系控制方程、流固交界面处线弹性有限元模型和结构的位移耦合条件、流固交界面处流体和结构的应力耦合条件进行空间离散，采用隐式有限差分方法对流固交界面处流体和结构的速度耦合条件进行时间离散。

本实施方式中，采用伽辽金有限元法对流体域控制方程进行空间离散，得到流体域控制方程的变分弱形式如式(18)所示：

$f_F={\∫}_{{\mathrm{\Ω}}^f}({\mathrm{\ω}}^{*}\▿\·v_f+-\mathrm{\ζ}\·{\mathrm{\ρ}}_f(\frac{{\∂v}_f}{\∂t}+v_f\·{\▿v}_f)-\▿\mathrm{\ζ}:{\mathrm{\σ}}_f+\mathrm{\ζ}\·f_f)\mathrm{d\Ω}=0---\left(18\right)$

其中，f_F为流体控制方程的变分弱形式，ω^*，ζ为检验函数，Ω^f为流体域。

采用伽辽金有限元法对结构域控制方程进行空间离散，得到结构域控制方程的变分弱形式如式(19)所示：

$f_{S1}={\∫}_{{\mathrm{\Ω}}_0^S}(\▿\mathrm{\η}:{\mathrm{\σ}}_l-\mathrm{\η}\·f_s)\mathrm{d\Ω}=0---\left(19\right)$

其中，f_S1为结构域控制方程的变分弱形式，η为检验函数，为结构域。

采用伽辽金有限元法对线弹性有限元模型控制方程进行空间离散，得到线弹性有限元模型控制方程的变分弱形式如式(20)所示：

$f_{\mathrm{LE}}={\∫}_{{\mathrm{\Ω}}_0^S}(\▿\mathrm{\η}:{\mathrm{\σ}}_l^{\mathrm{LE}})\mathrm{d\Ω}=0---\left(20\right)$

其中，f_LE为线弹性有限元模型控制方程的变分弱形式，为线弹性有限元模型Piola-Lagrange应力张量。

采用伽辽金有限元法对流固交界面处线弹性有限元模型和结构的位移耦合条件进行空间离散，得到流固交界面处线弹性有限元模型和结构的位移耦合条件的变分弱形式如式(21)所示：

$f_F^{\mathrm{FSI}}={\∫}_{{\mathrm{\Γ}}_N^f\∪{\mathrm{\Γ}}_{\mathrm{FSI}}}({\mathrm{\σ}}_f\·n_1^f)\mathrm{d\Γ}=0---\left(21\right)$

其中，为流固交界面处线弹性有限元模型和结构的位移耦合条件的变分弱形式，为流体边界面的外法线方向单位向量，为流体的Neumann边界，Γ为交界面。

采用伽辽金有限元法对流固交界面处流体和结构的应力耦合条件进行空间离散，得到流固交界面处流体和结构的应力耦合条件的变分弱形式如式(22)所示：

$f_{S1}^{\mathrm{FS}1}={\∫}_{{\mathrm{\Γ}}_{N0}^s\∪{\mathrm{\Γ}}_{\mathrm{FSI}0}}((\mathrm{\η}\·n_0^s):{\mathrm{\σ}}_l)\mathrm{d\Γ}+{\∫}_{{\mathrm{\Γ}}_{N0}^s\∪{\mathrm{\Γ}}_{\mathrm{FSI}0}}((\mathrm{\η}\·n_0^s):{\mathrm{\σ}}_l)\mathrm{d\Γ}=0---\left(22\right)$

其中，为流固交界面处流体和结构的应力耦合条件的变分弱形式，是流体边界面的初始外法线方向单位向量，为结构初始Neumann边界流固交界面，Γ_FSI0为初始流固交界面。

采用隐式有限差分方法对流固交界面处流体和结构的速度耦合条件进行时间离散，得到流固交界面处流体和结构的速度耦合条件的变分弱形式如式(23)所示：

$v_f^{n+1}=\frac{1}{\mathrm{\θ}}\frac{{\mathrm{dd}}_{{\mathrm{\Γ}}_{\mathrm{FSI}}}^{u,n+1}-{\mathrm{dd}}_{{\mathrm{\Γ}}_{\mathrm{FSI}}}^{u,n}}{\mathrm{dt}}-\frac{1-\mathrm{\θ}}{\mathrm{\θ}}{v}_f^n---\left(23\right)$

其中，为第(n+1)次迭代得到的流体速度，为第n次迭代得到的流体速度，θ为参数，取值范围为1/2≤θ≤1，为流固交界面上第(n+1)次迭代的结构位移，为流固交界面上第n次迭代的结构位移。

将流固耦合体系控制方程各方程的变分弱形式和流固交界面处的耦合条件的变分弱形式进行合并，得到流固耦合体系的强耦合整体方程，如式(24)所示：

$\left\{\begin{array}{c}{f}_F={\∫}_{{\mathrm{\Ω}}^f}(\mathrm{\ω}\▿\·v_f+-\mathrm{\ζ}\·{\mathrm{\ρ}}_f(\frac{{\∂v}_f}{\∂t}+v_f{\·\▿v}_f)-\▿\mathrm{\ζ}:{\mathrm{\σ}}_f+\mathrm{\ζ}\·f_f)\mathrm{d\Ω}=0\\ f_{S1}={\∫}_{{\mathrm{\Ω}}_0^S}(\▿\mathrm{\η}:{\mathrm{\σ}}_l-\mathrm{\η}\·f_s)\mathrm{d\Ω}=0\\ f_{\mathrm{LE}}={\∫}_{{\mathrm{\Ω}}_0^S}(\▿\mathrm{\η}:{\mathrm{\σ}}_l^{\mathrm{LE}})\mathrm{d\Ω}=0\\ f_F^{\mathrm{FSI}}={\∫}_{{\mathrm{\Γ}}_N^f\∪{\mathrm{\Γ}}_{\mathrm{FSI}}}\mathrm{\ζ}\·({\mathrm{\σ}}_f\·n_1^f)\mathrm{d\Γ}=0\\ f_{S1}^{\mathrm{FSI}}={\∫}_{{\mathrm{\Γ}}_{N0}^S\∪{\mathrm{\Γ}}_{\mathrm{FSI}0}}((\mathrm{\η}\·n_0^s):{\mathrm{\σ}}_l)\mathrm{d\Γ}+{\∫}_{{\mathrm{\Γ}}_{N0}^S\∪{\mathrm{\Γ}}_{\mathrm{FSI}0}}((\mathrm{\η}\·n_0^s):{\mathrm{\σ}}_1)\mathrm{d\Γ}=0\end{array}\right.---\left(24\right)$

步骤5：采用SST K-ω湍流模型模拟湍流，建立SST K-ω湍流模型控制方程，对SSTK-ω湍流模型控制方程进行空间离散，得到SST K-ω湍流模型控制方程的变分弱形式，用SST K-ω湍流模型的变分弱形式替换流固耦合体系的强耦合整体方程中的流体域控制方程的变分弱形式，得到考虑湍流模型的流固耦合体系的强耦合整体方程。

对SST K-ω湍流模型控制方程进行空间离散的具体方法为：采用伽辽金有限元法对SSTK-ω湍流模型的控制方程进行空间离散，即将滑移边界条件作为湍流的边界，将结构壁面的法线方向速度作为垂直于壁面的流体速度，得到离散后的SST K-ω湍流模型的控制方程。

SST K-ω湍流模型模拟湍流的公式如式(25)、(26)所示：

$\frac{\∂\left({\mathrm{\ρ}}_f\mathrm{\ω}\right)}{\∂t}+\frac{\∂}{{\∂x}_j}\left({\mathrm{\ρ}}_f{v}_{f}k\right)=P_k-{{\mathrm{\β}}^{*}\mathrm{\ρ}}_f\mathrm{\ωk}+\frac{\∂}{{\∂x}_j}\left[({\mathrm{\μ}}_p+{\mathrm{\σ}}_f{\mathrm{\μ}}_t)\frac{\∂k}{{\∂x}_j}\right]---\left(25\right)$

$\frac{\∂\left({\mathrm{\ρ}}_f\mathrm{\ω}\right)}{\∂t}+\frac{\∂}{{\∂x}_j}\left({\mathrm{\ρ}}_f{v}_f\mathrm{\ω}\right)=\mathrm{\α}\frac{\mathrm{\ω}}{k}{P}_k-{\mathrm{\β\ρ}}_f{\mathrm{\ω}}^2+\frac{\∂}{{\∂x}_j}\left[({\mathrm{\μ}}_p+{\mathrm{\σ}}_{\mathrm{\ω}}{\mathrm{\μ}}_t)\frac{\∂\mathrm{\ω}}{{\∂x}_j}\right]+{2\mathrm{\ρ}}_f(1-F_1)\frac{{\mathrm{\σ}}_{\mathrm{\ω}2}}{\mathrm{\ω}}\frac{\∂k}{{\∂x}_j}\frac{\∂\mathrm{\ω}}{{\∂x}_j}---\left(26\right)$

其中，k为湍动能，ε为湍动能的耗散率，ω=ε/k，v_fi为流体在水平方向x_i上的速度，v_fj为流体在竖直方向x_j上的速度，x_i为水平方向，x_j为竖直方向，μ_p为为分子扩散所造成的动力粘性，μ_t为流体脉动速度对应的动力粘性，β^*＝0.09。

为了避免在停滞区域中湍流的积累，SST模型引入一个产生限制器如式(27)所示：

${\tilde{P}}_k=\min (P_k,{10\mathrm{\β}}^{*}{\mathrm{\ρ}}_f\mathrm{k\ω})---\left(27\right)$

用φ₁＝{α₁，β₁，σ_k1，σ_ω1}代表原k-ω模型的常数，用φ₂＝{α₂，β₂，σ_k2，σ_ω2}代表变形后k-ω模型的常数，那么SST模型中的常数φ＝{α，β}的公式如式(28)所示：

φ＝F₁φ₁+(1-F₁)φ₂ (28)

本实施方式中，β₁＝0.075，β₂＝0.0828，α_k1＝0.5，α_k2＝1.0，α_ω1＝0.5，α_ω2＝0.856，α₁＝5/9，α₂＝0.44。

其中，F₁为第一个混合函数，

${\arg }_1=\min [\max (\frac{\sqrt{k}}{0.09\mathrm{\ωh}};\frac{{500v}_f}{h^2\mathrm{\ω}});\frac{{4\mathrm{\ρ\σ}}_{\mathrm{\ω}2}k}{{\mathrm{CD}}_{\mathrm{k\ω}}{h}^2}];$

${\mathrm{CD}}_{\mathrm{k\ω}}=\max ({2\mathrm{\ρ\σ}}_{\mathrm{\ω}2}\frac{1}{\mathrm{\ω}}\frac{\∂k}{\mathrm{\ω}{\∂x}_j}\frac{\∂\mathrm{\ω}}{{\∂x}_j},{10}^{-10});$

h是膜结构中心距离最远壁面的距离。

湍流涡粘性定义如式(29)所示：

${\mathrm{\υ}}_t=\frac{a_{1}k}{\max (a_1,\mathrm{\ω};{\mathrm{SF}}_1^{\′})}---\left(29\right)$

其中，υ_t为湍流涡粘性，S为应变率的不变量，F₁′为第二个混合函数： ${\arg }_2=\max (2\frac{\sqrt{k}}{0.09\mathrm{\ωy}};\frac{500v}{y^2\mathrm{\ω}}),a_1=0.31.$

采用伽辽金有限元法对SST K-ω湍流模型的控制方程进行空间离散，将滑移边界条件作为湍流的边界，垂直于壁面的流体速度由固体壁面的法线方向速度给定，结合的壁面摩擦边界条件如式(30)、(31)所示：

其中，k＝1，2，τ_k是流固边界上沿切线方向的正交单位向量，β为指定的常数，本实施方式中，β＝0，为流固交界面处的流体柯西(Cauchy)应力，为流固边界面的外法线方向单位向量。

得到SST K-ω湍流模型控制方程的变分弱形式如式(32)、(33)所示：

$f_{\mathrm{tur}1}={\∫}_{{\mathrm{\Ω}}^f}\left(\mathrm{\λ}\frac{\∂\left({\mathrm{\ρ}}_{f}k\right)}{\∂t}+\mathrm{\λ}\frac{\∂}{{\∂x}_j}\left({\mathrm{\ρ}}_f{v}_{f}k\right)-{\mathrm{\λP}}_k+{{\mathrm{\λ\β}}^{*}\mathrm{\ρ}}_f\mathrm{\ωk}-\mathrm{\λ}\frac{\∂}{{\∂x}_j}\left[({\mathrm{\μ}}_p+{\mathrm{\σ}}_f{\mathrm{\μ}}_t)\frac{\∂k}{{\∂x}_j}\right]\right)\mathrm{d\Γ}---\left(32\right)$

$f_{\mathrm{tur}2}={\∫}_{{\mathrm{\Ω}}^f}\left(\mathrm{\γ}\frac{\∂\left({\mathrm{\ρ}}_f\mathrm{\ω}\right)}{\∂t}+\mathrm{\γ}\frac{\∂}{{\∂x}_j}\left({\mathrm{\ρ}}_f{v}_f\mathrm{\ω}\right)-\mathrm{\γ}(\mathrm{\α}\frac{\mathrm{\ω}}{k}{P}_k-{\mathrm{\β\ρ}}_f{\mathrm{\ω}}^2)+\frac{\∂}{{\∂x}_j}\left[({\mathrm{\μ}}_p+{\mathrm{\σ}}_{\mathrm{\ω}}{\mathrm{\μ}}_t)\frac{\∂\mathrm{\ω}}{{\∂x}_j}\right]-{2\mathrm{\γ\ρ}}_f(1-F_1)\frac{{\mathrm{\σ}}_{\mathrm{\ω}2}}{\mathrm{\ω}}\frac{\∂k}{{\∂x}_j}\frac{\∂\mathrm{\ω}}{{\∂x}_j}\right)\mathrm{d\Γ}---\left(33\right)$

其中，λ和γ为检验函数。f_tur1和f_tur2分别表示SST K-ω湍流模型控制方程的变分弱形式。

用SST K-ω湍流模型的变分弱形式替换流固耦合体系的强耦合整体方程中的流体域控制方程的变分弱形式，即得到考虑湍流模型的流固耦合体系的强耦合整体方程，如式(34)所示：

$\left\{\begin{array}{c}{f}_{\mathrm{tur}1}={\∫}_{{\mathrm{\Ω}}^f}(\mathrm{\λ}\frac{\∂\left({\mathrm{\ρ}}_{f}k\right)}{\∂t}+\mathrm{\λ}\frac{\∂}{{\∂x}_j}\left({\mathrm{\ρ}}_f{v}_{f}k\right)-{\mathrm{\λP}}_k+{\mathrm{\λ\β}}^{*}{\mathrm{\ρ}}_f\mathrm{\ωk}-\mathrm{\λ}\frac{\∂}{{\∂x}_j}[({\mathrm{\μ}}_p+{\mathrm{\σ}}_f{\mathrm{\μ}}_t)\frac{\∂k}{{\∂x}_j}\left]\right)\mathrm{d\Γ}\\ f_{S1}={\∫}_{{\mathrm{\Ω}}_0^S}(\▿\mathrm{\η}{\mathrm{\σ}}_l-\mathrm{\η}\·f_s)\mathrm{d\Ω}=0\\ f_{\mathrm{LE}}={\∫}_{{\mathrm{\Ω}}_0^S}(\▿\mathrm{\η}{\mathrm{\σ}}_l^{\mathrm{LE}})\mathrm{d\Ω}=0\\ f_{\mathrm{tur}2}={\∫}_{{\mathrm{\Ω}}^f}(\mathrm{\γ}\frac{\∂\left({\mathrm{\ρ}}_f\mathrm{\ω}\right)}{\∂t}+\mathrm{\γ}\frac{\∂}{{\∂x}_j}\left({\mathrm{\ρ}}_f{v}_f\mathrm{\ω}\right)-\mathrm{\γ}(\mathrm{\α}\frac{\mathrm{\ω}}{k}{P}_k-{\mathrm{\β\ρ}}_f{\mathrm{\ω}}^2)+\mathrm{\γ}\frac{\∂}{{\∂x}_j}[({\mathrm{\μ}}_p+{\mathrm{\σ}}_{\mathrm{\ω}}{\mathrm{\μ}}_t)\frac{\∂\mathrm{\ω}}{{\∂x}_j}]-{2\mathrm{\γ\ρ}}_f(1-F_1)\frac{{\mathrm{\σ}}_{\mathrm{\ω}2}}{\mathrm{\ω}}\frac{\∂k}{{\∂x}_j}\frac{\∂\mathrm{\ω}}{{\∂x}_j})\mathrm{d\Γ}\\ f_{S1}^{\mathrm{FSI}}={\∫}_{{\mathrm{\Γ}}_{N0}^S\∪{\mathrm{\Γ}}_{\mathrm{FSI}0}}((\mathrm{\η}\·n_0^S):{\mathrm{\σ}}_l)\mathrm{d\Γ}+{\∫}_{{\mathrm{\Γ}}_{N0}^S\∪{\mathrm{\Γ}}_{\mathrm{FSI}0}}((\mathrm{\η}\·n_0^S):{\mathrm{\σ}}_l)\mathrm{d\Γ}=0\end{array}\right.---\left(34\right)$

步骤6：采用Newton方法将考虑湍流模型的流固耦合体系的强耦合整体方程进行线性化，得到线性化的流固耦合体系的强耦合整体方程。

得到线性化的流固耦合系统的强耦合整体方程如公式(35)所示：

$\frac{{\mathrm{df}}_{\mathrm{FSI}}}{\mathrm{dx}}{\mathrm{\Δz}}_i={-f}_{\mathrm{FSI}}^i---\left(35\right)$

其中，f_FSI表示强耦合整体方程，z_i为线性化的流固耦合系统的强耦合整体方程的所有未知量，是雅可比矩阵，i表示迭代时间步，为第i迭代步的风荷载作用力。

步骤7：构造线性化后的强耦合整体方程矩阵的预处理器，得到带有预处理器的流固耦合体系的强耦合整体方程。

构造线性化后的强耦合整体方程矩阵的预处理器的具体方法为：构造线性化后的强耦合整体方程矩阵的分块表达，在雅可比矩阵中使用块状矩阵结构，使得预处理矩阵模块化，得到线性化后的强耦合整体方程矩阵的预处理器。

构造的线性化后的带有预处理器的强耦合整体方程如式(36)所示：

$\left[\begin{array}{cccc}\frac{{\∂f}_{S1}}{{\∂u}_s}& \mathrm{\δ}\frac{{\∂f}_{S1}}{{\∂u}_s^{\mathrm{FSI}}}& & \\ \frac{{\∂f}_{S1}^{\mathrm{FSI}}}{{\∂u}_s}& \mathrm{\δ}(\frac{{\∂f}_{S1}^{\mathrm{FSI}}}{{\∂u}_s^{\mathrm{FSI}}}+\frac{{\∂f}_F^{\mathrm{FSI}}}{{\∂u}_{\mathrm{LE}}^{\mathrm{FSI}}})+\frac{{\∂f}_F^{\mathrm{FSI}}}{{\∂l}_f^{\mathrm{FSI}}}& \frac{{\∂f}_F^{\mathrm{FSI}}}{{\∂l}_f}& \frac{{\∂f}_F^{\mathrm{FSI}}}{{\∂u}_{\mathrm{LE}}}\\ & \mathrm{\δ}\frac{{\∂f}_F}{{\∂u}_{\mathrm{LE}}^{\mathrm{FSI}}}+\frac{{\∂f}_F}{{\∂l}_f^{\mathrm{FSI}}}& \frac{{\∂f}_F}{{\∂l}_f}& \frac{{\∂f}_F}{{\∂u}_{\mathrm{LE}}}\\ & \mathrm{\δ}\frac{{\∂f}_{\mathrm{LE}}}{{\∂u}_{\mathrm{LE}}^{\mathrm{FSI}}}& & \frac{{\∂f}_{\mathrm{LE}}}{{\∂u}_{\mathrm{LE}}}\end{array}\right]\left[\begin{array}{c}{\mathrm{\Δu}}_s^i\\ {\mathrm{\Δl}}_f^{\mathrm{FSI},i}\\ {\mathrm{\Δl}}_f^i\\ {\mathrm{\Δu}}_{\mathrm{LE}}^i\end{array}\right]={-f}_{\mathrm{FSI}}^i---\left(36\right)$

其中，为第i次迭代后流体域的所有未知量的增量，为第i次迭代后索膜结构位移的增量，δ＝θΔt，Δt为时间增量，为交界面上的位移和速度增量的关系，为第i迭代步的线弹性有限元模型的位移的增量，l_f表示流体域的所有未知量，包括流体速度和压力，u_s表示索膜结构位移，u_LE表示线弹性有限元模型的位移，为流固交界面上的流体速度和压力。

步骤8：采用Krylov子空间投影法求解带有预处理器的流固耦合体系的强耦合整体方程，得到强耦合整体方程流体压力、流体速度、索膜结构位移和线弹性有限元模型的位移。

步骤8.1利用第i迭代步得到的强耦合整体方程的流体压力和流体速度索膜结构位移和线弹性模型的位移采用基于线弹性有限元模型的网格更新方法进行流固耦合体系的网格，得到第(i+1)次迭代步计算时的网格体系。

步骤8.2采用Krylov子空间投影法求解带有预处理器的流固耦合体系的强耦合整体方程，得到增量得到第(i+1)次迭代步的强耦合整体方程的流体压力和流体速度索膜结构位移和线弹性模型的位移如式(37)、(38)(39)所示：

$u_s^{i+1}=u_s^i+{\mathrm{\Δu}}_s^i---\left(37\right)$

$l_f^{i+1}=l_f^i+{\mathrm{\Δl}}_f^i---\left(38\right)$

$u_{\mathrm{LE}}^{i+1}=u_{\mathrm{LE}}^i+{\mathrm{\Δu}}_{\mathrm{LE}}^i---\left(39\right)$

步骤9：利用流体压力、流体速度和索膜结构位移进行索膜结构的抗风设计，得到抗风设计参数：索膜结构风压系数、索膜结构风振响应和索膜结构气动特性。

本实施方式中，选用鞍型索膜结构作为索膜结构的典型代表，其几何尺寸示意图如图5(a)所示，基本参数如下：跨度L＝20m，高度h＝5m，矢跨比f/L＝1/8，预张力T＝2.0kN/m，薄膜厚度为1mm，单位面积的质量g＝1.25kg/m²，张拉刚度为Et＝8.0×10⁵N/m，剪切刚度Gt＝1.2×10⁴N/m，泊松比υ＝0.3。

入口平均风剖面服从指数变化规律，模拟C类地貌类型，α为0.22，入口处风速大小为10.32m/s。风向角为0°，定义为沿鞍形膜两高点对角线的方向90°为沿鞍形膜两低点对角线的方向。计算结构响应时，在模拟的初始阶段(t＝45s)考虑结构阻尼的作用以建立稳态流场。

计算域为240m×180m×40m，膜中心离入口的距离为80m，结点数为13.7万，网格数为27万，膜结构网格划分如图4所示，流域上壁面为滑移边界，下壁面和流固耦合边界均为无滑移边界。

先计算鞍型膜结构的风压系数，计算中时间步取为0.005s，这种小时间步会提高计算速度。由于膜结构的柔性较大，此处共进行了25次迭代。给定整个流体域的入口速度后，模拟计算开始，在经历了几个时间步的不稳定压力峰值后，速度场开始收敛，模型设置计算时长t＝200s。

如图5所示，图(a)为本发明具体实施方式中鞍型膜结构的几何尺寸图，图(b)为本发明具体实施方式中鞍型膜结构分区图，根据鞍形屋盖结构布置形式及风压分布特点，将鞍形屋盖分为角区(B1，B3，C1，C3)、边区(B2，C2，D1，E1)和中区(A1，A2，A3，A4)共12个区域，跨度L＝20m，本实施方式中，计算0°和90°风向角下该结构的考虑耦合和不考虑耦合作用时各分区的风压系数，其中0°是沿鞍形膜两高点对角线的方向，90°是沿鞍形膜两低点对角线的方向。

典型风向角下膜结构各分区风压系数如表1所示：

表1 典型风向角下膜结构各分区风压系数

计算表明，膜结构表面均以受风吸力为主，较大风吸力分布在迎风边缘位置且变化梯度较大；这是由于气流在结构迎风处前缘产生的较强分离造成的；不同风向角下鞍形膜面分区风压系数的分布有着比较明显的差异。0°风向角时，部分分区有负风压出现，高负压区出现在90°风向角时，此时所有分区均受风吸力作用。最大负压出现在90°风向角时，说明此时膜面前缘角部来流分离最为严重。

计算了考虑结构阻尼时不同风向角下考虑流固耦合作用和不考虑流固耦合作用(膜结构振动，但不考虑流固耦合效应)时膜结构中点的平均风压系数时程如图6所示，(a)为0°风向角下中点的平均风压系数时程图；(b)为90°风向角下中点的平均风压系数时程图，考虑流固耦合效应和不考虑流固耦合效应(膜结构振动，但不考虑流固耦合效应)的风压系数不同，但均受负风压。两种情况下风压系数的差异是由于考虑耦合效应时风场引起的结构变形改变了周围流场分布，从而进一步改变了结构上风压的分布。

在计算中发现，采用强耦合整体隐式程序计算考虑耦合作用时结构响应时，前50s时程计算大约耗时98小时，而计算不考虑耦合作用时，前50s时程计算大约耗时56小时。即强耦合整体技术比不考虑耦合作用多耗费约43％的机时。另外为了说明本方法的计算效率，采用弱耦合分区法的商业软件ANSYS-CFX对该鞍型结构进行了分析计算，前50s时程计算大约耗时80小时，本方法与其相比约多耗时18％，这是由于本实施方式采用整体隐式法，要计算负责的整体矩阵。

计算比较了更精细网格划分时的情况，结果表明当网格精度提高约20％时，计算精度仅提高约3％，稳定性基本不受影响，而耗费机时却增加了约30％。因此计算结果表明网格划分的精细程度对于计算结果的精度和稳定性影响是不大的。

膜结构抗风设计中非常重要的气动参数是附加质量。所谓附加质量是指结构在空气中发生振动时，将会带动部分周围空气一起运动，因此结构的有效质量除包含结构本身质量外还应包括因结构运动而产生的附加质量。对于膜结构而言，附加质量与薄膜质量可能在同一个数量级上，甚至更大，因此在讨论像薄膜结构这种轻型结构的动力特性时，附加质量的影响是不可以忽略的。

这里结合强耦合整体技术计算得到的膜结构位移，计算了开敞膜结构的附加空气质量在不同风向角下、不同矢跨比下随振型的变化规律，如图7所示。从图7看出，无论风向角或矢跨比如何变化，膜结构的附加空气质量都是随着振型数的增加而减小，且减小幅度较大，即振型阶数越低，附加空气质量越大。第一振型对应的附加空气质量最大，无量纲附加空气质量约在0.54～0.63之间；90°风向角时的附加空气质量小于0°风向角时的附加空气质量，0°风向角时前5阶振型对应的附加空气质量约在0.63～0.64之间，90°风向角时前5阶振型对应的附加空气质量约在0.53～0.60之间；附加空气质量随矢跨比增大而减小，且减小的幅度较大，这主要是由于矢跨比增大，膜结构刚性也增大，使膜结构和空气间的流固耦合效应降低的缘故。

Claims (7)

1.一种基于强耦合整体技术的索膜结构抗风设计方法，其特征在于，包括以下步骤：

步骤1：根据实际膜结构材料的物理性质和实际膜结构的初始预应力，通过找形方法建立索膜结构的初始形态；

步骤2：建立流固耦合体系控制方程和流固交界面处的耦合条件，流固耦合体系控制方程包括：流体域控制方程、结构域控制方程和线弹性有限元模型控制方程；流固交界面处的耦合条件包括流固交界面线弹性有限元模型和结构的位移耦合条件、流固交界面处流体和结构的应力耦合条件、流固交界面处流体和结构的速度耦合条件；

步骤3：采用线弹性有限元模型处理流体域变形，利用线弹性有限元模型控制方程对流固耦合体系网格进行实时更新；

步骤4：对流固耦合体系的控制方程和流固交界面处的耦合条件进行时空离散，得到流固耦合体系的各个控制方程的变分弱形式和流固交界面处的各耦合条件的变分弱形式，合并各个控制方程和耦合条件的变分弱形式，得到流固耦合体系的强耦合整体方程；

步骤5：采用SST K-ω湍流模型模拟湍流，建立SST K-ω湍流模型控制方程，对SSTK-ω湍流模型控制方程进行空间离散，得到SST K-ω湍流模型控制方程的变分弱形式，用SST K-ω湍流模型的变分弱形式替换流固耦合体系的强耦合整体方程中的流体域控制方程的变分弱形式，得到考虑湍流模型的流固耦合体系的强耦合整体方程；

步骤6：采用Newton方法将考虑湍流模型的流固耦合体系的强耦合整体方程进行线性化，得到线性化的流固耦合体系的强耦合整体方程；

步骤7：构造线性化后的强耦合整体方程矩阵的预处理器，得到带有预处理器的流固耦合体系的强耦合整体方程；

步骤8：采用Krylov子空间投影法求解带有预处理器的流固耦合体系的强耦合整体方程，得到强耦合整体方程流体压力、流体速度、索膜结构位移和线弹性有限元模型的位移；

步骤9：利用流体压力、流体速度和索膜结构位移进行索膜结构的抗风设计，得到抗风设计参数：索膜结构风压系数、索膜结构风振响应和索膜结构气动特性。

2.根据权利要求1所述的基于强耦合整体技术的索膜结构抗风设计方法，其特征在于，所述的通过找形方法建立索膜结构有限元模型具体方法为：采用小杨氏模量逐步刚度逼近的方法建立索膜结构的初始形态。

3.根据权利要求1所述的基于强耦合整体技术的索膜结构抗风设计方法，其特征在于，所述的建立流固耦合体系控制方程的具体方法为：采用纳维叶-斯托克斯方程建立流体域控制方程，采用拉格朗日方程建立结构域控制方程，采用动力平衡方程的半离散有限元方程的准静态形式建立线弹性有限元模型控制方程。

4.根据权利要求1所述的基于强耦合整体技术的索膜结构抗风设计方法，其特征在于，所述的采用线弹性有限元模型处理流体域变形的具体方法为：将流体网格视为线弹性固体，利用流体单元刚化矩阵中的刚化函数得到线弹性有限元模型单元的单元刚化矩阵，得到流固耦合体系的网格更新公式，对流固耦合体系网格进行实时更新，采用控制网格单元的纵横比控制更新后网格质量。

5.根据权利要求1所述的基于强耦合整体技术的索膜结构抗风设计方法，其特征在于，所述的对流固耦合体系的控制方程和流固交界面处的耦合条件进行时空离散的具体方法为：采用伽辽金有限元法对流固耦合体系控制方程、流固交界面处线弹性有限元模型和结构的位移耦合条件、流固交界面处流体和结构的应力耦合条件进行空间离散，采用隐式有限差分方法对流固交界面处流体和结构的速度耦合条件进行时间离散。

6.根据权利要求1所述的基于强耦合整体技术的索膜结构抗风设计方法，其特征在于，所述的对SST K-ω湍流模型控制方程进行空间离散的具体方法为：采用伽辽金有限元法对SST K-ω湍流模型的控制方程进行空间离散，即将滑移边界条件作为湍流的边界，将结构壁面的法线方向速度作为垂直于壁面的流体速度，得到离散后的SST K-ω湍流模型的控制方程。

7.根据权利要求1所述的基于强耦合整体技术的索膜结构抗风设计方法，其特征在于，所述的构造线性化后的强耦合整体方程矩阵的预处理器的具体方法为：构造线性化后的强耦合整体方程矩阵的分块表达，在雅可比矩阵中使用块状矩阵结构，使得预处理矩阵模块化，得到线性化后的强耦合整体方程矩阵的预处理器。