膜结构柔性大、非线性强, 其控制荷载是风荷载。当膜结构产生风致振动时, 会产生明显的空气流体和结构相耦合的现象, 即流固耦合效应。在过去的十几年里, 各国学者对流体-固体耦合计算问题作了大量的研究, 并发展了一些新的计算方法[1,2]。Peskin首先提出的浸入边界法[3,4] (IB) 就是其中之一, 它是将与流体具有同样密度的浸入弹性体 (纤维) 看成是Navier-Stokes方程中的一系列等价体来建模的。同时, 许多学者也对浸入结构形状、质量、柔性进行了大量的研究[5,6], 因此近几年浸入边界法已经被成功地百富策略白菜网于不同的领域[7]。这种在浸入边界法中引入虚拟域的方法, 通常被称为浸入物体法 (IOM) 。对浸入物体法进行求解时, 一般需要计算Jacobian矩阵, 浸入物体法中的Jacobian矩阵的存储量要求在n2数量级上, 而对于有上百万个自由度的大型体系来说, Jacobian 矩阵需要1012的存储量, 这远远超出了一般计算机的容量。

目前对于采用浸入物体法研究膜结构风致振动中的流固耦合效应, 国内外尚处于发展阶段, 而对于采用带有预定条件的无矩阵迭代法求解浸入物体法建立的膜结构和空气形成的流固耦合体系模型的研究还少见报道, 在建模和求解过程中还有很多问题需要进一步研究。

本文给出了采用浸入物体法 (IOM) 的膜结构和空气流体的有限元方程, 提出采用带有预定条件的无矩阵Newton-Krylov 迭代算法求解浸入物体法, 该方法可以避免计算Jacobian 矩阵, 并引入预定条件矩阵。同时, 将提出的方法百富策略白菜网于一双坡型膜结构的风振耦合分析中, 得到了比较满意的结果。

考虑典型的流固系统, 浸入变形体的固体域为Ωs (t) 。为解决流体域不断变化的问题, 这里引入虚拟域, 即引入对于固体域Ωs (t) 来说占据同样空间域的人造流体域。流固耦合系统的主要未知量是流体速度v、流体压力p以及固体压力ps。

定义Sobolev空间, 那么控制方程的弱形式可以定义为

∫ΩhwiINvIρv?hidΩ−∫ΓhfwiINvIfΓhfidΓ+∫Ωh(wiINvI,jτij−phwiINvI,i)dΩ+∫Ωhs[wsiJNuJ(ρs−ρ)(v?hi−gi)+wsiJNuJ,j(σsij−σfij)]dΩ−∫ΩhwiINvIρgidΩ+∫ΩhqlNpI(vhj,j+ps,hκ)dΩ+∫ΩhsqsJNpJ(J3−1+ps,hκs)⋅dΩ=0 (1)∫ΩhwiΙΝΙvρv?ihdΩ-∫ΓfhwiΙΝΙvfiΓfhdΓ+∫Ωh(wiΙΝΙ,jvτij-phwiΙΝΙ,iv)dΩ+∫Ωsh[wiJsΝJu(ρs-ρ)(v?ih-gi)+wiJsΝJ,ju(σijs-σijf)]dΩ-∫ΩhwiΙΝΙvρgidΩ+∫ΩhqlΝΙp(vj,jh+ps,hκ)dΩ+∫ΩshqJsΝJp(J3-1+ps,hκs)⋅dΩ=0(1)

式中 NvIΙv和NpIΙp分别为节点I的速度向量和压力的内插函数;

vI, wI, pI, qI分别为离散速度向量的节点值;

NuJ 和NpsJJps分别为节点J的唯一向量和未知压力;

ws, h, ps, h, qs, h分别为离散位移向量的节点值;

ρ为流体压力;

τij为剪应力;

fΓhfiiΓfh为作用在流体域上的外力;

ps为固体压力;

σsij 和σfij分别为固体域和流体域的Cauchy应力;

v?hiv?ih为I点的加速度;

gi为重力加速度;

κs为固体体积模量;

J3为变形斜率的行列式。

本文提出采用预定条件技术无矩阵Newton-Krylov 迭代算法来求解以上的浸入物体法。在无矩阵Newton-Krylov 迭代算法中, 无需形成Jacobian矩阵[8,9]。无矩阵Newton-Krylov迭代法中, 最重要的一步就是用有限差分计算代替w=Jqi,

Jqi≅r(Θm+1,k−1+eqi)−rΘm+1,k−1e (2)Jqi≅r(Θm+1,k-1+eqi)-rΘm+1,k-1e(2)

式中 J为Jacobian矩阵;

qi为n维Krylov子空间向量;

r为固体节点I的核函数;

Θm+1, k-1为第m+1时间步的第k-1次Newton-Raphson迭代的近似解;

e通常为计算机误差平方根值 (文献[8]) 。

在建立了n×n阶矩阵Hn和n× (n+1) 阶矩阵HˉˉˉˉnΗˉn后, 对于j=1～n, i=1～ (j-1) , 对Hn进行因式分解:

hij=cihij+sih(i+1)j,h(i+1)j=−sihij+cih(i+1)j,r=h2jj+h2(j+1)i−−−−−−−−−−√,cj=hjj/r,sj=h(j+1)j/r (3)hij=cihij+sih(i+1)j,h(i+1)j=-sihij+cih(i+1)j,r=hjj2+h(j+1)i2,cj=hjj/r,sj=h(j+1)j/r(3)

式中 hij为矩阵Hn中的元素;

r, c和s为Hn中元素计算得到的系数。

最后的解向量可以表达为,

ΔΘk,m=ΔΘk,0+∑i=1nyiqi (4)ΔΘk,m=ΔΘk,0+∑i=1nyiqi(4)

式中 ΔΘk, m为第m时间步的第k次迭代解;

ΔΘk, 0为初始向量;

yi为Krylov子空间内的列向量。

如果初始向量ΔΘk, 0在较小的Krylov子空间Kn内没有产生预期的结果, 那么初始向量可更新为ΔΘk, n, 然后继续迭代过程, 直到最后获得准确的解。

采用本文算法, 对一双坡型膜结构和空气间的流固耦合作用进行了计算模拟, 计算模型如图1所示。V10为10 m高处风速, H为建筑物檐口高度, 纵坐标y表示高度。模拟B类地貌, 流动雷诺数Re=1.75×107。边界条件为:进口风速25 m/s, 上方和右方出口压力为零, 下方边界速度为零。空气密度取为1.21 kg/m3, 空气粘性取为17.9×10-6 Pas。膜材厚度为3 mm, E=3.3×108 N/m2, 单位面积质量为1.25 kg/m2。流体域的时间步为ΔtF=1.25×10-3, ΔtS=1.25×10-4。湍流模型采用标准k-ε模型。

图2和图3分别给出了在不同时刻膜结构的风压分布和风速矢量图。

为了说明预定条件的高效性, 本文分别采用带有预定条件和不带预定条件的无矩阵Newton-Krylon迭代法计算了上述算例, 并对计算结果进行了对比。时间步的总数为50, 其中分别采用了3种衡量准则, 即每一时间步的Newton-Raphson迭代数、每一Newton-Raphson迭代的GMRES迭代, 以及在短时间内的总CPU时间。最后结果如表1所示。

从表中可以看出, 预定条件方法与无预定条件方法相比, 效率要高得多。但无论是否采用带有预定条件方法, 最后的计算结果都在允许的迭代误差范围之内。

对于膜结构风振中的流固耦合问题, 本文采用浸入物体法 (IOM) 对膜结构和空气流体进行建模, 并采用带有预定条件的无矩阵Newton-Krylov 迭代算法求解浸入物体法。利用上述方法对一双坡型膜结构与空气的耦合进行了计算, 得到以下结论:

(1) 浸入物体法可以作为膜结构这种高柔性结构与空气耦合作用问题的建模方法。

(2) 带有预定条件的无矩阵迭代算法, 可以用于求解浸入物体法, 并能取得准确结果。

(3) 预定条件的无矩阵迭代法在计算膜结构与风的流固耦合问题时, 可以得到准确结果, 且与无预定条件方法相比, 可以使计算效率大大提高。