膜结构是近年来国内应用较多的一种新型张力结构体系, 因其对风荷载的作用较为敏感, 所以如何准确地确定结构风振响应一直是膜结构抗风设计中的关键问题。膜结构在风荷载作用下的受激振动问题在理论上可归结为不可压缩粘性流体与几何非线性弹性体之间的非定常流固耦合问题[1]。目前, 在计算流体动力学领域 (computational fluid dynamics) , 对结构大变形下流固耦合问题的有效求解是十分具有挑战性的。由于对这一问题进行数值模拟的计算量十分巨大, 国际上求解类似问题通常采用具有矢量并行处理器 (vector parallel processor) 的超级计算机[2]。因此, 在现有情况下应用个人计算机对膜结构的风振响应进行简化的二维数值计算具有明显的现实意义。

本文结合膜结构的风振响应特点, 采用流固耦合问题数值计算方法, 通过一定的简化, 实现了对膜结构风场的二维数值模拟。通过数值计算, 可以对膜结构的风振响应的机理进行探讨。二维数值模拟可以在一定程度上反应出膜结构的风振响应的真实过程, 为膜结构的设计提供帮助。同时, 二维数值计算验证了本文对流固耦合问题的求解方法的可行性, 为实现对膜结构的三维风场数值模拟奠定基础。

在ALE (arbitrary lagrangian-eulerian) 描述中[3], 有限元的剖分是针对独立于结构和流体的参考坐标系进行的, 网格点即为参考点, 它可以在空间以任意形式运动。流体域和结构域的物理量可以通过Jacobi行列式映射到参考坐标系上, 从而实现运动描述方式的统一。

对于二维不可压缩流体而言, 控制方程包括质量守恒方程、动量守恒方程和本构方程, 当引入湍流模型时还要增加对流-扩散方程。

质量守恒方程

∂ui∂xi=0   (1)$\frac{\partial \text{u}{}_i}{\partial \text{x}{}_i}=0(1)$

动量守恒方程

∂ui∂t+(uj−uˆj)∂ui∂xj=fi+1ρ∂σij∂xj   (2)$\frac{\partial \text{u}{}_i}{\partial t}+(\text{u}{}_j-\widehat{\text{u}}{}_j)\frac{\partial \text{u}{}_i}{\partial \text{x}{}_j}=\text{f}{}_i+\frac{1}{\rho }\frac{\partial \text{σ}{}_{ij}}{\partial \text{x}{}_j}(2)$

本构方程

σij=−pδij+μ(∂ui∂xj+∂uj∂xi)   (3)$\text{σ}{}_{ij}=-\text{p}\text{δ}{}_{ij}+\mu (\frac{\partial \text{u}{}_i}{\partial \text{x}{}_j}+\frac{\partial \text{u}{}_j}{\partial \text{x}{}_i})(3)$

其中, ui (uj) 是流体速度;p是压力;uˆ$\widehat{\text{u}}$i (uˆ$\widehat{\text{u}}$j) 是网格速度;ρ是空气密度;μ是动力粘性系数;δij是三阶单位矩阵;fi是流体体积力;σij是柯西应力张量。将本构方程 (3) 代入动量守恒方程 (2) 中, 导出Navier-Stokes方程

∂ui∂t+(uj−uˆj)∂ui∂xj=fi−1ρ∂p∂xi+μρ∂2ui∂xj∂xi   (4)$\frac{\partial \text{u}{}_i}{\partial t}+(\text{u}{}_j-\widehat{\text{u}}{}_j)\frac{\partial \text{u}{}_i}{\partial \text{x}{}_j}=\text{f}{}_i-\frac{1}{\rho }\frac{\partial \text{p}}{\partial \text{x}{}_i}+\frac{\mu }{\rho }\frac{\partial {}^2\text{u}{}_i}{\partial \text{x}{}_j\partial \text{x}{}_i}(4)$

质量守恒方程与Navier-Stokes方程组成求解ui和p的基本方程。

对结构质点运动的描述采用Lagrange坐标系, 有限元的剖分是针对结构进行的, 网格点随结构一起运动。

动量守恒方程

ρˉ∂2dˉi∂t2=σˉij,i+ρˉfˉi   (5)$\bar{\rho }\frac{\partial {}^2\bar{\text{d}}{}_i}{\partial t{}^2}=\overline{\text{σ}}{}_{ij,i}+\bar{\rho }\bar{\text{f}}{}_i(5)$

平衡条件

σˉijnj=stˉi   (6)$\overline{\text{σ}}{}_{ij}\text{n}{}_j={}^{\text{s}}\bar{\text{t}}{}_i(6)$

本构方程

σˉij=Dijkl⋅εkl   (7)$\overline{\text{σ}}{}_{ij}=\text{D}{}_{ijkl}\cdot \text{ε}{}_{kl}(7)$

其中, dˉ$\bar{\text{d}}$i是结构位移;ρˉ$\bar{\rho }$是材料密度;Dijkl是材料弹性张量;εkl是无穷小应变张量;fˉi$\bar{\text{f}}{}_i$是结构体积力;stˉ$\bar{\text{t}}$i是外部施加的表面应力;nj是结构表面的法向量。

流固耦合问题是流体域和结构域之间的非线性动态耦合问题, 其基本原理是通过分别满足流体域和结构域之间的耦合边界上的运动学平衡方程和动力学平衡方程, 将流体域和固体域耦合起来[4]。

运动学条件

df=ds   (8)$d{}_{\text{f}}=d{}_{\text{s}}(8)$

动力学条件

n⋅τf=n⋅τs   (9)$\text{n}\cdot \tau {}_{\text{f}}=\text{n}\cdot \tau {}_{\text{s}}(9)$

对于膜结构, 流固耦合边界属于无滑移的情况, 耦合边界上流体速度条件可以从运动学条件导出

vf=ds   (10)$v{}_{\text{f}}=d{}_{\text{s}}(10)$

其中, df和ds分别为流体域和结构域的耦合边界位移;τf和τs分别为耦合边界上流体和结构应力;vf为耦合边界流体域速度;n为流固耦合边界的法向方向矢量。

结构域和流体域的有限元模型是分开建立的, 两个模型之间的耦合是通过分别在其实际的作用面定义流固耦合边界实现的, 两者的网格可以不一致, 但要满足一定的容差要求。耦合面上流体节点的位移通过附近结构节点的位移插值得到, 而结构节点受到的流体作用力通过附近流体应力的积分得到

F(t)=∫hdτfdS   (11)$F(t)={\displaystyle \int h{}_{\text{d}}\tau {}_{\text{f}}\text{d}S}(11)$

其中, hd是结构虚位移。

如图1, 在流固耦合边界上, Nf和Ns分别代表流体单元和结构单元节点; f和s分别表示流体域和结构域节点上物理量 (位移、应力) ;fˆ$\widehat{f}$和sˆ$\widehat{s}$则分别表示各物理量在流固耦合边界处结构域和流体域节点上的插值。它们的映射关系表示如下

从结构域到流体域

{sj|j∈Ns}→{sˆi|i∈Nf}   (12)$\{s{}_j|j\in {\rm N}{}_{\text{s}}\}\to \{\widehat{s}{}_i|i\in {\rm N}{}_{\text{f}}\}(12)$

从流体域到结构域

{fi|i∈Nf}→{fˆj|j∈Ns}   (13)$\{f{}_i|i\in {\rm N}{}_{\text{f}}\}\to \{\widehat{f}{}_{\text{j}}|\text{j}\in \text{Ν}{}_s\}(13)$

在二维情况下, 用施加预应力的索单元来模拟膜结构。结构域的耦合边界对应着流体域位于同一几何位置的两条边界, 并且由于它们被结构单元隔离开, 两条边界上的流体变量是不同的。因此, 流体域耦合边界两侧的流体单元节点应该交错分布。在耦合边界的角点处, 流体域的流体变量是相同的, 因而两条耦合边界的流体单元节点应在此重合[5]。膜结构二维流固耦合边界处的流体单元节点分布详见图2。

膜结构的绕流属于大雷诺数、低速流动。本文选用大涡模拟 (large eddy simulation, LES) 作为湍流的数值计算方法。把湍流运动看成是由许多大小不同的漩涡组成, 质量、动量和能量的输运主要通过大漩涡实现, 小漩涡的作用表现为耗散。大涡模拟将流动变量划分为大尺度量和小尺度量, 通过某种滤波方式将小尺度部分滤出;小尺度运动对大尺度运动的影响在运动方程中表现为亚格子 (SGS) 雷诺应力项[6,7]。本文采用高斯滤波和Smagorinsky亚格子模型。

大涡模拟的基本方程为经过滤波后的连续性方程和Navier-Stokes方程

∂uˉi∂xi=0   (14)∂uˉi∂t+∂(uˉiuˉj)∂xj=−1ρ∂pˉ∂xi+μρ∂2uˉi∂xj∂xj−∂τij∂xj   (15)τij=−2νT(Sˉij)   (16)νT=(cΔ)2[∂uˉj∂xi(∂uˉi∂xj+∂uˉj∂xi)]   (17)$\begin{array}{l}\frac{\partial \bar{\text{u}}{}_i}{\partial \text{x}{}_i}=0(14)\\ \frac{\partial \bar{\text{u}}{}_i}{\partial t}+\frac{\partial (\bar{\text{u}}{}_i\bar{\text{u}}{}_j)}{\partial \text{x}{}_j}=-\frac{1}{\rho }\frac{\partial \bar{\text{p}}}{\partial \text{x}{}_i}+\frac{\mu }{\rho }\frac{\partial {}^2\bar{\text{u}}{}_i}{\partial \text{x}{}_j\partial \text{x}{}_j}-\frac{\partial \text{τ}{}_{ij}}{\partial \text{x}{}_j}(15)\\ \text{τ}{}_{ij}=-2\nu {}_{\text{Τ}}(\bar{\text{S}}{}_{ij})(16)\\ \nu {}_{\text{Τ}}=(c\Delta ){}^2\left[\frac{\partial \bar{\text{u}}{}_j}{\partial \text{x}{}_i}\left(\frac{\partial \bar{\text{u}}{}_i}{\partial \text{x}{}_j}+\frac{\partial \bar{\text{u}}{}_j}{\partial \text{x}{}_i}\right)\right](17)\end{array}$

式 (14) ～ (17) 中, uˉ$\bar{\text{u}}$i (uˉ$\bar{\text{u}}$j) 为滤波后的流体速度, pˉ$\bar{\text{p}}$为滤波后的压力, τij为亚格子应力张量 (SGS) , Sˉ$\bar{\text{S}}$ij为滤波后的变形速率张量, νT为涡粘系数, c为亚格子应力常数, 取值范围为0.1～0.23, Δ为空间网格宽度。

在流体域, 本文采用隐式时间积分:TR-BDF法 (Trapezoidal Rule / Backward-Differentiation-Formula) 对Navier-Stokes方程进行时间域的离散求解[8]。对于瞬态方程:∂u∂t=f(u)$\frac{\partial u}{\partial t}=f(u)$, 假定已完成时间步t的求解, 在t+Δt时刻, TR-BDF方法将通过连续的两个时间子步对瞬态方程进行离散求解

u(t+γΔt)=u(t)+γΔtf[u(t+12γΔt)]   (18)u(t+Δt)=u(t+βγΔt)+(1−α)Δtf[u(t+Δt)]   (19)$\begin{array}{l}u(t+\gamma \text{Δ}t)=u(t)+\gamma \text{Δ}tf[u(t+\frac{1}{2}\gamma \text{Δ}t)](18)\\ u(t+\text{Δ}t)=u(t+\beta \gamma \text{Δ}t)+(1-\alpha )\text{Δ}tf[u(t+\text{Δ}t)](19)\end{array}$

式中, u (t+βγΔt) = (1-β) u (t) +βu (t+γΔt) , γ=2-1/α, β=α2/ (2α-1) 。当12<α<1$\frac{1}{2}<\alpha <1$时, TR-BDF法具有二阶精度并且是无条件稳定的。尽管把每一个时间步分为两个时间子步进行求解会使计算耗费的资源翻倍, 但由于TR-BDF方法计算的准确性较好, 可以采用较少的时间步数, 从而可使整体上花费的计算时间减少。

对于结构非线性动力响应分析, 本文采用Newmark直接积分法。这种方法具有二阶精度和无条件稳定性[9]。

对耦合系统的求解包括迭代求解和直接求解, 这两种方法都要保证动态分析时结构模型和流体模型之间时间积分点的一致性;由于耦合系统的非线性, 两种方法都要通过迭代才能得到有限元方程的解。迭代解法中, 流体域和结构域的求解变量是完全耦合的, 但是对流体域和结构域方程分别求解, 其中利用耦合系统求解另一部分的最新信息[10]。本文采用直接解法, 其有限元模型是将流体模型和结构体模型装配在同一矩阵中进行求解。

耦合系统的整体求解向量为:X= (Xf, Xs) , Xf和Xs分别为流体域和结构域的求解向量;耦合系统的有限元方程组为:F= (Ff, Fs) =0, Ff=0和Fs=0分别为流体域和结构域的有限元方程组。本文采用Newton-Raphson法求解耦合系统的有限元方程组

Xk+1=Xk−[∂F(Xk)∂X]−1F(Xk)   (20)∂F(X)∂X=⎡⎣∂Ff∂Xf ∂Ff∂Xs∂Fs∂Xf ∂Fs∂Xs⎤⎦=[Aff AfsAsf Ass]   (21)$\begin{array}{l}\text{X}{}^{k+1}=\text{X}{}^k-\left[\frac{\partial F(\text{X}{}^k)}{\partial \text{X}}\right]{}^{-1}F(\text{X}{}^k)(20)\\ \frac{\partial F(\text{X})}{\partial \text{X}}=\left[\begin{array}{c}\frac{\partial F{}_{\text{f}}}{\partial X{}_{\text{f}}}\frac{\partial F{}_{\text{f}}}{\partial X{}_{\text{s}}}\\ \frac{\partial F{}_{\text{s}}}{\partial X{}_{\text{f}}}\frac{\partial F{}_{\text{s}}}{\partial X{}_{\text{s}}}\end{array}\right]=\left[\begin{array}{c}A{}_{\text{f}\text{f}}A{}_{\text{f}\text{s}}\\ A{}_{\text{s}\text{f}}A{}_{\text{s}\text{s}}\end{array}\right](21)\end{array}$

其中, k为迭代次数。迭代求解过程中, 当计算残差小于预设的容差时, 计算结束, 所获得的求解向量为有限元方程组的解。

这里需要特别注意等效矩阵A的集成过程, 具体步骤:

(1) 采用单一的结构体的有限元模型的集成过程, 对Fs和Ass进行集成。

(2) Ff和Aff的集成过程采用流固耦合边界的速度条件, 见式 (10) 。流固耦合边界处的等效矩阵Afs同时进行集成。流固耦合边界的流体节点i的速度平衡方程和等效矩阵计算过程

[Ff]i=vi−[ds(t+Δt)−ds(t)Δt]i   (22)[Aff]i=I‚[Afs]i=[−1Δt]i   (23)$\begin{array}{l}\left[F{}_{\text{f}}\right]{}_i=v{}_i-\left[\frac{d{}_{\text{s}}(t+\text{Δ}t)-d{}_{\text{s}}(t)}{\text{Δ}t}\right]{}_i(22)\\ \left[A{}_{\text{f}\text{f}}\right]{}_i=\text{Ι}‚\left[A{}_{\text{f}\text{s}}\right]{}_i=\left[-\frac{1}{\text{Δ}t}\right]{}_i(23)\end{array}$

(3) 按式 (11) , 考虑流体应力对结构的作用, 更新结构平衡方程Fs=0, 并计算Asf。对于流固耦合边界的结构节点j, 流体的作用力添加到结构平衡方程和等效矩阵的过程如下

在数值计算中, 结构单元选用2节点索单元。假定索是理想柔性的, 只能承受拉力, 并引入大位移小应变假定。膜材材性:E=900MPa;ρ=1068kg/m3;ν=0.3。流体单元为3节点流体单元。由于结构域的耦合边界对应着流体域位于同一几何位置的两条边界, 流体域耦合边界两侧的流体单元节点交错分布。两条耦合边界的流体单元节点在耦合边界的角点处重合。单元剖分局部详图见图3, 上述特点在图3中均得到体现。流体域共含有22605个节点和44500个单元。流域上壁面为滑移边界, 下壁面和流固耦合边界均为无滑移边界。

在算例1中, 模拟张拉膜的拉索长1m, 索单元施加初应变后具有4.737MPa的预应力。计算过程中, 入流流速从零逐渐增加至10m/s;将流速达到10m/s时作为零时刻, 步长0.25s, 持时10s。所模拟流动雷诺数的数量级为10-6。整体流域范围和索单元的变形见图4。索单元最大水平位移是1.527mm。为方便观察, 图4中单元的变形放大了1500倍。在t=5s时, 流场的整体和局部流线图见图5和图6, 流域相对压力分布见图7。流线图显示出背风区由上下两个反向的大的涡系组成, 它们的旋转方向分别为顺时针和逆时针。主涡系在水平方向的尺度约为拉索长度的3倍。对应于静压力分布图, 两个主涡系的中心亦是负压区的两个中心, 并对应着极大负压值。最大正压值出现在迎风面拉索的中部, 与最大节点位移相对应。最大负压的绝对值明显大于最大正压。

大涡模拟可以较好的模拟流场的脉动特性, 在结构上表现为索单元的变形和应力表现出脉动变化。图8为索中点最大水平向位移的时程曲线, 图9为索单元应力的时程曲线。对比两图发现, 索单元的节点位移和内力基本上同步变化。索单元应力的脉动均值为4.742MPa, 比找形后索单元的应力值大0.12%。节点位移和索单元应力的脉动幅值均很小。

索单元的变形很小, 对流场的影响十分有限。为证明算法的有效性, 本文利用计算流体动力学软件FLUENT6.0 平台对相同的流体刚性模型进行了模拟, 流场布局和单元划分与算例1完全相同, 湍流模型仍选择大涡模拟。在t=5s时流线图见图10, 流域压力分布见图11。

通过对比发现, 流线和静压力分布的趋势十分相近, 说明本文的算法是有效的。定量的分析发现, 背风区FLUENT 6.0计算的静压力分布偏大约15%。

算例2模拟的结构体型为典型的全开敞张拉膜结果, 与算例1相比更具有一般性。流场布局和膜结构的几何尺寸见图12。张拉膜结构曲面的稳定是依靠经向和纬向两个主轴方向反向的曲率来保障的。为模拟二维情况下反向的张拉刚度, 对称布置了4条拉索, 但图中的虚线部分所代表的拉索不参与流固耦合计算, 只提供弹性的张拉刚度。计算过程中, 入流流速为12m/s, 步长0.2s, 持时为60s。t=30s时, 索单元节点变形见图13。

索单元在流场中的受力以上吸力为主, 变形主要表现为竖向变形, 这与我国荷载规范中的类似结构的体型系数的规定是相符合的。在流速为12m/s时, 最大竖向位移出现在背风面索单元的中部, 详见图13。由于迎风面高度很小, 体型接近流线型, 流场受到的影响很小, 在流线图上表现为流线间没有相互干扰。入流速度稳定时, 流场对结构的影响比较稳定, 脉动变化幅度很小。

与算例1相比, 索单元的应力变化十分显著。索单元找形后, 具有4.95MPa的预应力。在入流流速为12m/s的流场作用下, 索单元的应力范围是5.20～5.22MPa, 增幅达到5%。这证明流场对膜结构的影响在一些情况下比较明显。算例1中, 索单元应力变化较小的主要原因是拉索两端固定后, 很难产生弹性变形。

本文结合膜结构的风振响应特点, 详细阐述了膜结构流固耦合计算中的结构域和流体域的网格关系以及耦合系统的求解过程。采用流固耦合问题数值计算方法, 实现了二维情况下对膜结构风场的数值模拟。计算表明, 湍流模型选择大涡模拟可以较好的模拟二维情况下膜结构在流场中的脉动反应。对于一般性的算例, 瞬态流场中膜的应力增幅较大, 这说明考虑流固耦合对于膜结构的风振响应计算的必要性。膜结构的体型对风振响应的计算结果会产生显著的影响, 合适的体型可以减小膜结构风振响应的脉动幅度, 同时减少流场的扰动变化。二维的数值计算验证了本文对流固耦合问题的求解方法的可行性, 为实现膜结构的三维风场数值模拟打下基础。