平面梁单元MATLAB有限元程序.docx

1、平面梁单元MATLAB有限元程序平面梁单元MATLAB有限元程序function =beam2n() clear all close all clear close %- BEAM2 - disp(=); disp( PROGRAM BEAM2 ); disp( Beam Bending Analysis ); disp( T.R.Chandrupatla and A.D.Belegundu ); disp(=); InputData; Bandwidth; Stiffness; ModifyForBC; BandSolver; ReactionCalc; Output; %- functio

2、n InputData - function = InputData(); global NN NE NM NDIM NEN NDN global ND NL NCH NPR NMPC NBW global X NOC F AREA MAT SMI S global PM NU U MPC BT REACT global CNST global TITLE FILE1 FILE2 global LINP LOUT global NQ disp(blanks(1); FILE1 = input(Input Data File Name ,s); LINP = fopen(FILE1,r); FI

3、LE2 = input(Output Data File Name ,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,

4、NMPC= deal(TMP(1),TMP(2),TMP(3); NPR=1; % E %- Coordinates - DUMMY = fgets(LINP); for I=1:NN TMP = str2num(fgets(LINP); N, X(N,:)=deal(TMP(1),TMP(2:1+NDIM); end %- Connectivity - DUMMY = fgets(LINP); for I=1:NE TMP = str2num(fgets(LINP); N,NOC(N,:), MAT(N), SMI(N) = . deal(TMP(1),TMP(2:1+NEN), TMP(2

5、+NEN), TMP(3+NEN); end %- Specified Displacements - DUMMY = fgets(LINP); for I=1:ND TMP = str2num(fgets(LINP); NU(I,:),U(I,:) = deal(TMP(1), TMP(2); end %- Component Loads - DUMMY = fgets(LINP); F = zeros(NQ,1); for I=1:NL TMP = str2num(fgets(LINP); N,F(N)=deal(TMP(1),TMP(2); end %- Material Propert

6、ies - DUMMY = fgets(LINP); for I=1:NM TMP = str2num(fgets(LINP); N, PM(N,:) = deal(TMP(1), TMP(2:NPR+1); end %- Multi-point Constraints B1*Qi+B2*Qj=B0 if NMPC 0 DUMMY = fgets(LINP); for I=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),TM

7、P(4),TMP(5); end end fclose(LINP); %- function Bandwidth - function =Bandwidth(); global NN NE NM NDIM NEN NDN global ND NL NCH NPR NMPC NBW global X NOC F AREA MAT SMI S global PM NU U MPC BT STRESS REACT global CNST global TITLE FILE1 FILE2 global LINP LOUT %- Bandwidth Evaluation - NBW = 0; for N

8、=1:NE NABS = NDN*(abs(NOC(N, 1) - NOC(N, 2) + 1); if (NBW NABS) NBW = NABS; end end for I=1:NMPC NABS = abs(MPC(I, 1) - MPC(I, 2) + 1; if (NBW 0) S(NR, NC) = S(NR, NC) + SE(I, J); end end end end end end %- function ModifyForBC - function =ModifyForBC(); global NN NE NM NDIM NEN NDN global ND NL NCH

9、 NPR NMPC NBW global X NOC F AREA MAT SMI S global PM NU U MPC BT REACT global CNST global NQ %- Decide Penalty Parameter CNST - CNST = 0; for I = 1:NQ if CNST 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

10、) = F(I2) + CNST * BT(I, 2) * BT(I, 3); end %- function BandSolver - function =BandSolver(); global NN NE NM NDIM NEN NDN global ND NL NCH NPR NMPC NBW global X NOC F AREA MAT SMI S global PM NU U MPC BT REACT global CNST global NQ %- Equation Solving using Band Solver - disp(Solving using Band Solv

11、er(bansol.m); F = bansol(NQ,NBW,S,F); %- function ReactionCalc - function =ReactionCalc(); global NN NE NM NDIM NEN NDN global ND NL NCH NPR NMPC NBW global X NOC F AREA MAT SMI S global PM NU U MPC BT REACT global CNST for I = 1:ND N = NU(I); REACT(I) = CNST * (U(I) - F(N); end %- function Output -

12、 function =Output(); global NN NE NM NDIM NEN NDN global ND NL NCH NPR NMPC NBW global X NOC F AREA MAT SMI S global PM NU U MPC BT REACT global CNST global TITLE FILE1 FILE2 global LINP LOUT disp(sprintf(Output for Input Data from file %sn,FILE1); fprintf(LOUT,Output for Input Data from file %sn,FI

13、LE1); disp(TITLE); fprintf(LOUT,%sn,TITLE); disp( Node# X-Displ Rotation); fprintf(LOUT, Node# X-Displ Rotationn); I=1:NN; % print a matrix disp(sprintf( %4d %15.4E %15.4En,I,F(2*I-1),F(2*I); fprintf(LOUT, %4d %15.4E %15.4En,I,F(2*I-1),F(2*I); %- Reaction Calculation - disp(sprintf( DOF# Reaction); fprintf(LOUT, DOF# Reactionn); for I = 1:ND N = NU(I); R = CNST * (U(I) - F(N); disp(sprintf( %4d %15.4E,N,REACT(I); fprintf(LOUT, %4d %15.4En,N,REACT(I); end fclose(LOUT); disp(sprintf(The Results are available in the text file %s, FILE2);