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

有限元程序设计

-1-

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

程序设计

1、程序功能及特点

a.该程序采用四边形4节点等参单元，能解决弹性力学的平面应力应变问题。 b.前处理采用网格自动划分技术，自动生成单元及结点信息。 b.能计算受集中力、自重体力、分布面力和静水压力的作用。 c.计算结点的位移和单元中心点的应力分量及其主应力。 d.后处理采取整体应力磨平求得各个结点的应力分量。 e.算例计算结果与ANSYS计算结果比较，并给出误差分析。 f.程序采用Visual Fortran 5.0编制而成。

2、程序流程及图框

启动输入原始数据自动划分网格形成MA，计算N，NH，MX形成整体刚度矩阵K形成荷载列向量RLU分解K=LU回代并求得结点位移输入结点位移计算单元应力及主应力等整体应力磨平结点应力停机

图2-１ 程序流程图

-2-

MAINPROGRAMINPUTHUAFENCBANDSKOSTIFCONCRDECOPFOBASTRESSOUTDISTREBODYRSUMSTRSSUMSFACERFDNXFUN8GAUSSSTRESS

图2-２子程序框图

其中，各子程序的主要功能为： INPUT――输入原始数据

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

CONCR――计算集中力引起的等效结点荷载{R}e BODYR――计算自重体力引起的等效结点荷载{R} FACER――计算分布面力引起的等效结点荷载{R} DECOP――支配方程LU三角分解 FOBA――LU分解直接解法中的回代过程 OUTDISP――输出结点位移分量 STRESS――计算单元应力分量 OUTSTRE――输出单元应力分量 STIF――计算单元刚度矩阵

ee??NiFDNX――计算形函数对整体坐标的导数???xFUN8――计算形函数及雅可比矩阵[J] SFUN ――应力磨平-单元下的‘K’=NCN‘ SCN――应力磨平-单元下的右端项系数‘CN‘

?Ni?，i?1,2,3,4。 ?y??TSUMSKN――应力磨平-单元下的右端项集成到总体的‘P‘

-3-

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——面力批号

-4-

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的结果，写出所有插值函数：

111111N1?(1??)(1??)?N5?N8N2?(??1)(1??)?N5?N6422; 422 111111N3?(??1)(1??)?N6?N7N4?(1??)(??1)?N7?N8422 422 11N5?(1??2)(1??)N6?(1??2)(??1)22 11N7?(1??2)(??1)N8?(1??2)(1??)22

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

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

-5-

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

?v?移函数表示成结点位移的函数，即u??Nu ,v??Nviii?1i?188ii，写成矩阵的形式为：

?u??Nu?????1?v??00N1N200...N80N2...?u1??v??1??u2?0????v2???IN1N8???...????u8??v??8?IN2?a1??a?...IN8??2?

?...????a8?u??N1N2...N8?ae?Nae

eN称为插值函数矩阵或形函数矩阵，a为单元结点位移列阵。

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

??x???????y??Lu?LNae?L?N1??xy????????x?Bi?LNi??0??????y?0?????Ni?0?y???????x?N2...N8?ae??B1B2...B8?ae?Bae

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

??Ni???x0????0?Ni????Ni???y?0?????N1?0?y???????x??0???Ni???y??Ni???x?，

B??B1B2?????x?...B8???0??????y0N1N200...N80N2...0?N8??

-6-

??u???????x???x??x??????v??????y??????0?y????xy??????u?v?????????y?x???y

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

?0?????u???v? ?y???????x???x???????y??D??DBae?Sae??xy???

其中

S?DB?D?B1B...2B8???S1S...2S?8，S称为应力矩阵，B是应变矩

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

???1?0??E?D??10??1??2?1????00??2?。

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

4.3、等参元变换：

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

x?f????y?g????,

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

x??Ni'xii?0m,

y??Ni'yii?0m

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

?segdS???seg(x,y)dS

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

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

-7-

设x??N'x, y??N'yii88ii，xi，yi，zi是结点在总体坐标内的坐标值，N’i为

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

??Ni?Ni?x?Ni?y?????x.????y.???函数Ni对ξ的偏导数可以表示成： ?

??Ni??Ni.?x??Ni.?y??x???y???????Ni???x??????????集合成矩阵形式就是： ???Ni???x??????????8?y???Ni???Ni???x??????x??.???J?? ?y???Ni???Ni?????????y?????y?8iiii式中的J称为雅克比矩阵，利用x?函数，即：

?Nx, y??Nyi?1i?1，J可以表示为自然坐标的

?8???Nixi?i?1??x,y????J???8???,?????Nixi?i?1??????Ny?ii???Ni?1??1??????????N8?Niyi??1?????i?1????8?N2???N2???N8??x1...?????x2??N8??......??????x8y1?y2??...??y8?,

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

??Ni???Ni???????x??1????J.???Ni???Ni????????y????

其中J是J的逆矩阵,它可以按照下式计算得到：J?1?*Jij?11*J,J是J的伴随矩阵，J它的元素是J的元素

Jij的代数余子式。

4.4、总体应力磨平

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

M??e?1Ve?)TCN*dV?0,(i?1,2,...,N)(?*??i

其另一种矩阵形式为：

-8-

?)TCN(?*)TNTCN?(?

令

Ke??Ni*?CNi*T

?Pe??a??i*??

TCNi*

并利用以下表达式，

?x*?(N1N2N3N4)(?*1x?*2x?*3x?*4x)T

?y*?(N1?xy*?(N1和

N2N2N3N3N4)(?*1yN4)(?*1xy?*2y?*3y?*4y)T?*2xy?*3xy?*4xy)T1?v/(1?v)0????C??1?v2?/E??v/(1?v)10??002/(1?v)???

以及

?N1?Ni*????N2N1N1N2N2N3N3N3N4N4???N4??

，只不过，这里有Ke为

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

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

5、算例

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

2

kN/m，不考虑自重体力。

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

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

2

-9-

长度单位：

图5-1 算例模型

图5-2 ANSYS计算模型图

实际计算结果与ANSYS分析结果的比较（详细计算结果见后面）： （1） 位移比较 比较项目 本程序计算结果 ANSYS分析结果 误差 X方向最大位移 0.0863mm 0.1200mm 28.0% Y方向最大位移 5.890mm 6.000mm 1.83%

-10-

&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 READ(1,70) TL 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,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为原始数据输入子程序

-16-

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,70)TL

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

10 FORMAT(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 READ(1,70)TL 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 READ(1,70) TL 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

-17-

READ(1,70) TL 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 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

-18-

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)

-19-

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

-20-

YMAX=U(2) JNYMAX=I ENDIF

WRITE(*,'(5X,I5,10X,2E14.3)') I,U WRITE(2,'(5X,I5,10X,2E14.3)') I,U ENDDO

650 FORMAT(/25X,'结点位移计算结果'/8X,'结点',13X,'X向位移',8X,'Y向位移') WRITE(2,651)XMAX,JNXMAX,YMAX,JNYMAX

651 FORMAT('X向最大位移＝',E14.3,5X,'所在结点号＝',I4/'Y向最大位移＝

&',E14.3,5X,'所在结点号＝',I4)

RETURN END

!************************************************************** ！ !输出单元应力分量子程序

! ************************************************************* SUBROUTINE OUTSTRE(NE,STRE) DIMENSION STRE(3,*),ST(6) REAL STREMAX STREMAX=0.0 NESTREMAX=0 WRITE(*,700) WRITE(2,700)

DO IE=1,NE !对单元循环 SX=STRE(1,IE) SY=STRE(2,IE) SXY=STRE(3,IE) ST(1)=SX ST(2)=SY ST(3)=SXY H1=SX+SY

H2=SQRT((SX-SY)*(SX-SY)+4.0*SXY*SXY) ST(4)=(H1+H2)/2.0 !计算最大应力 ST(5)=(H1-H2)/2.0 !计算最小应力 IF(ABS(SXY).LT.1.0E-4)THEN

IF (SX.GT.SY) ST(6)=0.0 !计算最大最小应力夹角 IF (SX.LE.SY) ST(6)=90.0 ELSE

ST(6)=ATAN((ST(4)-SX)/SXY)*57.29578 ENDIF

WRITE(*,'(6X,I4,3X,6F11.3)') IE,ST WRITE(2,'(4X,I4,5X,6F11.3)') IE,ST IF(ABS(STREMAX).LT.ABS(ST(4)))THEN STREMAX=ST(4) NESTREMAX=IE

-26-

ENDIF ENDDO 700 FORMAT(/30X,'单元应力'/5X,'单元号',8X,'X向应力',3X,'Y向应力',2X,'XY剪

&切应力',1X,'最大应力',1X,'最小应力',8X,'角度')

WRITE(2,701)STREMAX,NESTREMAX

701 FORMAT('最大单元应力值＝',F11.3/'所在单元号＝',I5) RETURN END

!************************************************************** ！ 计算单元应力子程序

! ************************************************************* SUBROUTINE STRESS(COOR,MEL,JR,AE,R,STRE)

DIMENSION XYZ(2,4),DNX(2,4),AE(4,*),STRE(3,*),COOR(2,*),MEL(5,*),

&JR(2,*),RJAC(2,2),SIG(3),B(3,8),R(*),iven(4)

COMMON /CMN1/ NP,NE,NM,NR

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

COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2) DO 106 IE=1,NE !对每个单元循环

NME=MEL(5,IE) !取IE号单元的材料类型号 DO 300 K=1,4 !对单元四个结点循环 IEK=MEL(K,IE) !取单元的结点号

DO 310 M=1,2 !对结点的两个自由度循环 310 XYZ(M,K)=COOR(M,IEK) !取结点的坐标

DO 320 M=1,2 !对结点的两个自由度循环 JRR=2*(K-1)+M !取单元的自由度号

320 NN(JRR)=JR(M,IEK) !取该自由度的总体自由度号 300 CONTINUE

E=AE(1,NME) !取材料的弹性模量 U=AE(2,NME) !取材料的泊松比

D1=E*(1.00-U)/((1.00+U)*(1.00-2.00*U)) !求应力应变关系矩阵系数 D2=E*U/((1.00+U)*(1.00-2.00*U)) D3=0.50*E/(1.00+U) SS=0.0 RR=0.0

CALL FDNX (XYZ,DNX,DET,RR,SS,RJAC,IVEN,IE) !计算形函数对整体坐标的导

&数，存放在DNX数组中

DO 30 I=1,4 !对单元四个结点循环 II=2*(I-1) J1=II+1 J2=II+2

BI=DNX(1,I) CI=DNX(2,I)

B(1,J1)=BI !计算应变位移关系矩阵系数 B(2,J1)=0.

-27-

B(3,J1)=CI B(1,J2)=0. B(2,J2)=CI B(3,J2)=BI 30 CONTINUE

DO 55 II=1,3 SIG(II)=0.00 55 CONTINUE DO 70 K=1,8 NA=NN(K) !取总体自由度号

IF (NA.EQ.0) GOTO 70 !如果NA＝0，该结点自由度有约束，位移为0，跳过 DO 60 L=1,3

SIG(L)=SIG(L)+B(L,K)*R(NA) !计算应变 60 CONTINUE 70 CONTINUE

SX=D1*SIG(1)+D2*SIG(2) !计算应力 SY=D2*SIG(1)+D1*SIG(2) SXY=D3*SIG(3)

STRE(1,IE)=SX !将应力值赋值给应力数组 STRE(2,IE)=SY STRE(3,IE)=SXY 106 CONTINUE

CALL OUTSTRE(NE,STRE) !输出应力结果 RETURN END

!************************************************************** ！ 计算集中力引起的等效结点荷载子程序

! ************************************************************* SUBROUTINE CONCR(NCP,R,JR)

DIMENSION R(*),JR(2,*),XYZ(2)!定义荷载列阵，自由度集中力引起的荷载分量 DO 100 I=1,NCP

READ (1,*) IP,PX,PY !读入结点信息 XYZ(1)=PX XYZ(2)=PY DO 95 J=1,2

L=JR(J,IP) !判断结点自由度 IF(L.EQ.0) GO TO 95

R(L)=R(L)+XYZ(J) !自由度为1，将集中力分量叠加到荷载列阵 95 CONTINUE 100 CONTINUE

!WRITE(2,*)(R(I),I=1,NH) !90 FORMAT('JIEDIANLI',F10.2) RETURN END

-28-

!************************************************************** ！ 计算自重体力引起的等效结点荷载子程序

! ************************************************************* SUBROUTINE BODYR(R,MEL,COOR,JR,AE)

DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),AE(4,*),XYZ(2,4),IVEN(4) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH

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

COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2) COMMON /GAUSS/ RSTG(3),H(3) H(1)=1.0 !定义高斯积分权系数 H(2)=1.0

RSTG(1)=-0.5773502691896260 !定义高斯积分点坐标 RSTG(2)=-RSTG(1) DO 10 IE=1,NE

DO I=1,8 !定义单元结点体力分量 RF(I)=0.00 ENDDO NEE=IE

NME=MEL(5,NEE) !读入单元材料信息 GAMMA=AE(3,NME) !引入材料重度 DO 75 K=1,4

IEK=MEL(K,NEE) !引入单元结点号 IVEN(k)=IEK DO 95 M=1,2 JJ=2*(K-1)+M

NN(JJ)=JR(M,IEK) !将整体结点自由度值赋于单元局部结点自由度 95 XYZ(M,K)=COOR(M,IEK) !整体结点坐标赋于单元局部结点坐标 75 CONTINUE

DO 99 IS=1,2 !定义高斯积分点和权系数 S=RSTG(IS) SH=H(IS)

DO 98 IR=1,2 RR=RSTG(IR) RH=H(IR)

CALL FUN8 (XYZ,RR,S,DET) !调用函数FUN8计算形函数和雅各比行列式 DO 30 I=1,4 J=2*I

RF(J)=RF(J)-FUN(I)*RH*SH*DET*GAMMA !体积力计算公式计算单元结点y向分量 30 CONTINUE 98 CONTINUE 99 CONTINUE

CALL ASLOAD (R) !调用ASLOAD,将自重引起的荷载分量叠加到荷载列阵 10 CONTINUE

-29-

RETURN END

!************************************************************** ！ 计算面力引起的等效结点荷载子程序

! ************************************************************* SUBROUTINE FACER(IEW,NSE,R,MEL,COOR,JR,WG)

DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),wg(*), &XYZ(2,4),IEW(*),PR(2)

COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH

COMMON /CMN3/ RF(8),SKE(8,8),NN(8) COMMON /CMN4/ NEE,NME

COMMON /GAUSS/ RSTG(3),H(3) H(1)=1.0 !定义高斯积分权系数 H(2)=1.0

RSTG(1)=-0.5773502691896260 !定义高斯积分点坐标 RSTG(2)=-RSTG(1)

NWF=0 !引入参数判断面力类型 NNF=0

IR=WG(1)+0.1 !引入面力类型值

if(IR.EQ.1)NWF=1 !判断受静水压力的面力类型

if(IR.EQ.2)NNF=1 !判断受均布法向压力的面力类型 DO 510 IE=1,NSE DO I=1,8

RF(I)=0.00 !定义单元结点面力分量 ENDDO

NEE=IEW(IE) !受面力作用的单元的单元号 DO 575 K=1,4

IEK=MEL(K,NEE) !引入单元结点号 DO 595 M=1,2 JJ=2*(K-1)+M

NN(JJ)=JR(M,IEK) !将整体结点自由度值赋于单元局部结点自由度 595 XYZ(M,K)=COOR(M,IEK) !整体结点坐标赋于单元局部结点坐标 575 CONTINUE

IF(NWF.EQ.1) then !面力是静水压力类型 GAMA=WG(2) !水压密度 Z0=WG(3) !最高水位的y坐标 NSU=WG(4)+0.1 !

CALL SURLOD (NSU,XYZ,PR,Z0,GAMA,1) !调用SURLOD，计算结点面力分量 endif

IF(NNF.EQ.1) then !面力是均布法向压力类型 q=WG(2) !均布值 NSU=WG(4)+0.1 ! do j=1,2

-30-

本文来源：https://www.bwwdw.com/article/r2ih.html