(结构分析的有限元法课件)第三章杆件系统的有限单元法

第三章杆件系统的有限单元法有限元法与结构力学中刚架计算的矩阵位移法有密切关系，因此我们就从杆件系统的有限元法入手，来了解有限元法的基本概念和过程。第三章杆件系统的有限单元法有限元法与结构力学中刚架计算的矩阵1杆单元(TrussorBar)杆的定义：两端铰接，只受轴向力作用的基本构件ijxui(FNi)uj(FNj)杆单元的结点位移和结点力杆单元(TrussorBar)杆的定义：两端铰接，只受轴2根据材料力学的有关知识，我们可以立刻写出杆单元的结点位移与结点力之间的关系为为材料的弹性模量，为杆单元的长度

杆单元的刚度矩阵为杆单元的横截面面积写成矩阵形式就是根据材料力学的有关知识，我们可以立刻写出杆单元的结点位移与结3杆单元的刚度矩阵杆单元的刚度矩阵4单元的刚度矩阵单元的刚度矩阵5平面梁单元(Beam)平面梁单元(Beam)6平面梁单元的位移模式平面梁单元的位移模式7待定参数的确定待定参数的确定8MATLAB不仅可以进行数值运算，也能进行符号运算。如式(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不仅可以进行数值运算，也能进行符号运算。如式(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上述命令执行后，得到结果如下：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]上述命令执行后，得到结果如下：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);执行后得到的等效结点力公式是fe=[1/6*L*(p1+2*p2)][1/6*L*(2*p2+p3)]按这个思路可以推导更高阶次分布轴力的等效结点力。如果想推导其他分布形式的等效结点力，只要修改前面程序中的p的18分布横向力分布横向力19(结构分析的有限元法课件)第三章杆件系统的有限单元法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);执行后，得到等效结点力为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单元的坐标变换在前面几节中，推导单元刚度矩阵时采用的全部是局部坐标系，它的坐标轴方向是由单元即杆或梁的截面主方向确定的。采用这样的坐标系，可以得到具有统一形式的单元刚度矩阵。但是，实际的结构一般由具有不同方向和处于不同位置的杆或梁所构成的。由不同方向单元组成的结构，它的整体刚度矩阵并不能由局部坐标下的单元刚度矩阵简单地叠加而成，所以必须建立一个统一的整体坐标系。计算时先将单元上的结点力和位移转换到整体坐标系，单元刚度矩阵亦作坐标变换，才可按叠加规则直接相加组成整体刚度矩阵。单元的坐标变换在前面几节中，推导单元刚度矩阵时采用24坐标变换坐标变换25平面杆单元的坐标转换矩阵平面杆单元的坐标转换矩阵26整体坐标系下的单元刚度矩阵整体坐标系下的单元刚度矩阵27平面梁单元的转换矩阵平面梁单元的转换矩阵28平面梁单元的转换矩阵平面梁单元的转换矩阵29(结构分析的有限元法课件)第三章杆件系统的有限单元法30(结构分析的有限元法课件)第三章杆件系统的有限单元法31(结构分析的有限元法课件)第三章杆件系统的有限单元法32整体刚度矩阵1324[1][2][3][4]整体刚度矩阵1324[1][2][3][4]331324[1][2][3][4]1324[1][2][3][4]34(结构分析的有限元法课件)第三章杆件系统的有限单元法351324[1][2][3][4]集成后的整体刚度矩阵1324[1][2][3][4]集成后的整体刚度矩阵36(e)ij(e)ij37整体刚度矩阵整体结点力向量整体平衡方程整体刚度矩阵整体结点力向量整体平衡方程38(结构分析的有限元法课件)第三章杆件系统的有限单元法39整体刚度矩阵的特性整体刚度矩阵的特性40边界约束条件的处理边界约束条件的处理41划行划列法设有方程约束条件划行划列法设有方程约束条件42划行划列的程序划行划列的程序43乘大数法乘大数法44乘大数的程序乘大数的程序45第二类约束条件：例子1第二类约束条件：例子146第二类约束条件：例子2第二类约束条件：例子247第二类约束条件：例子3第二类约束条件：例子348第二类约束条件的处理办法第二类约束条件的处理办法49第二类约束条件程序实现第二类约束条件程序实现50单元结点反力1、结点力中包含了分布力的等效结点力2、工程中经常要用的剪力、轴力是在单元的局部坐标系下定义的单元结点反力1、结点力中包含了分布力的等效结点力2、工程中经51单元结点反力的程序单元结点反力的程序52温度应力温度变化结构变形结构受到约束产生温度应力温度应力温度变化结构变形结构受到约束产生温度应力53有限元中的温度应力处理有限元中的温度应力处理54梁单元的温度荷载梁单元的温度荷载55温度荷载的程序温度荷载的程序56算例一算例一57有限元程序框架有限元程序框架58算例二算例二59有限元离散模型有限元离散模型60第三章杆件系统的有限单元法有限元法与结构力学中刚架计算的矩阵位移法有密切关系，因此我们就从杆件系统的有限元法入手，来了解有限元法的基本概念和过程。第三章杆件系统的有限单元法有限元法与结构力学中刚架计算的矩阵61杆单元(TrussorBar)杆的定义：两端铰接，只受轴向力作用的基本构件ijxui(FNi)uj(FNj)杆单元的结点位移和结点力杆单元(TrussorBar)杆的定义：两端铰接，只受轴62根据材料力学的有关知识，我们可以立刻写出杆单元的结点位移与结点力之间的关系为为材料的弹性模量，为杆单元的长度

杆单元的刚度矩阵为杆单元的横截面面积写成矩阵形式就是根据材料力学的有关知识，我们可以立刻写出杆单元的结点位移与结63杆单元的刚度矩阵杆单元的刚度矩阵64单元的刚度矩阵单元的刚度矩阵65平面梁单元(Beam)平面梁单元(Beam)66平面梁单元的位移模式平面梁单元的位移模式67待定参数的确定待定参数的确定68MATLAB不仅可以进行数值运算，也能进行符号运算。如式(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不仅可以进行数值运算，也能进行符号运算。如式(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上述命令执行后，得到结果如下：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]上述命令执行后，得到结果如下：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);执行后得到的等效结点力公式是fe=[1/6*L*(p1+2*p2)][1/6*L*(2*p2+p3)]按这个思路可以推导更高阶次分布轴力的等效结点力。如果想推导其他分布形式的等效结点力，只要修改前面程序中的p的78分布横向力分布横向力79(结构分析的有限元法课件)第三章杆件系统的有限单元法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);执行后，得到等效结点力为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单元的坐标变换在前面几节中，推导单元刚度矩阵时采用的全部是局部坐标系，它的坐标轴方向是由单元即杆或梁的截面主方向确定的。采用这样的坐标系，可以得到具有统一形式的单元刚度矩阵。但是，实际的结构一般由具有不同方向和处于不同位置的杆或梁所构成的。由不同方向单元组成的结构，它的整体刚度矩阵并不能由局部坐标下的单元刚度矩阵简单地叠加而成，所以必须建立一个统一的整体坐标系。计算时先将单元上的结点力和位移转换到整体坐标系，单元刚度矩阵亦作坐标变换，才可按叠加规则直接相加组成整体刚度矩阵。单元的坐标变换在前面几节中，推导单元刚度矩阵时采用84坐标变换坐标变换85平面杆单元的坐标转换矩阵平面杆单元的坐标转换矩阵86整体坐标系下的单元刚度矩阵整体坐标系下的单元刚度矩阵87平面梁单元的转换矩阵平面梁单元的转换矩阵88平面