大跨度膜结构是20 世纪70 ~ 80 年代发展起来的一种新型张拉结构,由于质量轻、覆盖面积大、造型多变等优点膜结构已发展成为大跨度张拉结构的主流形式之一。膜结构柔性大、质量轻等特点导致了风荷载是膜结构设计中的控制荷载。由于膜结构较强的几何非线性,其抗风设计中需要重点解决也是尚未很好解决的问题就是空气和膜结构间的相互作用,即流固耦合效应。对于膜结构风振中流固耦合问题的研究方法主要以风洞试验和数值模拟为主。随着计算机软硬件技术的迅速提高,这一方法越来越显示出其优势。

目前流固耦合问题的数值模拟方法包括弱耦合分区法、强耦合分区法和强耦合整体法。目前对于膜结构风振中流固耦合效应的数值模拟研究国内外采用的基本是弱耦合分区法和强耦合分区法,其中以弱耦合分区法为主,而对于强耦合整体方法的研究文献还十分有限[1,2],但强耦合整体方法在稳定性和精度等方面却显示出其优势。文献[3]采用强耦合整体方法对风与柔性结构的流固耦合作用进行了研究,并得到了较理想的结果。但由于强耦合整体方法采用一个非线性方程组表示流固耦合系统并进行求解,该非线性方程组的求解计算量非常大,因此如何处理强耦合整体方程,提高其求解效率是提高强耦合整体方法适应性所需解决的关键问题之一。

对于预处理方法的大多数研究是将其百富策略白菜网于流体动力方程的求解中。Kim[4]等将预处理算法百富策略白菜网于边界有限元法中,对包括湍流模型在内的流体动力方程进行了求解;Swanson[5]等采用带有隐式预处理的三阶段Runge-Kutta (RK) 算法求解了基于RANS湍流模型和一方程湍流模型的流体动力学方程;Rossow[6]和Swanson[7]采用全隐式算子用于二维和三维的RANS方程中对流体动力方程进行了求解,并取得了良好结果。

本文在采用Newton-Raphson法线性化非线性强耦合整体方程后,引入预处理方法,在雅可比矩阵中使用特殊的块结构,使得预处理矩阵模块化,分别推导了指定边界条件的线弹性模型预处理矩阵以及强耦合整体方程的预处理矩阵。结合广义极小残余法(generalized minimum residual,简称GMRES)给出了采用预处理方法的强耦合整体方法的求解步骤,最后将该方法百富策略白菜网于风与膜结构的流固耦合分析。

系统的强耦合整体方程可以写为,具体形式如下[3]

方程中的未知量包括流体速度uf,压力pf和线弹性模型位移 χps,以及结构位移 χs。该强耦合整体方程采用Newton-Raphson法线性化,即得到线性方程组

由于A为大型稀疏线性方程组的系数矩阵。广义极小残余法(generalized minimum residual,简称GMRES) 采用预处理方法,是求解大型非对称线性方程组的Krylov子空间投影法[8],因其数值稳定性好,被认为是求解大型非对称线性方程组最适合的迭代方法之一。但一般由于A的条件数很大,计算中易发生方程组病态,为了提高Krylov子空间方法的稳定性,可以在求解过程中采用引入预处理的方法,以改善矩阵的条件数,提高计算精度和稳定性。这里则是采用预处理的GMRES方法求解线性化后的强耦合整体方程。

这里采用结构的实际变形作为线弹性模型的边界位移条件,这样离散的流体方程就不直接与实际结构变形相关,线弹性模型也不与流体直接相关,这样就可以在雅可比矩阵中使用特殊的块结构,使得预处理矩阵模块化,可以在计算中反复使用该预处理矩阵。

对式(6) 中的虚拟位移加入拉格朗日乘数项,即边界上的应变和约束应满足下列条件,

式中第一项表示体系应变,第二项便是体系的约束。其中

其中 ξ 是拉格朗日坐标,ζ 是边界坐标,即 );Λ 是拉格朗日乘数,可以看成是作用在边界 ΓI上的表面牵引力,目的是产生所需的边界变形。

对于公式(6) 中给定的边界位移us,将ups对于结点坐标进行离散,

其中N表示网格中结点的数量,Xij表示第j个结点的第i个坐标,φj表示结点j的有限元基函数。这里采用相同的基函数对拉格朗日坐标进行扩展,

其中 ξi是第i个拉格朗日坐标,这里将拉格朗日坐标的结点值设为结点的初始欧拉位置,即 Θij= X0ij,这样可以保证未变形的位形上没有应力作用。这里只在约束边界 ΓI上定义拉格朗日乘数场并扩展为,

其中?表示在约束边界上的结点, 是通过将内部基函数约束在边界 ΓI上获得的基函数。Lij是未知量拉格朗日算子的坐标分量。

那么带有边界条件(6) 的线弹性模型的预处理矩阵可以写作,

式中Rps表示线弹性模型的雅可比矩阵,均表示分块矩阵,其中E是为约束线弹性模型问题的切线刚度矩阵,非对角线上的两分块矩阵Cxl和Cl X互为转置,即Cxl= ClTX,由在位移约束上施加拉格朗日乘数产生;X和L表示未知量向量。

方程(13) 的预处理形式可以写为[9,10]

其中 Υps是高效预处理器,W是在保证扩展矩阵Eaug= E + CxlW-1Clx稀疏条件下的任意矩阵,选取原则是尽量是使采用Newton-Raphson方法线性化后的体系中的Eaug和W求解简便。这里选用

其中 δ = ‖E‖∞

对于采用Newton-Raphson方法线性化后的线性体系,有如下的块状矩阵结构,

其中矩阵Rs是雅可比矩阵,其对角线上的矩阵块包括N-S雅可比行列式F,受流体荷载作用的结构切线刚度矩阵行列式S,以及在式(13) 中定义的线弹性模型雅可比行列式Rps。向量F包含了流体域N-S方程中的速度和压力未知量;S表示受流体荷载作用发生变形的结构域未知量;X表示流体网格的结点位置,L表示加在边界上的离散拉格朗日乘数;Cfx表示线弹性模型未知量(流体网格的结点位置) 对离散的N-S方程的影响分块矩阵Csf表示受流体荷载作用的结构表面的流体牵引力影响,该牵引力也会受流体网格节点位置的影响,这种依赖关系用矩阵Csx表示。Cls表示方程(6) 中的流体网格指定边界位移目前是与受流体荷载作用的结构位移场相关。

将式(13) 中定义的线弹性模型预处理方程带入到式(16) 中,得到流固耦合问题的预处理矩阵为,

从上式可以看出,

对上式做如下矩阵变换,

此时在方程组(7) 的两端左乘预处理矩阵 Υps-1,得到

令 ,得到

再用GMRES法求解方程组(20)。

总结以上预处理模型和强耦合整体方法,基于循环型GMRES(m) 算法[11],给出了采用预处理的强耦合整体方程的求解步骤:

(1) 将流体域、结构域和线弹性模型的所有变量都以初始场为零开始迭代,并进行2 - 3 步连续迭代,线性化流体对流项和Euler-Lagrange应变张量非线性部分,所有流场变量均存储在网格结点上;隐式化处理整个流固耦合体系的边界条件和交界面处的荷载。

(2) 采用Newton-Raphson方法线性化强耦合整体方程,得到线性方程组A·x = b,左乘预处理矩阵,得到方程(20);采用GMRES法求解方程(20),方法如下述。

(3) 随机选取初始估计x0,确定Krylov子空间维数m,定义矩阵Hm∈ R(m+1) ×m,并初始化其所有元素hi为0,其中Hm∈ R(m+1) ×m是上Hessenberg矩阵。

(4)选择适当大小的m,开始执行Arnoldi过程:

(1)计算u0=b-Ax0,α=‖u0‖2,λ1=u0/α

(2)计算cj=Auj,q=γ-1pscj(j=1,2…,m)

(3)计算 {(i=1,2…,j),pj+1,j=‖c‖2,uj+1=q/pj+1,j,定义矩阵Cm=[c1,c2,…,cm]

(5)计算xm=x0+Cmym

其中ym= arg miny‖αη1- Hmym‖

(6) 当满足 ‖b - Axm‖ ≤ ε 时停止迭代(ε 是设定的误差限),否则令x0= xm,回到步骤(5) 重新计算。

(7) 将求解结果保存为多维动态数组,在所有的时间步循环过后,将体系的位移、压力、速度等最后结果存储、输出。

这里以典型的鞍型膜结构为例,百富策略白菜网以上的强耦合整体方法程序进行风振响应分析计算。膜结构的计算简图如图1 所示,其基本参数如下:跨度L = 20 m,高度H = 5 m,矢跨比f/L = 1 /8,预张力T = 2. 0 k N/m,薄膜厚度为1 mm,单位面积的质量g = 1. 25 kg /m2,张拉刚度为Et= 8. 0 × 105N / m,剪切刚度Gt= 1. 2 × 104N / m,泊松比" = 0. 3。入口风剖面服从指数变化规律,平均风速为12 m / s,模拟C类地貌类型,α 为0. 22,则入口出风速大小为10. 32 m/s。风向角为0°,定义为沿鞍形膜两高点对角线的方向,90°为沿鞍形膜两低点对角线的方向。计算时结构阻尼比 ξ = 0. 02[13]。另外,计算结构响应时,在模拟的初始阶段(t = 45 s)考虑结构阻尼的作用以建立稳态流场,而在流固耦合模拟的主要阶段(即当流场达到真实流态(旋涡充分发展、脱落)及膜结构进入真实的振动状态)。电脑硬件条件为:Intel i5 3450S;CPU频率,2 800 MHz;内存容量:4GB DDR3,1 600 MHz。

雷诺数Re = 1. 8 × 106,湍流强度I = 0. 025,湍动能k = 0. 300,湍流耗散率 ε = 5. 931,其中湍流强度I参照日本规范[12]中的第Ⅱ类地貌取值。湍流模型采用标准k-ε 模型。计算域为240 m × 180 m × 40 m,膜中心离入口的距离为80 m,结点数为13. 7 万,网格数为27 万,膜结构网格划分如图2 所示,流域上壁面为滑移边界,下壁面和流固耦合边界均为无滑移边界。

这里计算了膜结构的风压分布,如图3 所示,为验证本文算法的正确性,将数值模拟计算结果与风洞实验结果(图4)进行了对比,通过对比可以看出,数值模拟结果与风洞实验结果总体分布规律一致,但有一定的差异性,这是由于该风洞实验为刚性模型,并未考虑流固耦合作用。从数值模拟的整体分布规律来看数值模拟结果还是可靠的。

这里还计算了不同风向角下考虑流固耦合作用和不考虑流固耦合作用(膜结构振动,但不考虑流固耦合效应)时膜结构中点的风压时程,结果如图5 所示。

从图5 可以看出,考虑流固耦合效应和不考虑流固耦合效应(膜结构振动,但不考虑流固耦合效应)的风压系数不同,但均受负风压。两种情况下风压系数的差异是由于考虑耦合效应时风场引起的结构变形改变了周围流场分布,从而进一步改变了结构上风压的分布。这表明采用预处理后强耦合整体方法的计算结果仍是可靠的。

表1 给出了强耦合整体方程在预处理前后矩阵条件数和计算时间(采用Gauss-Jordan消去法)的比较。

表2 给出了采用预处理前后不同时间步和迭代次数时的收敛情况。

分析表1 和表2 中的结果,可以看出:

(1)当采用预处理后系数矩阵的条件数由处理前的4. 296 2 × 106降低到处理后的231. 560 9,矩阵条件数大大减小,此时再用GMRES算法,计算时间大大缩短:预处理后达到同样指定收敛值所需时间比预处理前减少约43% 的时间,;计算达到收敛指定值所需迭代次数也大大减少,计算效率大大提高。预处理前条件数很大的原因是由雷诺方程离散化得到的系数矩阵的条件数偏大的缘故。

(2)在同样迭代次数的情况下,采用预处理后的计算残差要比预处理前小得多,所需计算时间也比预处理前大大减少,约较少41% 的时间,再次证明采用预处理方法后可以大大提高计算精度,提高计算效率。

表3 给出了GMRES迭代次数随不同雷诺数不同时间的变化规律。

表4 给出了GMRES迭代次数随时间步和拉格朗日乘数场的变化规律。

分析表3 和表4,可以看出,

(1)GMRES迭代次数随雷诺数和时间的增加而增加,但值得注意的是,在雷诺数相差很大时,GMRES迭代次数增加的速度并不是很快,说明采用预处理方法可以有效地减少计算耗时,提高计算效率。

(2)GMRES迭代次数随时间步而增大,随拉格朗日乘数场的增大而减少,再次说明预处理器的引入提高了计算效率。

图6 和图7 分别给出了未采用预处理和采用预处理后计算残差随计算时间步和计算总时间的变化。

分析图6 和图7,可以看出:采用预处理方法后计算残差随时间步和计算总时间的增大而迅速降低,达到预定残差的计算时间与未采用预处理相比大大减少,计算效率得到了较大提高。

本文在采用Newton-Raphson法线性化非线性强耦合整体方程后,引入预处理方法,,并结合GMRES法总结了采用预处理方法的强耦合整体方法的求解步骤,将该方法百富策略白菜网于风与膜结构的流固耦合分析,得到的主要结论有:

(1)采用预处理方法后仍可准确计算膜结构的风振响应、风压系数等重要参数。

(2)采用预处理方法后计算精度和效率大大提高,计算时间节约40%以上,计算效率大大提高。

(3)预处理方法对于强耦合整体方法是一种高效的处理方式。