飞机主要采用由骨架和蒙皮组成的薄壁结构, 结构复杂、尺度较大, 特别是大型运输机等研制需求, 使航空领域面临着日益复杂化、精细化、超大规模化等严峻挑战。计算机技术的不断革新、集群技术的快速发展, 使得快速并行计算、精简存储等高效数值算法备受瞩目。有限元撕裂合并(finite element tearing and interconnecting, FETI)^[1]正是一种并行求解大规模结构的数值计算技术, 因具有良好的稳定性、收敛性等特点得到不断发展, 先后出现了A-FETI^[2]、FETI-DP^[3]、L-FETI^[4]等形式。

然而, 实际应用中, 对不同子域并行网格划分时通常会采用不同的网格生成标准或工具, 且部分独立子域也会有网格加密的需要, 往往就导致部分相邻子域公共界面上的网格不匹配。实现非匹配界面“粘接”的方法有很多^[5], 如多点约束(multi-point constraint, MPC)、耦合有限单元(coupling finite elements, CFEs)^[6]等, 插值法就是一类重要且实用的数据传递技术。一些常用插值方法有最邻近插值(nearest neighbour interpolation)^[7]、投影法(projection method)^[8]及无限平板样条、薄板样条、非均匀有理B样条^[9]等样条插值法。文献[10]对一些插值方法的精度和效率进行了研究, 其中, 径向基函数(radial basis function, RBF)插值方法得到普遍关注^[11-14]。Smith等^[15]对几种常用的径向基函数进行了比较, Beckert和Wendland^[16]指出具有紧支撑的径向基函数具有更优异的特性。

为实现非匹配多域有限元问题的快速并行求解, 本文将L-FETI方法中相邻子域间界面节点力细划为子域边界节点力和分区框架节点力; 通过建立非匹配分区框架上的载荷平衡条件和位移协调条件, 推导了非匹配多域问题中框架方程的组织形式; 同时, 基于RBF插值技术, 在局部坐标系下实现了非匹配界面上数据的有效传递。最后, 通过简单数值算例验证了本文算法的有效性。

1 非匹配多域FETI算法

L-FETI方法在子域界面引入三重变量:界面节点位移、界面节点力、分区框架上耦合节点位移, 并基于3个子域条件和2个框架条件, 得到一组近似解耦的框架方程, 使得各个子域的计算相对于经典FETI算法更加独立。然而, 由框架载荷平衡条件(1)式和位移协调条件(2)式可看出, L-FETI方法仅适用于匹配网格多域有限元模型的并行求解

(1)

(2)

式中, N_s为子域总个数, u_b⁽ⁱ⁾和λ_b⁽ⁱ⁾分别为子域Ω_i边界上节点位移和节点力, u⁽ⁱ⁾为子域Ω_i中全部节点位移, u_f为总框架上节点位移。B⁽ⁱ⁾和L⁽ⁱ⁾为由0, 1构成的装配算子, 表示待提取自由度必然是被提取自由度的子集, 若待提取自由度不在被提取自由度中, 则该方法将束手无策。针对非匹配多域有限元问题的并行求解, 本文将介绍一种改进的FETI方法。

将模型划分成N_s个子域, 并在模型切割处引入虚拟框架Γ。设子域Ω_i和Ω_j相邻, 记两子域之间的局部框架为Γ_ij, 与子域Ω_i相邻的局部框架界面为Γ_i。一般来说, Γ_ij上节点分布与相邻子域边界Γ_i上节点分布无关, 但为了提高存储效率及对计算精度的考量, 当Γ_ij为非匹配界面时, 可令Γ_ij上节点分布与细网格边界节点一致。在Γ_ij处引入四重变量:子域边界节点位移u_b⁽ⁱ⁾和节点力λ_b⁽ⁱ⁾, 虚拟框架节点位移u_f⁽ⁱ⁾和节点力λ_f⁽ⁱ⁾, 如图 1所示。

对于非匹配多域模型, 将框架载荷平衡条件(1)式和位移协调条件(2)式做如下修正

(3)

(4)

式中, u_f为总框架Γ上节点位移, u_f⁽ⁱ⁾为与子域Ω_i相关的局部框架Γ_i上节点位移, G⁽ⁱ⁾为局部框架Γ_i节点到子域Ω_i边界节点的位移插值矩阵, 其余与前文描述一致。特别地, 当框架界面节点与子域边界节点匹配时G⁽ⁱ⁾为单位矩阵I⁽ⁱ⁾。

首先, 由子域有限元方程解得各子域位移解, 并联合子域外载条件代入(4)式, 可得框架位移协调方程

(5)

式中, R⁽ⁱ⁾为漂浮子域刚体模态阵, 可由漂浮子域节点坐标计算得到; α⁽ⁱ⁾为漂浮子域刚体位移; F_(s)为子域柔度阵, 对于约束子域来说, 其为子域刚度矩阵的逆矩阵, 对于漂浮子域来说, 其为子域广义柔度阵, 可根据文献[17]计算得到。

其次, 根据力的传递关系, 原框架平衡方程(3)式可改写为

(6)

最后, 将(5)式、(6)式和漂浮子域自平衡条件按子域顺序装配在一起, 可得非匹配多域模型中框架方程组的矩阵表达形式

(7)

式中

本文算法可以方便地嵌入第三方软件中, 充分利用现有有限元分析软件中成熟的单元库来获得各分割子域的刚度矩阵, 进而通过框架上有序节点的位移插值来组装框架方程组。

对于(7)式所示框架方程组, 利用预处理共轭投影梯度(PCPG)等方法进行求解, 可得到各子域边界节点力及各漂浮子域的刚体位移, 再通过子域问题的回代运算即可得到各子域位移解。值得一提的是, 在算法编程中, 操作算子B⁽ⁱ⁾和L⁽ⁱ⁾均可声明为待提取自由度(如子域边界自由度)在总自由度(如子域自由度或框架自由度)中的位置序列所构成的一维数组。这样, 通过压缩0元素, 既能节约存储空间, 也能提高检索计算效率。

2 局部框架数据传递

前文描述了适用于非匹配多域模型的改进FETI算法, 本节将介绍采用RBF插值来实现子域边界与局部框架间的数据传递。

RBF插值以界面能量守恒原理为基础, 即子域Ω_i边界作用于框架上的力λ_f⁽ⁱ⁾在框架任意位移δu_f⁽ⁱ⁾上所做的虚功与框架作用于子域Ω_i边界上的力λ_b⁽ⁱ⁾在子域边界任意位移δu_b⁽ⁱ⁾上所做的虚功恒相等。因为节点位移和节点力存在多个分量, 且分区边界可任意切割, 为简化各边界插值矩阵中自由度的选取, 本文引入切割边界局部坐标系。设子域Ω_i分割边界由多个简单边界(线或面)构成, 每个简单边界对应的局部框架为Γ_i^k。以简单边界端点为原点构建局部坐标系, 则局部坐标系下框架Γ_i^k节点集和对应简单边界节点集的某一自由度(例如u)方向位移插值分别如(8)式和(9)式所示。

(8)

(9)

式中, a为RBF插值中一次多项式系数, 其项数与节点分布相关, 当节点分布在直线、平面或空间曲面中, 项数N_p分别取2, 3或4;α为基函数系数; P和Q是与多项式相关的插值矩阵系数; M和T是与基函数相关的插值矩阵系数, 基函数的选择可参照文献[15]。由RBF原理可知, 此时, 非匹配界面上u自由度方向上节点力插值满足力守恒和矩守恒。

由(8)式解出系数项并代入(9)式可得u向位移插值

(10)

式中, G_uLoc^(ik)为局部坐标系下子域Ω_i第k条边界上的u向插值算子。同理, 可得其余方向的插值算子:G_vLoc^(ik), G_wLoc^(ik), G_θxLoc^(ik), G_θyLoc^(ik), G_θzLoc^(ik)。考虑到有限元方程组中待求解向量一般是逐节点将节点自由度集排列, 因此需引入初等行变换将所有自由度级下的插值重新排列, 再将重排后的节点位移转换到总体坐标系下, 可得

(11)

式中

(12)

G^(ik)为子域Ω_i第k条边界与局部框架Γ_i^k的位移插值算子, C_b^-(ik)和C_f^(ik)为初等行变换矩阵, S为局部坐标系与总体坐标系〈i, j, k〉节点位移之间的转换矩阵。遍历子域Ω_i的所有分割边界并将得到的插值算子G^(ik)按自由度组装, 便可得到子域Ω_i分割边界到相邻框架的插值算子G⁽ⁱ⁾。最后, 根据能量守恒原理, 可得局部框架Γ_i与子域Ω_i边界之间节点力的传递公式为

(13)

3 数值算例

以板弯模型为例, 长200 mm, 宽100 mm, 厚2 mm, 材料弹性模量E=70 Gpa, 泊松比μ=0.33。平板一端固支, 对边等距节点处各施加集中载荷f=10 N, 有限元网格采用四边形壳单元, 并被分割成4个独立子域, 如图 2所示。各子域网格粒度及框架信息如图 3所示, 框架1和3为匹配界面, 框架2和4为非匹配界面。

为验证计算精度, 本文另外分别构建了相应粗粒度和细粒度的2种匹配网格模型, 将3种模型中各自框架上节点集取并集, 并对合并后节点进行统一编号, 如图 4所示, 其中含有‘○’, ‘●’和‘□’分别表示全粗网格、全细网格和非匹配网格中框架节点。从三者的数值结果中提取出各自框架上节点的数值解进行对比, 各模型框架节点y向位移如图 5所示。表 1为非匹配框架界面2和4上各节点y向数值解, 表 2为匹配框架界面1和3上各节点y向数值解。

由图 5可见, 粗网格、细网格和非匹配网格模型中框架各节点y向位移解基本一致, 表明本文提出的在虚拟分区框架上引入独立节点位移和节点力的策略是完全可行的。

在表 1~2中, v_f, v_c和v_n分别表示细粒度、粗粒度和非匹配网格模型中相应节点y向位移。观察发现, 在非匹配框架界面2上, 从多子域交点(12号节点)到靠近模型自由边界节点(8号节点)误差越来越大, 在非匹配框架界面4上也一样, 因此, 分割边界上靠近模型自由边界处的节点误差是重要关注点。

由表 1~2可以看出, 在整个分区框架上最大误差在24号节点处, 为1.975%。利用商业分析软件Abaqus的绑定功能, 采用同样的壳单元对此3种粒度有限元模型分别进行数值仿真求解, 经相同的分析和比对方法得出框架节点y向位移的最大误差也在24号节点, 为3.951%, 表明本文提出的算法有效。

4 结论

以上计算分析表明了本文提出方法的可行性, 现将主要工作归结如下,

1) 在L-FETI区域分裂算法的基础上, 引入独立的框架节点力, 使子域边界和分区框架上各自的节点力和节点位移充分局部化, 修改了框架位移协调条件和载荷平衡条件, 推导出了非匹配多域问题中有限元撕裂合并的并行求解算法。

2) 以RBF插值技术为基础, 在非匹配界面(线或面)上引入局部坐标系, 使任意切割方向而界面形式(线或面)相同的分区边界上的传递矩阵具有统一格式, 避免了RBF插值中为防止插值矩阵奇异而带来的坐标选取难题, 便于编程处理。

3) 在四节点机群环境中完成四区域非匹配板弯模型的并行求解, 数值结果表明本文提出的非匹配多域并行处理技术具有良好的精度。