加权余量法推导泊松方程有限元解

外太空土豆儿

2月前浏览46

问题描述

考虑泊松方程在域上，边界为：

其中是未知函数，是已知源项函数。假设齐次 Dirichlet 边界条件：

数值求解步骤

1. 近似解构造

取近似解为基函数的线性组合：

其中：

为待定系数
需满足（齐次边界条件）
基函数通常选用分段多项式（有限元法）或全局正交函数（谱方法）

2. 余量定义

将近似解代入原方程得余量（残差）：

（因是近似解，余量不为零）

3. 加权余量原理

强制余量在加权积分意义上为零：

其中是权函数。在 Galerkin 法中取：

4. 弱形式推导（关键步骤）

对二阶导数项应用分部积分（格林恒等式）：

代入原方程得弱形式：

5. 离散代数系统

将代入弱形式：

定义：

刚度矩阵：
载荷向量：

得到矩阵方程：

其中为待求系数向量。

6. 数值求解流程

计算矩阵元素：
计算载荷项：
求解线性系统：
重构解：

物理与数学意义

步骤	数学操作	物理意义	数值优势
近似解	基函数展开	将连续问题投影到有限维空间	降低问题维度
弱形式	分部积分	将通量守恒转化为能量形式	降低导数阶次，允许连续基函数
矩阵组装	计算双线性积分	构建系统的能量平衡关系	生成稀疏对称正定矩阵

重要说明

基函数选择：

有限元法：局部支撑的分段多项式
谱方法：全局光滑的正交函数
均需满足

积分计算：

实际计算采用数值积分（如高斯积分）
在单元上分段积分后组装全局矩阵

边界条件处理：

非齐次 Dirichlet 条件需特殊处理
Neumann 边界条件会出现在边界积分项中

误差来源：

截断误差：有限维近似
离散误差：数值积分近似
代数求解误差

补充：格林恒等式推导

格林恒等式是矢量分析中的核心定理，在偏微分方程理论中具有基础地位。下面逐步推导格林第一恒等式。

预备知识

散度定理（高斯定理）：

其中：

是光滑矢量场
是边界的单位外法向量

乘积法则：

步骤1：构造矢量场

考虑两个标量函数和（要求。构造矢量场：

步骤2：应用散度定理

对应用散度定理：

步骤3：展开被积函数

左边应用乘积法则：

代入得：

步骤4：处理边界项

右边边界项可改写为：

其中是在边界上的法向导数。

步骤5：得到最终形式

整理得格林第一恒等式：

分量形式（三维直角坐标系）

特殊情形

当时：

（拉普拉斯算子的积分等于法向导数的通量）

当时：

物理意义

能量解释：在泊松方程中表示系统的总能量（左）与边界能量通量（右）的平衡
守恒律：描述标量场 (v) 通过算子 (u) 加权后的守恒关系
降阶作用：将二阶导数（拉普拉斯项）转化为一阶导数（梯度项）+ 边界项

应用场景

有限元方法：推导泊松方程的弱形式

电磁理论：证明电势解的唯一性
流体力学：推导不可压缩流的能量守恒方程

来源：有限元先生

ABAQUS输出节点荷载列阵

概述自己手撸有限元代码的时候，需要将自己计算的一些数据与商软对比，如全局/单元刚矩阵、质量矩阵和荷载列阵等等，我在早期的帖子中详细说过刚度矩阵、质量矩阵和阻尼矩阵的导出方法，有兴趣的可以阅读早期的帖子。本次帖子主要讲解在ABAQUS平台如何修改INP文件导出节点的荷载列阵。模型介绍建立尺寸为2x2x10的长方体，设置单元尺寸为0.5x0.5x0.5，划分后的网格为边界条件为一端固定。另一端施加切向面力，面力大小为1，边界条件示意图为施加动荷载，幅值曲线表达式为，计算总时长为10，增量步长为0.01.总计增量步数目为1000，下图为荷载的幅值曲线导出荷载列阵 ABAQUS无法通过GUI界面输出节点荷载列阵，必须修改INP文件中的关键字。废话不多说，下面是分析步的所有关键字。*Step, name=Step-1, nlgeom=NO, inc=1000*Dynamic,direct0.01,10.,*Dsload, amplitude=Amp-1Surf-1, TRSHR, 1., 1., 0., 0.*File Format, ASCII*Element Matrix Output, Elset=Part-1-1.Set-1, File Name=EMass, Output File=User Defined, mass=yes*Element Matrix Output, Elset=Part-1-1.Set-1, File Name=EStiffness, Output File=User Defined, stiffness=yes*Element Matrix Output, Elset=Part-1-1.Set-1, File Name=dload, Output File=User Defined, dload=yes*Restart, write, frequency=0*Output, field, variable=PRESELECT, frequency=1*Output, history, variable=PRESELECT*End Step 其中，最重要的有三行*Element Matrix Output, Elset=Part-1-1.Set-1, File Name=EMass, Output File=User Defined, mass=yes*Element Matrix Output, Elset=Part-1-1.Set-1, File Name=EStiffness, Output File=User Defined, stiffness=yes*Element Matrix Output, Elset=Part-1-1.Set-1, File Name=dload, Output File=User Defined, dload=yes 前两行输出了刚度矩阵和质量矩阵，最后一行输出单元节点的荷载列阵。计算结束之后，INP文件所在的文件夹会出现“DLOAD.mtx”文件，这里面就是节点荷载信息，下面截取mtx文件的部分内容图中给出了3组数据，主要涉及到编号为1、2和3的单元，注意并不是这3个单元的所有节点都有节点荷载，下面以第一个单元为例讲解**** ELEMENT NUMBER 1 STEP NUMBER 1 INCREMENT NUMBER 1** ELEMENT TYPE C3D8R *USER ELEMENT, NODES= 8, LINEAR** ELEMENT NODES** 106, 107, 112, 111, 1, 2, 7, 6 1, 2, 3** LOAD VECTOR FROM ELEMENT DISTRIBUTED LOADS. LOAD CASE 1***CLOAD** 1, 1, 6.24990E-04** 1, 2, -6.41966E-37** 2, 1, 6.24990E-04** 2, 2, -1.72014E-37** 5, 1, 6.24990E-04** 5, 2, -2.39585E-36** 6, 1, 6.24990E-04** 6, 2, -6.41966E-37 上面的信息包括了当前分析步和增量步位置，以及当前单元的所有节点编号，关键字下面就是具体的节点荷载列阵，一共有3列数据，第1列为施加荷载的节点编号，第2列是施加荷载的自由度，第3列是荷载数值。注意，在这个帖子的算例中，第二个自由度没有施加动荷载，但是这里依然有极小量数值，我猜测可能是因为软件算法为了数值平衡设置的一些数值，我们不用管。并不是所有的单元都会有节点荷载，但是ABAQUS在没有节点荷载的位置依然输出单元的有些信息。如**** ELEMENT NUMBER 271 STEP NUMBER 1 INCREMENT NUMBER 354** ELEMENT TYPE C3D8R *USER ELEMENT, NODES= 8, LINEAR** ELEMENT NODES** 458, 459, 464, 463, 353, 354, 359, 358 1, 2, 3** LOAD VECTOR FROM ELEMENT DISTRIBUTED LOADS. LOAD CASE 1***CLOAD 可以发现，这里只输出了单元的一些基本信息，不包括荷载列阵。我猜测这是内部代码没有对节点类型进行分类，直接在一个循环中索引了所有的节点，当检测到荷载列阵为空的时候就直接进行下一个单元。总结帖子讲解了在ABAQUS平台输出节点荷载列阵的方法。因为ABAQUS的GUI界面没有提供键鼠操作，因此需要修改INP文件，上面给出了具体的修改方法。来源：有限元先生