Quad2D4Node有限元基础教程4.3.2曾攀

最新推荐文章于 2025-07-31 22:28:05 发布

及时搬运学习

最新推荐文章于 2025-07-31 22:28:05 发布

阅读量824

点赞数 8

CC 4.0 BY-SA版权

文章标签：算法人工智能

本文链接：https://blog.csdn.net/qq_40299971/article/details/142171250

% Quad2D4Node_Stiffness(E,NU,h,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 为平面应变)，输出单元刚度矩阵 k(8X8)。
% Quad2D4Node_Assembly(KK,k,i,j,m,p) 
% 该函数进行单元刚度矩阵的组装，输入单元刚度矩阵 k，单元的节点编号 i、j、m、p，输出整体刚度
% 矩阵 KK。
% Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID) 
% 该函数计算单元的应力，输入弹性模量 E，泊松比 NU，厚度 h，4 个节点 i、j、m、p 的坐标
% xi,yi,xj,yj,xm,ym,xp,yp,，平面问题性质指示参数 ID(1 为平面应力，2 为平面应变)，单元的位移列阵 u(8X1)，
% 输出单元的应力 stress(3X1)，由于它为常应力单元，则单元的应力分量为 Sx,Sy,Sxy

clc
clear all;
 %---物理参数------------------
E=1;              %弹性模量,单位Pa
h=1;              %单元厚度,单位m
NU=0.25;             %泊松比
%-------------------------------
x1=0;
y1=0;
x2=1;
y2=0;
x3=1;
y3=1;
x4=0;
y4=1;
Quad2D4Node=Quad2D4Node_Func
[K,Bfirst,D]= Quad2D4Node.Stiffness(E,NU,h,x1,y1,x2,y2,x3,y3,x4,y4,1) 
k=K([3 4 5 6 8],[3 4 5 6 8])
p=[-1; 0; 1; 0 ;0]  %5*1  力边界条件
u=k\p   %由位移边界条件解
q=[0,0,u(1),u(2),u(3),u(4),0,u(5)]' % 全部节点位移
P=K*q  %节点力
fp=[0 0 p(1) p(2) p(3) p(4) 0 p(5)]'
F=P-fp  %支反力
% 单元位移

%单元应变
epusilon=Bfirst*q
%单元应力
sigama1=D*Bfirst*q

%%%%%%%%%%%% Quad2D4Node %%% begin %%%%%%%%%%%%% 
function Quad2D4Node=Quad2D4Node_Func
Quad2D4Node.Stiffness=@Quad2D4Node_Stiffness;
Quad2D4Node.Assembly=@Quad2D4Node_Assembly;
Quad2D4Node.Stress=@Quad2D4Node_Stress;
end
%%
function [k,Bfirst,D]= Quad2D4Node_Stiffness(E,NU,h,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 为平面应变) 
%输出单元刚度矩阵 k(8X8) 
%--------------------------------------------------------------- 
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