最新平面四边形4结点等参有限单元法

有限元程序设计

平面四边形4结点等参有限单元法

程序设计

1、程序功能及特点

a.该程序采用四边形4节点等参单元，能解决弹性力学的平面应力应变问题。

b.前处理采用网格自动划分技术，自动生成单元及结点信息。

b.能计算受集中力、自重体力、分布面力和静水压力的作用。

c.计算结点的位移和单元中心点的应力分量及其主应力。

d.后处理采取整体应力磨平求得各个结点的应力分量。

e.算例计算结果与ANSYS计算结果比较，并给出误差分析。

f.程序采用Visual Fortran 5.0编制而成。

2、程序流程及图框

图2-１程序流程图

图2-２子程序框图

其中，各子程序的主要功能为：

INPUT――输入原始数据

HUAFEN――自动网格划分，形成COOR(2,NP),X,Y的坐标值与单元信息CBAND――形成主元素序号指示矩阵MA（*）

SKO――形成整体刚度矩阵[K]

CONCR――计算集中力引起的等效结点荷载{R}e

BODYR――计算自重体力引起的等效结点荷载{R}e

FACER――计算分布面力引起的等效结点荷载{R}e

DECOP――支配方程LU三角分解

FOBA――LU分解直接解法中的回代过程

OUTDISP――输出结点位移分量

STRESS――计算单元应力分量

OUTSTRE――输出单元应力分量

STIF――计算单元刚度矩阵

FDNX――计算形函数对整体坐标的导数

，=

i1,2,3,4。

FUN8――计算形函数及雅可比矩阵[J]

SFUN ――应力磨平-单元下的‘K’=NCN‘

SCN――应力磨平-单元下的右端项系数‘CN‘SUMSKN――应力磨平-单元下的右端项集成到总体的‘P‘

SUMSTRS――应力磨平-单元下的集成到总体的‘K‘

GAUSTRSS――高斯消元求磨平后的应力

3、输入数据及变量说明

当程序开始运行时，按屏幕提示，键入数据文件的名字。

在运行程序之前，根据程序中INPUT需要的数据输入建立一个存放原始数据的文件，这个文件的名字为INDAT.DAT。数据文件包括如下内容：

（1）总控信息。共一条，4个数据。

L，B，NNL，NNB，NM，NR

L——矩形体长度

B--矩形体宽度

NNL--L方向上划分的结点数

NNB--B方向上划分的结点数

NM—单元材料类型数

NR——约束结点总数

（2）结点约束信息。共NR条，每条依次输入：

IP，IX，IY

IP——结点号

IX、IY——分别为IP结点在x，y方向的约束情况，如果约束填0，如果自由填1。（3）材料信息。共NM条，每条依次输入：

JJ,（AE（I，JJ），I=1，4）

JJ――材料类型号，（AE（I，JJ），I=1，4）――该材料的材料参数，共4个参数，排列顺序为：弹性模量、泊松比、容重、单元厚度。

（4）荷载信息

a.荷载控制信息。共一条，3个数据

NCP，IZ

NCP——受集中力作用的结点数

IZ——面力批数

b.若NCP>0，输入

IP，PX，PY

IP——结点号

PX、PY——分别为IP结点x，y方向的集中力分量。

c.若IZ>0，输入面力荷载信息，共IZ批，按批输入：

JS，NSE，（WG（I）I=1，4）

JS——面力批号

NSE ——第JS 批面力受到面力作用的单元个数，

（WG （I ），I ＝1，4）——该面力的特征参数共4个数据，第1个数为面力类型，填1表示受静水压力作用，填2表示受均布法向压力作用；第2个数为水压密度，如果是均布压力情况，就填均布压力的集度；第3个数为最高水位的y 坐标，如果是均布压力情况，可以填任意数；第4个数为面力作用的单元面的面号，单元面号的规定见图2-3。

（IEW （M ），M ＝1，NSE ），IEW （*）为受面力作用的单元的单元号，共NSE 个。

图3-1 单元结点编码与面号

4、理论基础和求解过程 4.1、构造插值函数：

本有限元计算采取的是四边形八结点等参元进行插值计算的。直接调用教材115页3..3.21的结果，写出所有插值函数：

158111(1)(1)422N N N ξη=----; 256

111(1)(1)422N N N ξη=+---

367111(1)(1)422N N N ξη=++-- 478

111(1)(1)422N N N ξη=-+-- 251(1)(1)2N ξη=-- 261(1)(1)

2N ηξ=-+ 271(1)(1)2N ξη=-+ 281(1)(1)

2N ηξ=--

4.2位移插值函数及应变应力求解：

在有限元方法中单元的位移模式一般采用多项式作为近似函数，多项式的选取应由低次到高次的完备多项式。位移模式的选取一般为：u ＝φβ。

u u v ??

=????

，φ为位移模式，β为广义坐标向量。根据方程求解得出广义坐标，可将位

移函数表示成结点位移的函数，即8

i i

i u N u

∑ ,8

i i

i v N v

∑，写成矩阵的形式为：

[]11121

282212

2888800...0 0

...0

......u v a u N N N a u u v IN IN IN N N N v a u v ????????????

??????

??=== ?????????????

??????

[]128...e e u N N N a Na ==

N 称为插值函数矩阵或形函数矩阵，e

a 为单元结点位移列阵。

确定了单元位移后，可以很方便地利用几何方程和物理方程求得单元的应变和应力。单元应变为：

[][]1

28128......x e e e e

y xy Lu LNa L N N N a B B B a Ba εεεγ????

======??????

B 称为应变矩阵，L 是平面问题的微分算子，其中：

0000

00i i

i i i i i i N x x N N B LN N y y N N y x y x ??????????????????

?????===???????????

?????

??????????????????，

[]1

281

28000...0 0

...0

x N N N B B B B N N N y y

x ?????????

????==???????

????????????

00x y xy u x x u v v y y u v y x y

x εεεγ?????????

????

???????????

????===????????????

????????????

?+????????????

根据物理方程可以求得单元应力，

x e e

y xy D DBa Sa σσσετ????

====??????

其中

[][]

128128......S DB D B B B S S S ===，S 称为应力矩阵，B 是应变矩

阵。由于是平面应力问题，E0和v0取为E 和v ，所以弹性矩阵

2101011002E D υυυυ?? ?

?= ?

- ???。

这部分内容参考了教材第2.2节。

4.3、等参元变换：

为了将局部（自然）坐标中几何形状规则的单元转换成总体坐标中几何形状扭曲的单元，以满足对一般形状求解域进行离散化的需要，必须建立坐标转换：

()x f ξη=,

()

y g ξη=

最方便的方法是将上式中表示成插值函数的形式，

i i

i x N x ==∑,

i i

i y N y ==∑

在有限元分析中，为建立求解方程，需要进行各个单元面积内的积分，其一般形式为：

(,)e

e s s gdS g x y dS

而g 中包含场函数对于总体坐标x ，y 的导数。

由于场函数是用自然坐标表述的，又因为在自然坐标内的积分限是规格化的，因此希望能在自然坐标内按规格化的数值积分方法进行上述积分。为此需要建立两个坐标系内导数之间的变换关系。

设8

i x N x ==

∑, 8

i y N y

==∑，xi ，yi ，zi 是结点在总体坐标内的坐标值，N ’i 为

形状函数，实际上它也是用局部（自然）坐标表示的插值函数。对于等参变换，结点数、结点号和插值函数都不变。

函数Ni 对ξ的偏导数可以表示成： ....i i i i i i N N N x y x y N N N x y x y ξξξηηη??????=+???????

???????=+???????

集合成矩阵形式就是： .i i i i i i N x

y N N x x J N N N x y y y ξξ

ξηηη?????????????

??????

????????

???????==?????????????

???????????????

???

??? 式中的J 称为雅克比矩阵，利用8

i i

i x N x ==∑, 8

i i

i y N y

==∑，J 可以表示为自然坐标的

函数，即：

()()8

111

81211

88...,......,...i i i i i i i i i i i i N x N y x y N N N x y x y J N N N N x N y x y ξξξ

ξξξηηη

ηηη====??

????

????????

???????????????≡=

=???????????????

????????????

??????

∑∑∑∑,

Ni 对于x,y 的偏导数可以用自然坐标显示地表示为:

.i i i i N N x J N N y ξη-??????

????

???

???=??????

??????????

其中1

J -是J 的逆矩阵,它可以按照下式计算得到：1*

1J J J

,J 是J 的伴随矩阵，它的元素*

J 是J 的元素

J 的代数余子式。

4.4、总体应力磨平

根据《有限单元法》第5章中的（5.3.15），即

°**1

()0,(1,2,...,)M

T i

e CN dV i N σσ=-==∑?

其另一种矩阵形式为：

°*()()T T T N CN CN σσ=

令

()**

e i i K N CN =

°()

e i P CN σ=

i a σ=

并利用以下表达式，

*****1

41234()()T x x

x x x N N N N σσσσσ=

*****12341234()()T

y y y y y N N N N σσσσσ= *****12341234()()T

xy xy xy xy xy N N N N σσσσσ=

和

()21/(1)01//(1)10002/(1)v v C v E v v v --?? ?

=--- ?

?-??

以及

4i N N N N N N N N N N N N N ?? ?=

? ??

就可以像单元元刚度矩阵到总体刚度矩阵一样，求出

i a σ=，只不过，这里有Ke 为

12*12的矩阵，而“总刚”K 为3*NP 阶的矩阵，NP 为结点数。

5、算例

应用本程序计算一个矩形土体受到均布荷载时的位移和应力，如下图所示，三面约束的土体尺寸为40m*10m ，取一半为20m*10m ，E=1.0×104 kN/m 2

，0.3μ=，在右端受均布荷载10.0q =

kN/m 2

，不考虑自重体力。

给定网格大小为0.20.5m m ?，有限元网格自动划分如图3-2所示，单元总数2000，节点总数2121。由于对称性，右端约束为滑移支座，只限制x 方向位移。

虽然土体一般不为弹性，但是方法是一样的，只是刚度矩阵改动即可使用弹塑性模型。

长度单位：

图5-1 算例模型

图5-2 ANSYS计算模型图

实际计算结果与ANSYS分析结果的比较（详细计算结果见后面）：

比较项目X方向最大位移Y方向最大位移

本程序计算结果0.0863mm 5.890mm

ANSYS分析结果0.1200mm 6.000mm

误差28.0% 1.83%

（2）应力比较

误差分析：

本程序计算得到的Y方向最大位移和ansys计算结果很相近，由于x方向上的位移不占主要部分，因此，其误差虽然有些大，但对总体位移影响不大。

由于网格较密，因此，计算得到的单元应力值（未磨平）与ansys结果相近，误差小于1%，不必要应力磨平也可以达到较好的精度。

本程序可以进行总体应力磨平，但是由于网格数较多，因此，磨平矩阵阶数较大，可能计算时间也很长，甚至无法计算。

（3）实际计算结果MATLAB图与ANSYS图比较：

(a)

(b)

图5-3 X方向应力图，（a)ANSYS （b) MATLAB

图5-4 主应力计算结果的MATLAB 图

图5-5 主应力计算结果的ANSYS 图

图5-6 剪应力计算结果的MATLAB 图

图5-7 剪应力计算结果的ANSYS 图

图5-8 ANSYS总位移图

图5-9 ANSYS X方向位移图

图5-10 ANSYS Y方向位移图

图5-11 ANSYS Y方向应力图

附录：

（1）输入文件数据：

平面四边形四结点单元输入数据

L B NNL NNB NM NR

20.0 10.0 101 21 1 141

********************************************************* 约束信息数据

结点号 X向约束信息 Y向约束信息

1 0 0

2 0 0

3 0 0

4 0 0

5 0 0

6 0 0

7 0 0

8 0 0

9 0 0

10 0 0

11 0 0

12 0 0

13 0 0

14 0 0

15 0 0

16 0 0

17 0 0

18 0 0

19 0 0

20 0 0

21 0 0

22 0 0

43 0 0

64 0 0

85 0 0

106 0 0

127 0 0

148 0 0

169 0 0

190 0 0

211 0 0

232 0 0

253 0 0

274 0 0

295 0 0

316 0 0

337 0 0

358 0 0

379 0 0

400 0 0

421 0 0

442 0 0

463 0 0

484 0 0

505 0 0

526 0 0

547 0 0

568 0 0

589 0 0

610 0 0

631 0 0

652 0 0

673 0 0

694 0 0

715 0 0

736 0 0

757 0 0

778 0 0

799 0 0

820 0 0

841 0 0

862 0 0

883 0 0

904 0 0

925 0 0

946 0 0

967 0 0

988 0 0 1009 0 0 1030 0 0 1051 0 0 1072 0 0 1093 0 0 1114 0 0 1135 0 0 1156 0 0 1177 0 0 1198 0 0

1240 0 0 1261 0 0 1282 0 0 1303 0 0 1324 0 0 1345 0 0 1366 0 0 1387 0 0 1408 0 0 1429 0 0 1450 0 0 1471 0 0 1492 0 0 1513 0 0 1534 0 0 1555 0 0 1576 0 0 1597 0 0 1618 0 0 1639 0 0 1660 0 0 1681 0 0 1702 0 0 1723 0 0 1744 0 0 1765 0 0 1786 0 0 1807 0 0 1828 0 0 1849 0 0 1870 0 0 1891 0 0 1912 0 0 1933 0 0 1954 0 0 1975 0 0 1996 0 0 2017 0 0 2038 0 0 2059 0 0 2080 0 0 2101 0 0 2102 0 1

2104 0 1

2105 0 1

2106 0 1

2107 0 1

2108 0 1

2109 0 1

2110 0 1

2111 0 1

2112 0 1

2113 0 1

2114 0 1

2115 0 1

2116 0 1

2117 0 1

2118 0 1

2119 0 1

2120 0 1

2121 0 1

*********************************************************

材料信息数据

类型号弹性模量泊松比容重单元厚度

1 10000.0 0.3 0.0 1.0

*********************************************************

荷载信息数据

受集中力作用的结点数面力批数

0 1

集中力数据

受集中力作用的结点号 X向集中力分量 Y向集中力分量

*********************************************************

面力数据

面力批号受面力单元个数面力类型面力集度最大集度面力作用面号

1 15

2 -10 -10 4

*********************************************************

面力作用的单元的单元号

2000 1980 1960 1940 1920 1900 1880 1860 1840 1820 1800 1780 1760 1740 1720

（3）源程序：

PROGRAM FEM

DIMENSION SK(200000),COOR(2,10000),AE(4,100),MEL(5,10000),WG(4), &JR(2,10000),MA(10000),R(10000),IEW(10000),STRE(3,10000),

&TOALFUN(10000,10000),SUMS(10000),GOODSTR(10000) COMMON /CMN1/ NP,NE,NM,NR

COMMON /CMN2/ N,MX,NH

COMMON /CMN3/ RF(8),SKE(8,8),NN(8)

OPEN(1,FILE='INDAT.DAT',STATUS='OLD')

OPEN(2,FILE='OUT.DAT',STATUS='NEW')

CALL INPUT(XL,B,IX,IY,AE,NCP,IZ,NNL,NNB,JR) !输入数据 CALL HUAFEN(XL,B,NNL,NNB,COOR,MEL)

CALL CBAND (MA,JR,MEL) !形成主元素序号指示矩阵MA

DO I=1,NH

SK(I)=0.0

ENDDO

CALL SK0(SK,MEL,COOR,JR,MA,AE) !形成整体刚度矩阵[K]

DO I=1,N

R(I)=0.0

ENDDO

IF(NCP.GT.0) CALL CONCR(NCP,R,JR)

CALL BODYR(R,MEL,COOR,JR,AE)

READ(1,70) TL

70 FORMAT(A)

IF(iz.GT.0)then

DO jj=1,iz

READ (1,*)JS,NSE,(WG(I),I=1,4)

READ(1,70) TL

READ(1,*)(IEW(M),M=1,NSE)

CALL FACER(IEW,NSE,R,MEL,COOR,JR,WG)

ENDDO

ENDIF

CALL DECOP (SK,MA)

CALL FOBA (SK,MA,R)

CALL OUTDISP(NP,R,JR)

CALL STRESS (COOR,MEL,JR,AE,R,STRE)

CALL GAUSTRSS(MEL,COOR,AE,TOALFUN,SUMS,GOODSTR,STRE,JR,R) !应力磨平并输出磨平后的应力值

WRITE(2,'(A)')' PROGRAM SAFF HAS BEEN ENDED'

WRITE(*,'(A)')' PROGRAM SAFF HAS BEEN ENDED' CLOSE(1)

CLOSE(2)

END

!INPUT为原始数据输入子程序

SUBROUTINE INPUT(XL,B,IX,IY,AE,NCP,IZ,NNL,NNB,JR)

DIMENSION AE(4,*),JR(2,*)

CHARACTER*400 TL

COMMON /CMN1/ NP,NE,NM,NR

COMMON /CMN2/ N,MX,NH

READ(1,70)TL

READ(1,*)XL,B,NNL,NNB,NM,NR

WRITE(2,10)XL,B,NNL,NNB,NM,NR

10FORMAT(5X,'平面四边形四结点单元有限元程序'//5X,'L=',F5.2,4X,' &B=',F5.2,4X,'NNL=',I3,4X,'NNB=',I3,4X,'NM=',I2,4X,'NR=',I3)

NP=NNL*NNB

DO 15 I=1,NP

DO 15 J=1,2

15 JR(J,I)=1

READ(1,70)TL

WRITE(2,110)

110 FORMAT(/5X,'约束信息'/5X,'结点号',5X,'X向约束',5X,'Y向约束')

DO 20 I=1,NR

READ(1,*)IP,IX,IY

WRITE(2,100)IP,IX,IY

100 FORMAT(8X,I5,5X,I2,9X,I2)

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

READ(1,70) TL

DO 55 J=1,NM

READ (1,*) JJ,(AE(I,JJ),I=1,4)

WRITE(2,910) JJ,(AE(I,JJ),I=1,4)

55 CONTINUE

910 FORMAT (/5X,'材料信息'/5X,'类型号',5X,'弹性模量',5X,'泊松比',5X,'容重& ',8X,'单元厚度'/5X,I5,4(5x,E8.3))

READ(1,70) TL

READ(1,*)NCP,IZ

WRITE(2,920)NCP,IZ

920 FORMAT (/5X,'荷载信息'/5X,'受集中力作用的结点数',8X,'面力批数&'/(2(12X,I5)))

READ(1,70) TL

70 FORMAT(A)

RETURN

END

!自动划分网格子程序，计算单元数，结点数，结点坐标，每个单元包含的结点

!单元的4个结点号与材料类型号，从左下角逆时针编号顺序；

SUBROUTINE HUAFEN(XL,B,NNL,NNB,COOR,MEL)

DIMENSION COOR(2,*),MEL(5,*)

COMMON /CMN1/ NP,NE,NM,NR

NP=NNL*NNB

NE=(NNL-1)*(NNB-1)

WRITE(2,*)'总结点数 NP=',NP

WRITE(2,*)'总单元数 NE=',NE

DL=XL/(NNL-1)

DB=B/(NNB-1)

N=1

M=1

T=NNL*NNB

DO 10 K=NNB,T,NNB

DO 20 I=N,K

IP=I

COOR(1,I)=DL*(K-NNB)/NNB

20 CONTINUE

N=K+1

10 CONTINUE

DO 30 J=NNB,T,NNB

DO 40 I=M,J

COOR(2,I)=DB*(I-M)

40 CONTINUE

M=J+1

30 CONTINUE

WRITE(2,100) (J,(COOR(I,J),I=1,2),J=1,NP)

100 FORMAT(/7X,'结点号',10X,'X坐标',8X,'Y坐标'/(8X,I5,5X,2F12.6)) N=1

DO 50 K=1,(NNL-1)

DO 60 I=N,(NNB-1)*K

NNE=I

MEL(1,I)=I+K-1

MEL(2,I)=MEL(1,I)+NNB

MEL(3,I)=MEL(1,I)+NNB+1

MEL(4,I)=MEL(1,I)+1

60 CONTINUE

N=N+(NNB-1)

50 CONTINUE

WRITE(2,111)(J,(MEL(I,J),I=1,4),J=1,NE)

111 FORMAT(/7X,'单元号',5X,'结点1',5X,'结点2',5X,'结点3',5X,'结点&4'/(7X,I5,5X,I5,5X,I5,5X,I5,5X,I5))

DO 70 IE=1,NE

MEL(5,IE)=1

70 CONTINUE

END

!**************************************************************

! LU三角分解子程序

! ************************************************************* SUBROUTINE DECOP (SK,MA) !方程LU分解

DIMENSION SK(*),MA(*)

COMMON /CMN2/ N,MX,NH

DO 50 I=2,N

L=I-MA(I)+MA(I-1)+1

K=I-1

L1=L+1

IF (L1.GT.K) GOTO 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) GOTO 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

RETURN

END

!************************************************************** ! LU三角分解后的回代求解

! ************************************************************* SUBROUTINE FOBA (SK,MA,R) !回代求解

DIMENSION SK(*),MA(*),R(*)

COMMON /CMN2/ N,MX,NH

DO 10 I=2,N

L=I-MA(I)+MA(I-1)+1

K=I-1

IF (L.GT.K) GOTO 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

II=MA(I)

45 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

55 R(J)=R(J)-SK(IJ)*R(I)

25 CONTINUE

30 CONTINUE

RETURN

END

!************************************************************** ！计算形函数对整体坐标的导数

! ************************************************************* SUBROUTINE FDNX (XYZ,DNX,DET,R,S,RJAC,IVEN,NEE) DIMENSION XYZ(2,*),DNX(2,4),RJAC(2,2),IVEN(*)

COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)

CALL FUN8 (XYZ,R,S,DET) !求形函数

IF (DET.LT.1.0E-5)THEN !判断J的大小

WRITE(2,600) NEE,R,S,det