求解定常势流正命题的全时-空变分有限元法

第25卷第2期刊Vol 25 No. 22004年6月CHINESE QUARTERLY OF MECHANICSJune 2004求解定常势流正命题的全时空变分有限元法陈池,陶毅,刘高联(上海大学力学所,上海200072)摘要:非定常流动变分原理的建立使得用有限元法来求解多工况点的设计问题成为可能。本文在刘高联的非定常变分理论的基础上,对定常变分问题进行时间相关有限元求解。但由于可压缩非定常位势流动的控制方程是双曲型的简单地把时间当作同空间一样的物理维来求解是不可行的。而现有的时空有限元法极其复杂增加了计算复杂度,使其很难用于工程设计中。为此,文[2、3]提出了求解一维非定常问题的新型时空有限元法。本文把该方法推广到二维流动用它求解二维弯管内的流动和翼型绕流问题。计算结果与用定常方法求得的结果几乎重合,说明该方法可以用于多维时间相关求解关键词:非定常;变分;时空有限元;计算流体力学中图分类号:0355文献标识码:A文章编号:02540053(2004)021576A New Time/Space FEM for Steady Potential FlowCHEN Chi, TAO Yi, LIU Gao-LianShanghai University, Shanghai Inst. App. Math. Mech, Shanghai 200072, China)Abstract: On the basis of Liu Gaolian's unsteady variational theories, the time-marching method is used tosolve steady variational problems. Because the governing equations for compressible unsteady potentialflow is hyperbolic, looking time dimension as space dimension in the same way is never appropriate. Unfortunately, the existing time/ space FEM is quite complicated, and difficult to be put into engineeringuse. A new time/space FEM was presented to calculate the one-dimensional unsteady problems in papers[2,3]. The method is now extended to calculate two-dimensional steady flow in a time-marching way. Asexamples, one pipe flow and one airfoil rounding flow are considered the results agree well with those ofthe steady flows, which demonstrates the usefulness of this method in solving multidimensional problemsKey words: unsteady variational calculus; time/ space FEM; CFD非定常问题是当今研究的热门问题之一,但到目前为止,多数非定常问题都是用有限差分或有限体积法求解的。有限元法以其能适应复杂的几何区域而逐渐在计算流体力学领域中得到应用。但是,用有限元方法求解非定常问题的文献也主要是基于 Galerkin变分理论的基础上进行的),对时间偏导数的处理大多数是采用的差分方法。这一方面是因为非定常的变分原理的建立是极其困难的,因为非定常问题的控制方程是双曲线型的,在时间方向只能给初始条件,而不能给终止条件。刘高联根据 Hamilton原理成功地建立了非定常问题的变分原理,为非定常问题的变分求解建立了完整周密的理论基础。另一方面是因为差分法相对于时空有限元法要容易。但我们注意到,有很多与时间相关的问题是很难用差分法解决的,例如降落伞的降落问题但如果用时空有限元法就容易处理°。变分法以其有强大的变域变分工具,可以为叶轮机进行多工况点反设计而逐渐得到重视。遗憾的是,迄今为止,还没有发现非定常变分计收稿日期:2003409-23基金项目:国家自然科学基金项目(50136030,10372055)作者简介:陈池(1973-),男,四川广安人,博土力学季刊第25卷算的有效方法。文献[7对二维振荡翼型进行了求解,但是文中对时间维的处理方式完全和空间维的处理方式等同。我们知道,非定常绝热位势流动的控制方程是双曲型方程,在时间方向绝对不能采取同空间方向一样的处理。为此,我们必须采用时空有限元法,但是现有的时空有限元法使用很不方便,公式也很复杂,而且也不适合变分计算。为此,我们需要提出新的求解方法。受双曲型方程特征的启发,刘高联提出在处理单元内的插值函数时对时间方向上采用展开处理方式的设想。经过一维的计算试验,证明是有效的,本文将用它来对二维的绝热势流问题进行时间相关求解。1二维非定常势流的控制方程及其变分理论维非定常势流其相应的无量纲(以滞止参数无量纲化)气动方程为:8t+·()=0(1)中=A(y-1)(atP(3)上述公式中,、A、、P、y分别为无量纲密度速度、位势压力和绝热指数。有量纲量的定义见文献[7]计算区域如图1所示初边值条件为(1)在初始边给定势函数∮和密度值P作为初始条件,稳定解与初始值关系不大。即φ=f2P=f3。(2)远场边界:应用 Riesman不变量来计算。(3)物型面上:元=0(4)激波面上:以g表示其法向分速,则有 Rankine-Hugniot激波关系。+P(A -g)[lA=0|p(An-g)=0图1计算块[|A,]=0Fig 1 A slab[|H门]=(y-1)g[|An](5)尾流线上下:an=0(6)尾缘点: Kutta无载荷条件。用变分推导的系统性途径可以得到以下的变分原理:上述非定常问题的解将使a=0,其中(1-(2∵()+2(1ydadydt+‖(pAn),dtp阳dady(φ-f2)-2∫3∮dxc第2期陈池,等:求解定常势流正命题的全时空变分有限元法2控制泛函的离散对势函数在有限单元内作如下的离散:1-[中k(,n,5=22+a1(1+51+§1-m+a2(1+5]+251274+a2(1+1-1[+a4(1+5](6)引入二维空间单元的形函数替代上式中含、?的各项,即有(,,3)=N[+a:(1+5](遵守求和约定)其中N=-41=2,M=9+41=2,N2=9+9+2,M,=99+2定义从空间-时间坐标(x,y,t)到有限单元体积坐标(、7之间的 Jacobian转换矩阵如下atdx ay atasas as从而有单元内势函数对空间-时间的偏导数8冲冲其中J1为 Jacobian矩阵的逆矩阵,J1、J2、J32分别为由 Jacobian矩阵的逆矩阵的第一、二、三列元素组成的行向量。v绰为势函数对体积坐标的散度,分别为a285[+a,(1+5]a7a71[+a,(1+5]a=Na,(遵守求和约定),1=12,4(下面对六k同)。由此我们可得到单元内泛函对任意a;的偏导数为e(1+5e[+a(1+5](1+5)JAAn2[+a1(1+5)/、、、1+5+J;18N,(1+)[+a/(1+5(1+J 2anD[+a1(1+]·J2aNo(1+5)NN力学季刊第25卷+‖(aA,)N(1+5)]·E3dd+‖pN(1+5·E1dln(9)上列诸式中J|为 Jacobian矩阵的行列式,E12、E13为面积分的转换系数为了把(9)式中含a的项线性化处理我们把和中的a提出来,从而可以得到单元内写成矩阵形式的代数方程组:b2b2 ba b2|a2p2(10)b3 b32 b3 b3abb ba b其中的参数为(1+)aN(1+y)(1+ya(1+5aN(1+5+J,iaN(1+)aN(1+)PJ3 aN(1+5)+j,aNJ,1(1+5NN(1+)(1+JIdedrds0N+‖(a,)N(1+5]·E2d5d3+‖pN(1+t·E1dηp式中对k有求和约定合并各单元的上述方程可以得到整个网格上的关于a的方程。初始密度假设为来流的密度,初始势函数为来流速度与当地横坐标的乘积,即:=A2。由于我们的计算方法是隐式的,可以取较大的时间步长。另外,为了消除上述线性方程组的奇异性,我们必须给定某些点上的a值。由此我们以特征线(或面)公式得到周围边界r上的速度和密度,从而可以算出当地的值。应用上述变分泛函求解,不需要其它边界条件,因为在泛函中我们已经考虑了所有的边界条件3数值结果3.1一维问题求解为了验证这种方法我们首先对一维非定常的流动进行了求解。我们的例子选自文[5],即要计算如图2中的大容器当阀门打开后管内的流动状况。用本文的方法所得的结果和文献[5]用特征线所求得的结果见图3和4二者几乎完全重合,说明我们的方法用于计算一维非定常流动是很成功的3.2二维非定常求解用上述方法我们对二维的可压缩绕流进行了非定常计算。计算的处理方式同一维问题求解的过程。第2期陈池,等:求解定常势流正命题的全时空变分有限元法161计算管没图2计算管段及相关参数Fig 2 Calculated pipe and related parameters10000100000图3压力比较(用符号衰示的为文[5]的结果图4速度比较(用符号示的为文[5]的结果Fig 3 Pressure comparison (results denotedFig 4 Velocity comparison(results denotedby character taken from paper [5])by character taken from paper [5])在这里我们做了两个算例。其中一个算例选自NACA的研究报告,文中是用流函数和势函数计算的定常流动。计算区域及网格如图5所示。进口无量纲速度是0.3993。我们的计算经过12个周期达到稳定。得到的结果( Time-Marching)和研究报告的结果( Stream-Potential)见图6和图7所示。可以看出它们是吻合得相当好的。势流函数法“鲁一喜一着一着一一看图5二维管道网格图6下管壁面速度Fig 6 Velocity on the lower wall of pipe另一个算例是计算的NACA0012翼型的绕流计算区域和网格见图8。计算经过40多个时间步(单步步长为10个无量纲时间)计算基本达到稳定(势函数绝对误差总和小于105),所得到的结果如图9所示。为了比较我们用定常的公式计算了同样的问题把它们的结果画在了一起。可以明显看出,本方法所得的结果和定常计算的结果吻合得相当好。实际的计算显示它们之间势函数差值的总和在10以内。力学季刊第25卷特流函数法图7上管壁面速度图8翼型网格Fig. 7 Velocity on the upper wall of pipeFig8 Airfoil grid势流函数法图9压力系数Fig 9 Pressure coefficient4结论针对非定常绝热流动的控制方程的双曲型特点,本文采用时空变分有限元来计算定常的变分问题。从所得到的结果可以看出,本方法用于计算二维时间相关的定常流动是切实可行的。通过上面的推导过程还可以看出,本方法可以推广应用于求解多维问题。参考文献:[11 Liu G L. A vital innovation in Hamilton principle and its extension to initial-value problems[C]. Proc 4th Intl Conf Nonlinear MechShanghai Univ Press, 2002, 90-97.[2]刘高联陶毅,李孝伟陈池新型全时空有限元法[C]水动力学学术会议论文集深圳,2003,[3]陶毅刘高联李孝伟陈池计算非定常流动的全时空有限元法[C]水动力学学术会议论文集深圳,2003[4 Chung T J. Finite element analysis in fluid dynamics[M]. McGraw-Hill Int Book Company, 1978[5 Hirsch C. Numerical computation of internal external flows, vol 2[M). John-wiley, New York, 1990[6] Hughes T J R, Hulbert G M. Space-time finite element methods for elastodynamics, Formulations and error estimates[J]. Comp Meth[7]郭加宏刘高联二维振荡机翼跨声速非定常绕流的变域变分有限元解空气动力学学报[1996,14(4);217-222[8]刘高联辆合热弹性动力学的统一的变分原理族[].力学学报,1999,31(2):165172