等效结点荷载计算机语⾔，基于FORTRAN的3D等效结点荷载
计算
司马丹琪++李元松++姜成潼++何泉
摘要： 在有限单元分析中，需将作⽤于单元的外⼒按虚功等效原则，移置到结点上成为等效结点⼒。由于数学运算的困难，往往限制了⼯程技术⼈员运⽤有限单元技术分析实际具有分布⼒的⼯程问题。以常⽤空间20结点单元为例，详述荷载等效过程，附注⾃编FORTRAN程序。并利⽤FEA软件进⾏计算对⽐验证程序正确性，为学者学习温度、预应⼒等类似⼯程问题的等效结点荷载计算提供有益参考。
Abstract： Abstract： In the finite element analysis， the external force acting on the element needs to be transferred to the node as the equivalent node force according to the principle of virtual work equivalence. Due to the difficulty of mathematical operations， engineers and technicians are often limited to use finite element technique to analyze practical distributed engineering problems. Taking the commonly used space 20-node element as an example， the equivalent process of load is described in detail. The software FEA is used to calculate and compare the correctness of the program，which provides a useful reference for scholars to study the equivalent nodal load of temperature， prestress and other similar engineering problems.
关键词： 离散体；等效结点荷载；虚功原理；程序设计
Key words： discrete body；equivalent node load；virtual work principle；program design
中图分类号：TP311.1 ⽂献标识码：A ⽂章编号：1006-4311(2018)08-0233-02
0 引⾔
在弹性介质静⼒问题的计算中外加作⽤因素很多，如集中⼒、分布⼒(包括引⼒、斥⼒等场⼒、惯性⼒⼀类体积⼒和⾯积分布的由⼒边界条件给定的接触⼒)等直接载荷，也可能有因为温度改变、装配因素、预应⼒作⽤等其它⼲扰⼒[1]。这些外加作⽤因素都可遵循⼒学等效原则(如静⼒等效、位移模式下的虚功等效)处理成结点载荷。这种等效处理往往涉及较为复杂的坐标变换运算[2]，对于⾮⼒学专业的⼯程技术⼈员⽽⾔，存在⼀定困难。并且由于计算机语⾔的限制，传统的有限元程序设计课程只能以⼆维有限元问题为例介绍程序设计过程，极少涉及三维有限元的编程。⽽实际⼯程问题⽆⼀不是三维问题，因此编写三维程序更具有实际价值。现有商业计算软件中边界处理条件功能均很强⼤，但很难涵盖⼯程实际中遇到的各种边界条件问题，⼀旦遇到软件中没有对应的处理⽅法，仍需根据⼯程条件⾃⾏开发程序解决。笔者在有限元程序系统(FEM-PS)开发过程中，对3D等参单元等效结点荷载的计算公式进⾏详细推导，并利⽤FORTRAN程序设计成通⽤模块，以期对类似问题的解决有所启⽰。本⽂以三维空间⼆⼗结点等参单元为例，介绍FORTRAN语⾔在程序设计中的使⽤，并针对三维空间设计过程中的等效结点荷载计算环节做程序详述。
1 等效结点荷载计算原理
1.1 集中⼒
设在单元内部或边界上任意点C(x，y，z)作⽤有集中荷载{Q}=[Qx，Qy，Qz]T，转化成等效结点荷载列阵{PE}(e)。根据转换前后虚功相等原则有[3]：
{PE} (e)=[Nc]T{Q}(1)
式中[Nc]为形函数在C点的值。
1.2 体积⼒
设单元内作⽤有均布体积⼒{p}=[px，py，pz]T，则可将单元微分体积dV上的体⼒{p}dV视为集中荷载，利⽤式(1)，在单元的体积V内积分，即得到该分布体积⼒的等效结点荷载{PE }
1.3 表⾯⼒
设在单元的某边界⾯上作⽤表⾯⼒{ }={ ， ， }T，可将微分⾯积dA上的⾯荷载{p}dA视为集中荷载，利⽤式(1)在⾯积A上积分，即得到该分布⾯荷载的等效结点荷载{PE}
式中Ni表⽰单元的形函数，对于空间六⾯体20结点单元，其表达式为：
由于⼯程实际中给出的表⾯⼒往往是沿作⽤⾯的法向⼒和切向⼒，对于空间问题，外⼒往往垂直作⽤于边界表⾯，以？孜=1的⾯为例，载荷集度以 0表⽰，{PE} (e)表达式变为[4]
在计算结点荷载向量时由于被积函数⽐较复杂，常采⽤三维⾼斯积分公式[5]求解：
2 FORTRAN的程序实现
c语言程序设计教程李丽娟作⽤在单元上的⾯⼒在等效成结点荷载的过程中，由于坐标轴的确定，⼒所作⽤⾯的⽅向会发⽣相应的变化。本⽂对三维⼆⼗结点的六个空间⾯进⾏编号，新定义变量NOACE(*)来表⽰⼒作⽤的單元⾯；将⼒作⽤⾯固定，如EXISP=-1.0(ξ=-1.0)，再利⽤FORTRAN程序内部存储各⾯形函数顺序数组，如NFACE=(/1，9，5，20，8，12，4，16/)，确定空间六⾯体结点拓扑关系：
IF(NOACE.EQ.1) THEN ！⾯⼒作⽤单元⾯的选择
NFACE=(/1，9，5，20，8，12，4，16/)
EXISP= -1.0
最后调⽤其他公共⼦模块矩阵[6]，求解分布⾯⼒的等效结点荷载。利⽤断点法对程序进⾏调试过程中，单元在受均布⾯⼒作⽤下进⾏等效之后，四周⾓点的等效结点荷载变为-1/12q，边中点变为1/3q，所述规律与⽂献[7]中相符。⾯⼒等效程序如下：
LGASH=(NFACE(IODEG)-1)*NDOFN+1 ！确定结点⾃由度编码
MGASH=(NFACE(IODEG)-1)*NDOFN+2
NGASH=(NFACE(IODEG)-1)*NDOFN+3 ！荷载等效
ELOAD(NEASS，LGASH)=ELOAD(NEASS，LGASH)+SHAPEI(NFACE(IODEG))*PXCOM* WEIGP(IGAUS)*WEIGP(JGAUS)
ELOAD(NEASS，MGASH)=ELOAD(NEASS，MGASH)+SHAPEI(NFACE(IODEG))*PYCOM*WEIGP(IGAUS)*WEIGP(JGAUS)
ELOAD(NEASS，NGASH)=ELOAD(NEASS，NGASH)+SHAPEI(NFACE(IODEG))*PZCOM* WEIGP(IGAUS)*WEIGP(JGAUS)
3 實例及分析
某空间正⽅体压⼒盒各边长为2.0m，上端⾯受均布⾯⼒荷载q=10.0kN/m2，其余⾯均固定，E=2.0×105MPa，μ=0.2，体密度
ρ=50.0kN/m3，采⽤20节点等参单元，单元个数8，节点总数81。(图1)
为验证程序的正确性，在划分相同⽹格的条件下，取有限元软件FEA计算结果与程序计算结果做⽐对，如图2、图3所⽰。
由图2、图3的对⽐分析可知，⽤FORTRAN编制的程序计算结果与利⽤FEA模拟出来的结果基本吻合，说明了编制程序的正确性。
4 结语
有限元商业软件的核⼼仍旧是有限元算法，只有掌握相应的逻辑流程及算法原理才能熟练操作⼀款软件，解决遇到的棘⼿问题。⽂章基于虚功原理推导了弹性问题空间20结点单元等效结点荷载计算公式，利⽤FORTRAN语⾔，设计成通⽤模块程序。采⽤算例验证算法与程序的正确性，为类似连续介质体问题的有限元法求解提供有益的参考。
参考⽂献：
[1]徐芝纶.弹性⼒学简明教程[M].三版.北京：⾼等教育出版社，1983.
[2]王勖成.有限单元法[M].北京：清华⼤学出版社，2003.
[3]赵更新，等.⼟⽊⼯程结构分析程序设计[M].北京：中国⽔利⽔电出版社，2001.
[4]王元汉，李丽娟，李银平，等.有限元基础与程序设计[M].⼴州：华南理⼯⼤学出版社，2002.
[5]朱加铭，欧贵宝，何蕴增.有限元与边界元法[M].哈尔滨：哈尔滨⼯程⼤学出版社，2002.
[6]Hinton E and D R J. Finite Element Programming[M].London： Academic Press，1977.
[7]朱伯芳.有限单元法原理与应⽤[M].⼆版.北京：中国⽔利⽔电出版社，1998.