Blog Links

• DalNur | 博客总目录

• Abaqus 二次开发 简介

• Abaqus 二次开发 应用实例

• Python 语言创建 Abaqus inp 文件

一、前言

  Abaqus 的二次开发主要有两种：求解器层次的 Fortran 和 前后处理层次的 Python ，前者对应用户子程序开发，后者对应用户图形界面程序开发。Abaqus 用户子程序（User Subroutine）是一个极其强大和灵活的分析工具，通过它用户可以定制各种各样的 Abaqus 功能，如编写材料本构关系 (UMAT/VUMAT)、自定义单元（UEL）、并行计算和与外部程序进行数据交换等等。

Abaqus分析流程

  从某种程度上讲，Abaqus 是一款功能强大的数值法解偏微分方程的求解器。有限单元法和有限差分法都是求解偏微分方程的数值方法，它们的核心思想之一是分片插值，而构建插值函数是其中最重要的工作，创建用户单元的过程就是构建插值函数。理论上，任何偏微分方程（组）都可以利用 Abaqus 进行求解，只需建立起相应的用户子单元。计算结果的后处理不能直接利用 Abaqus/CAE 进行，需借助 Python 程序。

二、软件配置

  Abaqus 用户子程序可以采用 C、C++ 或 FORTRAN 代码编写。一般采用 FORTRAN 语言编写用户子程序。显然，在进行用户子程序开发前，除了安装 Abaqus 外，还需安装 Fortran 。安装 Abaqus 和安装 Fortran 是两个独立的过程，特别注意软件版本间的对应关系。

  Fortran（Formula Translation）语言是数值计算领域（尤其是高性能计算）所使用的主要语言。它是编译型语言，程序在执行前需要一个专门的编译过程，将源代码编译生成机器语言，再由机器运行机器码（二进制），C/C++ 等都是编译型语言。。Fortran 代码的执行最终是由计算机的 CPU (中央处理器) 来完成的。但 CPU 能且只能直接识别并执行机器指令，它的表现形式是二进制编码，形如：

二进制代码

  编译器（Compiler）就是把 Fortran、C 等高级语言翻译成机器码，从而使计算机能够执行并得出相应结果的软件。由 Intel 公司推出的 Fortran 编译器 ifort （Intel Fortran compiler）被广泛使用，过去它作为软件套包 Parallel Studio XE 的一部分来发布，其中还包含著名的数学函数库 MKL 。如今（2021）这套产品被更名为 oneAPI 并更改了发布形式，允许编译器和函数库作为独立组件下载和安装，还简化了授权机制，为初学者使用带来了不少方便。如果希望在可视化开发环境（IDE）下编写、调试 Fortran 程序，还需要安装 Visual Studio。一般，先安装 Visual Studio，后安装 ifort 。显然，Visual Studio 的版本不能高于 ifort ，否则容易出现不兼容问题。

软件名称	版本	下载链接
Abaqus	6.14-1	Abaqus 6.14-1
Visual Studio	2013	Visual Studio Ultimate 2013
Intel Visual Fortran	2013	Intel Visual Fortran 2013 SP1

  Abaqus 版本与 Fortran 编译器版本间的匹配关系可在达索公司官网中查到，操作方式如下：点击 SIMULIA Platforms & Configuration Support，根据 Abaqus 版本进入相应页面，本文以 Abaqus 6.14-1 为例，如下图所示。

  Abaqus 版本与 Fortran 编译器版本要严格对应，这两个软件的安装相互独立，互不影响。但在安装英特尔公司开发的 Fortran 编译器 IVF（Intel Visual Fortran）时，IVF 会自动和电脑中已存在的 Visual Studio 进行适配。因此，在安装 IVF 前，需安装 Visual Studio。那么，安装的 Visual Studio 版本的发布时间应早于 IVF 的发布时间，一般选择低于 IVF 版本的 Visual Studio，如 IVF 的版本是 2017，则 Visual Studio 的版本不应高于 2017 ，2015 及以下版本为宜。另外，在安装上述软件时，可以不安装到 C 盘，但必须保证安装路径为纯英文。上述三个软件全部安装完成后，便可将它们关联起来，具体操作见：

  - Abaqus 用户子程序配置. ldlib.

  - Abaqus6.14+VS2013+IVF2013安装教程. qintianhaohao.

环境变量的添加与.bat文件的修改

  执行 Abaqus Verification 后，关联结果记录在 verify.log 文件内，若子程序配置失败，则有如下提示，此时需打开 std_user.log 文件查看详细信息。一般，fatal error LINK ×××× 的报错来自 Visual Studio ，以 VS error LINK ×××× 为关键词，便可搜索到解决办法。

三、数学/力学理论

3.1 偏微分方程

  凡表示未知函数、未知函数的导数与自变量之间的关系的方程，叫做微分方程。能使微分方程横成立的函数叫做微分方程的解。当未知函数是多元函数时，若方程式中有未知数对自变量的偏微分，这样的微分方程被称作偏微分方程。在弹性力学中，两类平面问题的平衡微分方程为：

∂σx∂x+∂τyx∂y+fx=0∂τxy∂x+∂σy∂y+fy=0}(3-1)\left.\begin{array}{l} \frac{\partial \sigma_{x}}{\partial x}+\frac{\partial \tau_{yx}}{\partial y}+f_{x}=0 \\[8pt] \frac{\partial \tau_{xy}}{\partial x}+\frac{\partial \sigma_{y}}{\partial y}+f_{y}=0 \end{array}\right\} \tag{3-1} ∂x∂σx​​+∂y∂τyx​​+fx​=0∂x∂τxy​​+∂y∂σy​​+fy​=0​⎭⎬⎫​(3-1)

  人类目前理解的物理世界，即三维空间和一维时间，通常都可以用偏微分方程 (PDE) 来完整描述。然而，遗憾的是，数学中的解析方法往往不能解决我们现实世界中的复杂问题，这些复杂性通常源自几何结构的复杂和物理场的复杂 ，也就是说，难以求出实际问题的解析解。解析解是偏微分方程的精确解，精确解很难得到，我们退而求其次，追求精度高的近似解。满足一定精度的近似解被称为数值解。偏微分方程的常见数值解法有：有限差分法（FDM）、有限单元法（FEM）、边界元法（EEM）、离散元法（DEM） 。

3.2 有限单元法

  有限单元法是求解偏微分方程数值解（高精度近似解）的一种方法。有限单元法的诞生，将现实世界 “ 时空离散化、数字化 ” ，构造出这些偏微分方程组的近似数值方程，并用数值方法求解。从某种程度上讲，大型通用有限元分析软件 Abaqus 和 ANSYS 就是偏微分方程的求解器。什么是有限元?

二维多连通域的有限元离散

  有限单元法有严格的数学证明作为其近似的客观性和合理性的保证。有限单元法不是从问题的微分方程和相应的定解条件出发，而是从与其等效的积分形式出发。等效积分的一般形式是加权余量法，它适用于普遍的方程形式。利用加权余量法的原理，可建立多种近似方法，例如最小二乘法、伽辽金法等都属于这一类数值分析方法。如果原问题的微分方程具有某些特定的性质，则它的等效积分形式的伽辽金法可以归结为某个泛函的变分。相应的近似解法实际上是求解泛函的驻值问题。里兹法就是属于这一类近似解法。有限单元法区别于传统的加权余量法和求解泛函驻值的变分法，该法不是在整个求解域上假设近似函数，而是在各个单元上分片假设近似函数（形函数），这样就克服了在全域上假设近似函数所遇到的困难，是近代工程数值分析方法领域的重大突破。

  采用使余量的加权积分为 0 来求得微分方程近似解的方法称为加权余量法（weighted residual method, WRM）。加权余量法：基于等效积分形式的近似方法。任何独立的完全函数集都可以用来作为权函数，按照对权函数的不同选择就得到不同的加权余量的计算方法并赋予不同的名称，常用的权函数的选择有以下几种：配点法、子域法、最小二乘法、力矩法和伽辽金法（Galerkin method）。伽辽金法（Galerkin method）是求变分方程近似解的一种方法。简单的利用近似解的试探函数序列作为权函数。​ 采用伽辽金法得到的求解方程的系数矩阵是对称的，这是在用加权余量法建立有限元格式时毫无例外地采用伽辽金法的主要原因，而且当微分方程存在相应的泛函时，伽辽金法与变分法往往导致同样的结果。当方程阶数很高时，这种对称性将给计算带来很大的方便。

  等效积分弱形式 （由等效积分形式通过分部积分得到，通过提高权函数的连续性要求来降低待求场函数的连续性要求）被有限元法经常利用为理论基础的正是等效积分的伽辽金弱形式。（伽辽金方法是有限单元法中使用最普遍的方法）更多见： 等效积分的“弱”形式 。

3.3 弹性力学

  弹性力学是研究弹性体由于受外力、边界约束或温度改变等原因而发生的应力、形变和位移。(研究弹性体的力学响应)。理论力学研究刚体，材料力学研究杆件，结构力学研究杆系，弹性力学研究各种形状的弹性体(如梁、板、壳、体等)。材料力学引入平截面计算假定，弹性力学没有。弹性力学的五大假定：五大基本假定：连续、均匀、各项同性，线弹性、小变形。

  理论力学、材料力学和弹性力学谈平衡，对平衡的要求是不一样的。弹性力学要求微元体严格平衡。理论力学取整体为研究对象，考虑整体的平衡，只决定整体的运动状态，但整体的平衡不能保证局部的平衡。材料力学取微段为研究对象，考虑微段（有限体） 的平衡（近似），但局部的平衡不能保证处处平衡。弹性力学取微元体为研究对象，考虑微元体的平衡（精确）能保证处处都平衡。处处平衡，局部必平衡，每个局部都平衡，整体必平衡。故弹性力学谈平衡是处处严格精确的平衡，满足弹性力学的平衡，确保材料力学的平衡，确保了理论力学的平衡，反之不然。弹性力学要求处处严格满足平衡方程，材料力学可以不是处处，局部满足也可以。

  弹性力学三大基本规律是：变形连续规律、应力-应变关系和运动（或平衡）规律。弹性力学中许多定理、公式和结论等，都可以从三大基本规律推导出来。反映变形连续规律的数学方程有两类：几何方程和位移边界条件。若应力和应变呈线性关系，这个关系便叫作广义胡克定律。反映运动（或平衡）规律的数学方程有两类：运动（或平衡）微分方程和载荷边界条件。

  对于三维问题，弹性力学的基本方程共 15 个，其中，运动（或平衡）方程 3 个，几何方程 6 个，物理方程 6 个。这 15 个基本方程中包含 15 个位置函数，即：3 个位移分量，6 个应变分量，6 个应力分量，基本方程的数目等于未知函数的数目。弹性力学的任务，就是在适当的边界条件下，从基本方程中求解这些未知函数。然而遗憾的是，仅少数情况能求得解析解答。对于大多数的复杂情况，我们退而求其次，采用有限单法可求得高精度的近似解，即数值解。综上，求解弹性力学问题，最终都会归结到偏微分方程的求解上来，而有限元法是解偏微分方程的一种数值解法。有限单元法，是一种有效解决数学问题的解题方法。其理论基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式 ，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。复杂的弹性力学问题可借助有限单元法进行求解，首先将连续体转化为离散化结构，然后再利用分片插值技术与虚功原理或变分方法进行求解。在进行有限元分析时，由于涉及到多维广义坐标下的运算，各种方程多以矩阵/向量的形式表达，力求简化形式，突出重点。

  弹性力学中有限元分析的物理实质：以假设的单元位移场来近似代替真实位移，使单元在 ∫dW刚 =0\int \mathrm{d} W_{\text {刚 }}=0∫dW刚 ​=0 意义下“平衡” 。在假设的位移场中包含刚体位移部分时，则可保证单元上力系平衡。有限元分析的实质：以单元上全部外力平衡、各结点全部平衡、单元的 ∫dW刚 =0\int \mathrm{d} W_{\text {刚 }}=0∫dW刚 ​=0 的离散化结点位移解答来作为实际平衡问题的近似解。有限元的平衡并没有弹性力学定义中的那么严格 。

四、构造单元

  采用有限单元法求解弹性力学问题的主要步骤依次是：构造单元、结构离散、单元分析、整体分析、求解计算和结果处理。在Abaqus 进行用户单元开发时，主要涉及构造单元和结果的后处理两个部分，其余步骤均可由 Abaqus 自动完成。所谓的构造单元就是根据问题的类型和性质选定单元的形式并构造单元的插值函数。（创造母元）

Abaqus中几种常见的单元

4.1 形函数

  构造单元主要是确定单元维度与形状（二维、三维或2.5维度）、结点个数，每个结点自由度数和单元插值函数（形函数）。其中，最重要的工作就是确定单元的插值函数即形函数（根据各种理论/假设，经过一系列推导，得到形函数矩阵）。若无特殊说明，本文提及的有限元均为弹性力学中的有限位移元，单元以结点位移为基本未知量。

  有限单元法采用能量原理 (虚位移原理) 进行单元分析，因而必须提前给定位移函数的具体形式。一般而言，位移函数的选取影响计算结果的精度。在弹性力学中，恰当的选取位移函数是相当不容易的。有限单元法中，当单元划分的足够小时，把位移函数设定为简单的多项式形式也可以得到具有相当精度的结果。无论何种形式的位移函数，位移函数中必须包含能反应刚体位移和常应变状态的常数项和一次项，该条件构成了单元的完备性准则。此外，位移函数在单元内要连续，在相邻单元间应尽量保持协调 (即单元间位移函数也应尽量保持连续)，这是单元的协调性条件。理论和实践都已证明，完备性准则是有限元解收敛于真实解的必要条件。

  应用插值公式，可由单元结点位移 δe\boldsymbol {\delta_e}δe​，求得单元位移函数 d\boldsymbol {d}d，这个插值公式表示了单元中位移的分布形式，因此称为位移模式。有限位移元以结点位移为基本未知量，由单元的结点位移列阵 δe\boldsymbol {\delta_e}δe​ 和单元的形函数矩阵 N\boldsymbol {N}N 可得到单元位移场，即一旦结点位移确定后，单元内任意一点的位移可由下式求得：

d=N⋅δe(4-1)\boldsymbol {d}=\boldsymbol {N} \cdot \boldsymbol {\delta_e} \tag{4-1} d=N⋅δe​(4-1)

式中， N\boldsymbol {N}N 为形函数矩阵，对于非梁壳单元 (二维平面单元、三维体单元)，它是一个 单个结点自由度数 × 单元结点总自由度数 阶矩阵。

  一般情况下，单元内任一点的实际位移是坐标的很复杂的函数，仅利用单元的结点位移是不能精确表示的，由结点位移表示的单元内位移只是实际位移的一部分，或者说是实际位移的一个近似。实践表明，当合理地选择由结点位移可以确定的、用以代替单元的实际位移的位移形式，随着单元尺寸的减小，结果会收敛于实际位移。我们把在有限元分析中用来代替实际位移的位移形式称为位移函数。位移函数的选择直接关系到结果的收敛性，是很关键的一步。一般从泰勒级数展开的意义出发，选多项式作为位移函数。这不仅运算简单，并且可由项数的多少直接控制结果的精度。

  形函数是定义于单元内部的、坐标的连续函数。作为单元的插值函数，可将单元内任意一点的位移用结点位移表示。主要取决于单元形状、结点类型和结点数目。形函数主要有两种类型：Lagrange 型和 Hermite 型。其中，Lagrange 型不需要形函数在结点上斜率或曲率。Hermite 型需要形函数在结点上斜率或曲率。单元插值形函数一般采用不同阶次的幂函数多项式形式， 所采用的多项式的幂次即为形函数的幂次。形函数的物理意义是，结点 i 发生单位位移时，所引起单元内位移分布。 建立形函数主要有两种方法：广义坐标法和直接法（试凑法）。

  对于平面矩形双线性单元，其单元位移函数 d=(u,v)T\boldsymbol {d}=(u, v)^Td=(u,v)T，结点位移列阵 δe=(u1,v1,u2,v2,u3,v3,u4,v4)T\boldsymbol {\delta_e} = (u_1,v_1,u_2,v_2,u_3,v_3, u_4,v_4)^Tδe​=(u1​,v1​,u2​,v2​,u3​,v3​,u4​,v4​)T，根据有限元理论有：

d=(uv)=(N10N20N30N400N10N20N30N4)(u1v1u2v2u3v3u4v4)(4-2)\boldsymbol {d} = \begin{pmatrix}u\\[8pt]v\\\end{pmatrix} = \left (\begin{array}{cc|cc|cc|cc} N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4 & 0 \\[8pt] 0 & N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4 \end{array}\right) \begin{pmatrix}u_1\\v_1\\ \hline u_2\\v_2\\ \hline u_3\\v_3\\ \hline u_4\\v_4\\ \end{pmatrix} \tag{4-2} d=(uv​)=(N1​0​0N1​​N2​0​0N2​​N3​0​0N3​​N4​0​0N4​​)⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛​u1​v1​u2​v2​u3​v3​u4​v4​​​⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞​(4-2)

  将上式展开，有：

{u=N1u1+N2u2+N3u3+N4u4v=N1v1+N2v2+N3v3+N4v4(4-3)\left \{\begin{aligned} u & = N_1u_1+N_2u_2+N_3u_3+N_4u_4 \\[8pt] v & = N_1v_1+N_2v_2+N_3v_3+N_4v_4 \end{aligned} \right. \tag{4-3} ⎩⎨⎧​uv​=N1​u1​+N2​u2​+N3​u3​+N4​u4​=N1​v1​+N2​v2​+N3​v3​+N4​v4​​(4-3)

其中，u,v,wu,\,v,\,wu,v,w 和 N1,N2,N3,N4N_1,\,N_2,\,N_3,\,N_4\,N1​,N2​,N3​,N4​ 均为位置坐标 x,y,zx,\,y,\,zx,y,z 的函数。 在正则坐标系下，对于长为 2b 宽为 2a 的矩形单元有：

N=(N10N20N30N400N10N20N30N4)(4-4)\boldsymbol {N} = \left (\begin{array}{cc|cc|cc|cc} N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4 & 0 \\[8pt] 0 & N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4 \end{array}\right) \tag{4-4} N=(N1​0​0N1​​N2​0​0N2​​N3​0​0N3​​N4​0​0N4​​)(4-4)

{N1=+14(ξ−1)(η−1)N2=−14(ξ+1)(η−1)N3=−14(ξ+1)(η+1)N4=−14(ξ−1)(η+1)(4-5)\left \{\begin{aligned} N_{1}&=+\frac{1}{4}(\xi-1)(\eta-1) \\[8pt] N_{2}&=-\frac{1}{4}(\xi+1)(\eta-1) \\[8pt] N_{3}&=-\frac{1}{4}(\xi+1)(\eta+1) \\[8pt] N_{4}&=-\frac{1}{4}(\xi-1)(\eta+1) \end{aligned} \right. \tag{4-5} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧​N1​N2​N3​N4​​=+41​(ξ−1)(η−1)=−41​(ξ+1)(η−1)=−41​(ξ+1)(η+1)=−41​(ξ−1)(η+1)​(4-5)

ξ=xa,η=yb(4-6)\xi = \frac{x}{a} , \quad \eta = \frac{y}{b} \tag{4-6} ξ=ax​,η=by​(4-6)

矩形双线性单元（左：直角坐标系下的单元示意，右：正则坐标系下的单元示意）

4.2 应变矩阵

  显然，式（4-6）仅对矩形单元成立，对任意四边形单元式（4-4）、式（4-5）仍然成立，式（4-6）不再成立。将式（4-1）代入弹性力学几何方程，有：

ε=AT⋅d=AT⋅N⋅δe=B⋅δe(4-7)\boldsymbol\varepsilon = \boldsymbol A^T \cdot \boldsymbol d = \boldsymbol A^T \cdot \boldsymbol {N} \cdot \boldsymbol {\delta_e} = \boldsymbol B \cdot \boldsymbol {\delta_e} \tag{4-7} ε=AT⋅d=AT⋅N⋅δe​=B⋅δe​(4-7)

式中，A\boldsymbol AA 为微分算子矩阵，AT\boldsymbol A^TAT 为其转置矩阵，如式（4-8）所示；B\boldsymbol BB 为应变矩阵，单元结点位移列阵左乘应变矩阵变可得单元应变场，由式（4-7）可计算单元内任意一点处的应变。

A=[∂∂x0∂∂y0∂∂y∂∂x](4-8)\pmb{A}=\left[\begin{array}{ccc} \frac{\partial}{\partial x} & 0 & \frac{\partial}{\partial y} \\[8pt] 0 & \frac{\partial}{\partial y} & \frac{\partial}{\partial x} \end{array}\right] \tag{4-8} AAA=⎣⎡​∂x∂​0​0∂y∂​​∂y∂​∂x∂​​⎦⎤​(4-8)

  显然，应变矩阵 B\boldsymbol BB 是由形函数矩阵左乘微分算子矩阵的转置矩阵得到的（定义的）。对于矩形双线性单元，其应变矩阵如下所示。（可以看到微分算子矩阵定义了求导规则）

B=AT⋅N=[∂∂x00∂∂y∂∂y∂∂x](N10N20N30N400N10N20N30N4)(4-9)\pmb{B}= \pmb{A}^T \cdot \pmb{N} = \left[\begin{array}{cc} \frac{\partial\,}{\partial x} & 0 \\[8pt] 0 & \frac{\partial\,}{\partial y} \\[8pt] \frac{\partial\,}{\partial y} & \frac{\partial\,}{\partial x} \end{array}\right] \left (\begin{array}{cc|cc|cc|cc} N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4 & 0 \\[8pt] 0 & N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4 \end{array}\right) \tag{4-9} BBB=AAAT⋅NNN=⎣⎢⎢⎢⎡​∂x∂​0∂y∂​​0∂y∂​∂x∂​​⎦⎥⎥⎥⎤​(N1​0​0N1​​N2​0​0N2​​N3​0​0N3​​N4​0​0N4​​)(4-9)

  对于平面矩形双线性单元，其应变矩阵 BBB 为 3×8 阶矩阵，将上式展开，有：

B=AT⋅N=[∂N1∂x0∂N2∂x0∂N3∂x0∂N4∂x00∂N1∂y0∂N2∂y0∂N3∂y0∂N4∂y∂N1∂y∂N1∂x∂N2∂y∂N2∂x∂N3∂y∂N3∂x∂N4∂y∂N4∂x](4-10)\pmb{B}= \pmb{A}^T \cdot \pmb{N} = \left[\begin{array}{cccccccc} \frac{\partial{N_1}}{\partial x} & 0 & \frac{\partial{N_2}}{\partial x} & 0 & \frac{\partial{N_3}}{\partial x} & 0 & \frac{\partial{N_4}}{\partial x} & 0 \\[8pt] 0 & \frac{\partial{N_1}}{\partial y} & 0 & \frac{\partial{N_2}}{\partial y} & 0 & \frac{\partial{N_3}}{\partial y} & 0 & \frac{\partial{N_4}}{\partial y} \\[8pt] \frac{\partial{N_1}}{\partial y} & \frac{\partial{N_1}}{\partial x} & \frac{\partial{N_2}}{\partial y} & \frac{\partial{N_2}}{\partial x} & \frac{\partial{N_3}}{\partial y} & \frac{\partial{N_3}}{\partial x} & \frac{\partial{N_4}}{\partial y} & \frac{\partial{N_4}}{\partial x} \end{array}\right] \tag{4-10} BBB=AAAT⋅NNN=⎣⎢⎢⎢⎡​∂x∂N1​​0∂y∂N1​​​0∂y∂N1​​∂x∂N1​​​∂x∂N2​​0∂y∂N2​​​0∂y∂N2​​∂x∂N2​​​∂x∂N3​​0∂y∂N3​​​0∂y∂N3​​∂x∂N3​​​∂x∂N4​​0∂y∂N4​​​0∂y∂N4​​∂x∂N4​​​⎦⎥⎥⎥⎤​(4-10)

式中，NiN_iNi​ 按式（4-5）计算。

五、等参数单元

  对于矩形单元，它很难适应曲线边界和非正交的直线边界。同时，在划分单元时，改变单元的大小也很困难，很难在不同部分采用大小不同的单元。因此，矩形平面单元未能在实际中得到广泛应用。为此，我们引入等参元来解决这一问题。通过子单元与母单元间的映射，实现任意四边形单元的分析计算。为了将局部（自然）坐标中几何形状规则的单元转换成总体坐标中几何形状扭曲的单元，以满足对一般形状求解域进行离散化的需要，必须建立一个坐标变换（局部坐标系向整体坐标系变换），使得局部坐标系下的单元点与整体坐标系下的单元点有一一对应关系，即：

[xy]=T[ξη](5-1)\left[\begin{array}{l} x \\ y \end{array}\right]= \pmb{T} \left[\begin{array}{l} \xi \\ \eta \end{array}\right] \tag{5-1} [xy​]=TTT[ξη​](5-1)

式中，TTT 为坐标变换矩阵，矩阵左乘表示坐标变换。

  结构离散化时，建立了两种坐标系——结构整体坐标系和单元局部坐标系。两种坐标系下的同一物理量存在着相互转换关系，将局部坐标系下表达的量转换为整体坐标系下表达的量，或相反均称为坐标转换。

5.1 坐标变换

  在建立单元位移函数时，我们通过插值的方式用单元的结点位移来表示单元内任意一点处的位移。显然，也可以类似的方式，用单元上某些特征结点的坐标来表示单元内任意一点处的坐标值。这样就建立起单元形状在两个坐标系下的对应变换，为了建立相应的坐标变换，最方便的方法是：

{x=∑i=1mNi′xiy=∑i=1mNi′yi(5-2)\left \{\begin{aligned} x=\sum_{i=1}^{m} N_{i}^{\prime} x_{i} \\[14pt] y=\sum_{i=1}^{m} N_{i}^{\prime} y_{i} \end{aligned} \right. \tag{5-2} ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧​x=i=1∑m​Ni′​xi​y=i=1∑m​Ni′​yi​​(5-2)

式中，mmm 为用于描述单元形状/位置坐标而选用的特征结点个数，Ni′N_{i}^{\prime}Ni′​ 为单元形状插值函数，xix_{i}xi​、yiy_{i}yi​ 为特征结点坐标值。 其中，mmm 用以进行坐标变换的单元结点数；xi,yix_i, \ y_ixi​, yi​ 是这些结点在总体坐标系下的坐标值；Ni′N_i^′Ni′​ 称为形状函数，实际上它也是用局部（自然）坐标表示的插值函数。如果坐标变换和位移插值采用相同的结点，并且采用相同的插值函数，即 m=nm=nm=n，Ni′=NiN_i^′=N_iNi′​=Ni​，则称这种变换为等参变换。等参数变换？因为用于描述单元形状的结点、插值函数与描述单元位移场所采用的结点、插值函数完全相同，故名等参变换，和位移场参数相同。此外，还有超参元和亚参元。对于等参元，式（5-2）转化为：

{x=∑i=1nNi(ξ,η)xi=N1x1+N2x2+N3x3+N4x4y=∑i=1nNi(ξ,η)yi=N1y1+N2y2+N3y3+N4y4(5-3)\left \{\begin{aligned} x &=\sum_{i=1}^{n} N_{i}(\xi, \eta) x_{i} = N_{1}x_{1} + N_{2}x_{2} + N_{3}x_{3} + N_{4}x_{4} \\[14pt] y &=\sum_{i=1}^{n} N_{i}(\xi, \eta) y_{i} = N_{1}y_{1} + N_{2}y_{2} + N_{3}y_{3} + N_{4}y_{4} \end{aligned} \right. \tag{5-3} ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧​xy​=i=1∑n​Ni​(ξ,η)xi​=N1​x1​+N2​x2​+N3​x3​+N4​x4​=i=1∑n​Ni​(ξ,η)yi​=N1​y1​+N2​y2​+N3​y3​+N4​y4​​(5-3)

式中，NiN_iNi​ 为单元形函数，均为局部坐标 ξ,η\xi, \etaξ,η 的函数。

  通过上式可建立起两个坐标系之间的变换，从而将自然坐标系下形状规则的单元变换为总体笛卡尔坐标系下形状扭曲的单元。称前者为母单元，后者为子单元。只要给出母元中任意一点的坐标利用坐标变换关系式就可以找到子元中唯一一点与之对应。当母元确定后，即结点数与形函数确定后，坐标变换式只与结点坐标有关，而子元坐标是根据单元划分任意指定的，这导致一个母元与一族子元相对应，这也是分别称其为母元和子元的原因。为什么不直接在整体坐标系下构造扭曲单元的插值函数，而是在局部坐标系下构造插值函数后通过坐标变换得任意一点处的位移？即为什么要作等参变换？若直接构造扭曲单元的插值函数，由不同单元计算的交界上一点的位移将不一定相等，则位移函数在两单元公共边界上将不再连续，而连续性是弹性力学的基本假定。采用等参变化可以解决这一问题，能保证位移在单元间的协调连续。

线性单元的坐标变化（左：凸四边形子元，右：正方形母元）

  在全局直角坐标系/整体坐标系 xoyxoyxoy 中，任取一任意四边形单元，四边形的四个角点取为结点，各结点的坐标值为 (xi,yi)(x_i, \, y_i)(xi​,yi​)，对于这种任意四边形等参元，可令其实际形状所构成的单元为子单元，把子单元各边中点连线做一个局部坐标系/自然坐标系 ξo′η\xi o^{\prime} \etaξo′η，且令单元各结点的局部坐标分为为 1 或 -1，这样就把子单元映射到局部坐标系上，而称为标准的正方形单元，称此正方形单元为母单元。整体坐标系适用于所有单元，即适用于整个求解区域，而局部坐标系只适用于每一个特定单元。

  因此，我们就可以把矩形单元的位移插值函数式（4-3），单元内任意一点的坐标变换式（5-3），以及局部坐标变换式，应用在任意四边形等参元上，并重新写成 (n=4)(n=4)(n=4)：

{u=∑i=1nNi(ξ,η)ui=N1u1+N2u2+N3u3+N4u4v=∑i=1nNi(ξ,η)vi=N1v1+N2v2+N3v3+N4v4(5-4)\left \{\begin{aligned} u & =\sum_{i=1}^{n} N_{i}(\xi, \eta) u_{i} = N_1u_1+N_2u_2+N_3u_3+N_4u_4 \\[10pt] v & = \sum_{i=1}^{n} N_{i}(\xi, \eta) v_{i} = N_1v_1+N_2v_2+N_3v_3+N_4v_4 \end{aligned} \right. \tag{5-4} ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧​uv​=i=1∑n​Ni​(ξ,η)ui​=N1​u1​+N2​u2​+N3​u3​+N4​u4​=i=1∑n​Ni​(ξ,η)vi​=N1​v1​+N2​v2​+N3​v3​+N4​v4​​(5-4)

{ξ=∑i=1nNi(ξ,η)ξi=N1ξ1+N2ξ2+N3ξ3+N4ξ4η=∑i=1nNi(ξ,η)ηi=N1η1+N2η2+N3η3+N4η41=∑i=1nNi(ξ,η)(5-5)\left \{\begin{aligned} \xi & =\sum_{i=1}^{n} N_{i}(\xi, \eta) \xi_{i} = N_1\xi_1+N_2\xi_2+N_3\xi_3+N_4\xi_4 \\[10pt] \eta & = \sum_{i=1}^{n} N_{i}(\xi, \eta) \eta_{i} = N_1\eta_1+N_2\eta_2+N_3\eta_3+N_4\eta_4 \\[10pt] 1 & = \sum_{i=1}^{n} N_{i}(\xi, \eta) \end{aligned} \right. \tag{5-5} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧​ξη1​=i=1∑n​Ni​(ξ,η)ξi​=N1​ξ1​+N2​ξ2​+N3​ξ3​+N4​ξ4​=i=1∑n​Ni​(ξ,η)ηi​=N1​η1​+N2​η2​+N3​η3​+N4​η4​=i=1∑n​Ni​(ξ,η)​(5-5)

{x=∑i=1nNi(ξ,η)xi=N1x1+N2x2+N3x3+N4x4y=∑i=1nNi(ξ,η)yi=N1y1+N2y2+N3y3+N4y4(5-6)\left \{\begin{aligned} x &=\sum_{i=1}^{n} N_{i}(\xi, \eta) x_{i} = N_{1}x_{1} + N_{2}x_{2} + N_{3}x_{3} + N_{4}x_{4} \\[14pt] y &=\sum_{i=1}^{n} N_{i}(\xi, \eta) y_{i} = N_{1}y_{1} + N_{2}y_{2} + N_{3}y_{3} + N_{4}y_{4} \end{aligned} \right. \tag{5-6} ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧​xy​=i=1∑n​Ni​(ξ,η)xi​=N1​x1​+N2​x2​+N3​x3​+N4​x4​=i=1∑n​Ni​(ξ,η)yi​=N1​y1​+N2​y2​+N3​y3​+N4​y4​​(5-6)

  式（5-4）和式（5-5）是任意四边形在局部坐标系下的位移插值函数和单元内任意一点局部坐标插值公式，而式（5-6）是每个单元的局部坐标系与整体坐标系之间的坐标变换式。以上等参变换能保证在整体坐标系下满足相容条件，也就是说，在两相邻任意四边形单元公共边上的位移是连续的，坐标变换后仍然是连续的，两相邻单元公共边上的公共点在坐标变换后任为公共点，不会出现重叠和撕裂现象。

  在对等参元进行分析时，单元内仍然假定为连续的、均匀的、各项同性的完全弹性体。因为单元内部仍是连续体，弹性力学的变形连续规律（几何方程）仍然成立。在整体坐标系下，弹性力学的几何方程为：

ε=B⋅δe(5-7)\boldsymbol\varepsilon = \boldsymbol B \cdot \boldsymbol {\delta_e} \tag{5-7} ε=B⋅δe​(5-7)

B=AT⋅N=[∂N1∂x0∂N2∂x0∂N3∂x0∂N4∂x00∂N1∂y0∂N2∂y0∂N3∂y0∂N4∂y∂N1∂y∂N1∂x∂N2∂y∂N2∂x∂N3∂y∂N3∂x∂N4∂y∂N4∂x](5-8)\pmb{B}= \pmb{A}^T \cdot \pmb{N} = \left[\begin{array}{cccccccc} \frac{\partial{N_1}}{\partial x} & 0 & \frac{\partial{N_2}}{\partial x} & 0 & \frac{\partial{N_3}}{\partial x} & 0 & \frac{\partial{N_4}}{\partial x} & 0 \\[8pt] 0 & \frac{\partial{N_1}}{\partial y} & 0 & \frac{\partial{N_2}}{\partial y} & 0 & \frac{\partial{N_3}}{\partial y} & 0 & \frac{\partial{N_4}}{\partial y} \\[8pt] \frac{\partial{N_1}}{\partial y} & \frac{\partial{N_1}}{\partial x} & \frac{\partial{N_2}}{\partial y} & \frac{\partial{N_2}}{\partial x} & \frac{\partial{N_3}}{\partial y} & \frac{\partial{N_3}}{\partial x} & \frac{\partial{N_4}}{\partial y} & \frac{\partial{N_4}}{\partial x} \end{array}\right] \tag{5-8} BBB=AAAT⋅NNN=⎣⎢⎢⎢⎡​∂x∂N1​​0∂y∂N1​​​0∂y∂N1​​∂x∂N1​​​∂x∂N2​​0∂y∂N2​​​0∂y∂N2​​∂x∂N2​​​∂x∂N3​​0∂y∂N3​​​0∂y∂N3​​∂x∂N3​​​∂x∂N4​​0∂y∂N4​​​0∂y∂N4​​∂x∂N4​​​⎦⎥⎥⎥⎤​(5-8)

式（5-7）是在整体坐标系下推导得到的。其中，应变矩阵 BBB 中的每个元素都是单元形函数 NiN_iNi​ 对整体坐标 x,yx, \, yx,y 的偏导数。而任意四边形等参元的形函数 NiN_iNi​ 又是针对局部坐标的，即为局部坐标 ξ,η\xi, \etaξ,η 的函数。因此，需要进行坐标变换，统一变量。设任意四边形在整体坐标系下的位移插值函数式：‘

{u=u(x,y)v=v(x,y)(5-6)\left \{\begin{aligned} u= u(x, y) \\[14pt] v = v(x, y) \end{aligned} \right. \tag{5-6} ⎩⎪⎨⎪⎧​u=u(x,y)v=v(x,y)​(5-6)

  而该单元在局部坐标系下的位移插值函数式（5-4）可写成：

{u=∑i=1nNi(ξ,η)ui=u(ξ,η)v=∑i=1nNi(ξ,η)vi=v(ξ,η)(5-7)\left \{\begin{aligned} u & =\sum_{i=1}^{n} N_{i}(\xi, \eta) u_{i} = u(\xi, \eta) \\[10pt] v & = \sum_{i=1}^{n} N_{i}(\xi, \eta) v_{i} = v(\xi, \eta) \end{aligned} \right. \tag{5-7} ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧​uv​=i=1∑n​Ni​(ξ,η)ui​=u(ξ,η)=i=1∑n​Ni​(ξ,η)vi​=v(ξ,η)​(5-7)

  这两种形式的位移插值函数式通过坐标变换式（5-6）联系起来，为了方便把式（5-6）写成：

{x=∑i=1nNi(ξ,η)xi=x(ξ,η)y=∑i=1nNi(ξ,η)yi=y(ξ,η)(5-8)\left \{\begin{aligned} x &=\sum_{i=1}^{n} N_{i}(\xi, \eta) x_{i} = x(\xi, \eta) \\[14pt] y &=\sum_{i=1}^{n} N_{i}(\xi, \eta) y_{i} = y(\xi, \eta) \end{aligned} \right. \tag{5-8} ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧​xy​=i=1∑n​Ni​(ξ,η)xi​=x(ξ,η)=i=1∑n​Ni​(ξ,η)yi​=y(ξ,η)​(5-8)

5.2 雅可比矩阵

  根据复合函数求导法则有：

∂u(ξ,η)∂ξ=∂u(x,y)∂x⋅∂x∂ξ+∂u(x,y)∂y⋅∂y∂ξ=∂u∂x⋅∂x∂ξ+∂u∂y⋅∂y∂ξ(5-9)\begin{aligned} \frac{\partial u(\xi, \eta)}{\partial \xi} &=\frac{\partial u(x, y)}{\partial x} \cdot \frac{\partial x}{\partial \xi}+\frac{\partial u(x, y)}{\partial y} \cdot \frac{\partial y}{\partial \xi} \\[14pt] &=\frac{\partial u}{\partial x} \cdot \frac{\partial x}{\partial \xi}+\frac{\partial u}{\partial y} \cdot \frac{\partial y}{\partial \xi} \end{aligned} \tag{5-9} ∂ξ∂u(ξ,η)​​=∂x∂u(x,y)​⋅∂ξ∂x​+∂y∂u(x,y)​⋅∂ξ∂y​=∂x∂u​⋅∂ξ∂x​+∂y∂u​⋅∂ξ∂y​​(5-9)

同理，对于另一个坐标 η\etaη 也可写出类似的表达式，将它们集合成矩阵形式，则有：

{∂u∂ξ∂u∂η}=[∂x∂ξ∂y∂ξ∂x∂η∂y∂η]{∂u∂x∂u∂y}⇒{∂u∂ξ∂u∂η}=J{∂u∂x∂u∂y}(5-10)\left\{\begin{array}{l} \frac{\partial u}{\partial \xi} \\[8pt] \frac{\partial u}{\partial \eta} \end{array}\right\}=\left[\begin{array}{ll} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\[8pt] \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{array}\right]\left\{\begin{array}{l} \frac{\partial u}{\partial x} \\[8pt] \frac{\partial u}{\partial y} \end{array}\right\} \Rightarrow\left\{\begin{array}{l} \frac{\partial u}{\partial \xi} \\[8pt] \frac{\partial u}{\partial \eta} \end{array}\right\}= \pmb{J}\left\{\begin{array}{l} \frac{\partial u}{\partial x} \\[8pt] \frac{\partial u}{\partial y} \end{array}\right\} \tag{5-10} ⎩⎨⎧​∂ξ∂u​∂η∂u​​⎭⎬⎫​=⎣⎡​∂ξ∂x​∂η∂x​​∂ξ∂y​∂η∂y​​⎦⎤​⎩⎨⎧​∂x∂u​∂y∂u​​⎭⎬⎫​⇒⎩⎨⎧​∂ξ∂u​∂η∂u​​⎭⎬⎫​=JJJ⎩⎨⎧​∂x∂u​∂y∂u​​⎭⎬⎫​(5-10)

式中，J\pmb{J}JJJ 为雅可比矩阵或坐标变换矩阵，它是局部坐标的函数。对于任意四边形等参元，其具体形式为：

J(ξ,η)≡J=∂(x,y)∂(ξ,η)=[∂x∂ξ∂y∂ξ∂x∂η∂y∂η]=[∑∂Ni∂ξxi∑∂Ni∂ξyi∑∂Ni∂ηxi∑∂Ni∂ηyi](5-10)\pmb{J}(\xi, \eta) \equiv \pmb{J} = \frac{\partial(x, y)}{\partial(\xi, \eta)}=\left[\begin{array}{ll} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\[8pt] \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{array}\right]=\left[\begin{array}{cc} \sum \frac{\partial N_{i}}{\partial \xi} x_{i} & \sum \frac{\partial N_{i}}{\partial \xi} y_{i} \\[8pt] \sum \frac{\partial N_{i}}{\partial \eta} x_{i} & \sum \frac{\partial N_{i}}{\partial \eta} y_{i} \end{array}\right] \tag{5-10} JJJ(ξ,η)≡JJJ=∂(ξ,η)∂(x,y)​=⎣⎡​∂ξ∂x​∂η∂x​​∂ξ∂y​∂η∂y​​⎦⎤​=⎣⎡​∑∂ξ∂Ni​​xi​∑∂η∂Ni​​xi​​∑∂ξ∂Ni​​yi​∑∂η∂Ni​​yi​​⎦⎤​(5-10)

  若雅可比矩阵 J\pmb{J}JJJ 为可逆矩阵，式（5-9）等号两边同时左乘雅可比矩阵的可逆矩阵 J−1\pmb{J}^{-1}JJJ−1，有

{∂u∂x∂u∂y}=J−1{∂u∂ξ∂u∂η}(5-11)\left\{\begin{array}{l} \frac{\partial u}{\partial x} \\[8pt] \frac{\partial u}{\partial y} \end{array}\right\} = \pmb{J}^{-1}\left\{\begin{array}{l} \frac{\partial u}{\partial \xi} \\[8pt] \frac{\partial u}{\partial \eta} \end{array}\right\} \tag{5-11} ⎩⎨⎧​∂x∂u​∂y∂u​​⎭⎬⎫​=JJJ−1⎩⎨⎧​∂ξ∂u​∂η∂u​​⎭⎬⎫​(5-11)

  若雅可比矩阵 J\pmb{J}JJJ 为可逆矩阵，根据线性代数知识，J−1\pmb{J}^{-1}JJJ−1 可由下式计算：

J−1=J∗∣J∣=J∗detJ(5-12)\pmb{J}^{-1} = \frac{\pmb{J}^*}{|\pmb{J}|} = \frac{\pmb{J}^*}{{\rm det} \pmb{J}} \tag{5-12} JJJ−1=∣JJJ∣JJJ∗​=detJJJJJJ∗​(5-12)

式中，J∗\pmb{J}^*JJJ∗ 为雅可比矩阵 J\pmb{J}JJJ 的伴随矩阵；∣J∣|\pmb{J}|∣JJJ∣ 或 detJ{\rm det} \pmb{J}detJJJ 为雅可比矩阵 J\pmb{J}JJJ 的行列式。

  由微积分学知识，可知两个坐标之间一对一变换的条件是雅可比行列式 ∣J∣≠0|\pmb{J} |≠0∣JJJ∣​=0 ，等参变换作为一种坐标变换也必须服从此条件。若 ∣J∣=0|\pmb{J} |=0∣JJJ∣=0，则表明笛卡尔坐标中体积/面积微元为零，即在自然坐标中的体积微元对应笛卡尔坐标中的一点，这种变换显然不是一一对应。另外，若 ∣J∣=0|\pmb{J} |=0∣JJJ∣=0 ，则雅可比矩阵为不可逆矩阵，即 J−1\pmb{J}^{-1}JJJ−1 不存在，所以两个坐标之间的偏导数变换将不可能实现。所以，为了保证等参数变换得以进行，所谓的任意四边形所有内角都必须小于 180°，当内角接近 180° 时，∣J∣|\pmb{J}|∣JJJ∣ 接近于 0，因为 ∣J∣|\pmb{J}|∣JJJ∣ 出现在矩阵求逆过程中的分母上，这将导致计算误差迅速增大，因为在网格划分上不应使四边形单元过于歪斜。

单元划分的正常与不正常情况

  根据行列式的计算规则及式（4-5），对于任意 4 结点四边形等参元，其雅可比行列式 detJ{\rm det} \pmb{J}detJJJ 最终可写成如下形式：

detJ=∣J∣=aξ+bη+c(5-13){\rm det} \pmb{J}= |\pmb{J}| = a\xi+b\eta+c \tag{5-13} detJJJ=∣JJJ∣=aξ+bη+c(5-13)

式中，a,b,ca, \, b, \, ca,b,c 为常数，由以上相关公式求导、代入并展开即可得到。显然，任意 4 结点四边形等参元的雅可比行列式是 ξ,η\xi, \etaξ,η 的线性函数。

  有了上述坐标变换式，整体坐标系下的积分式可变换到自然坐标系下的规则化区域内进行，二维问题对面积微元 dAdAdA 进行积分，三维问题对体积微元进行积分 dV\mathrm{d} VdV，表示如下：

∫Af(x,y)dA=∬Af(x,y)dxdy=∫−11∫−11f(x(ξ,η),y(ξ,η))∣J∣dξdη(5-14)\int_{A} f(x, y) \, \, \mathrm{d} A=\iint_{A} f(x, y) \, \, \mathrm{d} x \mathrm{d} y=\int_{-1}^{1} \int_{-1}^{1} f(x(\xi, \eta), y(\xi, \eta))\, \,|\pmb{J} |\, \, \mathrm{d} \xi \mathrm{d} \eta \tag{5-14} ∫A​f(x,y)dA=∬A​f(x,y)dxdy=∫−11​∫−11​f(x(ξ,η),y(ξ,η))∣JJJ∣dξdη(5-14)

∫Af(x,y)dA=∫−11∫−11g(ξ,η)∣J∣dξdηg(ξ,η)=f(x(ξ,η),y(ξ,η))(5-15)\begin{aligned} \int_{A} f(x, y) \, \, \mathrm{d} A &=\int_{-1}^{1} \int_{-1}^{1} g(\xi, \eta) \,|\pmb{J} |\, \, \mathrm{d} \xi \mathrm{d} \eta \\[12pt] g(\xi, \eta) &= f(x(\xi, \eta), y(\xi, \eta)) \end{aligned} \tag{5-15} ∫A​f(x,y)dAg(ξ,η)​=∫−11​∫−11​g(ξ,η)∣JJJ∣dξdη=f(x(ξ,η),y(ξ,η))​(5-15)

六、单元分析

  在进行有限元分析时，单元内仍然假定为连续的、均匀的、各项同性的完全弹性体。因为单元内部仍是连续体，应按弹性力学方法进行分析。取单元各结点位移列阵： δe=(δ1,δ2,⋯,δi)T\boldsymbol {\delta_e} = (\delta_1,\;\delta_2,\;\cdots,\;\delta_i)^Tδe​=(δ1​,δ2​,⋯,δi​)T（i=i=i= 单个结点自由度数 ×单元结点个数）为基本未知量，则单元内任意一点处的位移可由结点位移表示。同样，单元内任意一点处的应变和应力均可由节点位移表示。

  任何变形连续体处于平衡状态的必要和充分条件是：对任意虚位移，外力所做的总虚功恒等于变形体所接受的总虚变形功，也即恒满足虚功方程。（随便找本有限元参考书都能查到）。在 FEM 中，用结点的平衡方程代替平衡微分方程，后者不在列出。有限元中，谈平衡为虚功方程的替换，以结点的平衡代替处处的平衡。几何方程与与物理方程严格遵守弹性力学，平衡方程不是。有限元分析削弱了弹性力学中的处处平衡，但又没有特别严重的削弱。

  有限单元法的概念，可以简述为：采用有限自由度的离散单元组合体模型去描述实际具有无限自由度的考察体，是一种在力学模型上进行近似的数值计算方法。其理论基础时分片插值技术与变分原理。有限元最一般的推导方式是等效积分的伽辽金弱形式。有限单元法的分析过程：结构离散 → 单元分析 → 整体分析。

  单元分析的主要内容：

    1. 应用插值公式，由单元结点位移 δe\boldsymbol {\delta_e}δe​，求单元的位移函数 d\boldsymbol {d}d，这个插值公式称为单元的位移模式，为： d=N⋅δe\boldsymbol {d}=\boldsymbol {N} \cdot \boldsymbol {\delta_e}d=N⋅δe​；

    2. 根据弹性力学几何方程，由单元的位移函数 d\boldsymbol {d}d，求出单元的应变，表示为 ε=B⋅δe\boldsymbol {\varepsilon}= \boldsymbol {B} \cdot \boldsymbol {\delta_e}ε=B⋅δe​，其中 B=AT⋅N\boldsymbol {B} = \boldsymbol {A}^T \cdot \boldsymbol {N}B=AT⋅N ；

    3. 根据弹性力学物理方程，由单元的应变 ε\boldsymbol {\varepsilon}ε，求出单元的应力，表示为 σ=D⋅B⋅δe=S⋅δe\boldsymbol {\sigma}= \boldsymbol {D} \cdot \boldsymbol {B} \cdot \boldsymbol {\delta_e}=\boldsymbol {S} \cdot \boldsymbol {\delta_e}σ=D⋅B⋅δe​=S⋅δe​ ；

    4. 根据虚功方程，由单元的应力 σ\boldsymbol {\sigma}σ， 求出单元的节点力，表示为 Fe=k⋅δe\boldsymbol {F_e} = \boldsymbol {k}\cdot \boldsymbol {\delta_e}Fe​=k⋅δe​ ；

    5. 按照虚功等效原则，将作用在单元上的非结点力力等效到单元结点上，得到单元等效结点荷载 FEe\boldsymbol {F_E^e}FEe​ 。

  以上各式中，A\pmb{A}AAA 为微分算子矩阵；B\pmb{B}BBB 为应变矩阵；D\pmb{D}DDD 为弹性矩阵；S\pmb{S}SSS 为应力转换矩阵。另外，弹性动力分析与弹性静力分析并无本质上的不同，只是在列单元平衡方程时，惯性力和阻尼力会作为体积力的一部分出现在平衡方程中。利用虚位移原理或势能原理的充分性，进行单元特性分析，可建立单元的运动方程/动平衡方程（列平衡方程的方法有很多，主要有：虚功原理法、能量变分原理法和伽辽金方法等），可得到单元的的质量矩阵和刚度矩阵，分别按式（6-1）、式（6-2）计算。

Me=∫VρNTNdV(6-1)\pmb {M}_e = \int_V \rho \, \pmb{N}^T \, \pmb{N} \, {\rm d}V \tag{6-1} MMMe​=∫V​ρNNNTNNNdV(6-1)

Ke=∫VBTDBdV(6-2)\pmb {K}_e = \int_V\pmb{B}^T \, \pmb{D} \, \pmb{B} \, {\rm d}V \tag{6-2} KKKe​=∫V​BBBTDDDBBBdV(6-2)

式中，Ke\pmb {K}_eKKKe​ 为单元在局部坐标系下的刚度矩阵；ρ\rhoρ 为材料得质量密度；Me\pmb {M}_eMMMe​ 称为协调质量矩阵或一致质量矩阵，这是因为导出它时，和导出刚度矩阵所根据得原理（伽辽金方法）及所采用位移插值函数是一致的，它是单元在局部坐标系下的质量矩阵。此外，在有限元法中还经常采用所谓集中质量矩阵。它假定单元的质量集中在结点上，这样得到的质量矩阵是对角矩阵。将单元协调质量矩阵 Me\pmb {M}_eMMMe​ 转换为单元集中质量矩阵，即对 Me\pmb {M}_eMMMe​ 进行对角化的方法，有多种方法可供选择。在实际分析中，协调质量矩阵和集中质量矩阵都有应用，一般情况下，两者给出的结果也相差不多。

  等参元通常也以位移作为基本未知量，因此在用最小势能原理变分得到的有限元一般格式对等参元同样适用。差别在于等参元的插值函数是用自然坐标给出的，等参元的一切计算都是在自然坐标系中形状规则的母单元内进行。单元矩阵计算时，可将式（6-1）和式（6-2）中的被积函数 N\pmb {N}NNN、B\pmb {B}BBB 等表示成自然坐标的函数，同时 dV{\rm d}VdV 及定积分上下限做相应变换。则对于任意 4 结点四边形等参元，有：

Me=t∫−11∫11ρNTN∣J∣dξdη(6-3)\pmb {M}_e= t\int_{-1}^{1} \int_{1}^{1} \rho \, \pmb{N}^T \, \pmb{N} \, |\pmb{J}| \, \mathrm{d} \xi \mathrm{d} \eta \tag{6-3} MMMe​=t∫−11​∫11​ρNNNTNNN∣JJJ∣dξdη(6-3)

Ke=t∫−11∫11BTDB∣J∣dξdη(6-4)\pmb {K}_e=t\int_{-1}^{1} \int_{1}^{1} \pmb{B}^T \, \pmb{D} \, \pmb{B} \, |\pmb{J}| \, \mathrm{d} \xi \mathrm{d} \eta \tag{6-4} KKKe​=t∫−11​∫11​BBBTDDDBBB∣JJJ∣dξdη(6-4)

∣J∣=∣∂x∂ξ∂y∂ξ∂x∂η∂y∂η∣=∣∑∂Ni∂ξxi∑∂Ni∂ξyi∑∂Ni∂ηxi∑∂Ni∂ηyi∣=aξ+bη+c(6-5)|\pmb{J}| = \left|\begin{array}{ll} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\[8pt] \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{array}\right|=\left|\begin{array}{cc} \sum \frac{\partial N_{i}}{\partial \xi} x_{i} & \sum \frac{\partial N_{i}}{\partial \xi} y_{i} \\[8pt] \sum \frac{\partial N_{i}}{\partial \eta} x_{i} & \sum \frac{\partial N_{i}}{\partial \eta} y_{i} \end{array}\right| = a\xi+b\eta+c \tag{6-5} ∣JJJ∣=∣∣∣∣∣∣​∂ξ∂x​∂η∂x​​∂ξ∂y​∂η∂y​​∣∣∣∣∣∣​=∣∣∣∣∣∣​∑∂ξ∂Ni​​xi​∑∂η∂Ni​​xi​​∑∂ξ∂Ni​​yi​∑∂η∂Ni​​yi​​∣∣∣∣∣∣​=aξ+bη+c(6-5)

七、数值积分

  等参元的单元质量矩阵、单元刚度矩阵和等效结点荷载列阵的计算都涉及定积分的计算，如式（6-3）、式（6-4）所示。由于一般情况下计算公式中被积函数十分复杂，很难用精确积分得到显式积分结果。因此，需要采用数值积分方法计算定积分，即在单元内选出某些点（称为积分点），算出被积函数在这些点处的值，再分别乘以权系数，然后以求其和作为定积分的近似值。

  数值积分（Numerical Integration）主要用于求定积分的近似值，用数值逼近的方法近似计算给定的定积分值。借助于电子计算设备，数值积分可以快速而有效地计算复杂的定积分。数值积分方法有很多，如梯形法、辛普森法和牛顿-柯斯特法等等。在有限元分析中，通常采用高斯积分法，因为它可以用较少的积分点达到较高的精度，从而可以节省计算时间。另外，在弹性力学有限元中，各种定积分的计算都有一个显著特点即被积函数均为多项式，这是有形函数的假定确定的。（形函数假定的就是多项式）

7.1 积分公式

  若函数 f(x)f(x)f(x) 在区间 [a,b][a, b][a,b] 上定积分存在，且定积分值为 III，即：

I=∫abf(x)dx(7-1)I = \int_{a}^{b} f(x) \mathrm{d} x \tag{7-1} I=∫ab​f(x)dx(7-1)

  若被积函数 f(x)f(x)f(x) 为初等函数，则容易用牛顿-莱布尼兹公式计算定积分 III 。当 f(x)f(x)f(x) 的形式非常复杂时，其原函数很难求得。此时，可采用如下所示的数值积分方法计算：

I=∫abf(x)dx=∑k=1nAkf(xk)+R(7-2)I = \int_{a}^{b} f(x) \mathrm{d} x = \sum_{k=1}^{n} A_{k} f\left(x_{k}\right) + R \tag{7-2} I=∫ab​f(x)dx=k=1∑n​Ak​f(xk​)+R(7-2)
I=∫abf(x)dx≈∑k=1nAkf(xk)(7-3)I = \int_{a}^{b} f(x) \mathrm{d} x \approx \sum_{k=1}^{n} A_{k} f(x_k) \tag{7-3} I=∫ab​f(x)dx≈k=1∑n​Ak​f(xk​)(7-3)

其中，RRR 为截断误差；nnn 为积分点个数，xkx_kxk​ 为积分点坐标，AkA_kAk​ 为求积系数，亦称伴随节点 xkx_kxk​ 的权，权 AkA_kAk​ 仅仅与积分点 xkx_kxk​ 的选取有关，而不依赖被积函数 f(x)f(x)f(x) 的具体形式。通常称式（7-3）为机械求积公式，xkx_kxk​ 的选取方式不同，造就了不同的积分公式。（插值型求积公式，用多项式来近似替代原函数）

  若某个求积公式对于所有次数不超过 mmm 的多项式都能准确成立，而对于某个 m+1m+1m+1 次的的多项式不准确成立，则称该求积公式具有 mmm 次代数精度，代数精度用衡量积分公式的精确性。通常，在积分点个数相等的情况下，代数精度愈高，求积公式愈精确。梯形公式具有一次代数精度；辛普森公式具有三次代数精度；具有 n+1n+1n+1 个求积节点的牛顿-柯斯特求积公式，至少具有 nnn 次代数精度，但当 nnn 较大时，由于 Runge 现象,收敛性也无法保证，故一般不采用高阶的牛顿-科特斯求积公式。

7.2 高斯积分

  定义：有 nnn 个节点的求积公式最高可具有 2n−12n-12n−1 次代数精度，称这种高精度的求积公式为高斯（Gauss）型求积公式，节点 xix_ixi​ 为 Gauss 点，伴随高斯点 xix_ixi​ 的加权系数为 WiW_iWi​。那么，根据高斯积分法，有：

I=∫abf(x)dx=∑i=1nWif(xi)+R(7-4)I = \int_{a}^{b} f(x) \mathrm{d} x = \sum_{i=1}^{n} W_{i} f\left(x_{i}\right) + R \tag{7-4} I=∫ab​f(x)dx=i=1∑n​Wi​f(xi​)+R(7-4)

  由于高斯积分点 xix_ixi​ 是勒让德（Legendre）多项式的根，因此对于各不同积分点数时，容易求得相应积分点的坐标，进而求得相对应的权重系数。在实际应用中，只需查表即可，下表列出了对于 (−1,1)(-1, 1)(−1,1) 积分域 nnn = 1 ~ 6 的积分点位置 ξi\xi_iξi​ 和权系数 HiH_iHi​ 的值。

高斯积分的积分点坐标和权系数

7.2 积分阶次

  完全积分与缩减积分。

八、UEL的实现

  通常，弹性力学问题最终都会归结为偏微分方程组的初边值问题。若偏微分方程组在给定初值条件和边值条件下得到解答，则与之对应的弹性力学问题亦得到解答。然而，遗憾是对于复杂问题，偏微分方程组很难得到精确的解析解答。我们退而求其次，希望能得到高精度的数值解答。

  有限单元法是求解偏微分方程数值解（高精度近似解）的一种方法。有限单元法有严格的数学证明作为其近似的客观性和合理性的保证。有限单元法不是从问题的微分方程和相应的定解条件出发，而是从与其等效的积分形式出发。等效积分的一般形式是加权余量法，它适用于普遍的方程形式。利用加权余量法的原理，可建立多种近似方法，如伽辽金法。如果原问题的微分方程具有某些特定的性质，则它的等效积分形式的伽辽金法可以归结为某个泛函的变分。相应的近似解法实际上是求解泛函的驻值问题。有限单元法区别于传统的加权余量法和求解泛函驻值的变分法，该法不是在整个求解域上假设近似函数，而是在各个单元上分片假设近似函数（形函数），这样就克服了在全域上假设近似函数所遇到的困难，是近代工程数值分析方法领域的重大突破。有的限元理论基础是分片插值技术与变分原理。有限元最一般的推导方式是等效积分的伽辽金弱形式。其在数学上坚实的理论分析与误差分析见本文第三章及相关拓展内容。

  在弹性力学领域，其经典的逆解法与半逆解法，具有很大的局限性，复杂问题需要借助有限单元法进行求解。弹性力学有限元的基本思想是将复杂的体系离散为有限个单元和结点，各单元在结点处连结。在外力作用下，保证每个单元受力平衡，每个结点处受力平衡，则体系一定受力平衡。单元分析保证每个单元都平衡，结点分析保证整个结构体系受力平衡。

核反应堆的高度复杂模型

  变形连续规律、应力-应变关系和运动（或平衡）规律是弹性的力学三大基本规律。弹性力学中的有限元分析，单元内及单元间仍然满足弹性力学的基本假定，即：连续、均匀、各项同性，线弹性、小变形。因此，几何方程（变形连续规律）和物理方程（应力-应变关系）对于有限元分析依旧精确成立。弹性力学问题之所以很难求解，一个方面是弹性力学对平衡要求的太过严格，弹性力学中的平衡指的是微元体的平衡，要求弹性体内任意一点处，处处平衡。而有限单元法放松了对平衡的要求。有限元的研究对象是单元和结点，它只要求单元上全部外力平衡和各结点全部平衡。

  对单元分析进行的方法有很多，如：牛顿第二定律、拉格朗日方程等。有限单元法采用能量原理 (虚位移原理) 进行单元分析。因而必须提前给定位移函数的具体形式。一般而言，位移函数的选取影响计算结果的精度。在弹性力学中，恰当的选取位移函数是相当不容易的。有限单元法中，当单元划分的足够小时，把位移函数设定为简单的多项式形式也可以得到具有相当精度的结果。无论何种形式的位移函数，位移函数中必须包含能反应刚体位移和常应变状态的常数项和一次项，该条件构成了单元的 完备性准则。此外，位移函数在单元内要连续，在相邻单元间应尽量保持协调 (即单元间位移函数也应尽量保持连续)，这是单元的 协调性条件。理论和实践都已证明，完备性准则是有限元解收敛于真实解的必要条件。

  在单元分析中，从单元的结点位移 → 求位移分布/位移函数/插值函数（形函数矩阵） → 求应变 (几何方程/应变矩阵) → 求应力 (物理方程/应力矩阵) → 求结点力 (单元间相互作用力)，为单元的内力分析；外荷载移置到结点荷载，为单元的外力分析。单元分析的主要目的是确定单元对结点的反作用力。有限元中所有的荷载统统处理到结点上去。 以单个结点为研究对象，则作用在结点上的荷载主要有：连接该结点的各单元作用在该结点上的反作用力 (包括单元结点力的反作用力和单元等效结点荷载) 和直接作用在该结点上的外载。作用在分析结点上的各力组成了 汇交力系。单元对结点的作用力与结点对单元的作用力为作用力与反作用力关系，根据牛顿第三定律，作用力与反作用力，等大反向属于同种性质的力。根据各结点的受力平衡，可得到整个系统的平衡方程，这是有限元分析的通俗易懂的思路。对于复杂结构体系，我们有更高级的方式对结构整体进行分析即创建整体的平衡方程，实际上是各结点的平衡方程，该方法是最小势能原理（势能驻值原理），值得说明的是势能原理与虚位移原理彼此是完全等价的。通过计算整个结构的势能，并使势能的一阶变分为 0 ，便可得到结构的整体质量矩阵和刚度矩阵。势能是泛函（函数的函数），泛函一阶变分（导数）为 0，则泛函取得驻值。

  用势能原理进行结构整体分析并不经济，目前在有限元程序中一般采用直接刚度法获得结构原始刚度矩阵 KKK 和质量矩阵 MMM，由单元刚度矩阵的贡献来获得结构原始刚度矩阵，一般采用先处理法集成结构刚度矩阵，结构边界条件（主要指零位移条件）在形成结构刚度方程前就进行处理。通过上述分析，我们便可得到（线弹性）多自由度系统的运动微分方程，统一表示为矩阵形式：

Mx¨+Cx˙+Kx=P(t)(8-1)\pmb{M}\pmb{\ddot{x}} ~~ + ~~ \pmb{C}\pmb{\dot{x}} ~~ + ~~ \pmb{K} \pmb{{x}} ~~ =~~ \pmb{P}(t) \tag{8-1} MMMx¨x¨x¨  +  CCCx˙x˙x˙  +  KKKxxx  =  PPP(t)(8-1)

式中，M\pmb{M}MMM、C\pmb{C}CCC 、 K\pmb{K}KKK 和 P(t)\pmb{P}(t)PPP(t) 分别为质量矩阵、阻尼矩阵、刚度矩阵和结点荷载列阵；x\pmb{x}xxx、x˙\pmb{\dot{x}}x˙x˙x˙ 和 x¨\pmb{\ddot{x}}x¨x¨x¨ 分别为位移向量、速度向量和加速度向量。式（8-1）为系统强迫振动的运动方程，由式（8-1）可以看出，作用在体系上的外荷载主要由虚加的惯性力、阻尼力和弹性恢复力来平衡。

  系统的质量矩阵 M\pmb{M}MMM、阻尼矩阵 C\pmb{C}CCC 、刚度矩阵 K\pmb{K}KKK 和结点荷载列阵 P(t)\pmb{P}(t)PPP(t) ，是分别由各自的单元矩阵和向量集成的，即：

M=∑MeC=∑CeK=∑KeP=∑Pe(8-2)\begin{array}{ll} \pmb{M}=\sum \pmb{M}_{e} & \pmb{C}=\sum \pmb{C}_{e} \\[12pt] \pmb{K}=\sum \pmb{K}_{e} & \pmb{P}=\sum \pmb{P}_{e} \end{array} \tag {8-2} MMM=∑MMMe​KKK=∑KKKe​​CCC=∑CCCe​PPP=∑PPPe​​(8-2)

其中

Me=∫VeρNTNdVCe=∫VeμNTNdVKe=∫VeBTDBdVPe=∫VeNTfdV+∫SσNTTdV(8-3)\begin{aligned} \pmb{M}_e &= \int_{\rm{V}_{e}} \rho \, \pmb{N}^{\rm{T}} \pmb{N} \, \rm{d} V \\[12pt] \pmb{C}_e &= \int_{\rm{V}_{e}} \mu \, \pmb{N}^{\rm{T}} \pmb{N} \, \rm{d} V \\[12pt] \pmb{K}_e &= \int_{\rm{V}_{e}} \pmb{B}^{\rm{T}} \pmb{D} \pmb{B} \, \rm{d} V \\[12pt] \pmb{P}_e &= \int_{\rm{V}_{e}} \pmb{N}^{\rm{T}} \pmb{f} \, \rm{d} V + \int_{\rm{S}_{\sigma}} \pmb{N}^{\rm{T}} \pmb{T} \, \rm{d} V \end{aligned} \tag{8-3} MMMe​CCCe​KKKe​PPPe​​=∫Ve​​ρNNNTNNNdV=∫Ve​​μNNNTNNNdV=∫Ve​​BBBTDDDBBBdV=∫Ve​​NNNTf​f​​fdV+∫Sσ​​NNNTTTTdV​(8-3)

  显然，Abaqus 软件最核心的功能就是对方程组（8-1）的求解，此处，我们仅讨论的是弹性力学有限元，其余类似。单元类型的不同直接影响的是形函数矩阵，不同的单元类型具有不同的形函数矩阵，如 C3D8 与 C3D20 本质上的差别就是对单元内位移场插值函数（形函数）的定义不同。若 Abaqus 自带的单元库无法满足需求，可根据需要开发自定义单元（UEL），即创造新单元。而创造新单元的过程就是确定形函数矩阵的过程。那么，在 Abaqus 中实现 UEL 首先要做的就是在理论层面把单元分析明白，把单元形状、结点数目、结点自由度、形函数矩阵确定下来，进而单元的应变矩阵便确定下来。所以本文绝大部分内容都在介绍有限元的基本理论。

  单元类型用于确定单元的形函数矩阵 N\pmb{N}NNN 和应变矩阵 B\pmb{B}BBB，材料属性用于确定单元的弹性矩阵 D\pmb{D}DDD，它们共同确定了单元的质量矩阵 Me\pmb{M}_{e}MMMe​、阻尼矩阵 Ce\pmb{C}_{e}CCCe​ 和刚度矩阵 Ke\pmb{K}_{e}KKKe​。在给定位移边界及荷载边界的条件下，还能确定单元的荷载列阵 Pe\pmb{P}_{e}PPPe​。Abaqus 根据用户提供的各种输入参数，并根据式（8-3）自动计算单元的质量矩阵、阻尼矩阵、刚度矩阵和结点荷载列阵。每个单元的质量矩阵、阻尼矩阵、刚度矩阵和结点荷载列阵确定后，Abaqus 将它们自动组装成系统的质量矩阵、阻尼矩阵、刚度矩阵和荷载列阵，进而方程组（8-1）的具体形式便确定下来。各种矩阵的计算与集成、方程组的求解都是由 Abaqus 求解器来完成的，而计算所需的各种输入参数，均记录在 Abaqus 的 inp 文件中。

  综上所述，若采用 Abaqus 实现 UEL，必须准备两个文件，一个是记录模型基本信息的 .inp 文件，另个一是描述单元性能的 .for 文件。.inp 文件可由 Abaqus/CAE 创建，创建后将其内默认的单元类型关键字修改为自定义单元，从而实现与自定义单元的关联。而记录单元性能的 .for（Fortran）文件，则根据官方提供的接口进行创建。

  更多关于 .for 文件的创建，见 Abaqus 软件帮助文档 Abaqus User Subroutines Reference Guide >> 1. User Subroutines >> 1.1 Abaqus/Standard subroutines >> 1.1.28 UEL (User subroutine to define an element.)，如下所示。同样，关于 .inp 文件单元关键字的定义与修改，见 Abaqus 软件帮助文档 Abaqus Keywords Reference Guide >> *USER ELEMENT 。特别提示：不用特意学习 Fortran ，用到哪里看哪里即可，也不用特意了解 inp 的书写规则，查帮助文档即可。

UEL帮助文档

8.1 问题描述

  本文以一平面应变问题为例，介绍 Abaqus 实现用户自定义单元的主要过程。待分析问题的基本信息如下，矩形平板的长和宽分别为 8 m 和 5 m 。材料为 Q235，其弹性模量为 2.06 × 10¹¹ Mpa，质量密度为 7850 kg/m³，泊松比为 0.3 。板的四边固定，对该板进行直接稳态动力学分析，分析频率为 10 Hz，荷载作用在板的中心点处，大小为 2022+1.12i2022 + 1.12 i2022+1.12i，作用方向为沿 xxx 轴正方向。inp 文件的创建由 Abaqus/CAE 完成，分析采用的 4 结点四边形平面应变等参元由 UEL 实现，具体操作如下。

8.2 inp文件的创建

直接稳态分析步的创建

网格的划分与单元属性的修改

边界条件的指定与荷载的施加

创建作业并写入inp文件（KT-CPE4R.inp）

8.3 inp文件的修改

  为了实现对用户自定义单元的引用，上述创建的 inp 文件（KT-CPE4R.inp）还需要做一定修改。首先，将系统默认的单元类型 CPE4R 修改为自定义单元类型 U2022，关键字为 *User element；其次，添加自定义单元的属性，如平面应变单元的厚度，弹性模量，泊松比等，关键字为 *Uel property，如下所示。显然，阻尼等参数可以从 *UEL PROPERTY 关键字中定义。另外，在 inp 文件中，两个 * 号开头的数据行为注释行。

单元类型的修改（*USER ELEMENT）

  关于关键字 *User element 后各参数的定义见帮助文档，其下的数据行的含义如下图中的红色方框，其实就是结点自由度。

单元属性的修改（*UEL PROPERTY）

8.4 for文件的创建

  用于定义自定义单元性能的 .for 文件的基本格式如下图所示，其中仅红色方框部分允许用户自由填充。用户在这部分定义形函数矩阵、高斯积分点、质量矩阵等等。另外，Fortran 语言中字符 C 为整行注释，! 为行内注释。

KT-CPE4R.for文件的创建

8.5 作业的提交计算

作业的提交

8.6 结果的后处理

计算结果（仅显示结点）

  Python语言创建 odb 文件，添加自定义场变量，主要思路：Displacements can be visualized, but only with pointers as can be seen above. The other variables like temperature or stresses can not be displayed. In the user manual of Abaqus is written that plotting User-elements is not supported, but, if the User-element contain displacement degrees of freedom, they can be overlaid with standard elements with no stiffness. With this method, only the displacements can be visualized. In order to visualize also stresses, temperatures, etc., use is made of a Python script. This script first reads the input file (.inp) and extracts there the information of the node-locations and which node-series form elements. Secondly it reads the result file (.dat) and extracts here the values for the different variables e.g. displacements, stresses, temperatures, etc. The third step is to create a standard Abaqus result file containingthe values of the User-element job. This new *.odb file will be created with standard

九、参考文献

[1]. ABAQUS非线性有限元分析与实例. 庄茁等.

[2]. ABAQUS的UEL开发技术在结合部特性分析中的应用. 刘超.

[3]. ABAQUS UEL 计算过程浅谈. 窗台上的叔本华.

[4]. 有限元分析——ANSYS理论与应用. Saeed Moaveni 等.

[5]. 等效积分的“弱”形式 . 啊门的博客.

[6]. 有限单元法 . 王勖成.

[7]. 有限单元法教程. 王焕定，王伟 编著.

[8]. 祝 KangTong 师弟毕业顺利！

若本文对你有所帮助，可微信扫码打赏
你的鼓励将是我创作的最大动力
如有问题，欢迎邮件交流
liyang@alu.hit.edu.cn
祝各位攻城狮们，天天开心，一切顺利！

Abaqus 用户子程序 UEL相关推荐

1. ABAQUS用户子程序一览表

   说明 ABAQUS用户子程序一览表 ABAQUSStandard subroutines Refence 说明 本系列文章本人基本没有原创贡献,都是在学习过程中找到的相关书籍和教程相关内容的汇总和梳理 ...

2. Abaqus用户子程序umat的学习

   Abaqus用户子程序umat的学习 说明:在文件中,!后面的内容为注释内容.本文为学习心得,很多注释是自己摸索得到.如有不正确的地方,敬请指正. ! ------------------------ ...

3. matlab有限元分析与应用_专栏 | UEL用户子程序开发步骤—有限元理论基础及Abaqus内部实现方式研究系列20...

   作者介绍 snowwave02 博士,高级工程师 snowwave02团队:设计仿真领域的软件开发团队,由软件.机械.物理等专业人员组成,10年以上CAD/CAE软件开发经验,精通Abaqus二次开发 ...

4. 【JY】 ABAQUS子程序UEL的有限元原理与应用

   不等待 即关注 [简述ABAQUS中UEL子程序] ABAQUS作为成熟的商用有限元软件,可为高级用户提供特定的分析需求.ABAQUS常见的二次开发子程序包括:UMAT.VUMAT.UGENS.UEL ...

5. ABAQUS学习记录1——用户子程序综述

   概述 ABAQUS提供了相当丰富的单元类型,材料属性等数据库可供用户选择,但是工程问题是千变万化的,为了满足用户的特殊工程要求,ABAQUS为用户提供了强大而又灵活的用户子程序接口(USER SUBR ...

6. abaqus6.13+vs2012+ivf2013用户子程序关联步骤

   abaqus6.13+vs2012+ivf2013用户子程序关联步骤 之前笔者在论坛上也是找了很多方法和关联步骤,许多步骤都十分繁琐,不容易懂,实在看不懂作者的心思.下面简要介绍怎么操作: **1.* ...

7. adams c语言,Adams2013编译C语言用户子程序

   1.操作系统:Windows xp 32位 2.软件版本:Adams 2013 32位.Visual Studio 2010专业版 32位 3.编译软件:Intel Visual Fortran 11 ...

8. ABAQUS材料子程序学习（20年12月2日）

   @ABAQUS材料子程序学习(20年12月2日) 前言 继续记录自己学习过程,本文针对<非线性本构关系在ABAQUS中的实现>第三章"黏弹性"本构的学习,UMAT子程序 ...

9. ABAQUS材料子程序学习（线性各向同性硬化塑性）

   ABAQUS材料子程序学习(线性各向同性硬化塑性) 前言 塑性力学增量形式实现 umat子程序 参数 计算结果 前言 记录自己学习abaqus软件umat子程序的t过程,本文主要参考了<非线性本 ...

10. c语言的子程序文件名,Adams2012编译C语言用户子程序

    1 概述 在Adams使用过程中,有些复杂的情况特别是涉及到一些逻辑表达,用函数表达式很难表达出来,这种情况需要使用用户子程序. Adams用户子程序支持C语言和Fortran语言,随着C语言的普及, ...

最新文章

1. 2014家电盘点：求变与创新
2. 【安全技术】关于几种dll注入方式的学习
3. C/C++实现快速排序
4. gmp会屏蔽掉base里的几个函数
5. 一些web开发中常用的、做成cs文件的js代码 - 搜刮来的
6. ws flv连接播放得测试代码
7. 【重磅】世界上最可信、最权威的人工智能数据和洞察来源：2021年人工智能指数报告...
8. resin session共享 redis_Spring Boot 利用Redis实现session共享
9. 文献综述_软件单元测试
10. Ionic 创建打包项目
11. EPICS -- asynRecord记录使用示例
12. n(n-1)表示什么？n(-n)表示什么？
13. 2020 年全国房价排行榜出炉
14. diagram使用(BLOCK DIAGRAM)
15. jsp简介及工作原理
16. node.jshe npm的区别
17. windows7英文版，变为中文版
18. P1489 猫狗大战
19. 年终总结——过去已逝，未来可期不可欺
20. html怎么做密码的判断,用户密码格式判断 .html

热门文章

1. jdk8安装和环境变量配置
2. android so文件解密器,【Android 原创】so文件动态加解密的CrackMe
3. netty 百度网盘 密码
4. steam遇到错误代码解决方案
5. python键值对是什么意思_python键值对
6. 轻松解决加密的PDF如何编辑简单技巧
7. QT5.1.0，QT4.8.0以及VC2010、VC2012的测试对比
8. c语言 error c2562,C语言之关键字（二） void,const
9. 结构体Sqlist L与Sqlist L的区别
10. iPad协议 一键转发 群发消息 获取群二维码 清理僵尸粉V:viplac