TJU计算流体力学 第3章

70第3章计算流体力学的有限元方法与差分法比较，有限元法具有适应复杂网格形状和复杂边界条件的特点。以下介绍有限元法的基本思想。3．1试函数方法设方程为)()(uGuF（3-1）其中，F、G均为微分算子。试解为iMiiMCu1（3-2）其中，),,(zyxii为线性无关的函数，称为基函数；iC为待定系数；在选定i后，要解决的关键问题是确定iC的方法。确定iC有两类方法可供选择：1）余量法；2）变分法。对试解的要求有3种选择：①满足边界条件，不满足方程，如古典的RITZ和GALERKIN方法；②不满足边界条件，满足方程，如有限基本解方法；③边界条件和方程均不满足，如有限差分和有限元方法。（此处满足指精确满足，计算流体力学中的离散方程是微分方程的近似方程，边界条件通常以后置修正法满足。）3．2加权余量法定解问题可以表示为fAu在中（在域内）gBu在上（在边界上）其中A，B为线性微分算子，为计算域，为计算域边界。设方程解为MiiiMCu1（3-3）使满足边界条件，与原方程解的误差为剩余量71fCAfAuRMiiiMe)(1（3-4）确定iC可以有以下几种方法。1）配置法在域上选取M个点，每个点表示所在区域的平均值。即iAiiudAAu1，Mi,,2,1分别代入方程（3-4）则有eMMMiMMMiiMMMeMiiiMRfzyxCAfAuRfzyxCAfAu111111111),,(),,(（3-5）使M个点上的余量eR等于零，形成M个联立方程，可以解出M个iC未知数。即MMMMMMMMfffCCCAAAAAAAAA2121212222111211（3-6）其中，第一个下标表示方程的序号；第二个下标表示基函数的序号。2）子区域法将域划分为M个子区域，eMeU1，kj，kj，即子区域覆盖整个区域，子区域无相交或交集为空集。将所有基函数在每一子区域上积分，（3-4）式积分后组成M个联立方程，并使每一子区域上余量eR积分为零。0)(eedRdfAueM，Me,2,1（3-7）即72MMMMfdfdfdCCCdAdAdAdAdAdAdAdAdAMMMM2122211121212121（3-8）形成M个联立方程，可以解出M个iC未知数。3）最小二乘法以子区域划分计算域，在域上剩余量的平方、加权、积分取极小值。)0min(2dwRe（3-9）其中，w为权重函数。上式极小化条件为（一阶导数为零）0])([21dfCAwCMiiij（3-10）积分边界确定，可以交换积分微分次序0])()[(1dfCAAwMiiij，Mj,2,1（3-11）即fdwAfdwAfdwACCCdAwAdAwAdAwAdAwAdAwAdAwAdAwAdAwAdAwAMMMMMMMM2121212221212111（3-12）734）矩量法以子区域划分计算域，选择完备函数序列中的一组线性无关函数12,,,1nxxx作为基函数i，权函数1w时，对剩余量fCARiiie)(1取不同阶次的静矩。对稳定平衡系统各阶合力矩为零，有0dRei，Mi,2,1（3-13）即fdxxfdfdCCCdAxxdAxxAdxdxAxxAxdxAddAxAxdAdnMnnnnnn1211121111（3-14）5）Galerkin法以子区域划分计算域，选基函数j与剩余量fCARiiie)(1正交，即j与eR内积为零，权函数1w时有0),(dRRejej，Mj,2,1（3-15）将fCARiiie)(1代入后展开为74fdfdfdCCCdAdAdAdAdAdAdAdAdAMMMMMMMM2121212221212111（3-16）除配置法外，以上几种方法均采用子区域划分计算域，子区域法是使各子域剩余量为零，最小二乘法是使剩余量的平方取极小值，矩量法是将j具体化的力、力矩平衡系统，Galerkin法将基函数j与剩余量表示为内积空间，内积中的函数正交，则内积为零0),(ejR，其物理意义是力在可能位移平面上的投影为零，即力所做的功为零。3．3变分法在理论研究中，有一些物理定律是通过“对某种积分取极值”的简明的数学原则推求出来的。当被积函数中含有未知函数时，有关确定积分极值的问题属于变分计算。已知)(xyy，在),(00yxA，),(11yxB两点间的积分值10),,(xxdxyyxFI（3-17）为最小，这是求)(xy的变分问题。这类问题有最速降线、曲面短程线、曲线等周最大面积、过两点最小旋转面、费马原理（光路最短时间）、悬索形状等问题。（见附录）设)(xy具有二阶导数，且)(0x=0，)(1x=0，α是微小参数（α相当于扰动量），则)(xy的近似函数为)(xy，)(x是)(xy趋近)(xy的竞争函数（优选函数），得)()()(xxyxy（3-18）)()()(xxyxy（3-19）75AB0x1x)(x)()()(xxyxy)(xyxy图3-1变分函数示意图表示是通过A、B，随α变化的一族曲线，使)(I有极值的必要条件是0)(II，是变分算子。10),,()(xxdxyyxFI（3-20）其中，y和y由（3-18）和（3-19）确定，与有关，则10][)(xxdxyyFyyFI（3-21）将（3-18）和（3-19）代入上式得10)]()([)(xxdxxyFxyFI（3-22）对积分号下的第二项做分部积分，并代入)(0x=0，)(1x=0条件得dxyFdxdxdxyFdxdxxyFdxxyFxxxxxxxx)()()()()()(10101010（3-23）76将（3-23）代入（3-22）整理后0)]()[(10xxdxyFdxdyFx（3-24）由于)(x和任意，则有0)(yFdxdyF（3-25）式（3-24）称为欧拉―拉格朗日方程。是函数极小化中的重要公式。对二维问题0),(),(yxfyxAu2),(Ryx0),(uyxu),(yx其中：为域内，为边界。求未知函数u的变分函数为dxdyyuxuuyxFI),,,,(（3-26）未知函数u的曲面簇如图3-1a所示，),(00yx),(10yx),(01yx),(11yxxyu),(yxu),(yx),(),(),(yxyxuyxu0),(cyx图3-1a二维变分函数示意图77边界条件表达为0),(00yx，0),(10yx，0),(01yx，0),(11yx；或0),(cyx，c表示边界曲线。),(),(),(yxyxuyxu在控制点上严格满足边界条件，且xxuxu，yyuyu变分函数写为dxdyyyuxxuuyxFI),,,,()(极值存在条件为0)(II，即0),,,,()(dxdyyyuxxuuyxFI其中),,,,(yyuxxuuyxF展开为)/()/()/()/(yuyuFxuxuFuuFFyyuFxxuFuF)/()/(由复合函数求导公式])/([)/(])/([xuFxxxuFxuFx])/([)/(])/([yuFyyyuFyuFy代入极值条件中78dxdyyyuFxxuFuFI])/()/([)(0]})/([])/([{}])/([])/([{dxdyyuFyxuFxdxdyyuFyxuFxuF由格林公式（GreenFormula）dsgfgdxfdydxdyygxfCCS)cossin()()(其中，在边界上逆时针为正（如图3-1b），dsdxcos，dsdysin，dndxsin，dndycos，ds为切线方向单位矢量，dn为法线方向单位矢量。oxySCdxdydsdn图3-1b边界微分变换示意图yxysyxsxssincosyxynyxnxncossin变分为790)cos)/(sin)/((}])/([])/([{)(dsyuFxuFdxdyyuFyxuFxuFIC注意到),(yx在边界上有0),(cyx，则上式为0}])/([])/([{dxdyyuFyxuFxuF泛函极值条件为0])/([])/([yuFyxuFxuF称为欧拉方程。设uyxfyuxuF),(])()[(2122其中，),(yxf是任意函数，代入欧拉方程得0)()(),(yuyxuxyxf或写为0),(2yxfu（3-27）说明未知函数u应满足泊松方程，这一原理称为变分原理，函数F称为泛函（数）。用0)],([2udxdyyxfuI表示最小能量条件，用dxdyuyxfyuxuI}),(])()[(21{22表示系统稳定条件。3．4边界条件二阶偏微分方程（如泊松方程）的定解问题可以表示为80上在上在上在上在内在432012000gnunuuuufu（3-28）1和2上的边界条件为狄里克雷即本质性边界条件。3和4上的边界条件为诺伊曼即自然边界条件。n为边界外法线方向。在数值解计算时，需要满足的边界条件是本质性边界条件，自然边界条件在方程离散过程中得到满足。3．5单元函数和形状函数3.5.1一维单元函数和形状函数对一维情况，将计算域分为结点1、2、3……n，n为节点总数。设解为njjjnxcu1)(（3-29）其中，jc为待定常数，)(xj