开放科学(资源服务)标识码(OSID): 
0 引言
针对有限元法[1]的网格剖分难题、计算稳定性和计算精度控制等问题,笔者借助数学流形思想,提出了基于连续介质的偏微分方程的数值计算新方法——独立覆盖流形法[2-3],在有限元法的插值方式之外,提出构造近似函数的新方式:求解域划分如图1所示,通过“分区”(图1(a)中各分块区域)降低区域内真实物理场分布的复杂度,在相邻分区之间自动生成窄“条带”的覆盖重叠区域(图1(b)中每个分区加上其周边的窄条带称为一个覆盖);以各覆盖的独立区域(称为独立覆盖,即覆盖中除去条带的部分)为主进行分析,采用完备级数(常用多项式级数)作为局部近似函数;在窄条带的覆盖重叠区域,用有限单元的线性形函数作为单位分解函数,将局部近似函数(级数)连接成连续的整体近似函数;采用伽辽金法的标准求解流程,获得各分区的级数系数,形成“分区级数解”作为偏微分方程的数值解答。
在各覆盖区域,由分区级数逼近真实解,通过覆盖重叠保证整体近似函数的连续性。由“分区”形成的覆盖网格具有任意形状、任意连接和任意加密的特性(图1(a)),有望从根本上解决数值计算的网格剖分难题,还能模拟实际求解域的精确几何边界并严格施加本质边界条件[4]。在“级数解”的选择上,可应用局部解析解模拟物理场的局部特性以加快收敛,比如用裂纹尖端(裂尖)级数解处理裂尖奇异性[2,5],或用完备级数的适当形式反映其整体特性,比如模拟梁板壳变形特征(无需推导梁板壳控制方程及相应的数值计算公式),避免了梁板壳有限单元的各种问题(如C1近似函数构造困难、自锁问题、不易与实体单元连接等),并实现了精确几何的曲梁和曲壳分析[6⇓-8]。
在各分区内采用完备级数,在理论上保证了随级数升阶而稳定逼近局部真实物理场,使该方法具有严格的收敛性(包括物理场导数的收敛)[3]。采用多项式级数时,各阶次的完全多项式便于计算精度的度量。无论采取何种方式,只要覆盖网格细分到适当大小,就一定能用适当阶次的完备级数逼近真实解,从而保证了计算结果的可靠性。根据以上收敛原理提出逐级一分为二的网格分裂算法[9],并初步推广到三维分析[10]。在二维结构线弹性静力分析中,配合误差估计进行h-p型混合自适应分析,初步实现了与CAD融合的自动计算[11]:只需在CAD中人工输入必要的几何信息和计算条件及参数,其他所有计算过程(包括后处理)无需人工参与。另外还尝试了材料非线性[12]和几何非线性分析(包括在固定网格中求解大变形问题[13]),初步研究了流体类型的微分方程求解[14-15],均取得了很好的效果。
独立覆盖流形法避免了有限元法常用的低阶插值方式所存在的过度刚性、各种自锁[1]等缺陷,以及采取相应措施后带来的非协调性、零能模式[1]等问题。另外,导数自然收敛,无需应力重构。该方法具备了常规数值流形方法[16]、广义有限元法[17]的高阶级数逼近性能,又避免了这些方法的方程组线性相关问题,相比低阶单元收敛快;采用完备的Williams级数分析裂纹问题,比扩展有限元[18]收敛更快,并能直接得到裂尖应力强度因子;在等几何分析方法[19]之外提出了几何保形性的新途径,可适用于任意几何外形;在网格划分和网格细化的便利性和灵活性上堪比无网格法[20],并能实现本质边界条件(含曲线、曲面边界)的严格施加;解决了间断伽辽金法[21]在相邻单元之间的不连续问题,推进了高阶级数单元以及传统解析级数的实际应用。但该方法高阶单元的计算复杂,带宽大,因此相比有限元法的计算量增大。
1 偏微分方程组及边界条件
1.1 偏微分方程组
$\boldsymbol{L} \boldsymbol{u}+\boldsymbol{f}=0 \text {, 在 } \boldsymbol{\Omega} \text { 上。 }$
(1)弹性力学运动微分方程
其中,u=[u v w]T为待求的位移场矢量,e= + + ,▽2u= + + ,f=[fx fy fz]T为体力荷载。
式中:E为弹性模量;μ为泊松比;ρ为密度; 分别为体力荷载的x、y、z方向的分量;t为时间。
(2)传导方程
。
式中:u为待求的标量物理场(如温度场、渗流场等);f为非齐次项;a2▽2u也可以换成kx +ky +kz 的形式,a、kx、ky、kz为参数;类比式(2)和式(3),可知式(3)中的时间导数项中ρ=1。
(3)波动方程
。
式中u为待求的标量物理场,可以是声场或电场、磁场中的某个标量。
上述方程若去掉时间导数项,则分别成为稳态求解的弹性力学平衡微分方程、泊松方程(去掉非齐次项为拉普拉斯方程)。
上述偏微分方程涉及3个控制参数:偏微分方程组的方程个数npde(即未知物理量的个数);计算维数ndim;关于坐标的2阶微分项的相关参数矩阵(如弹性矩阵)的维数nb。式(2)当npde>1时(二维或三维),npde=ndim:npde=2时分为平面应力和平面应变两种情况,nb=3;npde=3时,nb=6。对于式(3)和式(4),npde=1,还要定义ndim,且nb=ndim。
以下主要介绍偏微分方程式(2)的数值计算公式,很容易退化为偏微分方程式(3)和式(4)的计算公式。同时一般以ndim=3即三维形式表述,也容易退化到二维和一维形式。
限于篇幅,本文仅考虑了以上常见的三种偏微分方程。实际上很容易加入0阶微分项u,如赫姆霍兹方程。流体类的偏微分方程如对流-扩散方程、Burgers方程、Navier-Stokes方程,涉及坐标的一阶偏导数项,将另文介绍。
1.2 边界条件
(1)第一类边界条件
$u=u_{0} ; \quad v=v_{0} ; w=w_{0} \text {, 在 } \Gamma_{u} \text { 上 }$。
其中,Γu为位移边界,u0、v0、w0为给定值。式(3)和式(4)仅考虑u=u0。
(2)第二类边界条件
对于式(2)
$\sigma_{i j} n_{j}=q_{i} \text {, 在 } \Gamma_{\sigma} \text { 上 }$。
对于式(3)和式(4)
$\frac{\partial u}{\partial n}=q_{0} \text {, 在 } \Gamma_{q} \text { 上 。 }$
式中:Γσ为应力边界;σij为应力;nj为方向余弦;σijnj应用了爱因斯坦求和约定,对于ndim=3,i=1,2,3,j=1,2,3;qi为给定值;Γq为导数边界;q0为给定值。
(3)第三类边界条件(混合边界条件)
对于式(2)为弹性边界条件
其中,在Γσ上定义局部坐标系x'=(x',y',z'),具体表述为t-s-n(法向n(指向体内为正)和切向t、s),相应的位移为u'、v'、w'。
式中:τt、τs为定义在局部坐标系下的切向应力;σn为定义在局部坐标系下的法向应力;u'0、v'0、w'0为给定的边界位移值;kt、ks为切向刚度系数;kn为法向刚度系数。
对于式(3)和式(4),有
。
式中:k为参数;u0为给定值。
1.3 初值条件
u|t=0=bu ;v|t=0=bv ;w|t=0=bw ,
。
式中bu、bv、bw、hu、hv、hw为给定值。
式(3)和式(4)仅考虑u|t=0=bu和 =hu。
采用伽辽金法将以上的偏微分方程(组)及其边界条件转换到数值计算的具体算法和公式,此过程与有限元法类似,仅仅是局部近似函数不同,采用的是级数形式而不是插值函数。
2 局部近似函数及其偏导数
2.1 局部近似函数与形函数
本文涉及真实存在的各种复杂几何情况,需要定义多种坐标系,首先约定不同坐标系的表示方式为:整体直角坐标x,曲线或曲面参数坐标 ,积分坐标x',覆盖坐标 ,级数坐标 。
对于偏微分方程式(2),覆盖i的位移局部近似函数在覆盖坐标系下的通用表达式为
$\begin{aligned}\boldsymbol{\varphi}_{i} \tilde{\boldsymbol{u}}_{i}=\boldsymbol{\varphi}_{i}\left\{\begin{array}{l}\tilde{u} \\\tilde{v} \\\tilde{w}\end{array}\right\}_{i}= & \boldsymbol{\varphi}_{i} \sum_{m_{i}=1}^{n_{i}}\left[\begin{array}{lll}N_{11} & N_{12} & N_{13} \\N_{21} & N_{22} & N_{23} \\N_{31} & N_{32} & N_{33}\end{array}\right]_{m_{i}}\left\{\begin{array}{l}a \\b \\c\end{array}\right\}_{m_{i}}= \\& \sum_{m_{i}=1}^{n_{i}} \boldsymbol{\varphi}_{i} \boldsymbol{N}_{m_{i}} \boldsymbol{d}_{m_{i}}。\end{aligned}$
式中:φi为单位分解函数,独立覆盖的φi=1; 为覆盖坐标下的位移;ni为级数项数; =[a b c 为待求的第mi项级数的系数向量;Nij(i=1,2,3;j=1,2,3)为基函数矩阵中的元素; 为该级数项的基函数矩阵(以下简称级数矩阵),以多项式级数为例,最常见的是对角矩阵 =I (I为单位矩阵),
$N_{m_{i}}=\bar{x}^{n-l-k} \bar{y}^{l} \bar{z}^{k}$。
其中,级数公式用 =( , , )表达, 称为级数坐标,n=0,…,pi, k=0,…,n,l=0,…,(n-k),pi为多项式的最高阶次。式(13)去掉 相关项,取k=0退化为二维;再去掉 相关项,取l=0退化为一维。可见,级数矩阵 ( )由级数坐标 确定。
多项式级数的级数坐标定义为
。
式中:( , , )为覆盖的中心坐标; 、 和 分别表示该覆盖在各方向尺度的一半(对于各方向相差不大的情况,也可采用统一尺度)。
。
三维情况下,位移由覆盖坐标系转化到整体坐标系。
。
式中:Li( )为坐标转换矩阵,由覆盖坐标 确定;cos 、cos 、cos 分别是 轴关于整体坐标x、y、z的方向余弦,其他以此类推。
条带中的位移表达式为(仅考虑2个分区之间的窄条带,举例如图1(b)所示的红线两侧)
。
$\begin{array}\varphi_{1}=\frac{h_{2}+r_{1}}{h_{1}+h_{2}}, \quad \varphi_{2}=1-\frac{h_{2}+r_{1}}{h_{1}+h_{2}},\\（厚度为 h_{1} 的半条带）；\\\varphi_{1}=1-\frac{h_{1}+r_{2}}{h_{1}+h_{2}}, \quad \varphi_{2}=1-\frac{h_{1}+r_{2}}{h_{1}+h_{2}},\\（厚度为 h_{2} 的半条带）。\end{array}$
可见,φi(x')由条带积分坐标x'的法向坐标确定。
边界条带中的边界覆盖,其覆盖坐标 等同于上述积分坐标x',因此采用x'表达的多项式级数,级数坐标$\overline{\boldsymbol{x}}=(\bar{x}, \bar{y})$为
$\bar{x}=\frac{\tilde{x}-\tilde{x}_{i}}{l_{\tilde{x}_{i}}}=\frac{x^{\prime}-x_{i}^{\prime}}{l_{x^{\prime} i}}, \quad \bar{y}=\frac{\tilde{y}-\tilde{y}_{i}}{l_{\tilde{y}_{i}}}=\frac{y^{\prime}}{h}$。
式中x'i为边界覆盖的中心坐标。去掉不含y'的多项式级数项,且常数项等于边界条件的数值,则剩下的各项组合起来就会自动满足y'=0处(即边界上)的边界条件。
三维情况的条带和边界条带与二维情况类似,相邻分区的界面(或边界面)可以是平面或曲面。通过界面(或边界面)及其法向z'定义条带的积分坐标x'=(x',y',z')。曲面情况,曲面上的点(x0,y0,z0)由参数方程$\left\{\begin{array}{l}x_{0}=x_{0}(\hat{x}, \hat{y}) \\y_{0}=y_{0}(\hat{x}, \hat{y}) \\z_{0}=z_{0}(\hat{x}, \hat{y})\end{array}\right.$确定,$\hat{\boldsymbol{x}}=(\hat{x}, \hat{y})$为参数坐标,则积分坐标$\boldsymbol{x}^{\prime}=\left(\hat{x}, \hat{y}, z^{\prime}\right)$。
综上所述,对于第i个覆盖中的第mi级数项,局部近似函数中的形函数子矩阵 统一表达为
。
借用有限元法中的形函数概念,局部近似函数为各个形函数与相应系数的乘积之和。
对于式(3)和式(4),不存在Li,仅考虑φiui= φi 以及 =φi(x') ( )。
2.2 形函数的偏导数
关于第i个覆盖中的第mi级数项,导数子矩阵 为
$\boldsymbol{B}_{m_{i}}=\boldsymbol{A}\left(\boldsymbol{T}_{m_{i}}\right)=\boldsymbol{A}\left(\varphi_{i} \boldsymbol{L}_{i} \boldsymbol{N}_{m_{i}}\right)$。
式中A为微分算子矩阵,对于偏微分方程式(2),考虑应变子矩阵的微分算子
$\begin{align}\boldsymbol{A}=\left[\begin{array}{ccc}\frac{\partial}{\partial x} & 0 & 0 \\0 & \frac{\partial}{\partial y} & 0 \\0 & 0 & \frac{\partial}{\partial z} \\\frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 \\0 & \frac{\partial}{\partial z} & \frac{\partial}{\partial y} \\\frac{\partial}{\partial z} & 0 & \frac{\partial}{\partial x}\end{array}\right],(\text { npde }=3) ; \boldsymbol{A}=\left[\begin{array}{cc}\frac{\partial}{\partial x} & 0 \\0 & \frac{\partial}{\partial y} \\\frac{\partial}{\partial y} & \frac{\partial}{\partial x}\end{array}\right],(n p d e=2) ;\\\boldsymbol{A}=\left[\frac{\partial}{\partial x}\right],(\text { npde }=1)\end{align}$
对于偏微分方程式(3)和式(4)
可见,需要求出 ( )、 ( )和 ( ),然后放入矩阵 中的相应位置。令 = 。本文约定,用 ( )表示 ⊗( ),⊗为克罗内克积,即算子向量 的每一项作用于 矩阵的每一个元素(下文的 类似)。首先
。
对于 (Li ),Li分为常量矩阵和变量矩阵2种, 分为对角矩阵和非对角矩阵2种,组合成以下多种情况:
(1)最常见的简单情况,整体直角坐标定义的多项式级数,即覆盖坐标 =x。此时Li=I, =I ,则
。
。
。
。
以下具体讨论 的计算。以 为例,令 = ,根据级数表达式求得 之后,再由链式法则 = 得到 ,即
。
因此首先要求出级数坐标 关于整体坐标x的偏导数 。
对于情况(1), =x,仅需考虑 与 转换关系式,即式(14),
。
对于情况(2)和(3),覆盖坐标系 =( , , )为固定的坐标系,由于
。
。
对于情况(4),覆盖坐标系为定义在曲线或曲面边界条带上的随动坐标系,以曲面为例,边界覆盖的覆盖坐标 中的( , )=(x',y')=( , ),式(32)改写成
。
对于边界覆盖,通过式(32)或式(34)分别求得平面或曲面下的 ,则
。
同理,φi、Li各为坐标x'或 表述,可能是固定的直角坐标或随动的曲线曲面坐标(均在上文讨论过),相应地求得$\frac{\partial \varphi_{i}}{\partial \boldsymbol{x}^{\prime}} 、\frac{\partial \boldsymbol{L}_{i}}{\partial \tilde{\boldsymbol{x}}}$并按式(32)或式(34)找到x'或 与x关系后,即可得到 、 。
对于偏微分方程式(3)和式(4),式(25)不涉及Li,只需计算 和 ,但这2项计算中的坐标转换仍可能涉及Li。
2.3 小结
本节内容看似复杂,实际是为了应对各种情况(新方法的多种特色功能)进行的分门别类。只需掌握以下几条原则:
(1)在各覆盖中定义覆盖坐标系 ,最常见的是整体直角坐标系即 =x,也可能根据材料的方向、裂尖处的裂纹走向来定义,边界覆盖的覆盖坐标 等同于边界条带的积分坐标x'。级数坐标 与覆盖坐标 之间一般是归一化的关系式(14),裂尖所在覆盖是极(柱)坐标与 的关系。
(2)形函数子矩阵 =φi(x')Li( ) ( ),其中,φi(x')、Li( )、 ( )分别由条带积分坐标x'、覆盖坐标 和级数坐标 确定。若 由曲线或曲面描述,则Li( )为变量矩阵,否则为常量。式(25)中 的偏导数需要考虑Li(常量或变量)、 (对角或非对角)情况的各种组合。
(3)最终都要在整体坐标x下求解,因此需要根据不同形体(分区或条带)所采用的坐标情况(可能由整体坐标x或参数坐标 定义,见本系列文章中的第二篇文章)进行坐标转换,包括采用链式法则的微分运算。
3 单元矩阵
独立覆盖流形法的单元包括独立覆盖单元和条带单元(含边界条带)。质量子矩阵 为
。
式中:Ωe为单元的积分域。其中,对于独立覆盖,i=1,j=1;对于条带单元,两侧可连接不同类型的独立覆盖(如一般的独立覆盖与边界覆盖[4]),i=1,2,j=1,2;mi和mj依次对i覆盖和j覆盖的所有级数项循环,下同。
关于坐标的2阶微分项形成的刚度子矩阵 为
。
体荷载子向量 为
。
对于第二类边界条件式(6),面荷载或线荷载的子向量 为
。
式中:Γe为单元表面上荷载作用区的积分域,q= 。若面荷载按面的法向和切向分量定义为σ= ,即按面的积分坐标x'=(x',y',z')定义(具体表述为t-s-n),则q=Lσ,其中,L= ,cosαt、cosβt、cosγt分别是t轴关于整体坐标x、y、z的方向余弦,其他以此类推,法向n指向体内为正。此外,还可以考虑作用于点的集中力荷载 = f。
对于第三类边界条件式(8),分别形成刚度矩阵和荷载向量:
。
其中,Db= ;u'0= ;
Db,u'0定义于面的积分坐标系x'下。
以上各项单元子矩阵形成单元质量矩阵Me、刚度矩阵Ke、荷载向量fe以及待求未知量向量de(来自于 )。对于式(3)和式(4)的波动方程和传导方程,类似的M、K、F也分别借用质量矩阵、刚度矩阵和荷载向量的名词,只是将以上公式中的 替换为 ,f、q和u'0分别替换为f、q0和u0,D替换为a2,Db替换为k,且不涉及L矩阵。
4 时间积分
4.1 时间的一阶导数
式(3)含有关于时间的一阶导数。各单元矩阵组集K= Ke、M= Me、F= Fe和u= de后( 表示对单元进行组集),
。
令u、 、F分别为[24]
式中:u0、 、F0分别为前一步的物理量、时间导数和荷载;F1为当前步的荷载;u1、 分别为待求的当前步的物理量及其导数;θ为参数。则
$\begin{aligned}\left(\frac{\boldsymbol{M}}{\theta \Delta t}+\boldsymbol{K}\right) \boldsymbol{u}_{1} & =\left(\frac{\boldsymbol{M}}{\theta \Delta t}+\boldsymbol{K}-\frac{1}{\theta} \boldsymbol{K}\right) \boldsymbol{u}_{0}+ \\\boldsymbol{F}_{1} & +\frac{1-\theta}{\theta} \boldsymbol{F}_{0}。\end{aligned}$
式中Δt为时间步长。当θ≥1/2无条件稳定,θ=1为完全隐式方法。由式(10)的初值条件,逐步计算出每一步的物理场量u1。
4.2 时间的二阶导数
式(2)和式(4)中含有关于时间的二阶导数。以式(2)为例,组集矩阵后生成
。
考虑Rayleigh阻尼[1],令
。
式中α、β为参数。设u1、 分别为[24]
。
式中 和 分别为前一步和本步的二阶导数值,则
$\begin{array}{c}{\left[\left(\alpha+\frac{1}{\theta \Delta t}\right) \boldsymbol{M}+(\beta+\theta \Delta t) \boldsymbol{K}\right] \boldsymbol{u}_{1}=} \\\theta \Delta t \boldsymbol{F}_{1}+(1-\theta) \Delta t \boldsymbol{F}_{0}+\left(\alpha+\frac{1}{\theta \Delta t}\right) \boldsymbol{M} \boldsymbol{u}_{0}+ \\\frac{1}{\theta} \boldsymbol{M} \frac{\partial \boldsymbol{u}_{0}}{\partial t}+[\beta-(1-\theta) \Delta t] \boldsymbol{K} \boldsymbol{u}_{0}\end{array}$。
由式(10)和式(11)的初值条件,逐步计算出每一步的物理场量u1,则一阶导数和二阶导数值分别为:
。
5 结束语
本篇总结了独立覆盖流形法最基本的通用计算公式,将局部近似函数中的形函数表达为单位分解函数、坐标转换矩阵和级数矩阵的乘积,具体讨论了各种情况下的形函数及其求导方法,给出各种矩阵的表达式,以及时间积分方法,适用于弹性力学运动微分方程、传导方程、波动方程的一维至三维稳态和瞬态分析(含三类边界条件),涵盖了高阶级数、任意形状网格、精确几何边界模拟及本质边界条件的准确施加、裂纹尖端附近采用解析级数等独立覆盖流形法的特色功能(精确几何的梁板壳分析将另文讨论)。利用这些公式即可开展独立覆盖流形法通用程序的编写工作,系列文章中的第二篇文章将介绍通用程序设计方法。