等参元八节点matlab,四边形八节点等参元matlab程序
懸臂鋼梁，尺寸如圖一所示；v=0.3。h=1，E=2.1e11.
圖一 懸臂鋼梁
圖二 單元劃分與結點編號
Matlab 輸出結果
附錄Ⅰ：
有限元ANSYS分析結果
采用PLANE183單元(四邊形八節點)單元得出的結構Y向最大位移為-0.216E-04。約等于matlab平面四邊形八節點等參元結點Y向最大位移-2.4024E-5。
附錄Ⅱ：
%---------------四邊形八節點等參元 matlab計算程序----------------------------
%——— ———— ———— 主 程 序———— —————
%*******************************************************************%************************************
% 2012年
% 本程序只能處理集中荷載作用下的情況
% 只輸出了節點位移、單元中心點的應力
%*******************************************************************%***************
% 變量說明
% E v h
% 彈性模量 泊松比 厚度
% NPOIN NELEM NVFIX NNODE NFPOIN
% 總結點數 , 單元數, 約束結點個數, 單元節點數 ,受力結點數
% COORD LNODS
% 結構節點整體坐標數組, 單元定義數組,
% FPOIN FORCE FIXED
% 結點力數組， 總體荷載向量, 約束信息數組
% HK DISP
% 總體剛度矩陣,結點位移向量
%******************************
clear all
format short e
FP1=fopen(bjd.txt,rt); %打開數據文件
%%讀入控制數據
E=fscanf(FP1,%f,1); %彈性模量
v=fscanf(FP1,%f,1); % 泊松比
h=fscanf(FP1,%f,1); %厚度
NELEM=fscanf(FP1,%d,1); %單元數
NPOIN=fscanf(FP1,%d,1); % 總結點數
NNODE=fscanf(FP1,%d,1); %單元節點數
NFPOIN=fscanf(FP1,%d,1); %受力結點數
NVFIX=fscanf(FP1,%d,1); %約束結點個數
LNODS=fscanf(FP1,%f,[NNODE,NELEM]); % 單元定義： 單元結點號(逆時針)
COORD=fscanf(FP1,%f,[2,NPOIN]); % 結點號 x,y坐標(整體坐標下)
FPOIN=fscanf(FP1,%f,[3,NFPOIN]);
% 節點力：結點號、X方向力(向右正)，Y方向力(向上正)
FIXED=fscanf(FP1,%d,[3,NVFIX]);
%約束信息數組(n,3) n:受約束節點數目, (n,1):約束點號
%(n,2)與(n,3)分別為約束點x方向和y方向的約束情況,受約束為1否則為0
%*******************************************************************
%*******************************************************************
%========平面應力問題的求解==============
%
%*******************************************************************
%*******************************************************************
%— — — — —— — — —— — — —— — — — — — — —
% 剛度矩陣的生成
%計算剛度矩陣，并對約束條件進行處理
Ke=zeros(2*NNODE,2*NNODE); % 單元剛度矩陣并清零
HK=zeros(2*NPOIN,2*NPOIN); % 張成總剛矩陣并清零
%調用子程序 生成單元剛度矩陣
for m=1:NELEM %m為單元號
Ke=K(E,v,h,...
COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),...
COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),...
COORD(LNODS(m,5),1),COORD(LNODS(m,5),2),...
COORD(LNODS(m,7),1),COORD(LNODS(m,7),2)); %調用單元剛度矩陣
a=LNODS(m,:); %臨時向量,用來記錄當前單元的節點編號
%對總剛度矩陣的處理
for j=1:8
for k=1:8
HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+...
Ke(j*2-1:j*2,k*2-1:k*2);
end
end
end
%—————————————————————————————————
%對荷載向量進行處理
FORCE=zeros(2*NPOIN,1); % 張成總荷載向量并清零
for i=1:NFPOIN
b1=FPOIN(i,1)*2-1;b2=FPOIN(i,1)*2; %FPION(i,1)為作用點
FORCE(b1)=FPOIN(i,2); %FPION(i,2)為x方向的節點力
FORCE(b2)=FPOIN(i,3); %FPION(i,3)為y方向的節點力
end
%—————————————————————————————————
%將約束信息加入總剛,總荷載
for i=1:NVFIX
if FIXED(i,2)==1
c1=2*FIXED(i,1)-1;
HK(c1,:)=0; %將一約束序號處的總剛列向量清0
HK(:,c1)=0; %將一約束序號處的總剛行向量清0
HK(c1,c1)=1; %將行列交叉處的元素置為1
FORCE(c1)=0;
end
if FIXED(i,3)==1
c2=2*FIXED(i,1);
HK(c2,:)=0;
HK(:,c2)=0;
HK(c2,c2)=1;
FORCE(c2)=0;
end
end
%—————————————————————————————————
%===========================================================
%===========================================================
DISP=HK\FORCE %計算節點位移向量
%===========================================================
%===========================================================
%———————————求解單元應力————————————————
stress=zeros(3,NELEM);
for m=1:NELEM
u(1:16)=0;
d=LNODS(m,:); %臨時向量,用來記錄當前單元的節點編號
for i=1:NNODE
u(i*2-1:i*2)=DISP(d(i)*2-1:d(i)*2);
%從總位移向量中取出當前單元的節點位移
end
D=(E/(1-v*v))*[1 v 0;v 1 0;0 0 (1-v)/2];%彈性矩陣
%形成應變矩陣BM
BM=zeros(3,16);
for i=1:NNODE
J=Jacobi(COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),...
COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),...
COORD(LNODS(m,5),1),COORD(LNODS(m,5),2),...
COORD(LNODS(m,7),1),COORD(LNODS(m,7),2),0,0);
[N_s,N_t]=DHS(0,0);
B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);
B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);
BM(1:3,2*i-1:2*i)=[B1i 0;0 B2i;B2i B1i]/det(J);
end
stressm=D*BM*u;
stress(:,m)=stressm;
end
stress %輸出應力
function Ke=K(E,v,h,x1,y1,x3,y3,x5,y5,x7,y7)
%=========單元剛度矩陣===============
% E 彈性模量
% v 泊松比
% h 厚度
% x1,y1,x3,y3,x5,y5,x7,y7 為4個角結點的坐標
%矩陣尺寸為16 x 16
Ke=zeros(16,16);
D=(E/(1-v*v))*[1 v 0;v 1 0;0 0 (1-v)/2];%彈性矩陣
%高斯積分 采用 3 x 3 個積分點 書74頁
W1=5/9;W2=8/9;W3=5/9; %加權系數
W=[W1 W2 W3];
r=15^(1/2)/5;
x=[-r 0 r];%積分點
for i=1:3
for j=1:3
B=eleB(x1,y1,x3,y3,x5,y5,x7,y7,x(i),x(j));
J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,x(i),x(j));
Ke=Ke+W(i)*W(j)*B*D*B*det(J)*h;
end
end
function B=eleB(x1,y1,x3,y3,x5,y5,x7,y7,s,t)
%調用導函數
[N_s,N_t]=DHS(s,t);
%求Jacobi矩陣
J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t);
%求應變矩陣B
B=zeros(3,16);
for i=1:8
B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);
B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);
B(1:3,2*i-1:2*i)=[B1i 0;0 B2i;B2i B1i];
end
B=B/det(J);
function J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t)
%-------Jacobi-----------
%單元坐標
%2,4,6,8點的坐標
x2=(x1+x3)/2;y2=(y1+y3)/2;
x4=(x3+x5)/2;y4=(y3+y5)/2;
x6=(x5+x7)/2;y6=(y5+y7)/2;
x8=(x7+x1)/2;y8=(y7+y1)/2;
x=[x1 x2 x3 x4 x5 x6 x7 x8];
y=[y1 y2 y3 y4 y5 y6 y7 y8];
%%調用形函數對局部坐標的導數
[N_s,N_t]=DHS(s,t);
%求Jacobi矩陣的行列式的值
x_s=0;y_s=0;
x_t=0;y_t=0;
for i=1:8
x_s=x_s+N_s(i)*x(i);y_s=y_s+N_s(i)*y(i);
x_t=x_t+N_t(i)*x(i);y_t=y_t+N_t(i)*y(i);
end
J=[x_s y_s;x_t y_t];
function N=shape(s,t)
%ξ,η
N(1) = (1-s)*(1-t)*(-s-t-1)/4;
N(3) = (1+s)*(1-t)*(s-t-1)/4;
N(5) = (1+s)*(1+t)*(s+t-1)/4;
N(7) = (1-s)*(1+t)*(-s+t-1)/4;
N(2) = (1-t)*(1+s)*(1-s)/2;
N(4) = (1+s)*(1+t)*(1-t)/2;
N(6) = (1+t)*(1+s)*(1-s)/2;
N(8) = (1-s)*(1+t)*(1-t)/2;
function [N_s,N_t]=DHS(s,t)
%形函數求導
%ξ,η
N_s(1)=-1/4*(1-t)*(-s-t-1)-1/4*(1-s)*(1-t);
N_s(3)=1/4*(1-t)*(s-t-1)+1/4*(1+s)*(1-t);
N_s(5)=1/4*(1+t)*(s+t-1)+1/4*(1+s)*(1+t);
N_s(7)=-1/4*(1+t)*(-s+t-1)-1/4*(1-s)*(1+t);
N_s(2)=1/2*(1-s)*(1-t)-1/2*(1+s)*(1-t);
N_s(4)=1/2*(1+t)*(1-t);
N_s(6)=1/2*(1-s)*(1+t)-1/2*(1+s)*(1+t);
N_s(8)=-1/2*(1+t)*(1-t);
N_t(1)=-1/4*(1-s)*(-s-t-1)-1/4*(1-s)*(1-t);
N_t(3)=-1/4*(1+s)*(s-t-1)-1/4*(1+s)*(1-t);
N_t(5)=1/4*(1+s)*(s+t-1)+1/4*(1+s)*(1+t);
N_t(7)=1/4*(1-s)*(-s+t-1)+1/4*(1-s)*(1+t);
N_t(2)=-1/2*(1+s)*(1-s);
N_t(4)=1/2*(1+s)*(1-t)-1/2*(1+s)*(1+t);
N_t(6)=1/2*(1+s)*(1-s);
N_t(8)=1/2*(1-s)*(1-t)-1/2*(1-s)*(1+t);
bjd.txt 文件數據
2.1E11 0.3 1 5 28 8 1 3
1 2 3 13 20 19 18 12
3 4 5 14 22 21 20 13
5 6 7 15 24 23 22 14
7 8 9 16 26 25 24 15
9 10 11 17 28 27 26 16
0.0 0.0
0.5 0.0
1.0 0.0
1.5 0.0
2.0 0.0
2.5 0.0
3.0 0.0
3.5 0.0
4.0 0.0
4.5 0.0
5.0 0.0
0.0 0.5
1.0 0.5
2.0 0.5
3.0 0.5
4.0 0.5
5.0 0.5
0.0 1.0
0.5 1.0
1.0 1.0
1.5 1.0
2.0 1.0
2.5 1.0
3.0 1.0
3.5 1.0
4.0 1.0
4.5 1.0
5.0 1.0
17 0 -10000
1 1 1
12 1 1
18 1 1
展開閱讀全文
總結
以上是生活随笔為你收集整理的等参元八节点matlab,四边形八节点等参元matlab程序的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： centos7.0安装php,cento
- 下一篇： java 注释 连接,java – 如何