1、选择位移模式时,第选择位移模式时,第2章提到要考虑章提到要考虑解的收敛性解的收敛性,即要考虑,即要考虑到到位移模式的完备性和协调性位移模式的完备性和协调性。实际操作中,一般应考虑。实际操作中,一般应考虑位移模式的对称性。这是因为,有限元位移模式的选择实位移模式的对称性。这是因为,有限元位移模式的选择实际是以帕斯卡(际是以帕斯卡(Pascal)三解形基础上的(如图)三解形基础上的(如图7-2所示),所示),由低价至高阶,顺序选取,组成多项式。由低价至高阶,顺序选取,组成多项式。多项式中的项数多项式中的项数等于单元节点自由度数等于单元节点自由度数。如三节点三角形单元,位移。如三节点三角形单元,位移
2、模式取完全一次式,共模式取完全一次式,共3项。项。六节点三角形单元,位移模六节点三角形单元,位移模式取完全二次式共式取完全二次式共6项。如项。如果某一阶次不能全取,则果某一阶次不能全取,则应按对称性原则适当选取。应按对称性原则适当选取。 1 x y x2 xy y2 x3 x2y xy2 y3 x4 x3y x2y2 xy3 y4 图图7-2 多项式选择的多项式选择的 怕斯卡三角形怕斯卡三角形 1 x y x2 xy y2 x3 x2y xy2 y3 x4 x3y x2y2 xy3 y4 图图7-2 多项式选择的多项式选择的 怕斯卡三角形怕斯卡三角形例如在下节将要讨论的四例如在下节将要讨论的四
3、结点矩形单元中,位移模结点矩形单元中,位移模式不能取式不能取1,x,y,x2四项,四项,也不能取也不能取1,x,y, y2四项,四项,而应取而应取1,x,y,xy四项。四项。 7.2 四节点矩形单元四节点矩形单元 图图7-3示出的矩形单元,边长分别为示出的矩形单元,边长分别为2a和和2b。取。取4个角点为个角点为节点,编号为节点,编号为i,j,l,m。将。将x轴和轴和y轴置于单元的对称轴轴置于单元的对称轴上。上。单元的位移函数可取为:单元的位移函数可取为:1、位移函数、位移函数 在上式表示的位移模式中,在上式表示的位移模式中,a1, a2, a3, a5, a6, a7, a8反映了单元的刚体
4、位移和常反映了单元的刚体位移和常应变。在单元的边界应变。在单元的边界(x=a或或y =a)上(或),位移是上(或),位移是按线性分布的。因此,相邻按线性分布的。因此,相邻单元在公共边上的位移是连单元在公共边上的位移是连续的。这样,位移模式满足续的。这样,位移模式满足了解答收敛性的充分条件。了解答收敛性的充分条件。 ijlmxyaabb图图7-3在式(在式(7-1)中代入节点位移和节点坐标后,可解出)中代入节点位移和节点坐标后,可解出xyayaxaavxyayaxaau87654321(7-1)式中形函数为式中形函数为 :byaxNbyaxNbyaxNbyaxNmlji1141114111411
5、141(7-3)mmlljjiimmlljjiiuNuNuNuNvuNuNuNuNu(7-2)各待定系数各待定系数(a1 a8)。将这些系数再代入式(。将这些系数再代入式(7-1),可得:),可得:则式(则式(7-3)可简写为)可简写为),(4/ )1)(1 (mljiNiii(7-4)将位移函数写成矩阵形式,即有与式(将位移函数写成矩阵形式,即有与式(2-20)相同的形式)相同的形式 eNvuf(7-5)式中式中Tmmlljjiievuvuvuvu(7-6)令令 byax,在节点上的值为:在节点上的值为: ),(,mljibyaxiiiiININININNnlji(7-7)其中,其中,I为二
6、阶单位矩阵。为二阶单位矩阵。2、应变矩阵、应变矩阵根据几何方程,可得与式(根据几何方程,可得与式(2-25)同样的形式)同样的形式eB(7-8)把应变矩阵把应变矩阵B写成子矩阵形式写成子矩阵形式 mljiBBBBB (7-9) 其其 中中),()1 ()1 ()1 (00)1 (41mljibaababBiiiiiiiii(7-10) 由此可见,由此可见,B是是 、 的函数,即是的函数,即是x、y的函数。因的函数。因此此单元中的应变不再是常数单元中的应变不再是常数。 3、应力矩阵、应力矩阵根据应力根据应力-应变关系,可以计算单元中的应力,得到式(应变关系,可以计算单元中的应力,得到式(2-28
7、)同样形式)同样形式 eSD(7-11)应力矩阵应力矩阵S具有与式(具有与式(2-29)同样形式)同样形式BDS (7-12)将将S写成子矩阵形式写成子矩阵形式 (7-13)mljiSSSSS 其其 中中 )1 (21)1 (21)1 ()1 ()1 ()1 ()1 (42iiiiiiiiiiiiibaabababES),(mlji(7-14)上式对应平面应力情形。对于平面应变情形,只需将其中上式对应平面应力情形。对于平面应变情形,只需将其中的的E, 作相应的改变即可。作相应的改变即可。 4、单元刚度矩阵、单元刚度矩阵 单元刚度矩阵可采用式(单元刚度矩阵可采用式(2-33a)进行计算)进行计算
8、AThdxdyBDBk(2-33a)在四节点矩形单元中,在四节点矩形单元中,k是一个是一个88的矩阵。将的矩阵。将k写成分写成分块形式:块形式: 88mmmlmjmilmllljlijmjljjjiimilijiikkkkkkkkkkkkkkkkk(7-16)其中的子矩阵其中的子矩阵krs22可由下式计算可由下式计算 )1 (4222Ehkrssrsrsrsrsrsrsrsrsrsrsrsrbaabbaab31121311212131121311),(mljisr(7-17) 上式对应平面应力情形。对于平面应变情形,只须将上式上式对应平面应力情形。对于平面应变情形,只须将上式中的中的E、 作相
9、应的改变。作相应的改变。 5、等价节点力、等价节点力 单元体积力和表面力引起的节点力仍可用式(单元体积力和表面力引起的节点力仍可用式(2-45)和)和(2-46)进行计算。)进行计算。 对本问题给定的位移函数,若体积力是重力的情形(设对本问题给定的位移函数,若体积力是重力的情形(设重度为重度为 ),单元等价节点载荷列阵为:),单元等价节点载荷列阵为:AVTVhdxdyqNF(2-45)hdlqNFSTlS(2-46)TeVhabF10101010(7-18) 有了对单元的上述结果,便可应用第有了对单元的上述结果,便可应用第5章的方法组集结章的方法组集结构刚度矩阵和节点荷载向量;求解节点位移;计
10、算内力和构刚度矩阵和节点荷载向量;求解节点位移;计算内力和应力。应力。 四节点矩形单元采用较高阶的位移模式,具有比三节点四节点矩形单元采用较高阶的位移模式,具有比三节点三角形单元较高的计算精度。但矩形单元也有缺点三角形单元较高的计算精度。但矩形单元也有缺点, 在三角形单元在三角形单元i, j, m的各边中点增设一个节点,的各边中点增设一个节点,使每个单元具有使每个单元具有6个节点,个节点,得到图得到图7-4所示的六节点所示的六节点三角形单元。三角形单元。 这种单元具有这种单元具有12个个自由度,可以采用完全自由度,可以采用完全二次多项式的位移模式:二次多项式的位移模式: 一是不能适应斜线及曲线
11、边界,二是不便于采用大小不同一是不能适应斜线及曲线边界,二是不便于采用大小不同的单元的单元。7.3 六节点三角形单元六节点三角形单元 1、位移模式、位移模式 ijmi j m xy图图7-421211101098726524321yaxyaxayaxaavyaxyaxayaxaau(7-20)所取位移模式反映了单元的所取位移模式反映了单元的刚体位移和常应变刚体位移和常应变;单元内部;单元内部是是连续的连续的;在单元边界上位移分量按抛物线变化,而每条;在单元边界上位移分量按抛物线变化,而每条公共边界上有公共边界上有3个公共结点个公共结点,可以保证相邻两单元位移的,可以保证相邻两单元位移的连续性。
12、因此,上述位移模式满足连续性。因此,上述位移模式满足收敛的必要和充分条件收敛的必要和充分条件。 上述位移模式确定之后,可以用分析三节点三角形单元上述位移模式确定之后,可以用分析三节点三角形单元和四节点矩形单元相同的方法进行分析。得到形函数、应和四节点矩形单元相同的方法进行分析。得到形函数、应变矩阵、应力矩阵、单元刚度矩阵、等价节点力向量。但变矩阵、应力矩阵、单元刚度矩阵、等价节点力向量。但其过程十分繁复,采用其过程十分繁复,采用面积坐标面积坐标可以大大简化计算。可以大大简化计算。2、面积坐标、面积坐标 对于一个三角形对于一个三角形ijm(图图7-5),三角形内任一点),三角形内任一点P(x,y
13、)的)的位置,可以用如下的三个比值来确定:位置,可以用如下的三个比值来确定:ijmxy图图7-5PAALAALAALmmjjii(7-21)AiAjAm(1)定义)定义其中其中A为三角形为三角形ijm的面积,的面积,Ai, Aj, Am分别为三角形的分别为三角形的Pjm, Pmi, Pijd的面积。这三个比值的面积。这三个比值Li, Lj, Lm称为称为P点的面积坐标。点的面积坐标。 由于由于 AAAAmii则则 1mjiLLL(7-22)由此可见,由此可见,P点的三个面积坐标不是独立的。同时,面积点的三个面积坐标不是独立的。同时,面积坐标只是用以确定三角形内部某点的位置,因而是一种局坐标只是
14、用以确定三角形内部某点的位置,因而是一种局部坐标。下面进一步给出面积坐标的几个性质。部坐标。下面进一步给出面积坐标的几个性质。 (2)面积坐标与直角坐标的关系)面积坐标与直角坐标的关系在图在图7-5中,三角形中,三角形Pjm的面积为的面积为 )(2111121ycxbayxyxyxAiiimmjji(7-23)由式(由式(7-23),式(),式(7-21)化为)化为 ),()(21mjiycxbaALiiii(7-24)将式(将式(7-24)、()、(7-23a)和式()和式(2-18)、)、 (2-17)对比,)对比,可知,面积坐标就是三节点三角形单元的形函数可知,面积坐标就是三节点三角形单
15、元的形函数 jmmjiyxyxamjiyybmjixxc(7-23a)Ni、Nj、Nm。 将式(将式(7-24)的)的3个式子分别乘以个式子分别乘以xi, xj, xm,然后相加,然后相加,并利用关系式(并利用关系式(7-23a),有),有xLxLxLxmmjjiiyLyLyLymmjjii同理同理 (7-25)(3)面积坐标的导数公式)面积坐标的导数公式 根据面积坐标与直角坐标的关系,由复合函数的求导公根据面积坐标与直角坐标的关系,由复合函数的求导公式,有式,有mmjjiimmjjiimmjjiimmjjiiLAcLAcLAcLyLLyLLyLyLAbLAbLAbLxLLxLLxLx2222
16、22(7-26)(4)面积坐标的积分公式)面积坐标的积分公式下面给出面积坐标的幂函数积分公式。它们在计算单元下面给出面积坐标的幂函数积分公式。它们在计算单元刚度矩阵和等效结点载荷时有用。刚度矩阵和等效结点载荷时有用。 在三角形单元上进行积分时,有在三角形单元上进行积分时,有AcbacbadxdyLLLcmbjaiA2)!2(!(7-27)在三角形某一边(设在三角形某一边(设ij边,边长为边,边长为l)上进行积分时,有)上进行积分时,有 (7-28)lbabadlLLbjail)!1(!3、用面积坐标表示六节点三角形单元计算公式、用面积坐标表示六节点三角形单元计算公式 对应如图对应如图7-4所示
17、的六节点三角形单元,所示的六节点三角形单元,形函数可用面积形函数可用面积坐标表示为坐标表示为 (1)形函数和位移表达式)形函数和位移表达式ijmi j m xy图图7-6现利用形函数的性质检验式(现利用形函数的性质检验式(7-29)的正确性。先考虑三)的正确性。先考虑三角形的角点,例如图角形的角点,例如图7-6中的中的i点,有点,有 AAi0mjAA由式(由式(7-21)()(P16),有),有1iL0mjLL代入式(代入式(7-29),有),有 001mjimjiNNNNNN),;,(4:),)(12(:mjimjiLLNmjiLLNmjiiii边中结点角结点(7-29) 再考虑三角形的边中
18、点,例如再考虑三角形的边中点,例如i 点,面积划分如图点,面积划分如图7-7所示。显然有:所示。显然有: ijmi j m xy图图7-720AAAAmji由式(由式(7-21)(P16),有,有 210miiLLL代入式(代入式(7-29) (P16) ,001mjimjiNNNNNN进一步说明式(进一步说明式(7-29)所表示的)所表示的形函数的正确性形函数的正确性。 说明形函数说明形函数Ni在在i点等于点等于1,在其它节点等于,在其它节点等于0,因此是正确的,因此是正确的。形函数确定后,单元中任意一点的位移可以表示为:形函数确定后,单元中任意一点的位移可以表示为:eNf11212212(
19、7-30)其其 中中 Tmmjjiimmjjiievuvuvuvuvuvu112(7-31)122ININININININNmjimji(7-32)其中其中I为二阶单位阵,形函数由式(为二阶单位阵,形函数由式(7-29)确定。)确定。 (2)应变矩阵)应变矩阵单元中的应变仍可表示为:单元中的应变仍可表示为: eB11212313(7-33) 式中应变矩阵式中应变矩阵B为:为:(7-34)123mjimjiBBBBBBB其中其中),() 14() 14() 14(00) 14(2123mjiLbLcLcLbABiiiiiiiii(7-35)),;,()(4)(4)(400)(42123mjimj
20、ibLLbcLLccLLcbLLbABmjmjmjmjmjmjmjmji单元中的应力仍可表示为:单元中的应力仍可表示为: (3)应力矩阵)应力矩阵311212311212333133313SBDDe(7-36) 式中式中D是弹性矩阵,由式(是弹性矩阵,由式(2-9)确定;应变矩阵)确定;应变矩阵由式(由式(7-34)、()、(7-35)确定。根据矩阵乘法,可以给)确定。根据矩阵乘法,可以给出用面积坐标表示的应力矩阵出用面积坐标表示的应力矩阵S(4)单元刚度矩阵)单元刚度矩阵 单元刚度矩阵仍可表示为:单元刚度矩阵仍可表示为:hdxdyBDBkTA(7-37) 根据根据B、D的表达式以及面积坐标的
21、积分公式(的表达式以及面积坐标的积分公式(7-27),可以),可以求出求出k中元素的显式表示。由于较为繁复,这里就不列出详细结果。中元素的显式表示。由于较为繁复,这里就不列出详细结果。 (5)等价节点力向量)等价节点力向量 由于位移模式是非线性的,因此体积力和表面力引起由于位移模式是非线性的,因此体积力和表面力引起的节点力向量不能采用静力等效原理进行分配,而应采用的节点力向量不能采用静力等效原理进行分配,而应采用相应公式进行计算。相应公式进行计算。 单元体积力引起的等价节点力计算公式仍为:单元体积力引起的等价节点力计算公式仍为: AVTeVhdxdyPNF(7-38) 将由式(将由式(7-29
22、)、()、(7-32)表示的)表示的N代入,并应用积分式代入,并应用积分式(7-27),可以计算),可以计算 FV e。例如对于重力引起的。例如对于重力引起的 FV e ,有,有 TeVhAF 1010100000003它表示各边中点承担单元重力的它表示各边中点承担单元重力的1/3。 单元表面力引起的结点力计算公式仍为:单元表面力引起的结点力计算公式仍为: hdlPNFSTleS(7-39)设在设在ij边上受有边上受有x方向的均匀分布力方向的均匀分布力ps,对应的等价节点力,对应的等价节点力向量为(图向量为(图7-8) pslh/6pslh/64pslh/6ijmi j m xy图图7-9lp
23、sijmi j m xy图图7-8lpspslh/62pslh/6TseSlhPF040000000101 61TseSlhPF0320000000003121如在如在ij边上受到边上受到x方向的三角形分布面力,其集度在方向的三角形分布面力,其集度在i点为点为ps,在在j点为点为0。对应的节点力向量为(图。对应的节点力向量为(图3-9) 它表示边中点承担载荷的它表示边中点承担载荷的2/3,载荷集度大的角节点承担,载荷集度大的角节点承担1/3。 六结点三角形单元中的应变、应力不为常量,因此可以应用于六结点三角形单元中的应变、应力不为常量,因此可以应用于应力梯度较大的地方,精度较高。显然,其计算也
24、较复杂。应力梯度较大的地方,精度较高。显然,其计算也较复杂。 7.4 四节点四边形等参数单元四节点四边形等参数单元 1、等参数单元的概念、等参数单元的概念 现在,我们从任意四边形单元着手,介绍等参数单元现在,我们从任意四边形单元着手,介绍等参数单元的概念。的概念。1234xy图图7-10 任意四边形单元任意四边形单元 前面讲到的四节点矩形单元虽然比较简单,但难以应前面讲到的四节点矩形单元虽然比较简单,但难以应用于斜线边界。图用于斜线边界。图7-10所示四节点任意四边形单元容易适所示四节点任意四边形单元容易适应这种边界,应这种边界,但要在整体坐标系内,写出它的统一的形函但要在整体坐标系内,写出它
25、的统一的形函数又是相当复杂和困难的数又是相当复杂和困难的。 但是若能找到它与一个规则正方形的关系,就能写出它但是若能找到它与一个规则正方形的关系,就能写出它的统一的位移模式,这可以通过坐标变换来解决。在图的统一的位移模式,这可以通过坐标变换来解决。在图7-10所示四边形单元上,用等分四边的两族直线分割该四边所示四边形单元上,用等分四边的两族直线分割该四边形,以两族曲线的中心(形,以两族曲线的中心( =0、 =0 )为原点,沿)为原点,沿 、 增增大的方向作大的方向作 轴和轴和 轴,并令四边的轴,并令四边的 =1、 =1,就得出,就得出一组新坐标系一组新坐标系 (图(图7-11)。)。1234x
26、y图图11 实际单元实际单元 = -1 = 1 = 1 = -1 这里,这里, 、 是一种局部(单元)是一种局部(单元)坐标,它只应用于单元范围内。坐标,它只应用于单元范围内。而而x,y是整体(结构)坐标,是整体(结构)坐标,它适用于所有的单元。图中的它适用于所有的单元。图中的任意四边形单元是研究对象,任意四边形单元是研究对象,称为实际单元称为实际单元。参照式(参照式(7-2)和()和(7-3)P6,此基本单元位移函数可写为:,此基本单元位移函数可写为: ,41iiiuNu41iiivNv(7-40) 1234 = -1 = 1 = -1 = 1图图7-12 基本单元基本单元 为了得出实际单元
27、的位移模式和局部坐标与整体坐标为了得出实际单元的位移模式和局部坐标与整体坐标之间的变换关系,引入一个四节点的正方形单元,称之间的变换关系,引入一个四节点的正方形单元,称基本基本单元单元(图(图7-12) 。其中,形函数应为:其中,形函数应为:11411141114111414321NNNN引入新变量引入新变量 i、i (i=1,2,3,4)基本单元的形函数被写成:基本单元的形函数被写成:14, 1aa13,2aa12, 1bb14, 3bb)4 , 3 , 2 , 1()1)(1 (41iNiii(7-41) 现在,把基本单元的位移模式(现在,把基本单元的位移模式(7-40)和形函数式)和形函
28、数式(7-41)移用于图()移用于图(7-11)所示的实际单元)所示的实际单元,则实际单元,则实际单元的位移模式取为的位移模式取为: ,41iiiuNu41iiivNv(7-40)在结点处在结点处 :在其它结点处:在其它结点处:)4 , 3 , 2 , 1(1iNi)4 , 3 , 2 , 1(0iNi且,式中的形函数且,式中的形函数Ni仍由式(仍由式(7-41)确定。而把式()确定。而把式(7-41)中的中的 、 理解为理解为图图7-11所示所示实际单元的局部坐标,实际单元的局部坐标, i、 i便便是实际单元中节点是实际单元中节点i的局部坐标。的局部坐标。 41iixNx41iiyNy(7-
29、42)利用形函数的上述性质,可以将任意四边形的整体坐标写利用形函数的上述性质,可以将任意四边形的整体坐标写成:成:ixiy任意四边形单元中结点的整体坐标,任意四边形单元中结点的整体坐标,如果它已知,那么如果它已知,那么(7-42)表示了局部坐表示了局部坐标与整体坐标的变换标与整体坐标的变换 另一方面,式(另一方面,式(7-42)表明了实际单元中局部坐标()表明了实际单元中局部坐标( 、 )与整体坐标()与整体坐标(x、y)的一一对应关系,是一个坐标变)的一一对应关系,是一个坐标变换式。换式。 实际单元是任意四边形四节点单元,基本单元是正方实际单元是任意四边形四节点单元,基本单元是正方形单元,可
30、以认为:实际单元是对基本单元通过变换得来形单元,可以认为:实际单元是对基本单元通过变换得来的。的。由于实际单元的位移模式中采用了基本单元等同的形由于实际单元的位移模式中采用了基本单元等同的形函数,这个实际单元就称为等参数单元函数,这个实际单元就称为等参数单元。 类似于本章类似于本章3.2节进行的四结点矩形单元的特性分析,节进行的四结点矩形单元的特性分析,可以建立等参单元的应变矩阵、应力矩阵、刚度矩阵、节可以建立等参单元的应变矩阵、应力矩阵、刚度矩阵、节点力向量等的计算公式。与前面不同之处在于,在等参数点力向量等的计算公式。与前面不同之处在于,在等参数单元法中,要将对整体坐标单元法中,要将对整体
31、坐标x、y的导数计算和积分计算转的导数计算和积分计算转换为对局部坐标换为对局部坐标 、 的导数计算和积分计算。的导数计算和积分计算。例:实际单元的结点整体坐标如图例:实际单元的结点整体坐标如图(a)中括号内数字所示,中括号内数字所示,基本单元的结点局部坐标如图基本单元的结点局部坐标如图(b)中括号内数字所示。中括号内数字所示。图图(a) 实际单元实际单元1(0,0)(1) 试验证基本单元上的结点局部坐标与实际单元上对应点的整体坐试验证基本单元上的结点局部坐标与实际单元上对应点的整体坐标的对应关系。标的对应关系。(2) 求基本单元的局部坐标原点求基本单元的局部坐标原点( ),在实际单元上的整体坐
32、在实际单元上的整体坐标标(x,y)是多少?是多少?04(0,1)3(1,2)2(2,0)xy0(3/4,3/4)01(1,1)2(1,1)3(1,1)4(1,1)图图(b) 基本单元基本单元解解(1)以)以3结点为例,根据形函数的性质:结点为例,根据形函数的性质:在在3结点处结点处 :13N1240NNN应用式(应用式(7-42)3443322113xxNxNxNxNx3443322113yyNyNyNyNy说明了由基本单元上结点的局部坐标可映射出说明了由基本单元上结点的局部坐标可映射出实际单元上对应的结点整体坐标实际单元上对应的结点整体坐标(2)由于)由于 ,由(,由(7-40)P33 得:
33、得:041)1)(1(411N41)1)(1(412N41)1)(1(413N41)1)(1(414N代入(代入(7-42)43)0120(4141iixNx43)1200(4141iiyNy由上例可知:利用(由上例可知:利用(7-42)在基本单元上任意一点)在基本单元上任意一点 ,都可以在实际单元上找到一个对应点的坐标都可以在实际单元上找到一个对应点的坐标(x,y), 这样就把这样就把实际单元与基本单元紧密地联系起来。反之,则比较困难,实际单元与基本单元紧密地联系起来。反之,则比较困难,这是因为形函数是一个二次函数。为了避开这个困难,一这是因为形函数是一个二次函数。为了避开这个困难,一般都假
34、定基本单元上已知点去求实际单元上的对应点。般都假定基本单元上已知点去求实际单元上的对应点。、 2、应变矩阵、应变矩阵 单元的几何方程与式(单元的几何方程与式(7-8)、()、(7-9)相同,)相同,即:即:eB(7-8) mljiBBBBB (7-9)eeBBBBB4321(7-43)式中式中(7-44))4 , 3 , 2 , 1(00,iNNNNBxiyiyixii?这里采用记号这里采用记号 yNNxNNiyiixi, 由于形函数式(由于形函数式(7-41)是用局部坐标)是用局部坐标 、 给出的,给出的,将将 、 看作看作x、y的函数,根据的函数,根据复合函数的求导规则复合函数的求导规则,
35、有:有:yyNxxNNyyNxxNNiiiiii上式可记为:上式可记为: yixiiiNNyxyxNN,(7-45)上式右边第一个矩阵称为雅可比(上式右边第一个矩阵称为雅可比(Jacobi)矩阵:)矩阵: 其逆矩阵为:其逆矩阵为: 式中式中|J|为雅可比行列式为雅可比行列式 (7-48),xyyxyxyxJ由式(由式(7-42)p35,有,有 ,yxyxJ(7-46),11xyxyJJ(7-47)44,1144,11,iiiiiiiiiiiixNxxNxyNyyNy(7-49)由式(由式(7-41)p34,有,有(7-50)4/ )1 (4/ )1 (,iiiiiiNN由式(由式(7-45)p
36、41,有,有(7-51),1,iiyixiNNJNN式中式中分别由式(分别由式(7-47)和()和(7-50)确定。)确定。 从而由式(从而由式(3-43)、)、(3-44)确定出应变矩阵)确定出应变矩阵B。1J,iiNN和和3、应力矩阵、应力矩阵应力矩阵仍由下式得到应力矩阵仍由下式得到 BDS (7-52) 4、单元刚度矩阵、单元刚度矩阵单元刚度矩阵是一个单元刚度矩阵是一个88的矩阵,仍为的矩阵,仍为 AThdxdyBDBk 1111|ddJhBDBkT(7-53)由于由于B是用局部坐标系是用局部坐标系 、 给出的,坐标变换时有面积微给出的,坐标变换时有面积微元公式元公式 ddJdxdy|
37、因此,因此,k可由下式计算可由下式计算 5、等价节点力向量、等价节点力向量(1)体积力)体积力 设单元的体积力是设单元的体积力是 ,则等价节点力公式,则等价节点力公式为为TVyVxVppP,AVyVxiVyVxiVhdxdyppNppF)4 , 3 , 2 , 1(1111 iddJhppNVyVxi(7-54)(2) 表面力表面力 设单元的某边(如对应的设单元的某边(如对应的 =1)上作用有表面力)上作用有表面力 ps =psx psyT,则该边上节点的节点力为,则该边上节点的节点力为hdlppNppFSySxilSySxiSdyxhppNSySxi2211(7-55) 上式中,用到坐标变换
38、时线积分微元公式上式中,用到坐标变换时线积分微元公式(当当 =常量时常量时) p37dyxdl22对于对于 =1的边界表面力,等价节点力的计算过程完全一样。的边界表面力,等价节点力的计算过程完全一样。 6、关于高斯(、关于高斯(Gauss)积分)积分 在计算单元刚度矩阵式(在计算单元刚度矩阵式(7-53)和等价节点载荷向量式)和等价节点载荷向量式(7-54)、()、(7-55)时,由于被积函数比较复杂,通常可)时,由于被积函数比较复杂,通常可采用数值积分。即在单元上选择某些点,称为采用数值积分。即在单元上选择某些点,称为积分点,求出被积函数在这些积分点上的数值,再用一些积分点,求出被积函数在这
39、些积分点上的数值,再用一些权函数乘这些函数值后求和,就可得到权函数乘这些函数值后求和,就可得到近似积分值近似积分值。高斯。高斯积分法是数值积分法中具有较高精度的方法。现简要介绍积分法是数值积分法中具有较高精度的方法。现简要介绍如下:如下:(1)一维高斯积分公式)一维高斯积分公式 )()(111iinifHdf(7-56) 式中式中Hi对应于积分点对应于积分点 i的权函数,对于的权函数,对于n=24个积分个积分点,其坐标点,其坐标 i和权函数和权函数Hi列于表列于表7-1中。中。niHi20.57735026921.000000000030.77459666920.00000000000.555
40、55555560.888888888940.86113631160.33998104360.34785484510.6521451549表表3-1 高斯求积公式中的积分点坐标和权系数高斯求积公式中的积分点坐标和权系数(2)二维高斯积分公式)二维高斯积分公式 (7-57)),(),(111111ijjininifHHddf积分点积分点 j 、 i 和权函数和权函数Hi 、Hj同样可按表同样可按表3-1进行取值。进行取值。对式(对式(7-53)进行高斯积分,有)进行高斯积分,有 (7-58)hJBDBHHkTnmlnlm| 11对式(对式(7-54)进行高斯积分,有)进行高斯积分,有 (7-59)
41、|11JhppNHHFVyVxinmlnlmiV对式(对式(7-55)进行高斯积分,有)进行高斯积分,有221ddyddxhppNHFSySximlmiS(7-60)7.5 八节点四边形等参数单元八节点四边形等参数单元 上节讨论的四结点四边形等参数单元有时仍然不够理上节讨论的四结点四边形等参数单元有时仍然不够理想。想。一是其实际单元为直线边界,不能准确拟合物体的曲一是其实际单元为直线边界,不能准确拟合物体的曲线边界;二是位移模式的阶次还不够高,影响计算精度。线边界;二是位移模式的阶次还不够高,影响计算精度。为此,本节介绍一种精度更高、应用广泛的八节点四边形为此,本节介绍一种精度更高、应用广泛的
42、八节点四边形等参数单元。其实际单元和基本单元如图等参数单元。其实际单元和基本单元如图7-13、7-14所示。所示。图图7-14 基本单元基本单元 = -1 = 1 = -1 = 113567824图图7-13 实际单元实际单元xy12345678 = -1 = 1 = -1 = 1基本单元的位移模式可取为:基本单元的位移模式可取为: 2162152141321211109282726524321aaaaaaaavaaaaaaaau(7-61)采用形函数表示,将位移模式写成采用形函数表示,将位移模式写成iiiiiivNvuNu8181(7-62)式中式中2/ )1)(1 (4/ ) 1)(1)(
43、1 (2/ )1)(1 (4/ ) 1)(1)(1 (2/ )1)(1 (4/ ) 1)(1)(1 (2/ )1)(1 (4/ ) 1)(1)(1 (287265243221NNNNNNNN(7-63)仿照位移模式,将坐标变换式取为仿照位移模式,将坐标变换式取为 8181iiiiiiyNyxNx(7-64) 显然,该坐标变换式将显然,该坐标变换式将平面上的正方形映射为平面上的正方形映射为xy平面上的曲边四边形。平面上的曲边四边形。xy平面上每一条边都是一条二次平面上每一条边都是一条二次曲线,它完全由对应边上曲线,它完全由对应边上3个结点的坐标唯一确定。因此,个结点的坐标唯一确定。因此,单元是协
44、调的,同时也可证明,单元的位移函数反映刚单元是协调的,同时也可证明,单元的位移函数反映刚体位移和常应变,具有完备性。体位移和常应变,具有完备性。 有关八节点四边形单元的特性分析和等价节点力计算有关八节点四边形单元的特性分析和等价节点力计算过程与四节点四边形单元完全相同,具体公式形式也一致。过程与四节点四边形单元完全相同,具体公式形式也一致。区别仅在于两种单元有关矩阵的维数不同,具体是:区别仅在于两种单元有关矩阵的维数不同,具体是: (1)四节点单元的)四节点单元的e是一个是一个81列阵;而八结点列阵;而八结点单元的单元的e是一个是一个161列阵。列阵。 (2)四节点单元的)四节点单元的B是一个是一个38矩阵,由矩阵,由B1B4 4个子块个子块阵组成;而八节点单元的阵组成;而八节点单元的B是一个是一个316矩阵,由矩阵,由B1B8 8个个子块阵组成。个个子块阵组成。 (3)四节点单元的)四节点单元的k是一个是一个88矩阵,而八节点单元的矩阵,而八节点单元的k是一个是一个1616矩阵。矩阵。 (4)四结点单元某边表面力引起的等价节点力)四结点单元某边表面力引起的等价节点力 Fs i公式公式对应于该边上的对应于该边上的2个节点;而八节点单元对应于该边上的个节点;而八节点单元对应于该边上的3个节点。个节点。7.6 八节点四边形等参数单元计算程序八节点四边形等参数单元计算程序