平面梁单元MATLAB有限元程序.docx
《平面梁单元MATLAB有限元程序.docx》由会员分享，可在线阅读，更多相关《平面梁单元MATLAB有限元程序.docx（9页珍藏版）》请在冰豆网上搜索。
平面梁单元MATLAB有限元程序
平面梁单元MATLAB有限元程序
function[]=beam2n（）
clearall
closeall
clear
close
%------------------------BEAM2---------------------------
disp（'========================================'）;disp（'PROGRAMBEAM2'）;disp（'BeamBendingAnalysis'）;disp（'T.R.ChandrupatlaandA.D.Belegundu'）;disp（'========================================'）;
InputData;
Bandwidth;
Stiffness;
ModifyForBC;
BandSolver;
ReactionCalc;
Output;
%------------------------functionInputData---------------------------
function[]=InputData（）;
globalNNNENMNDIMNENNDN
globalNDNLNCHNPRNMPCNBW
globalXNOCFAREAMATSMIS
globalPMNUUMPCBTREACT
globalCNST
globalTITLEFILE1FILE2
globalLINPLOUT
globalNQ
disp（blanks
（1））;
FILE1=input（'InputDataFileName','s'）;LINP=fopen（FILE1,'r'）;
FILE2=input（'OutputDataFileName','s'）;LOUT=fopen（FILE2,'w'）;
DUMMY=fgets（LINP）;
TITLE=fgets（LINP）;
DUMMY=fgets（LINP）;
TMP=str2num（fgets（LINP））;
[NN,NE,NM,NDIM,NEN,NDN]=deal（TMP
（1）,TMP
（2）,TMP（3）,TMP（4）,TMP（5）,TMP（6））;
NQ=NDN*NN;
DUMMY=fgets（LINP）;
TMP=str2num（fgets（LINP））;
[ND,NL,NMPC]=deal（TMP
（1）,TMP
（2）,TMP（3））;
NPR=1;%E
%-----Coordinates-----
DUMMY=fgets（LINP）;
forI=1:
NN
TMP=str2num（fgets（LINP））;
[N,X（N,:
）]=deal（TMP
（1）,TMP（2:
1+NDIM））;end
%-----Connectivity-----
DUMMY=fgets（LINP）;
forI=1:
NE
TMP=str2num（fgets（LINP））;
[N,NOC（N,:
）,MAT（N）,SMI（N）]=...
deal（TMP
（1）,TMP（2:
1+NEN）,TMP（2+NEN）,TMP（3+NEN））;
end
%-----SpecifiedDisplacements-----DUMMY=fgets（LINP）;
forI=1:
ND
TMP=str2num（fgets（LINP））;
[NU（I,:
）,U（I,:
）]=deal（TMP
（1）,TMP
（2））;end
%-----ComponentLoads-----
DUMMY=fgets（LINP）;
F=zeros（NQ,1）;
forI=1:
NL
TMP=str2num（fgets（LINP））;
[N,F（N）]=deal（TMP
（1）,TMP
（2））;end
%-----MaterialProperties-----DUMMY=fgets（LINP）;
forI=1:
NM
TMP=str2num（fgets（LINP））;
[N,PM（N,:
）]=deal（TMP
（1）,TMP（2:
NPR+1））;
end
%-----Multi-pointConstraintsB1*Qi+B2*Qj=B0
ifNMPC>0
DUMMY=fgets（LINP）;
forI=1:
NMPC
TMP=str2num（fgets（LINP））;
[BT（I,1）,MPC（I,1）,BT（I,2）,MPC（I,2）,BT（I,3）]=...
deal（TMP
（1）,TMP
（2）,TMP（3）,TMP（4）,TMP（5））;
end
end
fclose（LINP）;
%------------------------functionBandwidth---------------------------
function[]=Bandwidth（）;
globalNNNENMNDIMNENNDN
globalNDNLNCHNPRNMPCNBWglobalXNOCFAREAMATSMISglobalPMNUUMPCBTSTRESSREACTglobalCNST
globalTITLEFILE1FILE2
globalLINPLOUT
%-----BandwidthEvaluation-----NBW=0;
forN=1:
NE
NABS=NDN*（abs（NOC（N,1）-NOC（N,2））+1）;
if（NBWNBW=NABS;
end
end
forI=1:
NMPC
NABS=abs（MPC（I,1）-MPC（I,2））+1;
if（NBWNBW=NABS;
end
end
disp（sprintf（'Bandwidth=%d',NBW））;
%------------------------functionStiffness---------------------------
function[]=Stiffness（）;
globalNNNENMNDIMNENNDN
globalNDNLNCHNPRNMPCNBWglobalXNOCFAREAMATSMISglobalPMNUUMPCBTREACT
globalCNST
globalNQ
%-----GlobalStiffnessMatrixS=zeros（NQ,NBW）;
forN=1:
NE
disp（sprintf（'FormingStiffnessMatrixofElement%d',N））;
%--------ElementStiffnessandTemperatureLoad-----
N1=NOC（N,1）;
N2=NOC（N,2）;
M=MAT（N）;
EL=abs（X（N1）-X（N2））;
EIL=PM（M,1）*SMI（N）/EL^3;
SE（1,1）=12*EIL;
SE（1,2）=EIL*6*EL;
SE（1,3）=-12*EIL;
SE（1,4）=EIL*6*EL;
SE（2,1）=SE（1,2）;
SE（2,2）=EIL*4*EL*EL;
SE（2,3）=-EIL*6*EL;
SE（2,4）=EIL*2*EL*EL;
SE（3,1）=SE（1,3）;
SE（3,2）=SE（2,3）;
SE（3,3）=EIL*12;
SE（3,4）=-EIL*6*EL;
SE（4,1）=SE（1,4）;
SE（4,2）=SE（2,4）;
SE（4,3）=SE（3,4）;
SE（4,4）=EIL*4*EL*EL;
disp（'....PlacinginGlobalLocations'）;
forII=1:
NEN
NRT=NDN*（NOC（N,II）-1）;
forIT=1:
NDN
NR=NRT+IT;
I=NDN*（II-1）+IT;
forJJ=1:
NEN
NCT=NDN*（NOC（N,JJ）-1）;
forJT=1:
NDN
J=NDN*（JJ-1）+JT;
NC=NCT+JT-NR+1;
if（NC>0）
S（NR,NC）=S（NR,NC）+SE（I,J）;
end
end
end
end
end
end
%------------------------functionModifyForBC---------------------------
function[]=ModifyForBC（）;
globalNNNENMNDIMNENNDN
globalNDNLNCHNPRNMPCNBW
globalXNOCFAREAMATSMIS
globalPMNUUMPCBTREACT
globalCNST
globalNQ
%-----DecidePenaltyParameterCNST-----CNST=0;
forI=1:
NQ
ifCNSTCNST=CNST*10000;
%-----ModifyforBoundaryConditions-----%---DisplacementBC---
forI=1:
ND
N=NU（I）;
S（N,1）=S（N,1）+CNST;
F（N）=F（N）+CNST*U（I）;
end
%---Multi-pointConstraints---
forI=1:
NMPC
I1=MPC（I,1）;
I2=MPC（I,2）;
S（I1,1）=S（I1,1）+CNST*BT（I,1）*BT（I,1）;
S（I2,1）=S（I2,1）+CNST*BT（I,2）*BT（I,2）;
IR=I1;
ifIR>I2;IR=I2;end
IC=abs（I2-I1）+1;
S（IR,IC）=S（IR,IC）+CNST*BT（I,1）*BT（I,2）;
F（I1）=F（I1）+CNST*BT（I,1）*BT（I,3）;
F（I2）=F（I2）+CNST*BT（I,2）*BT（I,3）;end
%------------------------functionBandSolver---------------------------
function[]=BandSolver（）;
globalNNNENMNDIMNENNDN
globalNDNLNCHNPRNMPCNBW
globalXNOCFAREAMATSMIS
globalPMNUUMPCBTREACT
globalCNST
globalNQ
%-----EquationSolvingusingBandSolver-----disp（'SolvingusingBandSolver（bansol.m）'）;
[F]=bansol（NQ,NBW,S,F）;
%------------------------functionReactionCalc---------------------------
function[]=ReactionCalc（）;
globalNNNENMNDIMNENNDN
globalNDNLNCHNPRNMPCNBW
globalXNOCFAREAMATSMIS
globalPMNUUMPCBTREACT
globalCNST
forI=1:
ND
N=NU（I）;
REACT（I）=CNST*（U（I）-F（N））;
end
%------------------------functionOutput---------------------------
function[]=Output（）;
globalNNNENMNDIMNENNDN
globalNDNLNCHNPRNMPCNBW
globalXNOCFAREAMATSMIS
globalPMNUUMPCBTREACT
globalCNST
globalTITLEFILE1FILE2
globalLINPLOUT
disp（sprintf（'OutputforInputDatafromfile%s\n',FILE1））;fprintf（LOUT,'OutputforInputDatafromfile%s\n',FILE1）;
disp（TITLE）;
fprintf（LOUT,'%s\n',TITLE）;
disp（'Node#X-DisplRotation'）;
fprintf（LOUT,'Node#X-DisplRotation\n'）;I=[1:
NN]';
%printamatrix
disp（sprintf（'%4d%15.4E%15.4E\n',[I,F（2*I-1）,F（2*I）]'））;fprintf（LOUT,'%4d%15.4E%15.4E\n',[I,F（2*I-1）,F（2*I）]'）;
%-----ReactionCalculation-----
disp（sprintf（'DOF#Reaction'））;
fprintf（LOUT,'DOF#Reaction\n'）;
forI=1:
ND
N=NU（I）;
R=CNST*（U（I）-F（N））;
disp（sprintf（'%4d%15.4E',N,REACT（I）））;
fprintf（LOUT,'%4d%15.4E\n',N,REACT（I））;end
fclose（LOUT）;
disp（sprintf（'TheResultsareavailableinthetextfile%s',FILE2））;
升级会员 