(結構分析的有限元法課件)第三章桿件系統(tǒng)的有限單元法

第三章桿件系統(tǒng)的有限單元法有限元法與結構力學中剛架計算的矩陣位移法有密切關系，因此我們就從桿件系統(tǒng)的有限元法入手，來了解有限元法的基本概念和過程。第三章桿件系統(tǒng)的有限單元法有限元法與結構力學中剛架計算的矩陣1桿單元(TrussorBar)桿的定義：兩端鉸接，只受軸向力作用的基本構件ijxui(FNi)uj(FNj)桿單元的結點位移和結點力桿單元(TrussorBar)桿的定義：兩端鉸接，只受軸2根據(jù)材料力學的有關知識，我們可以立刻寫出桿單元的結點位移與結點力之間的關系為為材料的彈性模量，為桿單元的長度

桿單元的剛度矩陣為桿單元的橫截面面積寫成矩陣形式就是根據(jù)材料力學的有關知識，我們可以立刻寫出桿單元的結點位移與結3桿單元的剛度矩陣桿單元的剛度矩陣4單元的剛度矩陣單元的剛度矩陣5平面梁單元(Beam)平面梁單元(Beam)6平面梁單元的位移模式平面梁單元的位移模式7待定參數(shù)的確定待定參數(shù)的確定8MATLAB不僅可以進行數(shù)值運算，也能進行符號運算。如式(3.20)中的矩陣Au和Av的求逆運算，我們可以在MATLAB的命令窗口下輸入>>symsL>>Au=[101L];>>Av=[100001001LL^2L^3012*L3*L^2];第一句是定義符號變量L，后面定義兩個矩陣Au和Av。然后我們再輸入下面求逆的命令>>inv(Au) ans=[1,0][-1/L,1/L]>>inv(Av)ans=[1,0,0,0][0,1,0,0][-3/L^2,-2/L,3/L^2,-1/L][2/L^3,1/L^2,-2/L^3,1/L^2]MATLAB不僅可以進行數(shù)值運算，也能進行符號運算。如式(39位移模式位移模式10應變矩陣B應變矩陣B11單元剛度矩陣單元剛度矩陣12上面剛度矩陣的推導也可以應用MATLAB的符號運算功能，計算的程序段如下>>symsELxy>>A=[100000010000001000100L0001L0L^2L^300102*L3*L^2];>>Hu=[100x00];>>Hv=[01x0x^2x^3];>>B=[diff(Hu,x,1);-y*diff(Hv,x,2)]*inv(A);>>Ke=int(E*transpose(B)*B,x,0,L)上面剛度矩陣的推導也可以應用MATLAB的符號運算功能，計算13上述命令執(zhí)行后，得到結果如下：Ke=[E/L,0,0,-E/L,0,0][0,12*E*y^2/L^3,6*E*y^2/L^2,0,-12*E*y^2/L^3,6*E*y^2/L^2][0,6*E*y^2/L^2,4*E*y^2/L,0,-6*E*y^2/L^2,2*E*y^2/L][-E/L,0,0,E/L,0,0][0,-12*E*y^2/L^3,-6*E*y^2/L^2,0,12*E*y^2/L^3,-6*E*y^2/L^2][0,6*E*y^2/L^2,2*E*y^2/L,0,-6*E*y^2/L^2,4*E*y^2/L]上述命令執(zhí)行后，得到結果如下：14等效結點力計算所謂等效結點力，是指非結點載荷按照虛功相等的原則分配到單元結點上的力。等效表面力的普遍公式對于作用在桿件單元上的分布力，我們可以簡寫為等效結點力計算所謂等效結點力，是指非結點載荷按照虛功相等的15分布軸向力分布軸向力16推導線性分布軸力的等效結點力公式的MATLAB程序如下symsLxx1x2p1p2Au=[101L];H1=[1x];Nu=H1*inv(Au);x1=0;x2=L;L1=(x-x2)/(x1-x2);L2=(x-x1)/(x2-x1);p=p1*L1+p2*L2;fe=int(transpose(Nu)*p,x,0,L);fe=simple(fe);MATLAB給出的等效結點力公式是fe=[1/6*L*(p2+2*p1)][1/6*L*(2*p2+p1)]推導線性分布軸力的等效結點力公式的MATLAB程序如下17如果想推導其他分布形式的等效結點力，只要修改前面程序中的p的形式。比如我們要推導二次拋物線分布軸力的等效結點力，我們只要把p寫成二次的Lagrange插值形式，即symsLxx1x2x3p1p2p3Au=[101L];H1=[1x];Nu=H1*inv(Au);x1=0;x2=L/2;x3=L;L1=(x-x2)*(x-x3)/(x1-x2)/(x1-x3);L2=(x-x1)*(x-x3)/(x2-x1)/(x2-x3);L3=(x-x1)*(x-x2)/(x3-x1)/(x3-x2);p=p1*L1+p2*L2+p3*L3;fe=int(transpose(Nu)*p,x,0,L);fe=simple(fe);執(zhí)行后得到的等效結點力公式是fe=[1/6*L*(p1+2*p2)][1/6*L*(2*p2+p3)]按這個思路可以推導更高階次分布軸力的等效結點力。如果想推導其他分布形式的等效結點力，只要修改前面程序中的p的18分布橫向力分布橫向力19(結構分析的有限元法課件)第三章桿件系統(tǒng)的有限單元法20推導線性分布橫向力等效結點力的MATLAB程序為symsLxx1x2p1p2Av=[100001001LL^2L^3012*L3*L^2];H2=[1xx^2x^3];Nv=H2*inv(Av);x1=0;x2=L;L1=(x-x2)/(x1-x2);L2=(x-x1)/(x2-x1);p=p1*L1+p2*L2;fe=int(transpose(Nv)*p,x,0,L);fe=simple(fe);執(zhí)行后，得到等效結點力為fe=[1/20*L*(7*p1+3*p2)][1/60*L^2*(3*p1+2*p2)][1/20*L*(3*p1+7*p2)][-1/60*L^2*(2*p1+3*p2)]推導線性分布橫向力等效結點力的MATLAB程序為21分布彎曲力矩分布彎曲力矩22symsLxx1x2m1m2Av=[100001001LL^2L^3012*L3*L^2];H2=[1xx^2x^3];Nv=H2*inv(Av);x1=0;x2=L;L1=(x-x2)/(x1-x2);L2=(x-x1)/(x2-x1);m=m1*L1+m2*L2;fe=int(transpose(diff(Nv,x,1))*m,x,0,L);fe=simple(fe);symsLxx1x2m1m223單元的坐標變換在前面幾節(jié)中，推導單元剛度矩陣時采用的全部是局部坐標系，它的坐標軸方向是由單元即桿或梁的截面主方向確定的。采用這樣的坐標系，可以得到具有統(tǒng)一形式的單元剛度矩陣。但是，實際的結構一般由具有不同方向和處于不同位置的桿或梁所構成的。由不同方向單元組成的結構，它的整體剛度矩陣并不能由局部坐標下的單元剛度矩陣簡單地疊加而成，所以必須建立一個統(tǒng)一的整體坐標系。計算時先將單元上的結點力和位移轉換到整體坐標系，單元剛度矩陣亦作坐標變換，才可按疊加規(guī)則直接相加組成整體剛度矩陣。單元的坐標變換在前面幾節(jié)中，推導單元剛度矩陣時采用24坐標變換坐標變換25平面桿單元的坐標轉換矩陣平面桿單元的坐標轉換矩陣26整體坐標系下的單元剛度矩陣整體坐標系下的單元剛度矩陣27平面梁單元的轉換矩陣平面梁單元的轉換矩陣28平面梁單元的轉換矩陣平面梁單元的轉換矩陣29(結構分析的有限元法課件)第三章桿件系統(tǒng)的有限單元法30(結構分析的有限元法課件)第三章桿件系統(tǒng)的有限單元法31(結構分析的有限元法課件)第三章桿件系統(tǒng)的有限單元法32整體剛度矩陣1324[1][2][3][4]整體剛度矩陣1324[1][2][3][4]331324[1][2][3][4]1324[1][2][3][4]34(結構分析的有限元法課件)第三章桿件系統(tǒng)的有限單元法351324[1][2][3][4]集成后的整體剛度矩陣1324[1][2][3][4]集成后的整體剛度矩陣36(e)ij(e)ij37整體剛度矩陣整體結點力向量整體平衡方程整體剛度矩陣整體結點力向量整體平衡方程38(結構分析的有限元法課件)第三章桿件系統(tǒng)的有限單元法39整體剛度矩陣的特性整體剛度矩陣的特性40邊界約束條件的處理邊界約束條件的處理41劃行劃列法設有方程約束條件劃行劃列法設有方程約束條件42劃行劃列的程序劃行劃列的程序43乘大數(shù)法乘大數(shù)法44乘大數(shù)的程序乘大數(shù)的程序45第二類約束條件：例子1第二類約束條件：例子146第二類約束條件：例子2第二類約束條件：例子247第二類約束條件：例子3第二類約束條件：例子348第二類約束條件的處理辦法第二類約束條件的處理辦法49第二類約束條件程序實現(xiàn)第二類約束條件程序實現(xiàn)50單元結點反力1、結點力中包含了分布力的等效結點力2、工程中經(jīng)常要用的剪力、軸力是在單元的局部坐標系下定義的單元結點反力1、結點力中包含了分布力的等效結點力2、工程中經(jīng)51單元結點反力的程序單元結點反力的程序52溫度應力溫度變化結構變形結構受到約束產(chǎn)生溫度應力溫度應力溫度變化結構變形結構受到約束產(chǎn)生溫度應力53有限元中的溫度應力處理有限元中的溫度應力處理54梁單元的溫度荷載梁單元的溫度荷載55溫度荷載的程序溫度荷載的程序56算例一算例一57有限元程序框架有限元程序框架58算例二算例二59有限元離散模型有限元離散模型60第三章桿件系統(tǒng)的有限單元法有限元法與結構力學中剛架計算的矩陣位移法有密切關系，因此我們就從桿件系統(tǒng)的有限元法入手，來了解有限元法的基本概念和過程。第三章桿件系統(tǒng)的有限單元法有限元法與結構力學中剛架計算的矩陣61桿單元(TrussorBar)桿的定義：兩端鉸接，只受軸向力作用的基本構件ijxui(FNi)uj(FNj)桿單元的結點位移和結點力桿單元(TrussorBar)桿的定義：兩端鉸接，只受軸62根據(jù)材料力學的有關知識，我們可以立刻寫出桿單元的結點位移與結點力之間的關系為為材料的彈性模量，為桿單元的長度

桿單元的剛度矩陣為桿單元的橫截面面積寫成矩陣形式就是根據(jù)材料力學的有關知識，我們可以立刻寫出桿單元的結點位移與結63桿單元的剛度矩陣桿單元的剛度矩陣64單元的剛度矩陣單元的剛度矩陣65平面梁單元(Beam)平面梁單元(Beam)66平面梁單元的位移模式平面梁單元的位移模式67待定參數(shù)的確定待定參數(shù)的確定68MATLAB不僅可以進行數(shù)值運算，也能進行符號運算。如式(3.20)中的矩陣Au和Av的求逆運算，我們可以在MATLAB的命令窗口下輸入>>symsL>>Au=[101L];>>Av=[100001001LL^2L^3012*L3*L^2];第一句是定義符號變量L，后面定義兩個矩陣Au和Av。然后我們再輸入下面求逆的命令>>inv(Au) ans=[1,0][-1/L,1/L]>>inv(Av)ans=[1,0,0,0][0,1,0,0][-3/L^2,-2/L,3/L^2,-1/L][2/L^3,1/L^2,-2/L^3,1/L^2]MATLAB不僅可以進行數(shù)值運算，也能進行符號運算。如式(369位移模式位移模式70應變矩陣B應變矩陣B71單元剛度矩陣單元剛度矩陣72上面剛度矩陣的推導也可以應用MATLAB的符號運算功能，計算的程序段如下>>symsELxy>>A=[100000010000001000100L0001L0L^2L^300102*L3*L^2];>>Hu=[100x00];>>Hv=[01x0x^2x^3];>>B=[diff(Hu,x,1);-y*diff(Hv,x,2)]*inv(A);>>Ke=int(E*transpose(B)*B,x,0,L)上面剛度矩陣的推導也可以應用MATLAB的符號運算功能，計算73上述命令執(zhí)行后，得到結果如下：Ke=[E/L,0,0,-E/L,0,0][0,12*E*y^2/L^3,6*E*y^2/L^2,0,-12*E*y^2/L^3,6*E*y^2/L^2][0,6*E*y^2/L^2,4*E*y^2/L,0,-6*E*y^2/L^2,2*E*y^2/L][-E/L,0,0,E/L,0,0][0,-12*E*y^2/L^3,-6*E*y^2/L^2,0,12*E*y^2/L^3,-6*E*y^2/L^2][0,6*E*y^2/L^2,2*E*y^2/L,0,-6*E*y^2/L^2,4*E*y^2/L]上述命令執(zhí)行后，得到結果如下：74等效結點力計算所謂等效結點力，是指非結點載荷按照虛功相等的原則分配到單元結點上的力。等效表面力的普遍公式對于作用在桿件單元上的分布力，我們可以簡寫為等效結點力計算所謂等效結點力，是指非結點載荷按照虛功相等的75分布軸向力分布軸向力76推導線性分布軸力的等效結點力公式的MATLAB程序如下symsLxx1x2p1p2Au=[101L];H1=[1x];Nu=H1*inv(Au);x1=0;x2=L;L1=(x-x2)/(x1-x2);L2=(x-x1)/(x2-x1);p=p1*L1+p2*L2;fe=int(transpose(Nu)*p,x,0,L);fe=simple(fe);MATLAB給出的等效結點力公式是fe=[1/6*L*(p2+2*p1)][1/6*L*(2*p2+p1)]推導線性分布軸力的等效結點力公式的MATLAB程序如下77如果想推導其他分布形式的等效結點力，只要修改前面程序中的p的形式。比如我們要推導二次拋物線分布軸力的等效結點力，我們只要把p寫成二次的Lagrange插值形式，即symsLxx1x2x3p1p2p3Au=[101L];H1=[1x];Nu=H1*inv(Au);x1=0;x2=L/2;x3=L;L1=(x-x2)*(x-x3)/(x1-x2)/(x1-x3);L2=(x-x1)*(x-x3)/(x2-x1)/(x2-x3);L3=(x-x1)*(x-x2)/(x3-x1)/(x3-x2);p=p1*L1+p2*L2+p3*L3;fe=int(transpose(Nu)*p,x,0,L);fe=simple(fe);執(zhí)行后得到的等效結點力公式是fe=[1/6*L*(p1+2*p2)][1/6*L*(2*p2+p3)]按這個思路可以推導更高階次分布軸力的等效結點力。如果想推導其他分布形式的等效結點力，只要修改前面程序中的p的78分布橫向力分布橫向力79(結構分析的有限元法課件)第三章桿件系統(tǒng)的有限單元法80推導線性分布橫向力等效結點力的MATLAB程序為symsLxx1x2p1p2Av=[100001001LL^2L^3012*L3*L^2];H2=[1xx^2x^3];Nv=H2*inv(Av);x1=0;x2=L;L1=(x-x2)/(x1-x2);L2=(x-x1)/(x2-x1);p=p1*L1+p2*L2;fe=int(transpose(Nv)*p,x,0,L);fe=simple(fe);執(zhí)行后，得到等效結點力為fe=[1/20*L*(7*p1+3*p2)][1/60*L^2*(3*p1+2*p2)][1/20*L*(3*p1+7*p2)][-1/60*L^2*(2*p1+3*p2)]推導線性分布橫向力等效結點力的MATLAB程序為81分布彎曲力矩分布彎曲力矩82symsLxx1x2m1m2Av=[100001001LL^2L^3012*L3*L^2];H2=[1xx^2x^3];Nv=H2*inv(Av);x1=0;x2=L;L1=(x-x2)/(x1-x2);L2=(x-x1)/(x2-x1);m=m1*L1+m2*L2;fe=int(transpose(diff(Nv,x,1))*m,x,0,L);fe=simple(fe);symsLxx1x2m1m283單元的坐標變換在前面幾節(jié)中，推導單元剛度矩陣時采用的全部是局部坐標系，它的坐標軸方向是由單元即桿或梁的截面主方向確定的。采用這樣的坐標系，可以得到具有統(tǒng)一形式的單元剛度矩陣。但是，實際的結構一般由具有不同方向和處于不同位置的桿或梁所構成的。由不同方向單元組成的結構，它的整體剛度矩陣并不能由局部坐標下的單元剛度矩陣簡單地疊加而成，所以必須建立一個統(tǒng)一的整體坐標系。計算時先將單元上的結點力和位移轉換到整體坐標系，單元剛度矩陣亦作坐標變換，才可按疊加規(guī)則直接相加組成整體剛度矩陣。單元的坐標變換在前面幾節(jié)中，推導單元剛度矩陣時采用84坐標變換坐標變換85平面桿單元的坐標轉換矩陣平面桿單元的坐標轉換矩陣86整體坐標系下的單元剛度矩陣整體坐標系下的單元剛度矩陣87平面梁單元的轉換矩陣平面梁單元的轉換矩陣88平面