二维Poisson方程有限元求解

1. 方程

二维Poisson方程的形式为：

$$ \frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2} = f(x, y) $$

这里不给出边界条件的原因是，采用数值方法（尤其是有限元方法）对于模型几何域的不敏感，任何形状都能够完成求解。边界条件的施加可以在控制方程得到离散后进行。

2. Galerkin法进行PDE方程的离散

Poisson方程属于边值问题，即物理场函数仅是空间坐标的函数，与时间无关。为了获得其有限元形式，需要首先构建出PDE对应的等效泛函形式。这里采用加权残余方法，当权函数取形函数时，该方法被称为Galerkin法。首先，有限元法应当先确定出单元的类型，通常为了计算方便选择低阶三角形单元，一方面在单元积分计算时能够直接进行显式表达，另一方面，三角形单元的保形能力较好，对于几何域的逼近效果要好。

2.1 单元类型

单元类型为三节点三角形单元，选定单元类型意味着形函数就会确定，这里本质要做的工作就是，在一个单元内，用这个单元三个节点处的场函数值，利用插值的思想去确定出单元内任意位置处的值。形函数的形式通常选择线性插值，利用单元三个节点的值，能够反向求解出插值系数，从而得到形函数的形式。即

$$ N_i = \frac{1}{2A} (\alpha_i + \beta_i x + \gamma_i y)\\ N_j = \frac{1}{2A} (\alpha_j + \beta_j x + \gamma_j y)\\ N_m = \frac{1}{2A} (\alpha_m + \beta_m x + \gamma_m y) $$

此时，单元内的任意点的场函数值可以表示为

$$ u(x, y) = N_i \cdot u_i + N_j \cdot u_j + N_m \cdot u_m $$

表示为矩阵形式为：

$$ u(x, y) = \begin{bmatrix} N_i & N_j & N_m \end{bmatrix} \cdot \begin{bmatrix} u_i \\ u_j\\ u_m \end{bmatrix}=[N] \cdot \{u_e\} $$

2.2 PDE等效积分形式

利用加权残余方法构建方程的积分形式。这里我们寻求的是数值解，是一个近似解答，如果记数值近似解为$\omega(x,y)$，既然为近似，那么代入原始的控制方程中必定是不满足的，因此，会存在一个残差$R$，为

$$ R = \frac{\partial ^2 \omega}{\partial x ^ 2} + \frac{\partial ^2 \omega}{\partial y^2} - f(x,y) $$

寻求近似解，自然要找到误差或者偏差最小的那一个，因此，为了找到误差最小的，强制使得的残差为0，应当是最好的。注意残差应当是相对于区域而言的，此时针对的是一个单元，所以需要函数内积的概念，加权函数选择形函数相同的，得到

$$ \iint_\Omega [N]^\mathrm{T} (\frac{\partial ^2 \omega}{\partial x^ 2} + \frac{\partial ^2 \omega}{\partial y^2} - f(x,y))\mathrm{d}x \mathrm{d}y $$

这个形式其实就是PDE所对应的等效积分强形式，其中函数仍然要求保证二阶导的连续，对于场函数的连续并没有降低，因此，称之为强形式。此时进行求解仍然具有较大的困难。并且，数学上二阶连续性的要求对于现实中的情况而言过高，降低连续性通常会得到更加真实的解答。因此，寻求弱形式似乎很有必要。此时利用高斯定理，对于二维问题高斯定律则退化为格林公式，本质上就是积分的分部积分。

$$ \iint_{\Omega}[N]^{\mathrm{T}} \frac{\partial ^ 2 \omega}{\partial x^2}\mathrm{d}x\mathrm{d}y= \int[\int[N]^\mathrm{T} \frac{\partial ^ 2 \omega}{\partial x ^ 2}\mathrm{d}x]\mathrm{d}y $$

进行分部积分，有

$$ \int[\int[N]^\mathrm{T} \frac{\partial ^ 2 \omega}{\partial x ^ 2}\mathrm{d}x]\mathrm{d}y= \int [{[N]^\mathrm{T}\frac{\partial \omega}{\partial x}|_{\partial \Omega}} - \int\frac{\partial [N]^{\mathrm{T}}}{\partial x}\frac{\partial \omega}{\partial x}\mathrm{d}x]\mathrm{d}y $$

注意到，在离散方程的时候，此时没有考虑边界条件，因此边界均是无约束的情况，此时分部积分的第一项的值应当为0，由此得到

$$ \iint_{\Omega}[N]^{\mathrm{T}} \frac{\partial ^ 2 \omega}{\partial x^2}\mathrm{d}x\mathrm{d}y= -\iint_{\Omega}\frac{\partial [N]^{\mathrm{T}}}{\partial x} \frac{\partial \omega}{\partial x}\mathrm{d}x\mathrm{d}y $$

关于$y$的偏导项进行同样的处理后，得到等效积分弱形式。

$$ \iint_{\Omega}[\frac{\partial [N]^{\mathrm{T}}}{\partial x}\frac{\partial \omega}{\partial x}+\frac{\partial [N]^{\mathrm{T}}}{\partial y}\frac{\partial \omega}{\partial y} - [N]^{\mathrm{T}}f(x,y)]\mathrm{d}x\mathrm{d}y=0 $$

上式即为等效积分弱形式，其中场函数要求一阶连续，连续性得到降低。为了得到有限元形式，在弱形式的基础上，利用形函数、节点处函数值表示的场函数代入后，得到

$$ \iint_{\Omega}\frac{\partial [N]^{\mathrm{T}}}{\partial x}\frac{\partial [N]}{\partial x} + \frac{\partial [N]^{\mathrm{T}}}{\partial y}\frac{\partial [N]}{\partial y}\mathrm{d}x\mathrm{d}y \cdot \{u_e\} = \iint_{\Omega}[N]^{\mathrm{T}}f(x,y)\mathrm{d}x\mathrm{d}y $$

此时，得到了离散后的有限元形式，左端将形函数矩阵以及形函数表达式代入后，可以发现，形函数的一阶导是域空间坐标无关的，因此，这个积分是直接可以进行计算的，积分就是这个单元的面积$A$。而形函数矩阵相乘后得到左端为

$$ [B] = \frac{1}{4A} \begin{bmatrix} \beta_i^2 + \gamma_i^2 & \beta_i \beta_j + \gamma_i \gamma_j & \beta_i \beta_m + \gamma_i \gamma_m\\ \beta_j \beta_i + \gamma_j \gamma_i & \beta_j^2 + \gamma_j^2 & \beta_j \beta_m + \gamma_j \gamma_m\\ \beta_m \beta_i + \gamma_j \gamma_i & \beta_m \beta_j + \gamma_m \gamma_j & \beta_m^2 + \gamma_m^2 \end{bmatrix} $$

现在只剩下右端的源项的积分了，因为形函数与$x,y$有关，并且源项$f(x,y)$也是$x,y$的函数，因此，这个积分不能直接进行计算，它是在三角形区域上的一个积分，采用数值积分应当是最好的形式。即需要用高斯积分进行计算。由于整个推导采用的都是基于全局坐标进行的，并没有引入等参描述，采用高斯积分时还需要利用等参描述，太过于麻烦。这样的情况下一般会采用近似积分来代替，即利用单元中心处的积分值乘以单元面积进行近似计算

$$ x_c = \frac{x_i + x_j + x_m}{3}\\ y_c = \frac{y_i + y_j + y_m}{3} $$

有：

$$ \iint_{\Omega}[N]^{\mathrm{T}}f(x,y)\mathrm{d}x\mathrm{d}y=[N(x_c, y_c)]^{\mathrm{T}} \cdot f(x_c, y_c) \cdot A $$

从而得到右侧的单元荷载向量，长度为3，记作$F$，则单元的线性方程组为

$$ [B]\{\omega_e\} = \{F\} $$

在获得单元刚度矩阵以及荷载向量后，通过直接集成的方法得到整个模型的刚度矩阵和向量，然后利用乘大数方法或者置1法进行边界条件的施加，然后通过线性方程组的求解，获得最终的解答。

本站文章如非特别说明，均为原创。未经本人同意，请勿转载！转载请务必注明出处！