平面三角形3节点有限元程序
平面三角形３结点有限元程序
1、程序名：FEMT3.FOR，FEMT3.EXE
2、程序功能
本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结点的支座反力。
程序采用Visual Fortran 5.0编制而成，输入数据全部采用自由格式。
3、程序流程及框图
图１-1 程序流程图
图１-2 程序框图
其中，各子程序的功能如下：
INPUT——输入结点坐标、单元信息和材料参数；
MR——形成结点自由度序号矩阵；
FORMMA——形成指标矩阵MA（N）并调用其他功能子程序，相当于主控程序；DIV——取出单元的3个结点号码和该单元的材料号并计算单元的b i，c i等；MGK——形成整体劲度矩阵并按一维存放在SK（NH）中；
LOAD——形成整体结点荷载列阵F；
OUTPUT——输出结点位移或结点荷载；
TREAT——由于有非零已知位移，对K和F进行处理；
DECOMP——整体劲度矩阵的分解运算；
FOBA——前代、回代求出未知结点位移δ；
ERFAC——计算约束结点的支座反力；
KRS——计算单元劲度矩阵中的子块K rs。
4、输入数据及变量说明
当程序开始运行时，按屏幕提示，键入数据文件的名字。
在运行程序之前，必须根据程序中输入要求建立一个存放原始数据的文件，这个文件的名字由少于8个字符或数字组成。数据文件包括如下内容：
⑴总控信息，共一条，9个数据
NP，NE，NM，NR，NI，NL，NG，ND，NC
NP——结点总数；
NE——单元总数；
NM——材料类型总数；
NR——约束结点总数；
NI ——问题类型标识，0为平面应力问题，1为平面应变问题；
NL ——受荷载作用的结点的数目；
NG ——考虑自重作用为1，不计自重为0；
ND ——非零已知位移结点的数目；
NC ——要计算支座约束反力的结点数目。
⑵ 材料信息，共NM 条，每条依次输入
EO ，VO ，W ，t
EO ——弹性模量（kN/m 2
）；
VO ——泊松比；
W ——材料容重（kN/m 3）；
t ——单元厚度（m ）。
这些信息都存放在数组AE （4，NM ）中。
⑶ 坐标信息，共NP 条，每条依次输入
IP ，X ，Y
IP ——结点号；
X ，Y ——分别为结点的x 坐标和y 坐标。
坐标信息存放在数组X(2，NP)中。
⑷ 单元信息，共NE 条，每条依次输入
JE ，L ，Io ，Jo ，Mo
JE ——单元号；
L ——为该单元的材料类型号。
Io ，Jo ，Mo ——分别为该单元i ，j ，m 的整体编码。
单元信息存放在数组MEO(4,NE)中。
⑸ 约束信息共NR 条，每条依次输入一个数 ??? ×× IP x I y I
IP ——结点号；
I x ，I y ——分别为该结点的约束情况，如果方向受约束时填0，如果自由则填1。
⑹ 荷载信息，共NL 条，每条依次输入y
IP ，F x ，F y
IP ——结点号；
F x ，F y ——分别为该结点的x ，y 方向的荷载分量（kN ）。
结点号存放在数组NF （NL ）中，结点荷载分量存放在数组FV （2，NL ）中。
⑺ 若ND ＞0，输入非零已知位移信息，共ND 条，每条依次输入
IP ，u x ，u y
IP ——结点号；
u x ，u y ——分别为该结点x ，y 方向已知位移分量（m ），若其中某方向为自由，则其相应分量为0。
结点号存放在数NDI （ND ）中，已知位移分量存放在数组DV(2，ND)中。
⑻ 支座反力信息，共NC 条，每条依次输入
IP ，M1，M2，M3，M4
IP ——支座结点号；
M1，M2，M3，M4 ——为与该支座结点相关的单元号，若不足4个，则用0补充。支座结点号存放在数组NCI(NC)中，相关单元号存放在数组NCE(4，NC)中。
以上数据须按如上顺序存放在数据文件中。除此之外，程序中还用到其他一些主要变量和数组，说明如下：
N ——结构自由度总数；
NH ——按一维存贮的整体劲度矩阵的总容量；
NX ——最大半带宽；
SK(10000)——维存贮的劲度矩阵；
R(1000)——开始存放等效结点荷载，求解方程以后，用来存放结点位移；
B(6)——存放单元应力σx ，σy ，τxy ，σ1，σ2，α；
MA(1000)——主元素序号指标矩阵；
JR(2，500)——结点自由度序号矩阵；
ME(3)——存放单元结点i ，j ，m 的整体编码；
NN(6)——单元结点自由度序号；
BI(3)，CI(3)——单元劲度矩阵计算公式中的b i ，b j ，b m 和c i ，c j ，c m ；
S ——三角形单元的面积；
H 11，H 12，H 21，H 22——单元劲度矩阵中子块K rs 的4个元素。
5、算例
一个正方形弹性体，厚度为1m ，四边受单位均布法向力作用，由于对称性，取其1/4进行计算，其有限元网格如图2-3所示，设2
7/1042m kn E ??=，167.0=μ，不考虑自重。该问题的精确解应力为x σ＝12kN /m ，y σ＝12kN /m ，xy τ＝0。
图１-3 有限元网格
（1）输入文件数据
6 4 1 5 0 3 0 0 5
2000.0 0.0 0.0 1.0
1 0.0 2.0
2 0.0 1.0
3 1.0 1.0
4 0.0 0.0
5 1.0 0.0
6 2.0 0.0
1 1 3 1 2
2 1 2 4 5
3 1 3 2 5
4 1
5
6 3
1 0 1
2 0 1
4 0 0
5 1 0
6 1 0
1 -0.5 -0.5
3 -1.0 -1.0
6 -0.5 -0.5
1 1 0 0 0
2 1 2
3 0
4 2 0 0 0
5 2 3 4 0
6 4 0 0 0
（2）输出文件结果
NODAL DISPLACEMENTS
NODE X-COMP. Y-COMP.
1 0.00000E+00 -0.10000E-02
2 0.00000E+00 -0.50000E-03
3 -0.50000E-03 -0.50000E-03
4 0.00000E+00 0.00000E+00
5 -0.50000E-03 0.00000E+00
6 -0.10000E-02 0.00000E+00
ELEMENT STRESSES
ELEMENT X-STRESS Y-STRESS XY-STRESS MAX-STRESS MIN-STRESS ANGLE
1 -1.000 -1.000 0.000 -1.000 -1.000 90.000
2 -1.000 -1.000 0.000 -1.000 -1.000 90.000
3 -1.000 -1.000 0.000 -1.000 -1.000 90.000
4 -1.000 -1.000 0.000 -1.000 -1.000 90.000
NODE STRESSES
NODE X-STRESS Y-STRESS XY-STRESS MAX-STRESS MIN-STRESS ANGLE
1 -1.000 -1.000 0.000 -1.000 -1.000 90.000
2 -1.000 -1.000 0.000 -1.000 -1.000 90.000
3 -1.000 -1.000 0.000 -1.000 -1.000 90.000
4 -1.000 -1.000 0.000 -1.000 -1.000 90.000
5 -1.000 -1.000 0.000 -1.000 -1.000 90.000
6 -1.000 -1.000 0.000 -1.000 -1.000 90.000
NODAL REACTIONS
NODE X-COMP Y-COMP
1 0.000 0.000
2 1.000 0.000
4 0.500 0.500
5 0.000 1.000
6 0.000 0.000
6、源程序
C FINITE ELEMENT PROGRAM FOR TWO DIMENSIONAL
C TRIANGLE ELEMENT
C
DIMENSION K(800000),COOR(2,3000),AE(4,11),
* MEL(5,2000),MA(6000)
CHARACTER*32 dat
COMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC WRITE(*,300)
300 FORMAT(///' '
* ':****'/'+PLEASE INPUT FILE NAME OF DATA')
READ(*,*)data
OPEN (4,FILE=data,STATUS='OLD')
OPEN (7,FILE='OUT',STATUS='UNKNOWN')
READ (4,*) NP,NE,NM,NR,NI,NL,NG,ND,NC
C WRITE (*,400) NP,NE,NM,NR,NI,NL,NG,ND,NC
C WRITE (7,400) NP,NE,NM,NR,NI,NL,NG,ND,NC
CALL INPUT (JR,COOR,MEL,AE)
CALL CBAND (MA,JR,MEL)
CALL SK0(SK,R,COOR,MEL,MA,JR,AE)
CALL LOAD (COOR,MEL,R,JR,AE) IF (ND.GT.0) CALL TREAT (SK,MA,JR,R) CALL DECOMP (SK,MA) CALL FOBA (SK,MA,R) WRITE(*,650) WRITE(7,650)
CALL OUTPUT(JR,R) WRITE(*,700) WRITE(7,700)
CALL CES (COOR,MEL,JR,R,AE) IF(NC.GT.0) call ERFAC (COOR,MEL,JR,R,AE)
400 FORMAT (/2X,'NP=',I3,2X,'NE=',I3,2X,'NM='
*,I3,2X,'NR=',I3,2X,'NI='I3,2X,'NL=',I3,2X,
*'NG=',I3,2X,'ND=',I3,2X,'NC=',I3)
500 FORMAT(1X,'TOTAL DEGREES OF FREEDOM N=',
*I4,1X,'TOTAL-STORAGE ','NH=',I5,1X,
*'MAX-SEMI-BANDWIDTH MX=',I3)
550 FORMAT(/20X,'TOTAL STORAGE IS
*GREATER THAN 50000')
600 FORMAT(30X,'NODAL FORCES'/8X,'NODE',
*11X,'X-COMP.',14X,'Y-COMP.')
650 FORMAT(/30X,'NODAL DISPLACEMENTS'/8X,
*'NODE',13X,'X-COMP.',12X,'Y-COMP.')
700 FORMAT(/30X,'ELEMENT STRESSES'/5X,
*'ELEMENT',5X,'X-STRESS',3X,'Y-STRESS',
*2X,'XY-STRESS',1X,'MAX-STRESS',1X,
*'MIN-STRESS',6X,'ANGLE'/)
STOP END C *********************************************
SUBROUTINE KRS (BR,BS,CR,CS)
COMMON /CB/ EO,VO,W,T,A,H11,H12,H21,H22
*,ME(3),BI(3),CI(3)
ET=EO*T/(1.0-VO*VO)/A/4.0 V=(1.0-VO)/2.0 H11=ET*(BR*BS+V*CR*CS) H12=ET*(VO*BR*CS+V*BS*CR)
H21=ET*(VO*CR*BS+V*BR*CS) H22=ET*(CR*CS+V*BR*BS) RETURN END
C *********************************************
SUBROUTINE INPUT (JR,COOR,MEL,AE)
DIMENSION JR(2,*),COOR(2,*),AE(4,*),MEL(3,*)
COMMON /CA/ NP,NE,NM,NR
COMMON /CC/ N,MX,NH
DO 70 I=1,NP
READ(4,*) IP,X,Y
COOR(1,IP)=X
COOR(2,IP)=Y
70 CONTINUE
DO 11 J=1,NE
READ(4,*)NEE,NME,(MEL(I,NEE),I=1,3)
MEL(3,NEE)=NME
11 CONTINUE
DO 10 I=1,NP
DO 10 J=1,2
10 JR(J,I)=1
DO 20 I=1,NR
READ(4,*) IP,IX,IY
JR(1,IP)=IX
JR(2,IP)=IY
20 CONTINUE
N=0
DO 30 I=1,NP
DO 30 J=1,2
IF (JR(J,I)) 30,30,25
25 N=N+1
JR(J,I)=N
30 CONTINUE
DO 55 J=1,NM
READ (4,*)jj,(AE(I,jj),I=1,4)
C WRITE(*,910)jj,(AE(I,jj),I=1,4)
IF(NI.eq.1)then
AE(1,jj)=AE(1,jj)/(1.0-AE(2,jj)*AE(2,jj)) AE(2,jj)=AE(2,jj)/(1.0-AE(2,jj))
endif
55 CONTINUE
910 FORMAT (/20X,'MATERIAL PROPERTIES'/(3X,I5,4(1x,E8.3)))
RETURN
END
C **********************************************
SUBROUTINE CBAND (MA,JR,MEL)
DIMENSION MA(*),JR(2,*),MEL(3,*),NN(6)
COMMON /CA/ NP,NE,NM,NR
COMMON /CC/ N,MX,NH
DO 65 I=1,N
65 MA(I)=0
DO 90 IE=1,NE
DO 75 K=1,3
IEK=MEL(K,IE)
DO 95 M=1,2
JJ=2*(K-1)+M
NN(JJ)=JR(M,IEK)
95 CONTINUE
75 CONTINUE
L=N
DO 80 I=1,6
NNI=NN(I)
IF(NNI.EQ.0) GO TO 80
IF(NNI.LT.L) L=NNI
80 CONTINUE
DO 85 M=1,6
JP=NN(M)
IF(JP.EQ.0) GO TO 85
JPL=JP-L+1
IF(JPL.GT.MA(JP)) MA(JP)=JPL
85 CONTINUE
90 CONTINUE
MX=0
MA(1)=1
DO 10 I=2,N
IF(MA(I).GT.MX) MX=MA(I)
MA(I)=MA(I)+MA(I-1)
10 CONTINUE
NH=MA(N)
WRITE (*,500) N,MX,NH
WRITE (7,500) N,MX,NH
500 FORMAT (/5X,'FREEDOM N='
*,I5,3X,'SEMI-BANDWI. MX=',I5,3X,
* 'STORAGE NH=',I7)
RETURN
END C***********************************************
SUBROUTINE SK0(SK,R,COOR,MEL,MA,JR,AE)
DIMENSION AE(4,*),COOR(2,*),MEL(3,*),JR(2,*),R(*),
*MA(*),SK(*),SKE(6,6),NN(6)
COMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC
COMMON /CB/ EO,VO,W,T,A,H11,H12,H21,H22,
*ME(3),BI(3),CI(3)
COMMON /CC/ N,NH
DO 10 I=1,NH 10 SK(I)=0.0 DO 70 IE=1,NE CALL DIV (IE,COOR,MEL,AE) DO 30 I=1,3
DO 30 J=1,3 CALL KRS (BI(I),BI(J),CI(I),CI(J))
SKE(2*I-1,2*J-1)=H11
SKE(2*I-1,2*J)=H12
SKE(2*I,2*J-1)=H21
SKE(2*I,2*J)=H22
30 CONTINUE
DO 40 I=1,3
J2=ME(I) DO 40 J=1,2 J3=2*(I-1)+J NN(J3)=JR(J,J2) 40 CONTINUE DO 60 I=1,6
DO 60 J=1,6 IF(NN(J).EQ.0.OR.NN(I).LT.NN(J)) GO TO 60
JJ=NN(I) JK=NN(J) JL=MA(JJ) JM=JJ-JK JN=JL-JM
SK(JN)=SK(JN)+SKE(I,J)
60 CONTINUE
70 CONTINUE
c WRITE(0,500) (SK(I),I=1,20)
500 FORMAT(/10X,'SK='/(6F12.5))
RETURN END
C **********************************************
SUBROUTINE LOAD (COOR,MEL,R,JR,AE) DIMENSION AE(4,*),COOR(2,*),MEL(3,*),R(*),JR(2,*),
* NF(50),FV(2,50)
COMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC
COMMON /CB/ EO,VO,W,T,A,A1(4),ME(3),BB(6)
COMMON /CC/ N,NH
DO 10 I=1,N 10 R(I)=0.0 IF(NG) 70,70,30
30 DO 60 IE=1,NE CALL DIV (IE,COOR,MEL,AE) DO 50 I=1,3 J2=ME(I) J3=JR(2,J2) IF(J3) 50,50,40 40 R(J3)=R(J3)-T*W*A/3.0
50 CONTINUE 60 CONTINUE 70 IF(NL) 110,110,80
80 READ(4,*) (NF(I),I=1,NL)
READ(4,*) ((FV(I,J),I=1,2),J=1,NL)
C WRITE(*,500) (NF(I),I=1,NL)
C WRITE(7,500) (NF(I),I=1,NL)
C WRITE(*,600) ((FV(I,J),I=1,2),J=1,NL)
C WRITE(7,600) ((FV(I,J),I=1,2),J=1,NL)
DO 100 I=1,NL JJ=NF(I) J=JR(1,JJ) M=JR(2,JJ) IF (J.GT.0) R(J)=R(J)+FV(1,I)
IF (M.GT.0) R(M)=R(M)+FV(2,I) 100 CONTINUE 110 RETURN 500 FORMAT(/20X,'NODES OF APPLIED LOAD***NF='
*/(1X,10I8))
600 FORMAT(/30X,'LUMPED-LOADS***FV='
*/(5X,5F15.3))
END C *********************************************
SUBROUTINE TREAT (SK,MA,JR,R) DIMENSION SK(*),MA(*),NDI(75),DV(2,75),JR(2,*),R(*) COMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC
COMMON /CC/ N,NH
READ(4,*) (NDI(J),J=1,ND) READ(4,*) ((DV(I,J),I=1,2),J=1,ND)
C WRITE(*,500) (NDI(J),J=1,ND)
C WRITE(7,500) (NDI(J),J=1,ND)
C WRITE(*,550) ((DV(I,J),I=1,2),J=1,ND) C WRITE(7,550) ((DV(I,J),I=1,2),J=1,ND)
DO 20 I=1,ND DO 20 J=1,2 IF(DV(J,I)) 10,20,10
10 JJ=NDI(I) L=JR(J,JJ) JN=MA(L) SK(JN)=1.0E30 R(L)=DV(J,I)*1.0E30 20 CONTINUE 500 FORMAT(/25X,'NODE NO.**NDI='/(1X,10I8))
550 FORMAT(/25X,'DISPLACEMENT-VALUES**DV='/
*(10X,6F10.6))
RETURN END C *********************************************
SUBROUTINE DECOMP (SK,MA)
DIMENSION SK(*),MA(*)
COMMON /CC/ N,NH
L=I-MA(I)+MA(I-1)+1
K=I-1 L1=L+1 IF (L1.GT.K) GO TO 30 DO 20 J=L1,K IJ=MA(I)-I+J M=J-MA(J)+MA(J-1)+1
IF (L.GT.M) M=L MP=J-1 IF (M.GT.MP) GO TO 20
DO 10 LP=M,MP
IP=MA(I)-I+LP JP=MA(J)-J+LP SK(IJ)=SK(IJ)-SK(IP)*SK(JP)
10 CONTINUE 20 CONTINUE 30 IF (L.GT.K) GO TO 50
DO 40 LP=L,K
IP=MA(I)-I+LP LPP=MA(LP) SK(IP)=SK(IP)/SK(LPP)
II=MA(I) SK(II)=SK(II)-SK(IP)*SK(IP)*SK(LPP)
40 CONTINUE
50 CONTINUE
C WRITE(0,500) (SK(I),I=1,20)
500 FORMAT(/10X,'SK='/(1X,6F12.4))
RETURN END C *********************************************
SUBROUTINE FOBA (SK,MA,R) DIMENSION SK(*),MA(*),R(*) COMMON /CC/ N,NH
DO 10 I=2,N L=I-MA(I)+MA(I-1)+1
K=I-1 IF (L.GT.K) GO TO 10
DO 5 LP=L,K IP=MA(I)-I+LP R(I)=R(I)-SK(IP)*R(LP)
5 CONTINUE
10 CONTINUE DO 20 I=1,N
R(I)=R(I)/SK(II)
20 CONTINUE DO 30 J1=2,N I=2+N-J1
L=I-MA(I)+MA(I-1)+1
K=I-1 IF (L.GT.K) GO TO 30
DO 25 J=L,K IJ=MA(I)-I+J R(J)=R(J)-SK(IJ)*R(I)
25 CONTINUE
30 CONTINUE RETURN END C *********************************************
SUBROUTINE CES (COOR,MEL,JR,R,AE) DIMENSION AE(4,*),COOR(2,*),MEL(3,*),JR(2,*),R(*),B(6) COMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC
COMMON /CB/ EO,VO,W,T,A,H11,H12,H21,H22,
*ME(3),BI(3),CI(3)
COMMON /CC/ N,NH
DO 100 IE=1,NE CALL DIV (IE,COOR,MEL,AE) ET=EO/(1.0-VO*VO)/A/2.0
DO 50 I=1,3 J2=ME(I) I2=JR(1,J2)
I3=JR(2,J2) IF(I2) 30,20,10
10 B(2*I-1)=R(I2)
GO TO 30 20 B(2*I-1)=0.0 30 IF(I3) 50,40,35
35 B(2*I)=R(I3)
GO TO 50 40 B(2*I)=0.0 50 CONTINUE H1=0.0 H2=0.0 H3=0.0 DO 60 I=1,3
H1=H1+BI(I)*B(2*I-1)
H2=H2+CI(I)*B(2*I)
H3=H3+BI(I)*B(2*I)+CI(I)*B(2*I-1)
60 CONTINUE
A1=ET*(H1+VO*H2)
A2=ET*(H2+VO*H1)
A3=ET*(1.0-VO)*H3/2.0
H1=A1+A2 H2=SQRT((A1-A2)*(A1-A2)+4.0*A3*A3)
B(4)=(H1+H2)/2.0
B(5)=(H1-H2)/2.0
IF (ABS(A3).GT.1E-4) GO TO 80
IF (A1.GT.A2) GO TO 70
B(6)=90.0 GO TO 90 70 B(6)=0.0 GO TO 90 80 B(6)=ATAN((B(4)-A1)/A3)*57.29578
90 B(1)=A1 B(2)=A2 B(3)=A3 WRITE(0,500) IE,B
WRITE(7,500) IE,B
100 CONTINUE 500 FORMAT(6X,I4,3X,6F11.3)
RETURN END C **********************************************
SUBROUTINE OUTPUT(JR,R) DIMENSION JR(2,*),R(*)
COMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC
COMMON /CC/ N,NH
DO 100 I=1,NP L=JR(1,I) IF(L) 30,20,10 10 S=R(L) GO TO 30 20 S=0.0
30 L=JR(2,I) IF(L) 60,50,40 40 SS=R(L) GO TO 60 50 SS=0.0 60 WRITE(*,500) I,S,SS
WRITE(7,500) I,S,SS
100 CONTINUE
500 FORMAT(5X,I5,2F20.5)
RETURN END C ********************************************
SUBROUTINE ERFAC (COOR,MEL,JR,R,AE) DIMENSION NCI(20),NCE(4,20),MEL(3,*),JR(2,*),R(*),
* AE(4,*),COOR(2,*)
COMMON /CA/ NP,NE,NM,NA(5),NC
COMMON /CB/ AB(5),H11,H12,H21,H22,
*ME(3),BI(3),CI(3)
COMMON /CC/ N,NH
READ(4,*) (NCI(J),J=1,NC)
READ(4,*) ((NCE(I,J),I=1,4),J=1,NC)
WRITE(*,500) (NCI(J),J=1,NC)
WRITE(*,600) ((NCE(I,J),I=1,4),J=1,NC)
WRITE(*,700) WRITE(7,500) (NCI(J),J=1,NC)
WRITE(7,600) ((NCE(I,J),I=1,4),J=1,NC)
WRITE(7,700)
DO 120 JJ=1,NC FX=0.0 FY=0.0
L=NCI(JJ) DO 110 M=1,4 IF(NCE(M,JJ)) 110,110,10
10 IE=NCE(M,JJ)
CALL DIV (IE,COOR,MEL,AE) DO 20 IM=1,3 K=IM
IF(L-ME(IM)) 20,30,20 20 CONTINUE WRITE(0,750) L WRITE(7,750) L
30 DO 100 IP=1,3 CALL KRS (BI(K),BI(IP),CI(K),CI(IP))
NL=ME(IP) JI=JR(1,NL)
JP=JR(2,NL) IF(JI) 60,40,50 40 S=0.0 GO TO 60 50 S=R(JI) 60 IF(JP) 70,70,80 70 SS=0.0
GO TO 90 80 SS=R(JP) 90 FX=FX+H11*S+H12*SS
FY=FY+H21*S+H22*SS 100 CONTINUE 110 CONTINUE WRITE(0,800) L,FX,FY
WRITE(7,800) L,FX,FY
120 CONTINUE 500 FORMAT(30X,'NODE NO.**NCI='/(1X,10I8))
600 FORMAT(30X,'ELEMENT-NO.**NCE='/
*(1X,10I8))
700 FORMAT(30X,'NODAL REACTIONS'/8X,
*'NODE',14X,'X-COMP',14X,'Y-COMP')
750 FORMAT(/10X,'ERROR OF ELEMENT MESSAGE'
*'****NODE NUMBE',I5)
800 FORMAT(6X,I5,2F20.3)
RETURN END
C **********************************************
SUBROUTINE DIV (IE,COOR,MEL,AE) DIMENSION COOR(2,*),AE(4,*),MEL(3,*)
COMMON /CB/ EO,VO,W,T,A,A1(4),ME(3),
*BI(3),CI(3) ME(1)=MEL(1,IE) ME(2)=MEL(2,IE) ME(3)=MEL(3,IE)
I=ME(1)
J=ME(2)
M=ME(3)
L=MEL(3,IE) BI(1)=COOR(2,J)-COOR(2,M) BI(2)=COOR(2,M)-COOR(2,I) BI(3)=COOR(2,I)-COOR(2,J)
CI(1)=COOR(1,M)-COOR(1,J) CI(2)=COOR(1,I)-COOR(1,M) CI(3)=COOR(1,J)-COOR(1,I) A=(BI(2)*CI(3)-CI(2)*BI(3))/2.0
EO=AE(1,L) VO=AE(2,L) W=AE(3,L) T=AE(4,L) RETURN
END
matlab有限元分析实例
MATLAB: MATLAB是美国MathWorks公司出品的商业数学软件，用于数据分析、无线通信、深度学习、图像处理与计算机视觉、信号处理、量化金融与风险管理、机器人，控制系统等领域。 MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室），软件主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式。 MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等。MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。 MATLAB有限元分析与应用:
《MATLAB有限元分析与应用》是2004年4月清华大学出版社出版的图书，作者是卡坦，译者是韩来彬。 内容简介： 《MATLAB有限元分析与应用》特别强调对MATLAB的交互应用，书中的每个示例都以交互的方式求解，使读者很容易就能把MATLAB用于有限分析和应用。另外，《MATLAB有限元分析与应用》还提供了大量免费资源。 《MATLAB有限元分析与应用》采用当今在工程和工程教育方面非常流行的数学软件MATLAB来进行有限元的分析和应用。《MATLAB有限元分析与应用》由简单到复杂，循序渐进地介绍了各种有限元及其分析与应用方法。书中提供了大量取自机械工程、土木工程、航空航天工程和材料科学的示例和习题，具有很高的工程应用价值。
平面三角形单元有限元程序设计
. 一、题目 如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。已知：P=150N/m ，E=200GPa ，=0.25，t=0.1m ，忽略自重。试计算薄板的位移及应力分布。 要求： 1. 编写有限元计算机程序，计算节点位移及单元应力。（划分三角形 单元，单元数不得少于30个）； 2. 采用有限元软件分析该问题（有限元软件网格与程序设计网格必 须一致），详细给出有限元软件每一步的操作过程，并将结果与程序计算结果进行对比（任选取三个点，对比位移值）； 3. 提交程序编写过程的详细报告及计算机程序； 4. 所有同学参加答辩，并演示有限元计算程序。 有限元法中三节点三角形分析结构的步骤如下： 1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。 2）单元分析，建立单元刚度矩阵。 3）整体分析，建立总刚矩阵。 4）建立整体结构的等效节点荷载和总荷载矩阵 5）边界条件处理。 6）解方程，求出节点位移。 7）求出各单元的单元应力。 8）计算结果整理。 一、程序设计 网格划分 如图，将薄板如图划分为6行，并建立坐标系，则
刚度矩阵的集成 建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。 通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种； （2）该单元对应总刚度矩阵的那几行哪几列 （3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列 循环又分为3层循环：（1）最外层：逐行计算 （2）中间层：该行逐个计算 （3）最里层：区分为第 奇/偶 数个计算 单元刚度的集成：[ ][][][][][]' '''''215656665656266256561661e Z e e e Z e Z e e e e k k k K k k k k k k +?++=? =?==?==?=?????? 边界约束的处理：划0置1法 X Y P X Y P
一三节点三角形单元.docx
有限元课程总结 一三节点三角形单元 1位移函数 移函数写成矩阵形式为： 确定六个待定系数 a4 v玉＞ 矩阵形式如下： J“= TV, 0 Nj 0 N m bJ _ 0 TV, 0 Nj 0 2单元刚度矩阵的计算 1）单元应变和节点位移的关系 由几何方程可以得到单元的应变表达式, 5 6 > = ----------------- b . 「2A ' 7 厂 f Mg Y — A ”——， Y Cd As _ u i 匕? 宀=[N]{5丫 V7 u i 8x dv du dv ----- 1 ---- dy dx J_ 2A C C i bj 0 0 Cj C J b J u j V J
2）单元应力与单元节点位移的关系 [KJ = [B r ]T [D][B s ] b r b s + — c r c s t s 2 * s “也+与仏 (T = i,jjn;s = i,jjn) 3) 单元刚度矩阵 卩心][K“] [K]J [K }i ] [K 〃] [心][K mj ] 3载荷移置 1）集中力的移置 图3 由虚功相等可得， （㈤丁附=（Q YJW ｛P ｝ 由于虚位移是任意的，则皿｝"=["卩｛鬥 2）体力的移置 [S M D I B .] = E 2A(1-Z /2) Mi Ci % 2 z 如图3所示, 令单元所受的均匀分布体力为｛〃｝= Et 4(1 —“2)A 地C$ + [DfB i % [K 加 [K 如 6
由虚功相等可得， （{J*r）r{7?r =^}>f[N]r{p}tdxdy {R}e =\\[N]r{p}tdxdy 3）分布面力的移置 设在单元的边上分布有面力{可二[片了r,同样可以得到结点载荷, {R}e=\[N]T{P}tds 4.引入约束条件，修改刚度方程并求解 1）乘大数法处理边界条件 图3?4所示的结构的约束和载荷情况，如图3?7所示。结点1、4上有水平 方向的位移约束，结点4、6上有垂直方向的约束，结点3上作用有集中力（'， 匕）。 整体刚度矩阵[K]求出后，结构上的结点力可以表示为: {F} = [K]{5} 根据力的平衡，结点上的结点力与结点载荷或约束反力平衡。用{?}表示结 点载荷和支杆反力，则可以得到结点的平衡方程： [K]0}={P} （3.4） 这样构成的结点平衡方程组，在右边向量{P}中存在未知量，因此在求解平衡
有限元2-弹性力学平面问题有限单元法(2.1三角形单元,2.2几个问题的讨论)综述
第2章 弹性力学平面问题有限单元法 2.1 三角形单元(triangular Element) 三角形单元是有限元分析中的常见单元形式之一，它的优点是: ①对边界形状的适应性较好，②单刚形式及其推导比较简单,故首先介绍之。 一、结点位移和结点力列阵 设右图为从某一结构中取出的一典型三角形单元。 在平面应力问题中，单元的每个结点上有沿x 、y 两个方向的力和位移，单元的结点位移列阵规定为： 相应结点力列阵为: (式2-1-1) 二、单元位移函数和形状函数 前已述及，有限单元法是一种近似方法，在单元分析中，首先要求假定（构 造）一组在单元内有定义的位移函数作为近似计算的基础。即以结点位移为已知量，假定一个能表示单元内部（包括边界）任意点位移变化规律的函数。 构造位移函数的方法是：以结点(i,j,m)为定点。以位移(u i ,v i ,…u m v m )为定点上的函数值，利用普通的函数插值法构造出一个单元位移函数。 在平面应力问题中，有u,v 两个方向的位移,若假定单元位移函数是线性的,则可表示成: (,)123 u u x y x y ααα==++ 546(,)v v x y x y ααα==++ (2-1-2)a 式中的6个待定常数α1 ,…, α6 可由已知的6个结点位移分量(3个结点的坐标) {}??? ?? ?????=????? ???? ?????????????=m j i m e d d d d m j j i v u v u v u i {} i i j j m X Y X (2-1-1)Y X Y i e j m m F F F F ?? ?? ???? ???? ??==??????????????????
有限元的MATLAB解法
有限元的MATLAB解法 1.打开MATLAB。 2.输入“pdetool”再回车，会跳出PDE Toolbox的窗口（PDE意为偏微分方程，是partial differential equations的缩写），需要的话可点击Options菜单下Grid命令，打开栅格。 3.完成平面几何模型：在PDE Toolbox的窗口中，点击工具栏下的矩形几何模型进行制作模型，可画矩形R，椭圆E，圆C，然后在Set formula栏进行编辑并（如双脊波导R1+R2+R3改为RI-R2-R3,设定a、b、s/a、d/b的值从而方便下步设定坐标） 用算术运算符将图形对象名称连接起来，若还需要，可进行储存，形成M文件。 4.用左键双击矩形进行坐标设置：将大的矩形left和bottom都设为0，width是矩形波导的X轴的长度，height是矩形波导的y轴的长度，以大的矩形左下角点为原点坐标为参考设置其他矩形坐标。 5.进行边界设置：点击“Boundary”中的“Boundary Mode”,再点击
“Boundary”中的“Specify Boundary Conditions”,选择符合的边界条件，Neumann为诺曼条件,Dirichlet为狄利克雷条件，边界颜色显示为红色。 6.进入PDE模式：点击"PDE"菜单下“PDE Mode”命令，进入PDE 模式，单击“PDE Specification”,设置方程类型，“Elliptic”为椭圆型，“Parabolic”为抛物型，“Hyperbolic”为双曲型，“Eigenmodes”为特征值问题。 7.对模型进行剖分：点击“Mesh”中“Initialize Mesh”进行初次剖分，若要剖的更细，再点击“Refine Mesh”进行网格加密。 8.进行计算：点击“Solve”中“Solve PDE”,解偏微分方程并显示图形解，u值即为Hz或者Ez。 9.单击“Plot”菜单下“Parameters”选项，打开“Plot Selection”对话框。选中Color,Height(3-D plot)和Show mesh三项，然后单击“Plot”按钮，显示三维图形解。 10.如果要画等值线图和矢量场图，单击“Plot”菜单下“Parameters”选项,打开“Plot Selection”对话框。选中Contour和Arrows两项，然后单击Plot按钮，可显示解的等值线图和矢量场图。 11.将计算结果条件和边界导入MATLAB中：点击“Export Solution”,再点击“Mesh”中“Export Mesh”。
有限元习题3
第一章 1.有限单元法求得的解为：[ ]3 A.精确解 B.解析解 C.近似解 D.整数解 2.弹性力学问题的基本解法有：[ ] ABD A.按位移求解 B.按应力求解 C.按单元刚度求解 D. 混合求解 E.按整体刚度求解 23.弹性力学问题的基本解法有：按位移求解，按应力求解和[ ]3 A. 按单元刚度求解 B. 按整体刚度求解 C. 混合求解 D.按平衡方程求解 24.弹性力学问题的基本解法有：按位移求解，混合求解和[]4 A. 按平衡方程求解 B. 按单元刚度求解 C. 按整体刚度求解 D. 按应力求解 25.弹性力学问题的基本解法有：按应力求解，混合求解和[ ]2 A. 按整体刚度求解 B. 按位移求解 C. 按单元刚度求解 D. 按平衡方程求解 3.用弹性力学经典解法解决实际问题的主要困难在于：[ ]4 A.对弹性体离散化的复杂性 B.刚度矩阵求解的困难性 C.受力分析的复杂性 D.求解偏微分方程的复杂性 4.用三角形单元的节点位移，可以表示单元中的：[ ]BDE A.弯矩 B.应变 C.扭矩 D.应力 E.结点力 26.用三角形单元的节点位移，可以表示单元中的应变，应力和[ ]3 A. 扭矩 B. 弯矩 C. 结点力 D.外力 27.用三角形单元的节点位移，可以表示单元中的应变，结点力和[ ]4 A.弯矩 B. 外力 C. 扭矩 D. 应力 28. 用三角形单元的节点位移，可以表示单元中的应力，结点力和[ ]4， A. 外力 B. 扭矩 C. 弯矩 D. 应变 5.将各个单元集合成离散化的结构模型进行整体分析，问题最后归结为求解[ ]。2 A.结点位移 B. 以结点位移为未知量的线性方程组 C.整体刚度矩阵 D.单元刚度矩阵 6.对于三角形三结点单元，其结点按照[]顺序进行排列。3 A.从左至右 B. 顺时针 C. 逆时针 D.以上均可 7.对于三角形三结点单元，每个结点位移在单元平面内有[ ]个分量 2 A.1 B.2 C.3 D.4 8.对于三角形三结点单元，共有[ ]个位移分量。4 A.3 B.4 C.5 D.6 9.形函数 N在结点i上的值等于[ ]。2 i A.0 B.1 C. -1 D.2 10.在单元中任意一点，三个形函数之和等于[ ]2 A.0 B.1 C.2 D.3 11.有了单元的位移模式，就可以应用[ ]求得单元的应变3 A.平衡微分方程 B.物理方程 C. 几何方程 D.积分方程 12.单元应力矩阵[S]与弹性矩阵[D]和单元应变矩阵[B]的关系是：[ ]C A. [S]= [D]+ [B] B. [S]= [D]—[B] C. [S]= [D] [B] D. [S]= [D]/ [B] 13.三角形三结点单元中，单元应力矩阵[S]是一个[ ]4 A.对称矩阵 B.零矩阵 C.非常数矩阵 D.常数矩阵 14.三角形三结点单元的应力分量为[ ]1 A.常量 B.变量 C.零 D.不确定
基于Matlab语言的按平面三角形单元划分的结构有限元程序设计模板
基于Matlab语言的按平面三角形单元划分的结构有限元程序设计 专业：建筑与土木工程 班级：建工研12-2 姓名：韩志强 学号： 471220580
基于Matlab语言的按平面三角形单元划分 结构有限元程序设计 一、有限单元发及Matlab语言概述 1. 有限单元法 随着现代工业、生产技术的发展，不断要求设计高质量、高水平的大型、复杂和精密的机械及工程结构。为此目的，人们必须预先通过有效的计算手段，确切的预测即将诞生的机械和工程结构，在未来工作时所发生的应力、应变和位移因此，需要寻求一种简单而又精确的数值分析方法。有限单元法正是适应这种要求而产生和发展起来的一种十分有效的数值计算方法。 有限元法把一个复杂的结构分解成相对简单的“单元”，各单元之间通过结点相互连接。单元内的物理量由单元结点上的物理量按一定的假设内插得到，这样就把一个复杂结构从无限多个自由度简化为有限个单元组成的结构。我们只要分析每个单元的力学特性，然后按照有限元法的规则把这些单元“拼装”成整体，就能够得到整体结构的力学特性。 有限单元法基本步骤如下： (1)结构离散：结构离散就是建立结构的有限元模型，又称为网格划分或单元划分，即将结构离散为由有限个单元组成的有限元模型。在该步骤中，需要根据结构的几何特性、载荷情况等确定单元体内任意一点的位移插值函数。 (2)单元分析：根据弹性力学的几何方程以及物理方程确定单元的刚度矩阵。 (3)整体分析:把各个单元按原来的结构重新连接起来，并在单元刚度矩阵的基础上确定结构的总刚度矩阵，形成如下式所示的整体有限元线性方程： {}[]{}δ F=① K 式中，{}F是载荷矩阵,[]K是整体结构的刚度矩阵,{}δ是节点位移矩阵。 (4)载荷移置:根据静力等效原理，将载荷移置到相应的节点上，形成节点载荷矩阵。 (5)边界条件处理:对式①所示的有限元线性方程进行边界条件处理。 (6)求解线性方程:求解式①所示的有限元线性方程，得到节点的位移。在该步骤中，若有限元模型的节点越多，则线性方程的数量就越多，随之有限元分析的计算量也将越大。 (7)求解单元应力及应变根据求出的节点位移求解单元的应力和应变。
平面三角形单元
第八章 平面问题的有限元分析及三角形单元的应用 第一节 概述 分析弹性力学平面问题时，最简单的单元式由三个结点组成的三角形单元。当用以分析平面应力问题时，可将其视为三角板；当用以分析平面应变问题时，则可式为三棱柱。各单元在结点处为铰结。图8-1所示位移悬臂梁离散为三角形单元的组合体 以矩阵形式列出弹性力学平面问题的基本量和基本方程。 谈形体所受体力分量可表示为 [ ] T y x y x p p p p p =??? ? ????= (8-1) 所受面力分量可表示为 [ ] T y x y x p p p p p =??? ? ????= （8-2） 体内任一点应力分量可表示为 []T xy y x τδδδ= （8-3） 任一点的应变分量可表示为 []T xy y x γεεε= （8-4） 任一点的位移分量可表示为 []T v u =δ （8-5） 弹性力学平面问题的几何方程的矩阵表达式为 ?? ???? ???????? ??? ???+??????=????????????=x u y v y v x u xy y x εεεε （8-6） 平面应力问题的物理方程的矩阵表达式为 ? ???? ? ?????????? ????????? ?--= ????? ? ??????xy y x xy y x E γεεμμμ μτσσ210 0010112 （8-7） 或简写成为 εσD = （8-8） 式中
???? ?? ? ?????? ?--=210 0010112μμμ μ E D （8-9） 称为弹性矩阵。 平面应变问题的物理方程也可写成式（8-8），但须将式（8-9）中的E 换成 2 1μ -E ，μ换成 2 1μμ -，因此得出 ???? ?? ????????? ?? ?-----+-= )1(2210 00110 11)21)(1()1(2 2 μμμμμμ μμμE D （8-10） 平衡微分方程及边界条件也可以用矩阵表示，但弹性力学有限元位移法中，通常用虚功 方程代替平衡微分方程和应力边界条件。虚功方程的矩阵表达式为 ?????***=+tdxdy tds p f ptdxdy f T T σε （8-11） 式中：[ ] T v u f ** * =，表示虚位移； []T xy x x * ***=γεεε，表示与虚位移相对应的虚应变。 为了便于计算，作用于弹性体上的体力和面力替换为作用在结点上的集中力，即等效结 点荷载。设作用于各个结点上的外力分量用如下列阵来表示 []T n n V U V U V U F ?=2211 与这些结点外力分量相对应得结点虚位移分量列阵为 []T n n v u v u v u * ******?=2211δ 则外力在虚位移上做的虚功为 F v V u U v V u U v V u U T n n n n ** *****=++?++++δ22221111 如平面弹性体的厚度为t ，该虚功除以t ，即可得出单位厚度薄板上的外力虚功。于是，式（8-11）所示虚功方程可写成 ??**=tdxdy F T T σεδ （8-11） 虚功方程不仅仅应用于弹性力学，也可用于塑性力学。其应用条件是：只要变形体的全部外力和应力满足平衡方程；位移是微小的，并满足边界条件，位移与应变满足几何方程。
《有限元基础教程》_【ANSYS算例】4.7.1(3) 基于3节点三角形单元的矩形薄板分析(GUI)及命令流
【ANSYS 算例】4.7.1(3) 基于3节点三角形单元的矩形薄板分析 如图4-20所示为一矩形薄平板，在右端部受集中力100 000N F =作用，材料常数为：弹性模量7110Pa E =?、泊松比1/3μ=，板的厚度为0.1m t =，在ANSYS 平台上，按平面应力问题完成相应的力学分析。 (a) 问题描述 (a) 有限元分析模型 图4–20 右端部受集中力作用的平面问题(高深梁) 解答 在ANSYS 平台上，完成的分析如下。 1. 基于图形界面的交互式操作(step by step) (1) 进入ANSYS(设定工作目录和工作文件) 程序 → ANSYS Interactive →Working directory (设置工作目录) →Initial jobname (设置工作文件名): 2D3Node →Run → OK (2) 设置计算类型 ANSYS Main Menu : Preferences… → Structural → OK (3) 选择单元类型 ANSYS Main Menu : Preprocessor →Element Type →Add/Edit/Delete… →Add… →Solid ：Quad 4node 42 →OK (返回到Element Types 窗口) → Options… →K3: Plane Strs w/thk(带厚度的平面应力问题) →OK →Close (4) 定义材料参数 ANSYS Main Menu : Preprocessor →Material Props →Material Models →Structural →Linear →Elastic → Isotropic: EX:1.0e7 (弹性模量)，PRXY: 0.33333333 (泊松比) → OK → 鼠标点击该窗口右上角的“ ”来关闭该窗口 (5) 定义实常数以确定平面问题的厚度 ANSYS Main Menu: Preprocessor →Real Constant s… →Add/Edit/Delete →Add →Type 1→ OK →Real Constant Set No: 1 (第1号实常数), THK: 0.1 (平面问题的厚度) →OK →Close (6) 生成单元模型 生成4个节点 ANSYS Main Menu: Preprocessor →Modeling → Create → Nodes → On Working Plane →输入节点1的x,y,z 坐标(2,1,0)，回车→输入节点2的x,y,z 坐标(2,0,0)，回车→输入节点3的x,y,z 坐标(0,1,0)，回车→输入节点4的x,y,z 坐标(0,0,0)，回车→OK 定义单元属性 ANSYS Main Menu: Preprocessor →Modeling → Create → Elements → Elem Attributes →Element type number:1 →Material number:1→Real constant set number:1 →OK 生成单元 ANSYS Main Menu: Preprocessor →Modeling → Create → Elements → User Numbered → Thru Nodes →Number to assign to element:1→Pick nodes:2,3,4→OK →Number to assign to element:2→Pick nodes:3,2,1→OK (7) 模型施加约束和外载 左边两个节点施加X,Y 方向的位移约束 ANSYS Main Menu: Solution → Define Loads → Apply →Structural → Displacement → On
单元节点和积分点有什么区别
单元节点和积分节点的联系和区别 有限元方法的实质： 通过变分原理极小值转化为矩阵的极小值（(变分原理)→ (最小势能原理) (虚功原理) 变分原理：把一个物理学问题（或其他学科的问题）用变分法化为求泛函极值（或驻值）的问题，就称为该物理问题 （或其他学科的问题）的变分原理。如果建立了一个新的变分原理，它解除了原有的某问题变分原理的某些约束条件，就称为该问题的广义变分原理；如果解除了所有的约束条件，就称为无条件广义变分原理，或称为完全的广义变分原理。1964年，钱伟长教授明确提出了引进拉格朗日乘子（Lagrange multiplier ）把有约束条件的变分原理化为较少（或没有）约束条件的变分原理的方法。 最小势能原理：最小势能原理就是说当一个体系的势能最小时，系统会处于稳定平衡状态。举个例子来说，一个小球在曲面上运动，当到达曲面的最低点位置时，系统就会趋向于稳定平衡。势能最小原理与虚功原理本质上是一致的。宇宙万物，如果其势能未达到“最小”（局部概念），它总要设法变化到其“相对”最小的势能位置。举个例子：一个物体置于高山上，它相对于地面来说有正的势能（非最小），因而它总有向地面运动的“能力”（向地面“跃迁”）（其力学本质是其处于一种不稳平衡状态）。因此，它试图（也只有）向下运动，才能保证其达到一个相对平稳的状态。λ 最小势能原理是势能驻值原理在线弹性范围里的特殊情况。对于一般性问题：真实位移状态使结构的势能取驻值（一阶变分为零），在线弹性问题中取最小值。形象的说，当你在一百米高的钢丝绳上走的时候你总是希望尽早回到地上，但其实只要你不动你也是平衡的，因为驻值也可以是极大值（此时称为随遇平衡）。而当你在一百米高的大楼里的办公室里时，你并不害怕，因为周围的物体的势能均不比你小，此时驻值取的是极小值而不是最小值。在有限元的理论中，最小势能原理是在所有满足给定边界条件的位移时，满足平衡微分方程的位移使得势能取得最小值。公式如下： 在江见鲸的有限元就很清晰的介绍了有限元的原理，再参考汪新老师的ppt ，理解就十分透彻。以下是我在网上搜索到的关于有限元的书，很值得一看。另外cook 的书也很好，是陈贡发老师介绍的，很值得一看。 1《The Finite Element Method 》O.C.Zienkiewicz ，R.Taylor 著，第五版，三卷本，有中文译本 《有限元法》(英)监凯维奇著（第四版中译本1985年出版，上下册，尹泽勇等译，权威著作，有限元研究者必读；第五版译名改成了《有限单元法》，曾攀译，2008年出版~~） 2.《Nonlinear Finite Element for Continua and Structures 》T .Belytschko 等著，有中文译本 《连续体和结构的非线性有限元》庄茁译，清华大学出版社，固体力学非线性有限元的集大成之作~~ ：外力对系统作的功势能：W U 21???? ??--+??? ??=+=∏???dS dV dV W U S T V T V T σf u f u εσ
有限元三角形等参单元
北方工业大学 高等有限元课程总结 姓名:韩双鹏学号: 2011303310105 专业班级:结构研-11 系(部、院):建筑工程学院 2012 年5 月25 日
高等有限元学习总结——六节点三角形等参数单元 1 概述 从弹性力学基本方程到有限元原理再到最新进展，经过本课程的学习，比较系统的掌握了有限元相关内容，更学习到了一种方法、一些生活中的哲理。首先从大方向掌握所学内容，避免迷失在局部造成一叶遮目不见泰山之悲剧，比如弹性力学原理从大方向说就是三类方程，以及其在各类问题中的应用；其次了解了科研的相关过程及创新之处，从已知的东西到无知的领域，正如老师所说，能成功地把某一领域的东西搬到相关领域，这就是一大创造，比如有限元中将梁弯曲的理论研究厚板弯曲问题，由有限元标准单元到等参元的研究等；再有，我们生活中的常识、学习中的某些东西值得我们细细品味，也许这就是平时所说的小事反应大道理，老师的理论：“很多想法都是错误的”“很好想到的方法也许很难走通”“有缺陷的东西才更体现出美”“平衡的理论，吃点亏也许是福”等等，受益匪浅。不再一一赘述，本文将取其中的一个知识点，总结六节点三角形等参单元的相关内容。 我们知道，无论三节点或者六节点三角形单元还是四节点或者八节点矩形单元，它们形状简单、规则但计算精度低，且对于复杂边界的适应性差，难以很好的拟合曲边边界，解决这一问题的通用方法是细分边界，以直代曲，利用更多的简单单元去拟合边界复杂的区域。但这样处理仍存在折线代替曲线所带来的误差，且这种误差不能通过提高单元位移函数的精度来补偿。那么能否构造出单元形状任意、边界适应性好、计算精度高的曲边单元，以便在给定的精度下用较少数目的单元去解决实际问题？这就是有限元中一类重要的单元——等参数单元。本文将总结等参数单元的基本概念，并以六节点三角形单元为例讲述等参元实现过程中的三种变换，以及该等参元的收敛性等问题。 2 等参数单元及实现过程 2.1 等参数单元概念 由于实际问题的复杂性，通常需要使用一些形状不规整和形状复杂的单元来离散边界形状复杂的原问题。如下图所示（a）中为常见的几何形状不规整的实际单元，称为实际单元，也称为参数单元。（b）中为对应的形状规整的单元，称为标准单元。对于形状复杂的实际单元的单元分析，若仍采用前面介绍的方法进行，则在单元位移函数的建立和单元刚度矩阵计算方面会遇到许多困难。由此可考虑利用前面介绍过的形状规整的标准单元的单元分析来研究实际单元，几何形状的不同可认为是坐标变换的结果。
有限元八种三维单元介绍
有限元八种三维单元介绍 有限元三维体单元常见单元有四面体4、10节点单元、六面体8、20、27节点单元、三棱柱6、15节点单元。我们在2000年新问世的四面体20节点单元。下面分别介绍如下： 1 四面体4节点单元（常应变单元、一次单元），见图一。 单元内部的位移插值函数为一次多项式，即只含常数项和Z Y X ,,四项。应变是位移的偏导数，故在单元内部，应力和应变为常数，位移和应力收敛速度都很慢，是非常落后的单元。 图一 四面体4节点单元（常应变单元） 2 四面体10节点单元（二次单元），见图二。 用体积坐标定义的单元：单元内位移插值函数为二次完全多项式，即含常数项和Z Y X ,,,YZ XZ XY Z Y X ,,,,,222十项，在单元内部，应力和应变为一次完全多项式，位移收敛速度很快，但应力收敛速度仍较慢。由于整体加密使用的节点数太多，而局部加密生成的单元奇异，刚度阵病态，故应力集中问题中很难得到精度较高的解，在不考虑应力集中、疲劳寿命的问题中，由于该单元使用节点较少、几何适应性强，被人们经常使用。 用直角坐标定义的单元：由六面体20节点单元通过节点重合退化得到。这种单元误差较大，无法求节点应力,只能求出 GAUSS 积分点的应力值，不推荐使用。 3 四面体20节点单元（三次单元），见图三。 用体积坐标定义的单元，单元内位移插值函数为完全三次多项式，即含常数项和Z Y X ,,, YZ XZ XY Z Y X ,,,,,222，XYZ Y Z X Z Z Y X Y Z X Y X Z Y X ,,,,,,,,,222222333二十项， 在单元内部，应力和应变为完全二次多项式，位移和应力收敛速度都很快，精度最高、几何适应性强，在应力集中、疲劳寿命分析问题中使用是非常有用和令人放心的单元。 4 三棱柱6节点单元（一次单元），见图四。 与四面体4 节点单元类似。
Matlab-PDE工具箱有限元法求解偏微分方程
在科学技术各领域中，有很多问题都可以归结为偏微分方程问题。在物理专业的力学、热学、电学、光学、近代物理课程中都可遇见偏微分方程。 偏微分方程，再加上边界条件、初始条件构成的数学模型，只有在很特殊情况下才可求得解析解。随着计算机技术的发展，采用数值计算方法，可以得到其数值解。 偏微分方程基本形式 而以上的偏微分方程都能利用PDE工具箱求解。 PDE工具箱 PDE工具箱的使用步骤体现了有限元法求解问题的基本思路，包括如下基本步骤： 1) 建立几何模型 2) 定义边界条件 3) 定义PDE类型和PDE系数 4) 三角形网格划分
5) 有限元求解 6) 解的图形表达 以上步骤充分体现在PDE工具箱的菜单栏和工具栏顺序上，如下 具体实现如下。 打开工具箱 输入pdetool可以打开偏微分方程求解工具箱，如下 首先需要选择应用模式，工具箱根据实际问题的不同提供了很多应用模式，用户可以基于适
当的模式进行建模和分析。 在Options菜单的Application菜单项下可以做选择，如下 或者直接在工具栏上选择，如下 列表框中各应用模式的意义为： ① Generic Scalar：一般标量模式（为默认选项）。 ② Generic System：一般系统模式。 ③ Structural Mech.，Plane Stress：结构力学平面应力。 ④ Structural Mech.，Plane Strain：结构力学平面应变。
⑤ Electrostatics：静电学。 ⑥ Magnetostatics：电磁学。 ⑦ Ac Power Electromagnetics：交流电电磁学。 ⑧ Conductive Media DC：直流导电介质。 ⑨ Heat Tranfer：热传导。 ⑩ Diffusion：扩散。 可以根据自己的具体问题做相应的选择，这里要求解偏微分方程，故使用默认值。此外，对于其他具体的工程应用模式，此工具箱已经发展到了Comsol Multiphysics软件，它提供了更强大的建模、求解功能。 另外，可以在菜单Options下做一些全局的设置，如下 l Grid：显示网格 l Grid Spacing…：控制网格的显示位置 l Snap：建模时捕捉网格节点，建模时可以打开 l Axes Limits…：设置坐标系围 l Axes Equal：同Matlab的命令axes equal命令 建立几何模型 使用菜单Draw的命令或使用工具箱命令可以实现简单几何模型的建立，如下 各项代表的意义分别为
平面三角形单元有限元程序设计
平面三角形单元有限元程序设计 P 9 m 9 m 一、题目 如图1所示,一个厚度均匀的三角形薄板,在顶点作用沿板厚方向均匀分布的竖向载荷。已知:P=150N/m,E=200GPa,=0、25,t=0、1m,忽略自重。试计算薄板的位移及应力分布。 要求: 1.编写有限元计算机程序,计算节 点位移及单元应力。(划分三角形单元,单元数不得少于30个); 2.采用有限元软件分析该问题(有 限元软件网格与程序设计网格必须一致),详细给出有限元软件每一步的操作过程,并将结果与程序计算结果进行对比(任选取三个点,对比位移值); 3.提交程序编写过程的详细报告及计算机程序; 4.所有同学参加答辩,并演示有限元计算程序。 有限元法中三节点三角形分析结构的步骤如下: 1)整理原始数据,如材料性质、荷载条件、约束条件等,离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。 2)单元分析,建立单元刚度矩阵。 3)整体分析,建立总刚矩阵。 4)建立整体结构的等效节点荷载与总荷载矩阵 5)边界条件处理。 6)解方程,求出节点位移。 7)求出各单元的单元应力。 8)计算结果整理。
一、程序设计 网格划分 如图,将薄板如图划分为6行,并建立坐标系,则 刚度矩阵的集成 建立与总刚度矩阵等维数的空矩阵,已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律,继而通过循环计算求得每个单元对应的节点序号。 通过循环逐个计算:(1)每个单元对应2种单元刚度矩阵中的哪一种; (2)该单元对应总刚度矩阵的那几行哪几列 (3)将该单元的单元刚度矩阵加入总刚度矩阵的对应行列 循环又分为3层循环:(1)最外层:逐行计算 (2)中间层:该行逐个计算 (3)最里层:区分为第 奇/偶 数个计算 X Y P X Y P 节点编号 单元编号
一用三结点三角形平面单元计算平面结构的应力和位移讲解
一：用三结点三角形平面单元计算平面结构的应力和位移。 1，设计说明书 计算简图，网格划分，单元及结点的编号如下图所示。由于结构对称，去四分之一结构分析。其中 E=2e10pa,mu=0.167,h=1m. 变量注释： Node ------- 节点定义 gElement ---- 单元定义
gMaterial --- 材料定义，包括弹性模量，泊松比和厚度 gBC1 -------- 约束条件 gNF --------- 集中力 gk------------总刚 gDelta-------结点位移 子程序注释： PlaneStructualModel ———定义有限元模型 SolveModel ———————求解有限元模型 DisplayResults ——————显示计算结果 k = StiffnessMatrix( ie )———计算单元刚度 AssembleStiffnessMatrix( ie, k )—形成总刚 es = ElementStress( ie )————计算单元应力 function exam1 % 输入参数：无 % 输出结果：节点位移和单元应力 PlaneStructualModel ; % 定义有限元模型 SolveModel ; % 求解有限元模型 DisplayResults ; % 显示计算结果 return ; function PlaneStructualModel % 定义平面结构的有限元模型 % 输入参数：无 % 说明： % 该函数定义平面结构的有限元模型数据： % gNode ------- 节点定义 % gElement ---- 单元定义 % gMaterial --- 材料定义，包括弹性模量，泊松比和厚度% gBC1 -------- 约束条件 % gNF --------- 集中力 global gNode gElement gMaterial gBC1 gNF % 节点坐标 % x y gNode = [0.0, 2.0 % 节点1 0.0, 1.0 % 节点2 1.0, 1.0 % 节点3 0.0, 0.0 % 节点4 1.0, 0.0 % 节点5 2.0, 0.0] ; % 节点6 % 单元定义 % 节点1 节点2 节点3 材料号 gElement = [3, 1, 2, 1 % 单元1
平面三角形单元有限元程序设计
平面三角形单元有限元程序设计
P 9 m 9 m 一、题目 如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。已知：P=150N/m，E=200GPa，=0.25，t=0.1m，忽略自重。试计算薄板的位移及应力分布。 要求： 1.编写有限元计算机程序，计算节点位移及单元应力。（划分三角形 单元，单元数不得少于30个）； 2.采用有限元软件分析该问题（有限元软件网格与程序设计网格必 须一致），详细给出有限元软件每一步的操作过程，并将结果与程序计算结果进行对比（任选取三个点，对比位移值）； 3.提交程序编写过程的详细报告及计算机程序； 4.所有同学参加答辩，并演示有限元计算程序。 有限元法中三节点三角形分析结构的步骤如下： 1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析，建立单元刚度矩阵。 3）整体分析，建立总刚矩阵。 4）建立整体结构的等效节点荷载和总荷载矩阵 5）边界条件处理。 6）解方程，求出节点位移。 7）求出各单元的单元应力。 8）计算结果整理。 一、程序设计 网格划分 如图，将薄板如图划分为6行，并建立坐标系，则 X Y P X Y P
边界约束的处理：划0置1法 适用：这种方法适用于边界节点位移分量为已知(含为0)的各种约束。 做法： （1）将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1，其他元素改为零；同 时将载荷列阵｛R｝中相应元素用已知位移置换。 ◎这样，由该方程求得的此位移值一定等于已知量。 （2）将〔K〕中已知位移相应的列的非主对角成元素也置0，以保持〔K〕的对称性。 ◎当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为 保证方程成立，须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。◎若约束为零位移约束时，此步则可省去。 特点： （1）经以上处理同样可以消除刚性位移(约束足够的前提下)，去掉未知约束反力。 （2）但这种方法不改变方程阶数，利于存贮。 （3）不过，若是要求出约束反力，仍要重新计算各个划去的总刚元素。 程序如下： 变量说明 NNODE 单元节点数 NPION 总结点数 NELEM 单元数 NVFIX 受约束边界点数 FIXED 约束信息数组 NFORCE 节点力数 FORCE 节点力数组
7a 六节点三角形单元
六节点三角形单元 单元节点位移向量： T e v u v u v u v u v u v u } {}{665544332211=δ 单元内任意一点的位移是单元节点位移的插值函数 已知u(x,y)和v(x,y)在六个顶点的函数值，可分别设u(x,y)和v(x,y)为具有六个待定系数的 插值多项式 2 652 432126524321y B xy B x B y B x B B v y A xy A x A y A x A A u +++++=+++++= 将已知条件代入，可解得其中的待定系数 u(x,y)和v(x,y)也可用Lagrenge 插值公式表示为： ∑=+++++==6 1665544332211k k k u N u N u N u N u N u N u N u ∑=+++++==61 665544332211k k k v N v N v N v N v N v N v N v 写成矩阵形式： ???????? ? ?????????? ?????????????????????????=??????66554433221165 4 3 2 1 654321v u v u v u v u v u v u 0 000000N N N N N N N N N N N N v u
缩写为： {}e ][δδN v u =? ?? ???= 单元内任意一点的应变 x v y u y v x u xy y x ??+ ??=??=??= γεε 矩阵形式： {}{}e N H v u H x y y x δε]][[][v u 0 0=??????=??? ???????? ?? ? ???????????? ???? = }[B]{}{ e δε=∴ []12 365 4 3 2 1 65 4 3 2 1 654321 0000000 0 0]][[][?=?? ??????????? ? ???????????????? ==B B B B B B N N N N N N N N N N N N x y y x N H B []?? ?? ??? ?????????????????=x N y N y N x N B k k k k k 0 0 k = 1,2,3...6 N k (k=1,2,3...6)为六节点三角形单元的形函数，N k 的表达式： 26524321),(y a xy a x a y a x a a y x N k k k k k k k +++++= k =1,2,3 (6) k i y x N y x N i i k k k k ≠== 0),(1 ),( y a x a a y N y a x a a x N k k k k k k k k 6535422 2 ++=??++=??∴ 单元内任意一点的应力：}]{][[}]{[}{e B D D δεσ==