机械工程用有限元法学习笔记（二）

梁单元的有限元法

梁类结构分为平面梁类结构和空间梁类结构。

平面梁类结构，是指结构物本身是一个直梁，所受的外力垂直于直梁的纵轴线，并且所有的外力或外力矩都作用在通过直梁纵轴线的一个平面内的结构。对于平面梁单元，每个节点只考虑挠度和转角两个自由度。

空间梁类结构，是指结构物本身是一个空间物（包括平面刚梁），所受的外力或外力矩是一组空间力系。对于空间梁单元，每个节点一般要考虑包括沿三个坐标轴的移动和绕三个坐标轴的转动在内的六个自由度。

一、 平面梁类结构的有限元法

平面梁单元的节点力和节点位移之间的关系

 式2-1

或简写为   

  式2-2

式中:

{F}—平面梁单元的载荷列阵

 式2-2a

[K]—平面梁单元的刚阵

 式2-2b

{q}—平面梁单元的节点位移列阵

 式2-2c

（式2-2b）的一般形式

[K]刚阵的物理意义是：[K]中第i行各元素就分别代表，在i号自由度上产生单位位移v_i、θ_zi、v_j、θ_zj时，对其他各个自由度（包括i号自由度）引起的节点力F_yi、M_zi、F_yj、M_zj的贡献。

式2-2表明，在节点力｛F｝的作用下，单元变形达到新的平衡状态后，节点力｛F｝与单元节点的弹性内力相平衡。

求解平面梁类结构的步骤：

⑴单元剖分。按单元剖分原则，把梁剖分成若干个梁单元；

⑵单元分析。写出单元i的刚阵[K]_i，根据式2-2写出各单元节点力和节点位移的关系

；

⑶单元综合。把各单元组合起来，形成原结构的整体，求出结构的总刚阵；

⑷利用边界条件，减缩方程组；

⑸求解变形；

⑹求解支反力；

⑺求解应力

题5：设有图3所示的某轴类结构，其材料的弹性模量为E，两段抗弯惯性矩的比值I_z2/I_z1=0.5，试用有限元求解⑴A点的挠度、转角以及B点的转角，⑵1、2、3处的支反力。

解：1.单元剖分

图3

根据单元剖分原则，把结构剖分为两个单元：①和②并标出各节点和各自由度的序号，如图3b所示。图中轴上的数字1，2，3是节点号，轴下方的数字1，2，3，4，5，6是节点自由度（位移）的序号（它们分别代表v₁,θ_Z1，v₂,θ_Z2，v₃,θ_Z3）

2. 单元分析

单元⑴的刚阵为

单元⑵的刚阵为

3.单元综合

按节点位移序号组成结构的总刚阵

总外力列阵

总位移列阵

4.利用边界条件，减缩方程组

由于v₁=v₃=0,则得减缩刚阵为

减缩后的刚阵的逆阵为

对应的减缩后的节点位移列阵和力列阵为

5.求解变形

由

，可解得各节点位移为

6.求解支反力

由

，可解得各节点支反力为

7. Matlab解

syms E I L P;

K1=E*I/L^3*[12 6*L -12 6*L;6*L 4*L^2 -6*L 2*L^2;-12 -6*L 12 -6*L;6*L 2*L^2 -6*L 4*L^2]

K2=0.5*K1

A=zeros(4,2)

B=zeros(2,6)

K10=[K1 A; B]

K20=[B; A K2]

K=K10+K20

A1=K(2:4,[2,3,4])

A2=K(2:4,[6])

A3=K(6,[2,3,4])

A4=K(6,[6])

Kr=[A1 A2;A3 A4]

Fr=[0;-P;0;0]

Krinv=inv(Kr)

qr=Krinv*Fr

q=[0;qr(1);qr(2);qr(3);0;qr(4)]

F=K*q

>> qr =

  -1/3/E/I*L^2*P

  -1/4/E/I*L^3*P

 -1/12/E/I*L^2*P

  5/12/E/I*L^2*P

>> F =

 1/2*P

  0

    -P

  0

 1/2*P

  0

题6：图4所示的某轴类结构，其材料的弹性模量为E，抗弯惯性矩的比值I₂/I₁=5。，试用有限元求解A点的挠度和转角以及B点的转角。

解：首先，将（图4-a）作用在变截面梁中点处的集中力转移到变截面处，图4-b所示。此时有

其次，把该轴剖分为三个单元：①、②和③并标出各节点和各自由度的序号，如图4-b所示。它们的刚阵分别为

单元⑴的刚阵为

单元⑵的刚阵为

单元⑶的刚阵为

再次，把各单元刚阵按节点序号对号入座，得结构总刚阵

总外力列阵

总位移列阵

由于v₁=v₄=0,则得减缩刚阵为

对应的减缩后的节点位移列阵和力列阵为

最后，将

值代入

，得

解方程组，得各节点位移为

二、 坐标变换

1.   坐标变换的目的

用单元的局部坐标系Oxyz定义的位移和力转换成用整个结构系统的统一坐标系──Oxyz来定义。

2.   坐标系的方向余弦

l₁＝cos（ _xx）＝λ _xx；l₂＝cos（ _yx）＝λ _yx；l₃＝cos（ _zx）＝λ _zx

m₁＝cos（ _xy）＝λ _xy；m₂＝cos（ _yy）＝λ _yy；m₃＝cos（ _zy）＝λ _zy

n₁＝cos（ _xz）＝λ _xz；n₂＝cos（ _yz）＝λ _yz；n₃＝cos（ _zz）＝λ _zz

λ _pq是统一坐标系的_q轴至局部坐标系的p轴的夹角的余弦。

正负号由右手螺旋规则确定。

3.   坐标系的变换

（1）  对于二维坐标系的变换

二维坐标系的变换关系

｛q｝＝［λ］｛_q｝ 式2-3a

（2）  对于三维坐标系的变换

三维坐标系的变换关系

｛q｝＝［λ］｛_q｝  式2-3b

4.   节点位移的坐标变换

当某一个节点i具有三个挠度u_i，v_i，w_i和三个转角θ_xi，θ_yi，θ_zi，这六个自由度时，有

  式2-4

[Λ]—节点坐标变换矩阵

    式2-5

对于二维

5.   单元节点位移的坐标变换

  式2-6

[T]—单元坐标变换矩阵

6.   单元节点外力的坐标变换

  式2-7

7.   单元刚阵外力的坐标变换

  式2-8

8.   空间坐标系的二次变换

局部坐标系Oxyz与统一坐标系──Oxyz之间的坐标系方向余弦矩阵

   

三、 空间梁类结构的有限元法

1.   拉—压变形的单元梁刚阵

  1    7

2.   扭转变形的单元梁刚阵

3.   弯曲变形的单元梁刚阵

4.   单元刚阵的组合

对于空间梁单元，每个节点有六个自由度，其矩阵表示形式为6×6阶的子方阵。而每个单元共12个自由度，其矩阵表示形式为12×12阶的方阵。

5.   单元刚阵的坐标变换

由式2-9，空间梁单元刚阵划分为4个6×6阶的子方阵。即

又由式2-9，得到单元刚阵的坐标变换式

由此可见，单元刚阵的坐标变换可以按子矩阵分块进行，使得运算简化。

题7：图5所示为一个由三根直杆组成的平面刚架，设各杆的尺寸和几何特性相同；长度l＝5m截面积A＝0.1m²，抗弯惯性矩Iz＝8×10^-4m⁴，弹性模量E＝20×10¹⁰N/m²。试用有限元求解⑴点2处的挠度、转角；⑵1、2、3处的支反力。

解：1.单元剖分

图5  平面刚架

根据单元剖分原则，把结构剖分为三个单元：①、②和③；四个节点：1、2、3、4；12个自由度：1、2、3、4、5、6、7、8、9、10、11、12（它们分别代表u₁，v₁，θ_Z1；u₂，v₂，θ_Z2；u₃，v₃，θ_Z3；u₄，v₄，θ_Z1）

2. 单元分析

单元⑴的刚阵为

单元⑵的刚阵为

单元⑶的刚阵为

把各局部坐标系定义的单元刚阵转换为以统一坐标系定义的单元刚阵。

由于单元⑴的局部坐标系和统一坐标系是一致的，即

考虑单元⑵时，假设它的局部坐标系和统一坐标系相差-90°，如图5所示。把cosψ＝cos(-90°)＝0，sinψ＝sin(-90°)＝-1，代入节点坐标变换式，对于平面刚架，则得每个节点具有三个自由度(u，v，θ_z)的节点坐标变换矩阵

 

单元坐标变换矩阵为

 

根据

，则有

考虑单元⑶时，假设它的局部坐标系和统一坐标系相差30°，如图5所示。把cosψ＝cos(30°)＝，sinψ＝sin(30°)＝，代入节点坐标变换式，对于平面刚架，则得每个节点具有三个自由度(u，v，θ_z)的节点坐标变换矩阵

 

单元坐标变换矩阵为

 

根据

，则有

3.单元综合

按节点自由度序号对号入座的方法，将相对于统一坐标系下各单元刚阵

进行组集，则得结构的总刚阵

为

结构的总外力列阵F和总位移列阵q为

F=｛0，0，0，0，0，100，0，0，0，0，0，0｝^T

q=｛u₁，v₁，θ_z1，u₂，v₂，θ_z2，u₃，v₃，θ_z3，u₄，v₄，θ_z4｝^T

4.利用边界条件，减缩方程组

由于节点1，3，4是刚性固定的。因此，u₁=v₁=θ_z1=u₃=v₃=θ_z3= u₄=v₄=θ_z4=0,将它们作为边界条件，对总刚阵K、总外力列阵F和总位移列阵q进行减缩处理，可得减缩刚阵为

×10⁷

减缩后的刚阵的逆阵为

对应的减缩后的节点位移列阵和力列阵为

5.求解变形

由

，可解得各节点位移为

 

6.求解支反力

由

，可解得各节点支反力为

F_x1=3.3987 N；F_y1=9.9933 N；M_z1=16.648 N·m

F_x2=0 N；F_y2=0 N；M_z2=100 N·m

F_x3=－9.9889 N；F_y3=－2.2328 N；M_z3=16.637 N·m

F_x4=6.5901 N；F_y4=－7.7605 N；M_z4=16.705 N·m

ΣX=3.3987+0-9.9889+6.5901=0，平衡。

ΣY=9.9933+0-2.2328-7.7605=0，平衡。

7. Matlab解

clc

E=20e10;A=0.1;I=8e-4;L=5;

format short g;

K1=[E*A/L 0 0 -E*A/L 0 0; 0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2; 0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L;

  -E*A/L 0 0 E*A/L 0 0;0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2;0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L];

 K2=[E*A/L 0 0 -E*A/L 0 0; 0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2; 0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L;

  -E*A/L 0 0 E*A/L 0 0;0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2;0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L];

 K3=[E*A/L 0 0 -E*A/L 0 0; 0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2; 0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L;

  -E*A/L 0 0 E*A/L 0 0;0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2;0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L];

K12=K1;

A2=[cos(-pi/2) sin(-pi/2) 0;-sin(-pi/2) cos(-pi/2) 0; 0 0 1];

T2=[A2 zeros(3); zeros(3) A2];

K23=T2'*K2*T2;

A3=[cos(pi/6) sin(pi/6) 0;-sin(pi/6) cos(pi/6) 0; 0 0 1];

T3=[A3 zeros(3); zeros(3) A3];

K24=T3'*K3*T3;

%单元⑴1，2，3；4，5，6

K11=[K12 zeros(6,6);zeros(6,12)];

%单元⑵4，5，6；7，8，9

K22=[zeros(3,12);zeros(6,3) K23 zeros(6,3);zeros(3,12)];

%单元⑶4，5，6；10，11，12

K33=[zeros(3,12);zeros(3,3) K24(1:3,[1:3]) zeros(3,3) K24(1:3,[4:6]);zeros(3,12); zeros(3,3) K24(4:6,[1:3]) zeros(3,3) K24(4:6,[4:6])];

K=K11+K22+K33

Kr=K(4:6,[4:6])

% F=[0 0 0 0 0 100 0 0 0 0 0 0]';

 % q=[u1 v1 angle1 u2 v2 angle2 u3 v3 angle3 u4 v4 angle4]';

 Fr=[0;0;100];

 Krinv=inv(Kr);

 qr=Krinv*Fr

q=[0;0;0;qr(1);qr(2);qr(3);0;0;0;0;0;0]

 F=K*q

>>Kr =

  7.0192e+009  1.7254e+009    1.92e+007

  1.7254e+009  5.0269e+009 -5.1446e+006

    1.92e+007 -5.1446e+006    3.84e+008

>>qr =

 -8.4968e-010

  5.5821e-010

  2.6047e-007

>>F =

    3.3987

    9.9933

    16.648

   0

   0

    100

   -9.9889

   -2.2328

    16.637

    6.5901

   -7.7605

    16.705

将E、A、I、L等数值代入上式，有

弹性力学平面问题的有限元法

一、    平面问题的基本方程

平面问题是指弹性体内一点的应力、应变或位移状态只与两个坐标方向的变量有关的二维问题。

1.  平面应力和平面应变

应力分量：

 式3-1

应变分量：

   式3-2

应力分量与应变分量间的关系：

   式3-3

或

   式3-4

式中：[D]—弹性矩阵

对于平面应力问题

   式3-5a

对于平面应变问题（E换成，μ换成）

   式3-5b

2.  平面问题的基本方程

平面问题的总位能表达式

当δΠ=0，可以得到用位移u和v表示的基本方程

   式3-6

如采用应力函数

，可得到平面问题的基本微分方程（双调和方程）

   式3-7

用有限元法求解平面问题的思路：

①   剖分和插值

把整个平面区域S用三角形板单元或矩形板单元等进行剖分并在单元内进行位移函数（形状函数）的插值。

   式3-8

②   单元分析

把形状函数代入位能泛函式Π_i，并按单元进行计算。

   式3-9

③   单元组集

把各单元重新组集起来。

   式3-10

说明：{q}是单元各节点的位移列阵；[K]是单元刚阵；{F}是单元的广义载荷列阵；q是整个区域上各单元节点的位移总和的列阵；K是总刚阵；F是整个区域上的广义载荷列阵，S是单元总数。

二、    位移函数

①   设定的位移函数是泛函的极限条件，即控制方程的近似解。

②   选择位移函数的阶次应考虑下列因素：

ⅰ、满足完备性和协调性。一般采用一个由低阶算起完全的多项式表示位移函数。如u=a₁+a₂x+a₃y+a₄x²+a₅xy+a₆y²……

ⅱ、对称性即该多项式位移函数与局部坐标系的方位无关。

ⅲ、多项式的项数与节点自由度相等。

ⅳ、收敛性

③   设定位移函数时应符合

ⅰ、在单元内部和边界上（包括节点处）处处都能满足力的平衡条件和变形协调条件。

ⅱ、在单元内部要求应变或应力最少应是常值（或线性变化的）。

ⅲ、包含有代表刚体运动的项。

一个单元内各点的位移实际上包含着两部分，一是单元本身变形引起的部分；二是刚体位移部分。位移函数必须能反映这两种位移。

三、    三角形单元分析

1.  位移函数{d}

图3-1 三角形板单元

设三角形板单元单元内某一点的位移函数为

代数形式：

    式3-11a

矩阵形式：

    式3-11b

那么单元的三个节点i、j、k可以写成

    式3-12

由式3-12解得，

    式3-13

式中：

为了不使A为负值，i、j、k的顺序必须是逆时针方向，如图3-1。

将

代人

中，整理得

    式3-14

式中：[N]形状函数

    式3-15

2.  应变{ε}

由

，可得

    式3-16

式中：[B] 应变矩阵

    式3-17

3.  应力{σ}

根据虎克定律，得

    式3-18

式中：[R] 应力矩阵

对于平面应力问题

 式3-19a

对于平面应变问题

 式3-19b

4.  刚阵[K]

根据虚位移原理，可推得

    式3-20

式中：[K] 刚度矩阵

    式3-21

式中：

对于平面应力问题

 

对于平面应变问题

 

四、    矩形单元分析

1.  位移函数{d}

矩形板单元

2.  应变{ε}

3.  应力{σ}

4.  刚阵[K]

5.