流固耦合界面信息传递方法研究

流固耦合界面信息传递方法研究

随着计算力学的发展和计算机计算计算的提高，多物理场的集成分析开始广泛应用于不同的领域，如航空、医学、磁学等。在结构工程领域，结构风梯度的破坏日益严重。为了揭示风和结构的耦合效应，有必要根据流固耦合理论进行分析，并从物理本质上揭示风和结构的耦合效应规律。此外，该项目结构也日益与其他学科相结合。例如，关于太阳能热态屋顶结构的研究包括分析压力场、位移场、磁体场和电炉等物理和化学性质的耦合效应。为能对复杂的耦合场进行分析,需要在不同领域的计算程序间搭建信息传递平台,因此,合理有效的界面信息传递方法是实现耦合分析的关键环节之一.本文针对流固耦合分析中的CFD、CSD间的界面信息传递进行系统的阐述和研究,旨在为实现流固耦合分析提供参考,文中多数方法可拓展至其他领域.1cfd和csd流的基本原理流固耦合的分析过程一般可以描述为:在流场(如风、水流等)作用下,结构特别是柔性结构(如大跨、高耸、膜结构等)会产生较大的变形和振动,因此会对周围流场会产生较大影响,而流场的改变会进一步改变作用在结构表面上的压力大小,从而形成流体与结构的相互耦合作用,图1.要合理有效地在CFD和CSD流之间进行数据交换,必须遵循一定的原理和法则.结合文献,概括为以下4个基本原理.1.1流体耦合边界位移分量的uf、vf、wf运动学连续性条件即在流固界面上流体、固体对应点位移、速度的一致性.位移连续条件可表示为:df=ds(1)写成分量的形式:式中,df为流体在耦合边界的位移,dfx、dfy、dfz分别为df在x、y、z方向的位移分量;ds为结构在耦合边界的位移,dsx、dsy、dsz分别为ds在x、y、z方向的位移分量.速度连续条件可表示为:uf=us(3)写成分量的形式:式中,uf为流体在耦合边界的速度,uf、vf、wf为uf在x、y、z方向的速度分量;us为结构在耦合边界的速度,us、vs、wf为us在x、y、z方向的速度分量.1.2耦合边界上的应力动力连续性条件即接触界面要满足力的守恒.根据界面上任一点的力的平衡:σs·ns=σf·nf(5)式中,σs、σf分别为结构、流体在耦合边界上的柯西应力张量;ns、nf分别为结构、流体在耦合边界上的外法线方向矢量.对固体域而言,耦合边界作为结构的纽曼边界条件,式(5)可写成分量的形式:式中,fx、fy、fz为固体在耦合边界上任一点沿x、y、z方向的应力分量;Px、Py、Pz为流体在耦合边界上任一点沿x、y、z方向的压强分量,其值由CFD计算给出.1.3耦合界面上固体、流体的变化耦合界面的能量守恒原理是指在耦合作用过程中,耦合界面上流体荷载(外力)、固体力(内力)在界面位移上所做的虚功相等:δW=δuΤs·fs=δuΤf·ff(7)式中,δus、δuf分别为耦合界面上固体、流体虚位移;fs、ff分别为耦合界面上固体、流体表面力.耦合界面上流体、固体的虚位移之间关系可以用下式表示:δuf=Hδus(8)H为由不同耦合方法得到的传递矩阵.将式(8)代入式(7)可得:fs=HTff(9)根据式(7)～式(9),可知传递矩阵H的求解致关重要.文献证明当转换矩阵H中行元素之和等于1时,可从能量守恒出发,导出界面力的守恒.1.4能量守恒原理效率的高低可以根据求解精度和求解时间的比值进行衡量.一般情况下,同时满足上述四个条件的理想的数据传输比较困难,而能量守恒原理作为自然界的一般法则,近些年来得到广泛重视和应用.2网格密度函数问题如图2所示,流体域和固体域在离散后,根据网格节点间的对应关系可分为一致网格和非匹配网格.通常情况下CFD计算要求的网格密度比CSD要密得多,由此产生耦合界面上两套非匹配网格之间的数据传递问题,在数学上这是一个双向插值问题,对该问题的研究也成为流固耦合的热点和难点之一.3局部插值法耦合界面信息传递方法从总体上可以分为局部插值法和整体插值法.局部插值法指在进行信息传递时,流体或固体界面网格点的未知物理量可根据固体或流体界面上的部分单元或网格点的已知物理量求得,而整体插值法需要根据所有单元或网格点上的已知物理量求得.以下介绍局部插值法中的不同方法.3.1固体主单元节点数计算映射点插值法是最直接、最简单的传递方法.其思想是:流固耦合边界上流体节点的位移和其对应的固体主单元中的映射点位移相等,而映射点位移由固体主单元中根据形函数插值得到,对力的传递则是一个相反的过程.关于主单元及映射点的概念和算法在3.4节给出.如图2b所示,流体边界单元中的节点A,其在结构中的主单元为abc,映射点为A′,根据映射点插值法:dfA=dfA′=ns∑i=1Νsidsi=Νsa(xA′)dsa+Νsb(xA′)dsa+Νsc(xA′)dsc(10)式中,ns为固体主单元节点数;dfA、dfA′分别为流体节点和映射点的位移矢量.写成分量的形式:dfAx=dfA′x=ns∑i=1Νsixdsix=Νsa(xA′)dsax+Νsb(xA′)dsax+Νsc(xA′)dscx(11a)dfAy=dsA′y=ns∑i=1Νsiydsiy=Νsa(xA′)dsay+Νsb(xA′)dsay+Νsc(xA′)dscy(11b)dfAz=dsA′z=ns∑i=1Νsizdsiz=Νsa(xA′)dsaz+Νsb(xA′)dsaz+Νsc(xA′)dscz(11c)采用上述方法,流体界面所有节点位移可由固体界面节点位移表示出来,写成矩阵的形式:式中,Hnf×ns为位移转换矩阵,根据映射点插值法组装得到,ns、nf分别为固体、流体耦合边界上的节点数.映射点插值法容易理解和应用,但缺点是没有考虑界面整体位移协调和力的平衡,使其应用范围受到很大的局限.3.2界面位移nisStein等根据任一点位移相等的概念,采用加权余量法进行信息传递.界面上位移守恒,其连续形式:ds(x)=df(x)(13)ds、df可近似表示为:ds=ns∑i=1Νis(x)dsi,df=nf∑j=1Νif(x)dfi(14)式中,ds、df分别为界面上结构、流体网格点上的位移矢量,Nis、Nif为结构、流体的基函数.根据迦辽金方法,对式(13)同乘以权函数ϕ(x),得到:∫Γϕ(x)ns∑i=1Νis(x)usidx=∫Γϕ(x)nf∑j=1Νif(x)dfidx(15)对于ϕ(x),可有两种选择,即可以为流体或结构的基函数.ns∑i=1[∫ΓΝkfΝisdx]dsi=nf∑j=1[∫ΓΝkfΝifdx]dfifork=1,...nf(16)可写成矩阵的形式:Afsus=Affdf(17)式中,Afs为nf×ns的矩阵,Aff为nf×nf的矩阵.由式(17),流体界面位移可由界面结构上已知位移求得:df=A-1ffAfsds(18)因此位移转换矩阵写为:H=A-1ffAfs(19)采用加权余量法时,会碰到形如∫ΓsΝfΝsdS的积分问题,即积分式里有关于流体和固体单元插值函数,但积分域为固体界面单元或流体界面单元.至今,这一混合积分未得到有效处理.3.3变密度插值Goura等人提出了一种常体积转换方法(CVT),该方法是一种与结构模态无关的局部插值方法.基本思想为:气动点的正交投影在结构三角形内保持不变,与固体表面单元的线性弹性特性一致.通过确保由结构点和气动点组成的四面体的体积守恒,气动点超出平板的分量就可以计算出来.徐敏、陈士橹在CVT基础上提出了一种改进的常体积转换法,即引入面积限制值来保证网格插值的质量.改进的CVT插值方法能避免原CVT方法可能出现的异常情况,且大大提高了插值精度.CVT法也需计算节点间的映射关系,且关于常体积的假定增加了数据传递的复杂性.3.4寻找主单元和映射点的算法应用局部插值法不可避免地涉及到映射点的搜寻问题,这需要单独的算法和计算程序.3.4.1单元内部映射流体或固体界面上的某一节点向固体或流体界面上单元所在平面作正投影,当映射点在单元内部时,称此单元为节点的主单元,其垂足为映射点,图3.3.4.2映射点位置的确定假设单元ijk为节点p(xp)的一个备选主单元,见图4:节点p(xp)在单元ijk中映射点p′的坐标xp′:xp′=2∑i=1[(xp-x0)⋅gi]gi(20)垂足:xn=2∑i=1[(xp-x0)⋅g3]g3(21)3.4.3单元特征特性考虑映射点p′(xp′),确定p′点是否在给定单元中的直接方法是根据p′点坐标值计算的其形函数值.根据有限元概念:xp′=∑iΝixi(22)形函数之和:∑iΝi=1(23)由式(22)、式(23),对于三维四面体单元:式(24)可简写为:xp=xN(25)因此,可以求得N:N=x-1xp(26)对于三维四面单元:由此,判断点xp′是否在单元中的条件为:min(Ni,1-Ni)≥0∀i(28)在某些情况下,按式(28)得到的主单元可能不止一个,在此情况下,需要附加的判定标准为:dn=|xn|≤δ(29)式中,dn为法向距离;xn为依式(21)计算所得的垂足;δ为指定容差.3.4.4求解算法2:完全搜寻法根据上述原理,可采用完全搜寻法(bruteforcesearch)进行主单元及映射点的搜寻.完全搜寻法首先假定映射点p′及其坐标xp′,然后对边界上所有的单元作一次循环计算,计算与xp′相关的形函数.其计算步骤可描述为:(1)根据式(20)、(21)计算映射点和垂足的坐标.(2)根据式(27)计算Ni,并按式(28)进行判断,如符合条件的主单元不止一个,进行步骤(3).(3)根据式(29)作进一步判定.回到步骤1,直到完成对所有边界单元节点扫描和计算.完全搜寻法算法,需要对所有单元进行循环,一般情况下要一半单元需要参与循环.为提高计算效率,文献还给出一些改进方法:Octreesearch,NeighborttoNeighborsearch,Advancingsearch方法等.这些方法都是基于不同的假定,采用不同数据结构类型,因此计算量比完全搜寻法算法有不同程度降低.如Octreesearch方法对计算域就有不同的数据结构类型,图5.需要指出的是Octreesearch方法通过分解域搜寻法可以提高计算效率,但由于单元组织相对复杂,可能在计算过程中失效.其他改进的方法亦可能出现类似的问题,此时需重新启用完全搜寻法.4一般插值方法4.1有限表面插值全域插值法最早来源于航空领域中对平面机翼变形的计算.Harder(1971)提出了无限平板样条法(IPS),该方法基于一无限板微分平衡方程的叠加求解.Yu对IPS作了进一步研究,并编制了相应的计算程序.IPS是程序中非常普遍采用的一种插值方法,如MSC/NASTRAN.Appa(1989)认为采用IPS法从结构内部网格点外插得到平板边网格点的值不可靠,并提出有限表面插值(FPS).FPS引入了一系列的约束条件,使得变形后的板经过给定数据点.由此,建立了结构和气动网格点位移映射,得到传递矩阵.根据传递矩阵可以对流固界面间变形、压强、温度等插值计算.此外,还有InverseIsoparametricMapping(IIM),方法以有限元中的等参元为基础.但此方法只适用二维问题,且无法进行外推.徐敏认为上述样条插值仅适合于薄板处于最小弯曲能(平衡位置)所确定的位置并且应在满足流体表面和结构表面一致的条件下才能得到理想的结果.在IPS、FPS之外,出现了旨对不规则曲面的插值计算的样条函数:薄板样条(TPS),Multiquadric-Biharmonic(MQ)等.Smith对上述方法作了技术上的全面阐述,并指出了每种方法适用范围.4.2距离平均加权和距离平均加权Shepard方法首先由气象学家与地质工作者提出,后来由于Shepard的工作称为Shepard方法,其基本思想是将插值函数定义为各数据点函数值的距离倒数加权和:s(x)=n∑j=0fjlj(x)=n∑j=0fj∏k≠jrk(x)∑nm=0∏k≠mrk(x)(30)式中,fj为x点处的函数值;lj(x)为权函数;rk(x)=|x-xj|22.Shepard插值虽然是光滑的,但缺点是在数据点附近表现为平点(其各阶导数为零),这样限制了Shepard插值的精度.如果已经知道了函数在数据点的高阶导数,可以通过构造改进的Shepard插值或带导数条件的Shepard插值来解决此问题.4.3rbf插值函数近些年来,数学界对径向基函数(radialbasisfunction,RBF)为基础的近似方法做了大量的研究工作.其中,德国Justus-Liebig大学的Martin.D.Buhmann教授和国内复旦大学数学系的吴宗敏教授对RBF函数的基础理论和应用研究起了重要推动作用.径向基函数插值方法可用于大型散乱数据的处理而被广泛应用于地形学、医学、计算机图形学等领域.Frank曾做了大量的各种散乱数据插值方法的实例比较,得到的结论是径向基函数插值的结果最能使人满意.在d维欧几里德空间里给定一组位置不同的点X={x1[WW(]...[WW)]xΝ}⊆Rd,称之为“中心点”.进一步,如果知道在这些中心点上的标量值g1,g2,...,gN,可以由此确定一个连续函数,使函数通过这些中心点.当连续函数采用式(1)表示的形式时,称之为RBF插值函数:s(x)=i=Ν∑i=1αiϕ(|x-xi|)(31)式中,s(x)为x点处的未知函数值;x为未知点坐标x(x,y,z);xi为第i已知数据点坐标xi(x,y,z);ϕ为选用的径向基函数;αi为相应于第i个数据点的待求系数;函数|x-xi|为欧几里得距离,对于三维空间可用r直接表示为:r=|x-xi|=√(x-xi)2+(y-yi)2+(z-zi)2(32)根据RBF插值算法,流体表面x、y、z方向的节点位移dfx、,dfy、dfz也由固体表面x、y、z方向的节点位移dsx、dsy、dsz左乘以传递矩阵H求得.得到传递矩阵H后,根据能量守恒原理即可编制相应的计算程序,对流固耦合界面位移、速度、压力等物理量的进行传递.值得注意的是对自变量空间,可以找到一个截断多项式,使得由其产生的径向函数在给定的维数空间正定且满足给定的连续性条件,即采用紧支径向基函数(compactlysupportedradialbasisfunction,简称CSRBF).这类函数适用于处理大规模计算,因为由此导出的线性方程组是一个稀疏矩阵.RBF插值方法允许CSD、CFD采用任意网格形式,因此在CSD和CFD计算程序之间可以开发独立的界面信息传递的接口程序,很容易实现弱耦合分析,