6-轴对称问题有限元分析

航空工程先进数值计算技术轴对称问题的有限元分析王晓军航空科学与工程学院固体力学研究所1/54轴对称问题的有限元分析轴对称问题是弹性空间问题的一个特殊问题，这类问题的特点是物体为某一平面绕其中心轴旋转而形成的回转体。由于一般形状的轴对称物体，用弹性力学的解析方法进行应力计算，很难得到精确解，因此采用有限元法进行应力分析，在工程上十分需要，同时用有限元计算得到的数值解，近似程度也比较高。空间轴对称问题，一般来说是三维问题，但是由于对称性，轴对称平面中的两个位移分量可以确定物体的应变和应力状态。采用柱坐标系，用z和r分别表示一点处的轴向和径向坐标，w和u为相应位移。2/54单元位移函数单元在子午面分割，绕对称轴z旋转后形成环单元（见图5-1），而在子午面形成的网格形同平面问题一样。单元界面可以是三角形，也可以是矩形、任意四边形，我们选三角形单元，如图5-2所示。图5-1轴对称物体图5-2三角形环单元3/54单元位移函数由于对称性，在对称面上，单元结点位移为TTeTTTijmiijjmmuwuwuw(5.1)仿平面问题，选线性位移模式123456,wurzrz(5.2)同平面问题处理方法一样，有,,,,iijjmmiiijmiijjmmiiijmuNuNuNuNuwNwNwNwNw(5.3)4/54单元位移函数其中()(,,)2iiiiabrczNijm(5.4)11121iijjmmrzrzrz(5.5)1,11,(,,)1jjjijmmjijmmmmjimjmrzzarzrzbzzrzzrcrrijmr(5.6)5/54单元位移函数表达成矩阵形式000000iijmejijmmeijmNNNuNNNNwNININI(5.7)其中1001I6/54单元应变与应力单元应变根据弹性力学对轴对称体进行应力分析。在弹性力学中，空间轴对称的几何方程为rzrzururwzwurz(5.8)7/54单元应变与应力与平面问题相比，多了一项。这项应变的产生主要是由于径向位移u引起的，因为径向方向有了位移u以后，原来的周长2r发生了改变，因而产生了环向应变，这样也就产生了环向应力。所以，轴对称问题的应力分量就不是3个，而是4个。将式(5.3)、式(5.4)代入式(5.8)求解后得到00000010002iijmiijmjijmjiijjmmmmubbbwfffucccwcbcbcbuw(5.9)8/54单元应变与应力(,,)iiiiaczfbijmrr(5.10)式(5.9)可改写为eeijmBBBB(5.11)其中00102iiiiiibfBccb(5.12)9/54单元应变与应力由此可见，轴对称问题的几何方程式(5.11)，在形式上和平面问题是一样的，但是轴对称问题中的B和并不完全是常量元素，其中各点的应变将随r、z的函数，故B是r、z的函数。由于B是r、z的函数，所以单元中各点的应变将随r、z而变化，即单元中各点的应变不同。为了简化计算，通常用单元形心坐标(,)zr近似代替if中的r、z值，即用单元形心(,)zr处的应变作为单元的平均应变，变成常应变单元，即10/54单元应变与应力1()31()3ijmijmiiiiiizzzzzrrrrraczffbrr(5.13)11/54单元应变与应力单元应力根据物理方程推导的单元应力公式为10111011(1)(1)(12)1011120002(1)rrzzrzrzE(5.14)12/54单元应变与应力eeeijmDDBssss(5.15)111131222()iiiiiiiiiiiibAfAcAbfAcAsAbfcAcAb(5.16)其中12312(1),,12(1)4(1)(12)EAAA13/54单元应变与应力为了简化计算，消除对称轴上由于0r引起的麻烦，用z、r近似代替r、z。由式(5.12)和式(5.16)计算所得的应变、应力则是单元形心处的应变、应力的近似值。14/54单元刚度矩阵同平面问题一样，用虚位移原理推导单元刚度矩阵。在轴对称情况下，单元的虚位移方程为**()()TeTeFrdrddz(5.17)等式左边为单元等效结点力eF所做的虚功，等式右边为整个三角形环单元中应力的虚功。假设单元的虚位移为**eN15/54单元刚度矩阵由于对称性202d**eB(5.18)代入式(5.17)，上式可得**()()22eTeeTTeeTeFBDBrdrdzFBDBrdrdz(5.19)单元刚度矩阵2TkBDBrdrdz(5.20)16/54单元刚度矩阵写成分快形式iiijikjijjjkkikjkkkkkkkkkkkk(5.21)其中子矩阵2(,,;,,)TststkBDBrdrdzsijmtijm(5.22)取单元形心坐标z、r代替B中的r、z，则2TststkBDBr(5.23)17/54单元刚度矩阵于是有112123122()()()2()sttsttsttssstststtstststbbAfffAbAccAcbfAccrAkAcbfAbbccAbb(5.24)18/54整体刚度矩阵单元数en，结点数n，把单元的e、eF和k扩大到整个结构的自由度是数，叠加得到112eenneTeeFBDBrdrdzd(5.25)整体载荷阵列1eneeFF(5.26)19/54整体刚度矩阵整体刚度矩阵为112eenneTeeKkBDBrdrdz(5.27)故有KF(5.28)整体刚度矩阵K写成分快形式20/54整体刚度矩阵1111111111ijmniiiijiminjjijjjmjnmmimjmmmnnninjnmnnKKKKKKKKKKKKKKKKKKKKKKKKKK(5.29)21/54整体刚度矩阵其中子矩阵1(1,2,,;1,2,,)enststeKksntn22/54等效结点力现在讨论等效结点力的计算。KF的右边载荷列阵121enTeTTTneFFFFF(5.30)其中(1,2,,)TiriziFFFin等效结点力位移按照虚功原理有****()()2()()eTeiTTTcgsvFrffrddsfrdrddz(5.31)23/54等效结点力式中cr——集中载荷gf作用点的径向坐标sf——表面力vf——体积力单元中各点虚位移为**eN单元中各点的等效结点力为eeeegsvFFFF(5.32)24/54等效结点力其中，集中力的等效结点力为2eTggFrNf(5.33)表面力的等效结点力为2eTsrFNfrds(5.34)体积力的等效结点力为2eTsvFNfrdrdz(5.35)25/54等效结点力对整体而言，载荷为1()eneeeeeeegsvgsveFFFFFFF(5.36)由于积分号后的被积函数有变量r，无法采用采用静力等效原理。26/54等效结点力体积力对三角形单元来说000000ivrivzTijmvrjvreTvvijmvzjvzmvrmvzNfNfNNNfNfFNfNNNfNfNfNf(5.37)27/54等效结点力(1)自重0vrf，vzf，为重度，于是单元的自重移置到结点,,ijm等效结点力为02eevriviievziFFNrdrdzF(,,)ijm(5.38)记iijjmmrrNrNrN，有()iiiijjmmNrdrdzNrNrNrNdrdz(5.39)28/54等效结点力由积分公式!!!2(2)!ijmNNNdxdy得()(3)(,,)6121212jimiirrrNrdrdzrrijm0==(,,)-(3)6eevrivieivziFFijmrrF(5.40)29/54等效结点力(2)离心力此时2,0rzPrP(5.41)式中——角速度于是单元的离心力移置到结点,,ijm的等效结点力为2==2(,,)0eevriviievziFrFNrdrdzijmF(5.42)30/54等效结点力对其积分得22()iiiijjmmNrdrdzNrNrNrNdrdz利用积分公式2222(6)30iijmijmNrdrdzrrrrrrr于是有222(92)==(,,)150eevriijmvievziFrrrrFijmF(5.43)31/54等效结点力表面力设rOz平面上单元ij边上受有线性分布的径向表面力，如图5-3所示。图5-3受表面力的三角形单元32/54等效结点力在结点i的集度为iq，在结点j的集度为jq，ij边长为l。在此情况下有,0riijjzqqNqNq于是结点i集中等效结点力为==20eiijjesrisiiesziqNqNFFNrdsF(5.44)33/54等效结点力在ij边上形函数0mN，利用积分公式!!(,,)(1)!ijlNNdslijm，有2()(3)12()()12iiiiijjijijijiijjijlNNrdsNrNrNdsrrlNNrdsNNrNrNdsrr11(3)()(3)()21212060iijjijeiijjijsiqrrqrrqrrqrrlF(5.45)34/54等效结点力类似有()(3)06iijjijesjqrrqrrlF(5.46)00esmF(5.47)35/54精确刚度矩阵的计算精确刚度矩阵的计算单元刚度矩阵中，有9个22子矩阵2TststkBDBrdrdz(,,;,,)sijmtijm(5.48)由于sB、tB与坐标r、z有关，积分复杂。前面用单元形心坐标z、r代替r、z得到了近似刚度矩阵。下面求解精确刚度矩阵记tsssBBB(5.49)36/54精确刚度矩阵的计算其中，sB为用z、r代入sB得到，则有00120ssssssssbaczbBrrccb(5.50)00120ssssssssbaczbBrrccb(5.51)37/54精确刚度矩阵的计算00012000000101()00200tsssssssssssBBBaczaczrraczaczrr(5.52)38/54精确刚度矩阵的计算把式(5.49)代入式(5.48)得''''''2()()2