黃科4節(jié)點程序

1、四邊形四節(jié)點等參元 matlab計算程序姓名：黃科 學(xué)號：% *變量說明*% NNODE： 單元節(jié)點數(shù) NFORCE： 受力結(jié)點數(shù) % LNODS： 單元定義數(shù)組 COORD： 結(jié)構(gòu)節(jié)點整體坐標數(shù)組 % ASLOAD： 總體荷載向量 DISP： 結(jié)點位移向量% HK： 總體剛度矩陣 Ke： 單元剛度矩陣% TJLOAD：體積力% FIXED：約束信息數(shù)組， (n,1)：約束點號；(n,2)與(n,3)：分別為約束點x % 方向和y方向的約束情況,受約束為1否則為0；n：受約束節(jié)點數(shù)目；% FORCE：結(jié)點力數(shù)組；(i,1)為作用點；(i,2)為x方向的節(jié)點力；(i,3)為y方向% 的節(jié)點力；n

2、：節(jié)點力數(shù)目；% E：彈性模量 NPOIN： 總結(jié)點數(shù)% h：厚度 NELEM： 單元數(shù)% v：泊松比 NVFIX： 約束結(jié)點個數(shù)%讀取初始數(shù)據(jù)clearclcformat short e FP1=fopen(4jd.txt,rt); % 打開數(shù)據(jù)文件%*讀入控制數(shù)據(jù)*NPOIN=fscanf(FP1,%d,1); %總結(jié)點數(shù)NELEM=fscanf(FP1,%d,1); %單元數(shù)NFORCE=fscanf(FP1,%d,1); %結(jié)點力數(shù)NVFIX=fscanf(FP1,%d,1); %受約束結(jié)點數(shù)NNODE=fscanf(FP1,%d,1); %單元結(jié)點數(shù)E=fscanf(FP1,%e,1

3、); %彈性模量v=fscanf(FP1,%f,1); %泊松比h=fscanf(FP1,%f,1); %厚度%*讀取結(jié)構(gòu)信息*LNODS=fscanf(FP1,%f,NNODE,NELEM); % 單元定義：單元結(jié)點編碼（逆時針）；共NELEM行, NFPOIN列COORD=fscanf(FP1,%f,2,NPOIN); % 結(jié)點號x坐標，y坐標(整體坐標下)；共NPOIN行,2列FORCE=fscanf(FP1,%f,3,NFPOIN); % 節(jié)點力：結(jié)點號、X方向力(向右為正)，Y方向力(向上為正)FIXED=fscanf(FP1,%d,3,NVFIX); % 約束信息數(shù)組：受約束節(jié)點號

4、，x，y方向受約束情況（受約束為1）；共NVFIX% 行, 3列%平面應(yīng)力問題求解%形成剛度矩陣%計算剛度矩陣，并對約束條件進行處理HK=zeros(2*NPOIN,2*NPOIN); %張成總剛度矩陣并清零%調(diào)用子程序生成單元剛度矩陣for ie=1:NELEM %ie為單元號Ke= StiffnessMatrix(ie,E,v,COORD,LNODS,h); %調(diào)用單元剛度矩陣a=LNODS(ie,:); %臨時向量，用來記錄當(dāng)前單元的節(jié)點編號for j=1:4 %對行按結(jié)點號進行循環(huán) for k=1:4 %對得按結(jié)點號進行循環(huán)HK(a(j)*2-1):a(j)*2,(a(k)*2-1):

5、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 endend%對荷載向量進行處理ASLOAD=zeros(2*NPOIN,1); %張成總荷載向量并清0for i=1:NFORCE %對每個節(jié)點力進行循環(huán)B1=FORCE(i,1)*2-1;B2=FORCE(i,1)*2; %FORCE（i,1）為作用點ASLOAD(B1)=FORCE(i,2); %FORCE（i,2）為x方向的節(jié)點力ASLOAD(B2)=FORCE(i,3); %FORCE（i,3）為y方向的節(jié)點力en

6、dTJLOAD=ele_LOAD(NELEM,COORD,LNODS,R,h);GTJLOAD(2*NPOIN,1)=0; %張成整體等效節(jié)點荷載向量并清0for ie=1:NELEM %對單元的個數(shù)進行循環(huán) for j=1:NNODE %對單元的結(jié)點個數(shù)進行循環(huán) M1=LNODS(ie,j)*2; %對節(jié)點的第二個位移進行編號 GTJLOAD(M1-1):M1,1)=GTJLOAD(M1-1):M1,1)+TJLOAD(j*2-1:j*2,1);% 集成總體的等效節(jié)點荷載向量 endendASLOAD=ASLOAD+GTJLOAD; %求出結(jié)構(gòu)的節(jié)點荷載向量%總剛、總荷載進行邊界條件處理%將

7、約束信息加入總剛,總荷載（邊界條件處理）for i=1:NVFIX if FIXED(i,2)=1 N1=FIXED(i,1)*2-1; HK(N1,;)=0; %將一約束號處的總剛列向量清0 HK(:,N1)=0; %將一約束號處的總剛行向量清0 HK(N1,N1)=1; %將行列交叉處的元素置為1 ASLOAD(N1)=0;end if FIXED(i,3)=1 N2=FIXED(i,1)*2; HK(N2,;)=0; %將一約束號處的總剛列向量清0 HK(:,N2)=0; %將一約束號處的總剛行向量清0 HK(N2,N2)=1; %將行列交叉處的元素置為1 ASLOAD(N2)=0;en

8、dend%DISP=HKASLOAD %計算節(jié)點位移向量%求解單元應(yīng)力stress=zeros(3,NELEM); %張成應(yīng)力矩陣賦空值strain=zeros(3,NELEM); %張成應(yīng)變矩陣賦空值for ie=1:NELEM %對每個單元進行循環(huán) EDISP(8,1)=0; %張成單元的位移矩陣并清0 d=LNODS(ie,:); for j=1:NNODE %對每個單元節(jié)點進行循環(huán) EDISP(j*2-1:j*2)=DISP(d(j)*2-1:d(j)*2);%從總位移向量中取出當(dāng)前單元的節(jié)點位移endD= (E/(1-v*v)*1 v 0;v 1 0;0 0 (1-v)/2; %彈性

9、矩陣DB= MatrixB( ie, 1,1,COORD,LNODS ); %調(diào)用計算應(yīng)變矩陣B函數(shù),賦值（1，1）處S=D*B; %計算應(yīng)力矩陣Sstress(:,ie)= S*EDISP; %將每個單元的應(yīng)力矩陣集合在一起strain(:,ie)=B*EDISP; %將每個單元的應(yīng)變矩陣集合在一起endstress %輸出應(yīng)力矩陣 strain %輸出應(yīng)變矩陣%子程序function Ke=StiffnessMatrix(ie,E,v,COORD,LNODS,h)% 計算平面應(yīng)變等參數(shù)單元的剛度矩陣 % 說明：用高斯積分法求解平面等參數(shù)單元的剛度矩陣 Ke=zeros(8,8); %張成單

10、元剛度矩陣并清零D= (E/(1-v*v)*1 v 0;v 1 0;0 0 (1-v)/2; %彈性矩陣D% 3 x 3 高斯積分點和權(quán)系數(shù) xi = -0.1483, 0,0.1483; %高斯點w = 0.5555,0.88889,0.5555; %權(quán)系數(shù)for i=1:length(xi) for j=1:length(xi) B = MatrixB( ie, xi(i), xi(j) ,COORD,LNODS) ; J = Jacobi( ie, xi(i), xi(j) ,COORD,LNODS) ; Ke = Ke + w(i)*w(j)*transpose(B)*D*B*det(

11、J) *h; end endreturn%function TJLOAD=ele_LOAD(ie,COORD,LNODS,R,h)% 計算單元等效節(jié)點體積力% 輸入?yún)?shù)：ie單元號;R密度;h厚度% 返回值：TJLOAD結(jié)構(gòu)的等效節(jié)點荷載向量TJLOAD(8,1)=0;x=COORD(LNODS(ie,1:4),1);%取出每個單元結(jié)點在整體坐標系下的x方向的坐標y= COORD(LNODS(ie,1:4),2);%取出每個單元結(jié)點在整體坐標系下的y方向的坐標%利用三次高斯積分集成單元的體積力的等效節(jié)點力矩陣% 3 x 3 高斯積分點和權(quán)系數(shù) xi = -0.1483, 0,0.1483; %

12、取高斯點w = 0.5555,0.88889,0.5555; %權(quán)系數(shù)for i=1:length(xi) for j=1:length(xi)N_kc, N_eta = N_kceta( ie, xi(i), xi(j); J = Jacobi( ie, xi(i),xi(j),COORD,LNODS );N=Shape(xi(i),xi(j);TJLOAD(1:2:7)=TJLOAD(1:2:7)+N*N*x*det(J)*w(i)*w(j)*h;TJLOAD(2:2:8)=TJLOAD(2:2:8)+N*N*y*det(J)*w(i)*w(j)*h;%單元的體積力等效成節(jié)點的等效節(jié)點力e

13、ndendTJLOAD=TJLOAD*R;return%function N=Shape(kc,eta)% 計算形函數(shù)的值% 輸入?yún)?shù)：ie單元號; kc,eta單元內(nèi)局部坐標% 返回值：N形函數(shù)的值N1 = ( 1 - kc ) * ( 1 - eta ) / 4;N2 = ( 1 + kc ) * ( 1 - eta ) / 4 ;N3 = ( 1 + kc ) * ( 1 + eta ) / 4 ;N4 = ( 1 - kc ) * ( 1 + eta ) / 4 ;N = N1 N2 N3 N4; return%function B = MatrixB( ie, kc, eta,COO

14、RD,LNODS ) % 計算單元的應(yīng)變矩陣 B % 輸入?yún)?shù)：ie 單元號; kc,eta 局部坐標% 返回值：B在局部坐標處的應(yīng)變矩陣 BN_kc,N_eta = N_kceta( ie,kc, eta ) ; %調(diào)用形函數(shù)對局部坐標求導(dǎo)的函數(shù)J = Jacobi( ie, kc, eta ,COORD,LNODS) ; %調(diào)用雅克比矩陣函數(shù)J1=J(2,2) -J(1,2);-J(1,2) J(1,1)/J(1,1)*J(2,2)-J(1,2)*J(2,1);B = zeros( 3, 8 ) ; %張成應(yīng)變矩陣并清0for i=1:4 B1=J1(1,1)*N_kc(i)+J1(1,2

15、)*N_eta(i);B2=J1(2,1)*N_kc(i)+J1(2,2)*N_eta(i); B(1:3,(2*i-1):2*i) = B1 0;0 B2;B2 B1; end return %function N_kc, N_eta = N_kceta( ie, kc, eta ) % 計算形函數(shù)對局部坐標的導(dǎo)數(shù) % 輸入?yún)?shù)： ie單元號; kc,eta局部坐標% 返回值： N_kc 在局部坐標處的形函數(shù)對 kc坐標的導(dǎo)數(shù)% N_eta 在局部坐標處的形函數(shù)對 eta 坐標的導(dǎo)數(shù) N_kc = zeros( 1, 4 ) ; N_eta = zeros( 1, 4 ) ; %張成形函數(shù)對

16、局部坐標導(dǎo)數(shù)的矩陣并清0N_kc = zeros( 1, 4 ) ; N_eta = zeros( 1, 4) ; N_kc(1) = -(1-eta)/4; N_eta(1) = -(1-kc)/4 ; N_kc(2) = (1-eta)/4; N_eta(2) = -(1+ kc)/4 ; N_kc(3) = (1+ eta)/4; N_eta(3) = (1+ kc)/4 ; N_kc(4) = -(1+ eta)/4; N_eta(4) = (1-kc)/4 ;return %function J = Jacobi( ie, kc, eta,COORD,LNODS ) % 計算雅克比矩陣 % 輸入?yún)?shù)： ie單元號; kc,eta 局部坐標% 返回值： J在局部坐標(kc,eta)處的雅克比矩陣x = COORD(LNODS(ie,1:4),1) ; %取出每個單元結(jié)點在整體坐標系下的x方向的坐標y = COORD(LNODS(ie,1:4),2) ; %取出每個單元結(jié)點在整體坐標系下的y方向的坐標N_kc,N_eta = N_kceta( ie, kc, eta ) ; %調(diào)用形函數(shù)對局部坐標的導(dǎo)數(shù)x_kc = N_kc