有限元基础教程（曾攀）算例4.6.3（2）Matlab程序实现

---

算例4.6.3（2）

如图所示为一个平面4节点四边形等参元，试采用4点Gauss积分计算该单元的刚度矩阵。材料的弹性模量为E=30×106MPa，泊松比为μ=0.3，厚度为t=0.1m。

解题步骤

1. 输入单元信息

>>X0=[1;2;2.25;1.25];
>>Y0=[0;0;1.5;1];
>>E=3E13;nv=0.3;t0=0.1;

2. 计算雅可比矩阵及其相关参量

>> [N_cxi,N_eta,JA,dt_JA,inv_JA]=Jacobian(X0,Y0);

输出如下：
N_cxi =[eta/4 - 1/4, 1/4 - eta/4, eta/4 + 1/4, - eta/4 - 1/4]
N_eta =[cxi/4 - 1/4, - cxi/4 - 1/4, cxi/4 + 1/4, 1/4 - cxi/4]
JA =[1/2, eta/8 + 1/8]
[1/8, cxi/8 + 5/8]
dt_JA =cxi/16 - eta/64 + 19/64
inv_JA =[(8*(cxi + 5))/(4cxi - eta + 19), -(8(eta + 1))/(4cxi - eta + 19)]
[-8/(4cxi - eta + 19), 32/(4*cxi - eta + 19)]

3. 计算弹性系数矩阵

>>D=Quad2D4Node_D(E,nv,t0,X0(1),Y0(1),X0(2),Y0(2),X0(3),Y0(3),X0(4),Y0(4),1)

输出如下：
D =
1.0e+13 *
3.2967 0.9890 0
0.9890 3.2967 0
0 0 1.1538

4. 求解单元刚度矩阵

>>KE=Quad2D4Node_KE(D,X0,Y0,t0)

输出如下：

函数文件

主函数（main）

clc;clear;
X0=[1;2;2.25;1.25];
Y0=[0;0;1.5;1];
E=3E13;nv=0.3;t0=0.1;
[N_cxi,N_eta,JA,dt_JA,inv_JA]=Jacobian(X0,Y0);
D=Quad2D4Node_D(E,nv,t0,X0(1),Y0(1),X0(2),Y0(2),X0(3),Y0(3),X0(4),Y0(4),1)
KE=Quad2D4Node_KE(D,X0,Y0,t0)

函数文件

1. Jacobian.m

function [N_cxi,N_eta,JA,dt_JA,inv_JA]=Jacobian(X0,Y0)
%根据四个结点坐标，推算雅可比矩阵、雅可比行列式、雅可比逆阵
%输入参数：4个结点坐标值,列向量
%返回值：带有符号的表达式——形函数对局部坐标导数、雅可比矩阵、雅可比行列式、雅可比逆阵
syms cxi eta
%定义符号参量
%形函数对无量纲局部坐标偏导数
N_cxi=[-(1-eta),(1-eta),(1+eta),-(1+eta)]/4
N_eta=[-(1-cxi),-(1+cxi),(1+cxi),(1-cxi)]/4
%计算雅可比矩阵
JA=[N_cxi*X0,N_cxi*Y0;N_eta*X0,N_eta*Y0]
%计算雅可比行列式
dt_JA=det(JA)
%计算雅可比逆阵
inv_JA=inv(JA)
return 

2. Quad2D4Node_D.m

function D=Quad2D4Node_D(E,NU,~,xi,yi,xj,yj,xm,ym,xp,yp,ID)
%该函数计算单元的刚度矩阵
%输入弹性模量E、泊松比NU和厚度h
%输入4个节点i，j，m，p的坐标xi，yi，xj，yj，xm，ym，xp，yp
%输入平面问题性质指示参数ID（1为平面应力，2为平面应变）
%输出单元弹性矩阵D（8*8）
%-----------------------------
syms s t;
a=(yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s))/4;
b=(yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t))/4;
c=(xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t))/4;
d=(xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s))/4;
B1=[a*(t-1)/4-b*(s-1)/4 0;0 c*(s-1)/4-d*(t-1)/4;c*(s-1)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4];
B2=[a*(1-t)/4-b*(-1-s)/4 0;0 c*(-1-s)/4-d*(1-t)/4;c*(-1-s)/4-d*(1-t)/4 a*(1-t)/4-b*(-1-s)/4];
B3=[a*(t+1)/4-b*(s+1)/4 0;0 c*(s+1)/4-d*(t+1)/4;c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4];
B4=[a*(-1-t)/4-b*(1-s)/4 0;0 c*(1-s)/4-d*(-1-t)/4;c*(1-s)/4-d*(-1-t)/4 a*(-1-t)/4-b*(1-s)/4];
Bfirst=[B1 B2 B3 B4];
Jfirst=[0 1-t t-s s-1;t-1 0 s+1 -s-t;s-t -s-1 0 t+1;1-s s+t -t-1 0];
J=[xi xj xm xp]*Jfirst*[yi;yj;ym;yp]/8;
B=Bfirst/J;
if ID==1
    D=(E/(1-NU*NU))*[1 NU 0;NU 1 0;0 0 (1-NU)/2];
elseif ID==2
    D=(E/(1+NU)/(1-2*NU))*[1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2];
end

3. Quad2D4Node_B.m

function [B,dt]=Quad2D4Node_B(cxi,eta,X0,Y0)
%定义符号参量
%计算形函数对局部坐标的导数
N_cxi=[-(1-eta),(1-eta),(1+eta),-(1+eta)]/4;
%形函数对无量纲局部坐标偏导数
N_eta=[-(1-cxi),-(1+cxi),(1+cxi),(1-cxi)]/4;
%计算雅可比矩阵
JA=[N_cxi*X0,N_cxi*Y0;N_eta*X0,N_eta*Y0];
%计算行列式
dt=det(JA);
%计算H
H=[JA(2,2),-JA(1,2),0,0;
    0,0,-JA(2,1),JA(1,1);
    -JA(2,1),JA(1,1),JA(2,2),-JA(2,1)]/dt;
%计算Q
for i=1:1:4
    Q(1:4,(2*i-1):2*i) = [ N_cxi(i),0;
                            N_eta(i),0;
                            0,N_cxi(i);
                            0,N_eta(i)];
end
B=H*Q;
return

4. Quad2D4Node_KE.m

function KE=Quad2D4Node_KE(D,X0,Y0,t0)
%t0为厚度
KE=zeros(8);
gx=[-0.577350269189626,0.577350269189626];
%2*2高斯积分点和权系数
w=[1,1];
%将gs、w输入不同阶次数据，适用相应阶次高斯积分
n0=length(w);
%确定积分阶数，可根据w的元素个数，自行调整积分阶次
for i=1:n0
    for j=1:n0
        %调用函数计算几何矩阵、雅可比行列式
        [B,dt]=Quad2D4Node_B(gx(i),gx(j),X0,Y0);
        KE=KE+t0*transpose(B)*D*B*dt*w(i)*w(j);
        %一个高斯点的单刚
    end
end
return

参考文献

有限元基础教程（曾攀）

总结

这里对文章进行总结：
书中关于雅可比行列式的计算似乎有误，故计算结果与书中相比有出入，如有错误，还请指正！