虚拟区域法在流固耦合问题中的应用1.doc

《虚拟区域法在流固耦合问题中的应用1.doc》由会员分享，可在线阅读，更多相关《虚拟区域法在流固耦合问题中的应用1.doc（15页珍藏版）》请在三一文库上搜索。

1、精品论文虚拟区域法在流固耦合问题中的应用1及春宁 天津大学建筑工程学院天津市港口与海洋工程重点实验室，天津 (300072) E-mail：摘要：本文将多相流领域内的虚拟区域法引入到流固耦合问题的分析中，将固体视为应变 率 为零 的虚拟 流体 ，对流 体和 虚拟流 体均 以速度 和压 强作为 基本 变量， 均采 用 Navier-Stokes 方程作为控制方程，同时求解流体域和虚拟流体域得到整个计算域的流场分 布，应用分布式拉格朗日乘子法在虚拟流体域上施加刚体约束以保持虚拟流体的刚体外形和 运动形式，最终建立一套流固耦合模型及其数值求解方法。通过对粒子流问题和流固耦合问 题进行数值模拟，验证了

2、此模型的正确性和求解大变形/运动流固耦合问题的有效性。 关键词：虚拟区域法，流固耦合，联合控制方程，Navier-Stokes 方程，有限元 中图分类号：O357.11引言流固耦合问题涵盖广泛，应用范围涉及航空、航天、水利、建筑、石油、化工、海洋 以及生物等领域。如液体晃动对火箭飞行稳定性的影响，大型贮液罐在地震激励作用下受力， 液体湍振对输液管道的影响，以及波浪与结构物的相互作用等。流固耦合的重要特征为两相 介质之间的相互作用。即：固体在流体力的作用下产生变形和运动和固体的变形和运动使流 场发生变化。两种作用相互耦合，同时进行。流固耦合问题的求解方法主要分为交替模式求解方法和耦合模式求解方法

3、。 交替模式求解方法1-4通过将耦合系统分割为两个独立的子域（流体域和固体域），根据各个子域的物理特性分别应用适用的数值方法交替对子域进行求解，并通过流固接触边界 上的速度、压强相容条件达到两个子域的耦合。因此，其具有以下优点：无需重新建立联合 控制方程，对子域的求解可以利用非常完善的现有计算模型和程序，降低了耦合系统的求解 难度；针对各个子域的物理特性分别采用相应的基本变量（对流体采用速度和压强为基本变 量，对固体采用位移为基本变量），使子域的控制方程得以简化。但是，弱耦合方法也存在 自身难以弥足的缺点。如：有条件收敛，过大的时间步长会导致求解发散，因此对时间步长 的限制非常严格；即便子域模

4、型具有高阶的时间精度，整个耦合系统的时间精度却仅为一阶。耦合模式求解方法5-10同步求解流体域和固体域，将流体和固体的控制方程进行联立， 流体响应和固体响应同时由一套联合控制方程中解得，因而是控制方程层面上的耦合，即完 全意义上的流固耦合。耦合方法克服了交替方法在收敛性和求解精度方面的不足，放宽了对 时间步长的限制，可以应用于求解大变形、强耦合的流固耦合问题，更适于从本质上揭示复 杂的流固耦合过程。但由于流体和固体在物理性质和动力特征上存在的明显差异，两者的控 制方程在形式上截然不同，建立一套流固耦合联合控制方程仍存在诸多难点。虚拟区域法11-13是近年来发展起来的多相流动数值模拟领域内的一种

5、数值方法，主要用 于模拟刚体粒子在流场中的运动问题。虚拟区域法通过将刚体视为虚拟流体，对流体域进行 扩展，将包含固体（动边界）的随时间变化的流体域转化为不随时间变化的扩展流体域（真 实流体域+虚拟流体域）；在整个扩展域上以速度和压强为基本变量，以 Navier-Stokes 方程 为控制方程，计算并得到扩展域的流场；应用分布式拉格朗日乘子法在虚拟流体域上施加刚 体约束，使虚拟流体保持刚体外形和运动形式，同时使虚拟流体对周围流场的影响与刚体粒 子相同；最终得到刚体粒子在流场中的运动过程。虚拟区域法的主要特点是：此模型是无条 件收敛的，具有较高的求解精度和计算效率；刚体粒子与流体的接触边界被转化为

6、流体之间1本课题得到高等学校博士学科点专项科研基金（项目编号：20070056120）的资助。-15-（真实流体和虚拟流体）的内部边界，其上的速度、压强相容条件在控制方程中自动满足，大大简化了求解过程；虚拟流体的刚体约束条件通过分布式拉格朗日乘子施加，该乘子实际 上代表了单位体积的虚拟流体保持刚体外形和运动形式所需要的附加体积力；采用固定的欧 拉网格对扩展流体域进行离散，避免了采用拉格朗日网格所导致的网格畸变以及消除畸变所 必需而又复杂的网格重构算法。由以上可知，虚拟区域法对流体和刚体粒子的基本变量和控制方程进行了统一描述，建 立了流体和刚体粒子的联合控制方程，同步求解刚体粒子的运动和周围流场

7、的分布。因此， 其与流固耦合问题中的耦合模式求解方法类似，两者之间可以相互借鉴。本文在求解流固耦合问题时借鉴虚拟区域法，将刚体视为应变率为零的虚拟流体，对流 体和虚拟流体均以速度和压强作为基本变量；根据虚拟流体域上的刚体约束条件构造虚拟流 体质点的运动控制方程，并与真实流体的控制方程进行合并，得到流固耦合联合控制方程； 同时求解流体域和虚拟流体域得到整个计算域的流场分布，并应用分布式拉格朗日乘子法在 虚拟流体域上施加刚体约束以保持虚拟流体的刚体外形和运动形式；最终建立完全意义上的 流固耦合模型及其数值求解方法。2流固耦合控制方程2.1 控制方程f设整个区域由 Rd (d = 2, 3) 表示，

8、被密度为 、粘性为 的粘性不可压缩流体充满， 其中有可移动的刚体 B，如图 1 所示。图中， 为整个计算域的外边界， B 为流固接触边界。图 1 二维流固耦合问题示意图Fig. 1. A schematic diagram of 2-D FSI problems流体的控制方程为 Navier-Stokes 方程，in B(t )(1)(2)with u0 = 0(3) ( u + (u )u) = g + f tf u = 0in B(t)u( x, 0) = u0 ( x),x B(0),u = g0on ,with g0 nds = 0(4)控制方程在 B(t ) 上的边界条件由公式(13

9、)确定。对于牛顿流体，应力张量 定义为： = p + (u + ut )(5)d d du(= ui i =1 ) 和 p 分别表示流体的速度和压强， x(= xi i =1 ) 表示 R 中的空间点， n 表示外边 界 的单位外法线向量， g0 表示外边界 上的本质边界条件， 为单位对角矩阵， 为梯 度算子， g 为体积力向量， 表示流体的动力粘性系数， f 表示流体的密度，上标 t 表示矩阵的转秩。刚体的控制方程为欧拉方程（牛顿运动定律），M dV = Mg + F i + F e (G, )dtI d I = T i + T e (G, )dt(6) (7)dG = V(8)dtd =

10、(9)dtG(0) = G0 ,(0) = 0 ,V (0) = V 0 ,(0) = 0(10)其中，G 和 为刚体的重心坐标和转动角度，V 和 分别为刚体的平动速度和角速度，M和 I 分别为刚体的质量和转动惯量。 F i 为流体作用在刚体上的力，表示为：F i = Bnds(11)T i 为流体作用在刚体重心 G 上的力矩，表示为：JJJKT i = BGx (n)ds(12)F e 、T e 分别为流固耦合系统的外界（如计算域的外边界 ）作用在刚体上的力和力矩，与刚体的形状、位置以及刚体的重心坐标和转动角度有关，视具体情况确定。由于流固接触边界为无滑移边界，流固边界上流体质点的速度与刚体

11、质点在其上的速度 相同，即：JJJJJJKu( x, t) = V (t) + (t) G(t) x,x B(t) .(13)1.2 基于虚拟区域法的联合控制方程引入泛函空间，gV ( B(t) = u | u (H 1 ( B(t)d ,0JJJJJJK(14)u = g0 (t) on ,u = V (t) + (t) G(t) xon B(t)V ( B(t) = v,Y , | v (H 1 ( B(t)d ,0v = 0 on ,JJJJJJKv = Y + G(t) xon B(t), with Y Rd , R3 (15)0L2 ( B(t ) = q | q L2 ( B(t

12、), B (t ) qd = 0(16)其中，v 、Y 和 分别为流体质点速度、刚体平动速度和刚体角速度的变分， q 为压强的变 分。对公式(1)和(2)、公式(6)和(7)应用虚功原理，可得：u f B (t )(+ (u )u) vd + 2t B (t )D(u) : D(v)d B (t )p vd (17) g vd +(n) vds,= f B (t ) Bv V0 ( B(t ) B (t )q ud = 0,q L2 ( (M dVdt Mg Fi F e) Y = 0,( I ddt I Ti T e) = 0,B(t )(18)Y V0 ( B(t )(19) V0 ( B

13、(t )(20)JJJK其中， D(u) = (u + ut ) / 2 。根据变分 v 的定义，在边界 B 上 v = Y + Gx 。由此，可将公式(17)中的边界积分项改写为：B(n) vds = BJJJK(n) (Y + Gx)dsJJJK(21)= B(n)ds Y + BGx (n)ds = F i Y T i 由此可知，边界积分项即为流体作用在刚体上的力和力矩。回代，公式(17)可改写为u= f f B (t )(+ (u )u) vd + 2t B (t )i iD(u) : D(v)d B (t )p vd (22) g vd F B (t ) Y T ,v V0 ( B(

14、t )累加公式(22)、(19)和(20)，可得流固耦合控制方程的变分形式为：对于任意 t 0 ，求解g 0u(t ), p(t ), G(t), (t), V (t), (t) ，使得 u(t) V ( B(t) ， p(t) L2 ( B(t) ， G(t) Rd ，0(t) R3 ，V (t) Rd ， (t) R3 ，并且满足：u f B (t ) ( t+ (u )u) vd + 2 B (t )D(u) : D(v)d B (t )p vd + (M dV Mg F e ) Y + ( I d I T e ) dtdt(23)= f g vd , B (t )v,Y , V0 (

15、B(t ) B (t )q ud = 0,q L2 ( B(t )(24)以及刚体运动方程公式(8)和(9)和初始条件公式(3)和(10)。对于整个流固耦合系统， F i 、T i为内部作用，在公式(23)中与流体边界积分项相互抵消。此外，在流固耦合边界上，流体速g度自动满足无滑移边界条件(隐含在泛函空间V ( / B(t ) 的定义中)。由此可知，在流固耦0合边界上，速度、压强相容条件自动得到满足，无需显式体现在求解过程中，大大降低了耦合系统的求解难度。 以上流固耦合控制方程仅是对流体和刚体的控制方程进行简单累加，并未得到流固耦合的联合控制方程，流体和刚体仍是割裂开的（公式(23)和(24)

16、仅在 / B(t) 上积分）。为此，本文对以上变分问题应用虚拟区域法。采用以下两个步骤。(1) 将刚体代之以“虚拟”流体，并在其上施加刚体约束。(2) 应用分布式拉格朗日乘子法，“放松”刚体约束，使问题可解，得到基于虚拟区域 法的流固耦合联合控制方程。步骤(1)：用虚拟流体代替刚体 B(t) 。在刚体域内引入速度场 u ，且满足刚体约束条件， 即虚拟流体内部质点的运动控制方程为：JJJJJJKu( x, t) = V (t) + (t) G(t) x,x B(t )(25)为了建立流固耦合联合控制方程，须使虚拟流体内部的控制方程与真实流体控制方程在 形式上相似。为此，将公式(25)对时间求全导

17、数，并在 B(t) 内积分，可得：u f (+ (u )u) vd tB (t ) = fsM dVdt Y +fs( I d0dt I ) ,v,Y , V ()(26)其中， 为刚体的密度，空间V () 定义为：s0V0 () = v,Y , | v,Y , V0 ( B(t),JJJJJJKx B(t),(27)v( x, t) = Y + G(t) x, x B(t)与空间V0 ( / B(t) 相比，空间V0 () 增加了在区域 B(t ) 内部的定义。出于同样目的，对体积力向量 g 在 B(t) 内积分，有： fv gg Yv Y (28) f B (t ) d =M s, , ,

18、 V0 ()0此外，根据空间V () (公式(27)中 v 在 B(t ) 内的定义（刚体运动条件），有： 因此， v = 0in B(t),v,Y , V0 () (29)B (t ) p vd = 0,同理，根据公式(25)，有：v,Y , V0 ()(30)可得，D(u) = 0 u = 0in B(t )in B(t )(31) (32)D uD vv Y (33)B (t ) ( ) :( )d = 0, , , V0 ()q ud = 0,q L2 (B(t)(34)B (t ) 将公式(26)、(28)、(33)和(34)累加合并，可得形如真实流体域上的流固耦合控制方程的虚拟 流

19、体域上的流固耦合控制方程为：u f (+ (u )u) vd + 2 D(u) : D(v)d p vd B (t ) tB (t ) B (t ) = dVdv gd + f M Y +f ( I I ) (35)f B (t ) dtdt fsMg Y ,ss0v,Y , V () q ud = 0,q L2 ()(36)将公式(35)和(36)与公式(23)和(24)合并，可得流固耦合的联合控制方程为：对于任意g 0t 0 ，求解u(t), p(t ), G(t), (t), V (t), (t) ，使得 u(t) W () ， p(t) L2 () ，G(t) Rd ，0(t) R3

20、，V (t) Rd ， (t) R3 ，满足：u f (+ (u )u) vd + 2 D(u) : D(v)d p vd t dVd+ (1 f )(M Y + ( I I ) ) F e Y T e (37)sdtdt f= f g vd + (1 s)Mg Y ,v,Y , V0 () q ud = 0,q L2 ()JJJK 0 0 0(38)u( x, 0) = u0 ( x),x B(0),u( x, 0) = V+ G x,x B(0)(39)且满足刚体约束条件(25)、刚体运动方程公式(8)和(9)和初始条件公式(10)。上式中，泛函空g间W () 定义为：01dgW () =

21、 v | v (H () ,0v = g0 (t )on 步骤(2)：在以上的联合控制方程中，公式(25)表示刚体约束条件，即刚体域 B(t) 内的任 意一点的速度必须满足此条件，由此得到的解为控制方程的“强解”。为了得到其有限元“弱 解 ”，需 要“放松 ”刚体约 束条件。 为此，引 入分布式 拉格朗日 乘 子(t), (t) (t) = (H 1 (B(t)d ，对刚体约束条件在刚体域 B(t ) 上进行内积，并使刚体域内部的速度满足积分刚体约束条件，得到了基于虚拟区域法的流固耦合联合控制方程的变分形式为：对于 任意 t 0 ，求 解 u(t), p(t), G(t), (t), V (t

22、), (t), (t) ，使得 u(t ) Wg0() ，0p(t) L2 () ， G(t) Rd， (t) R3，V (t) Rd， (t) R3， (t ) (t ) ，并且满足： f ( u+ (u )u) vd + 2 D(u) : D(v)d p vd t JJJK f ( , v Y Gx) (t ) + (1 )(M dV Y + ( Id I ) )sdtdt(40) F e Y T e = f g vd + (1 f )Mg Y ,0v (H 1 ()d ,Y Rd ,s R3 q ud = 0,q L2 ()JJJJJJK(41)( , u(t ) V (t ) (t )

23、 G (t ) x) (t ) = 0, (t )(42)和刚体的运动方程(8)和(9)、刚体的初始条件公式(3)和流体的初始条件(39)。上式中，内积(, ) (t ) 定义为：或者，( , v ) (t ) =( v + 2 : v )d ,B (t ) , v (t )(43)( , v ) (t ) =( v + 2 D( ) : D(v )d ,B (t ) , v (t )(44)其中， 为特征长度，可取值为刚体的最大直径。与公式(37)相比，公式(40)增加了内积项JJJK( , v Y Gx) (t ) ，此项代表了虚拟流体保持刚体外形和运动形式所必需的附加体积力。1.3 联合

24、控制方程的时间离散格式对流固耦合联合控制方程应用破开算子法，可得联合控制方程的时间离散格式为： 给定初始时刻的 u0 , p0 , 0 , G0 , 0 , V 0 , 0 。步骤 1：已知第 n (n 0) 时间步的 un , pn , n , G n , n , V n , n ，求解 un +1 4 ，满足：un +1 4 unt vd +(un+1 4 )un+1 4 vd = 0,v W0 h(45)un +1 4 W n +1g0 h步骤 2：求解 un + 2 4 和 pn + 2 4 ，满足： un + 2 4 un +1 4 vd + 2 D(un + 2 4 ) : D(v

25、)d pn + 2 4 vd f t h= f g vd ,v W0 hq un + 2 4 d = 0, q L2(46)un + 2 4 W n +1 ,pn + 2 4 L2g0 h0 h步骤 3：预估刚体的位置、角度、速度和角速度。赋初值 G n + 2 4,0 = G n ，n + 2 4,0 = n ，V n + 2 4,0 = V n ， n + 2 4,0 = n 。令 k = 1, N ，执行以下循环。V n + 2 4,k = V n + 2 4,k 1 + (t / N ) g + (t / 2N )(1 )1 M 1 F e (G n + 2 4,k 1 , n + 2

26、 4,k 1 )f sf s n + 2 4,k = n + 2 4, k 1 + (t / 2N )(1 )1 I 1T e (G n + 2 4, k 1 , n + 2 4,k 1 )G n + 2 4,k = G n + 2 4,k 1 + (t / 4N )(V n + 2 4,k + V n + 2 4, k 1 ) n + 2 4, k = n + 2 4,k 1 + (t / 4N )( n + 2 4,k + n + 2 4, k 1 ) n + 2 4,k n + 2 4,k 1 1 1 en + 2 4,kn + 2 4, kf sV= V+ (t / N ) g + (

27、t / 4N )(1 ) M F (G, )f s+ (t / 4N )(1 )1 M 1F e (G n + 2 4,k 1 , n + 2 4,k 1 )(47)f sn + 2 4,k = n + 2 4, k 1 + (t / 4N )(1 )1 I 1T e (G n + 2 4, k , n + 2 4,k )f sG+ (t / 4N )(1 )1 I 1T e (G n + 2 4,k 1 , n + 2 4,k 1 ) n + 2 4,k = G n + 2 4,k 1 + (t / 4N )(V n + 2 4,k + V n + 2 4, k 1 )n + 2 4, k

28、= n + 2 4,k 1 + (t / 4N )(n + 2 4,k + n + 2 4, k 1 )令 G n + 2 4 = G n + 2 4, N ， n + 2 4 = n + 2 4, N ，V n + 2 4 = V n + 2 4, N ， n + 2 4 = n + 2 4, N 。步骤 4：求解 un +3 4 ， n + 3 4 ，V n + 3 4 ， n + 3 4 ，满足： f un + 3 4 un + 2 4 vd + (1 f )(MV n +3 4 V n + 2 4 Y+ ( Itn + 3 4 n + 2 4tst+ In + 2 4 n + 2 4

29、) )JJJJJJJKBn+2 4 = ( n + 3 4 , v Y G n + 2 4 x),hJJJJJJJKv W0 h ,Y Rd , R3(48)h( , un +3 4 V n + 3 4 n +3 4 G n + 2 4 x)n+2 4 Bhh= 0, n + 2 4gun +3 4 W n +1 ,0 hn + 3 4n + 2 4 ,V n + 3 4 Rd ,n +3 4 R3步骤 5：校正刚性固体的位置、速度和角速度。令 G n +1,0 = G n + 2 4 ， n +1,0 = n + 2 4 ，V n +1,0 = V n +3 4 ， n +1,0 = n +

30、3 4 。令 k = 1, N ，执行以下循环。V n +1, k = V n +1,k 1 + (t / 2N )(1 )1 M 1F e (G n +1,k 1 , n +1, k 1 )f sf s n +1,k = n +1,k 1 + (t / 2N )(1 )1 I 1T e (G n +1,k 1 , n +1,k 1 )G n +1,k = G n +1,k 1 + (t / 4N )(V n +1,k + V n +1,k 1 ) n +1,k = n +1, k 1 + (t / 4N )( n +1,k + n +1,k 1 ) n +1, k n +1,k 1 1 1

31、en +1,kn +1,kf sV= V+ (t / 4N )(1 ) M F (G, )f s+ (t / 4N )(1 )1 M 1F e (G n +1,k 1 , n +1,k 1 )(49)f sn +1,k = n +1,k 1 + (t / 4N )(1 )1 I 1T e (G n +1,k , n +1,k )f s+ (t / 4N )(1 )1 I 1T e (G n +1,k 1 , n +1, k 1 )= G+ V+ VG n +1,kn +1,k 1( t / 4N )(n +1,kn +1,k 1 )n +1,k = n +1, k 1 + (t / 4N )(

32、n +1,k + n +1,k 1 )令 G n +1 = G n +1, N ， n +1 = n +1, N ，V n +1 = V n +1, N ， n +1 = n +1, N 。步骤 6：令 un +1 = un + 3 4 ， pn+1 = pn+ 2 4 ， n+1 = n +3 4 。在公式(47)和公式(49)中，N 表示刚体运动预测和校正的分步数，本文取 N=10。显而易见，步骤 1 为纯对流问题，可以通过 Characteristic 法求解。步骤 2 和步骤 4 为一般化的 Stokes问题，可以通过 Uzawa 共轭梯度法求解。步骤 3 和步骤 5 为代数迭代算法。

33、3流固耦合模型正确性验证为了检验以上流固耦合模型，本文将其应用于 Glowinski11-13的二维粒子流算例，以期 对其正确性进行验证。此算例中，一个四周封闭的矩形区域内充满了粘性不可压缩流体，初 始时刻在流体中放置一个圆形的刚体。在重力的作用下，刚体克服流体的阻力开始自由下落。 在下落的过程中，刚体由于受到的流体阻力不断增大，其下落速度逐渐趋于稳定；流体则由 于受到刚体运动的影响，在刚体的周围产生速度场。与 Glowinski 的算例相同，本文取矩形 区域的高为 5.0 cm 宽为 2.0 cm，圆形刚体的直径为 0.25 cm，其初始位置为(1.0, 4.0) cm，流fs体的密度 =

34、1.0 g / cm3 ，流体的粘性系数 = 0.01 g(cm s) ，刚体的密度 = 1.5g / cm3 ，重力加速度 g = 981.0 cm / s2 ，如图 2 所示。图 3 给出了刚体粒子下落过程中的流函数等值线图，可以看出粒子后部出现了类似卡门涡街的漩涡，并随时间脱落。图 4 给出了刚体粒子的速度和位移时程曲线，并与 Glowinski 的结果进行了对比。可以看出两者吻合较好。图 2 粒子流的平面布置和边界条件Fig. 2. The layout and the boundary conditions of the particulate flow(a) 刚体粒子竖向速度时程曲

35、线(a) Time history of the rigid bodys vertical velocity图 3 粒子下落过程中的流函数等值线图Fig. 3. The contour of the stream function during the falling of the rigid body(b) 刚体粒子竖向位移时程曲线(b) Time history of the rigid bodys vertical displacement图 4 刚体粒子下落过程中的竖向速度和竖向位移时程曲线Fig. 4. Time histories of the rigid bodys vertic

36、al velocity and vertical displacement4流固耦合问题算例及分析本文设计了以下两个流固耦合算例。算例中，一个边长为 0.25 m 的正方形结构物放置 在一个长 5.0 m、高 1.0 m 充满粘性流体的长方形域内，如图 5 所示。流体域的左边界为均匀来流边界，速度分别为V0= 1.0m / s 和V0= 2.0m / s ，右边界为出流边界，上下边界为不可滑移边界。结构物的密度 = 2000kg / m3 ，流体的密度 = 1000kg / m3 ，流体的粘性系数s ff = 0.001kg /(m s) ，重力加速度 g = 9.81m / s2 ，结构物与

37、边界的滑动摩擦系数 = 0.1 。图 5 流固耦合问题平面布置及边界条件Fig. 5. The layout and the boundary conditions of the FSIF + Fproblem图 6 结构物的力学简化模型Fig. 6. The simplified mechanical model of the structure为了对结构物的运动进行模拟，需要对结构物进行简化。本文将刚性结构物简化为质量弹簧摩擦系统，如图 6 所示。图中，G 为重力，Fcl 和 Fcr 分别为左右弹簧反力，Ffl 和 Ffr分别为左右摩擦力，表示为： FFcl ,V ,=Fflhx clcrV(50) F x ,