X.4 二维平面应力
前言
嗯！
背景
目前为止，我们已经学习了一维梁的应力。
接下来，我们考虑一个二维的膜，它在遭受z轴方向作用时产生的应力和应变。
架子鼓的鼓膜就是一个很好的参考。
控制方程
考虑对如图所示情况的控制方程，又名2D泊松方程（Poisson Equation)：
\[\begin{equation}T\left(\frac{\partial^2w}{dx^2}+\frac{\partial^2w}{dy^2}\right)+f(x,y)=0
\end{equation}
\]
我们定义：w为膜的挠度，x、y为膜的坐标，T是膜每单位长度受到的张力，f表示膜每单位面积受到的横向力。
该控制方程的边界条件类似1D，为：在边界\(\Gamma\)上，
- 固定端（dirichlet boundary condition）：\(w=0\)，表示固定膜的边界；
- 自由端（neumann boundary condition）：\(\frac{dw}{dn}=0\)，表示膜在边界上的法向斜率为0。
若拓展到三维，控制方程变成：
\[\begin{equation}\left(\frac{\partial^2w}{dx^2}+\frac{\partial^2w}{dy^2}+\frac{\partial^2w}{dz^2}\right)=-\frac{1}{T}f(x,y,z)=F(x,y,z)
\end{equation}
\]
下面逐步分析2D泊松方程的推导。其核心思想为：通过考虑微小的单元，在垂直方向上的力平衡，来构造方程。
-
对于位移\(w(x+dx,y)\)（x方向的位移增加量），使用一阶泰勒展开，结果为
\(w(x,y)+\frac{\partial w(x,y)}{\partial x}dx\)。
这表示点\((x+dx,y)\)处的位移，可以近似等于\((x,y)\)点的位移，加上一个与dx成正比的增量，该增量由x方向的偏导数\(\frac{\partial w}{\partial x}\)决定；
对位移\(w(x,y+dy)\)（即y方向的位移增加量），一阶泰勒展开的结果为
\(w(x,y)+\frac{\partial w(x,y)}{\partial y}dy\)。
-
对上述泰勒展开的结果分别求一阶导：
- \(\frac{\partial w(x+dx,y)}{\partial x}=\frac{\partial w(x,y)}{\partial x}+\frac{\partial^2w(x,y)}{\partial x^2}dx\)
- \(\frac{\partial w(x,y+dy)}{\partial y}=\frac{\partial w(x,y)}{\partial y}+\frac{\partial^2w(x,y)}{\partial y^2}dy\)
- 这里需要注意，膜的张力产生的垂直分力，和膜的斜率有关。对于dxdy微小单元，如果左右两边的位移不同，膜会发生倾斜，斜率为\(\frac{\partial w}{\partial x}\)，而x方向斜率和y方向的斜率分别由上述两个方程表示。
膜的张力T将始终沿着膜的切线方向，当膜发生倾斜时，张力T将产生一个垂直于原平面的分力，该分力的大小会和张力T本身的大小，和膜的斜率有关，因为斜率意味着有角度存在。
（这里可以参考一维梁控制方程的推导，本节不再重复）
-
构建力平衡方程：
\[\begin{equation}
\begin{aligned}&-T\frac{\partial w}{\partial x}dy
+T(\frac{\partial w}{\partial x}+\frac{\partial^2w}{\partial x^2}dx)dy
-T\frac{\partial w}{\partial y}dx
\\&+T(\frac{\partial w}{\partial y}
+\frac{\partial^2w}{\partial y^2}dy)dx+f(x,y)dxdy
=0
\end{aligned}
\end{equation}
\]
这个方程代表了微小单元dxdy在垂直方向的力平衡，其中：
- \(-T\frac{\partial w}{\partial x}dy\)表示该单元左侧边在x方向上张力的垂直分量。
我们假设挠度很小，斜率很小，因而\(\sin\theta\approx\tan\theta\approx\theta\)，其中\(tan \theta = \frac{\partial w}{\partial x}\)就是斜率，\(\theta\)为膜的倾角。
这里负号表示方向向下，原式应为\(-Tsin(\theta)\cdot dy\)，这里乘以dy是因为T是膜每单位长度受到的张力，然后代入上述结论可得。
- \(T(\frac{\partial w}{\partial x}+\frac{\partial^2w}{\partial x^2}dx)dy\)表示右侧边在x方向上张力的垂直分量，一样是
\(T\times 斜率\times 单元边长\)。
- \(f(x,y)dxdy\)表示作用在微元体上，每单位体积的外力大小f，乘以微元体面积。
- 根据牛顿第二定律，微元体在垂直方向上受到的合力为零，把所有力相加得式3。
-
对上述方程化简，并除以dxdy，可得控制方程。
热的传递：流动的冷却
如上图所示，考虑一个管道流动的一小段体积，已知管壁温度\(T_w\)，对上述系统建立物理模型。
为简化问题，做出如下假设：
- 稳态流动，\(\frac{\partial T}{\partial t}=0\)；
- 流体的速度分布（velocity profile）是平推流（plug flow）或理想流动，\(U\neq f(r,x)\)：
表示流体在管道的任何横截面上，都以相同的速度流动。
这忽略了黏性作用导致的边界层效应，没有速度梯度。
- 流体的物理性质（密度、比热容、导热系数等）不变；
- 管道壁面的温度保持恒定；
- 在任何一个横截面上，温度都是均匀的，\(T\neq f(r)\)；
- 沿管道轴向的热传导，和热对流相比，可以忽略不计。
构建方程：
\[热量积聚速度=产生热量速度+热流入的速度-热流出速度
\]
下面一项项来看。
- 热量积聚速度：由于稳态流动，\(A\Delta x \rho C_P \frac{d\bar{T}}{dt}=0\)
- 热产生的速度：\(\dot{q}=0\)
- 热流入的速度：\(AU_0\rho C_PT(x)\)
- 热流出的速度：\(AU_0\rho C_PT(x+\Delta x)+2\pi R\Delta xh(T-T_w)\)，，其中第一项表示随着流体流动而流出的热，第二项表示通过管壁流出的热。
据此，我们可以构建方程如下：
\[\begin{equation}AU_0\rho C_PT(x)-AU_0\rho C_PT(x+\Delta x)-2\pi R\Delta xh(T-T_w)=0
\end{equation}
\]
把上式除以\(AU_0\rho C_P \Delta x\)，并求当\(\Delta x\)趋近于0时的极限：
\[\begin{equation}\lim_{\Delta x\to0}\left\{\frac{\left(T(x+\Delta x)-T(x)\right)}{\Delta x}+\frac{2\pi Rh}{AU_0\rho C_P}(T-T_w)\right\}=0
\end{equation}
\]
令\(\frac{2\pi Rh}{AU_0\rho C_P}=\lambda\)，简化可得最后的公式：
\[\begin{equation}\frac{dT(x)}{dx}+\lambda(T-T_w)=0
\end{equation}
\]
移项并积分，得：
\[\begin{equation}\int\frac{dT(x)}{T-T_w}=\int-\lambda dx
\end{equation}
\]
两边积分得：
\[\begin{equation}\ln(T-T_w)=-\lambda x+C
\end{equation}
\]
代入初始条件：当\(x=0\)时，\(T(0)=T_0\)：
\[\begin{equation}\begin{array}{c}\ln(T_0-T_w)=-\lambda(0)+C\end{array},~C=\ln(T_0-T_w)=\ln(K)
\end{equation}
\]
回代：
\[\begin{equation}\begin{array}{l}\ln(T-T_w)=-\lambda x+\ln(T_0-T_w)\\\end{array}
\end{equation}
\]
把带ln的移到同一边，然后两边同时取以e为底的指数：
\[\begin{equation}\frac{T-T_w}{T_0-T_w}=e^{-\lambda x}=\Theta
\end{equation}
\]
该函数图像如下图所示：
上述推导是基于流体流动速度u在管道各处都是恒定的。
如果我们引入一个“抛物线”形的速度分布，如下图：
该速度分布可用如下公式描述：
\[\begin{equation}u=2u_0(1-(\frac{r}{R})^2)
\end{equation}
\]
考虑上述速度分布后，新的能量守恒方程变成：
\[\begin{equation}k\frac{\partial^2T}{\partial x^2}+k\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial T}{\partial r}\right)=2u_0\rho C_p\left(1-\left(\frac{r}{R}\right)^2\right)\frac{\partial T}{\partial x}
\end{equation}
\]
接下来，基于上图中的灰色微元（内半径为r，厚度为\(\Delta r\)，长度为\(\Delta x\)），详细解释上述能量守恒方程的推导过程。
考虑如下三种能量项：
-
轴向热传导：
-
流入x侧面的热量：
\[\begin{equation}q_x=-kA_x\frac{\partial T}{\partial x}|_x=-k(2\pi r\Delta r)\frac{\partial T}{\partial x}|_x
\end{equation}
\]
其中，负号表示热量流入为正，k是热导率，单位W/(mK)，\(A_x\)是x侧面的面积，\(\frac{\partial T}{\partial x}|_x\)表示在x方向处，温度沿x方向的梯度。
-
流出\(x+\Delta x\)侧面的热量：
\[\begin{equation}q_{x+\Delta x}=-kA_{x+\Delta x}\frac{\partial T}{\partial x}|_{x+\Delta x}=-k(2\pi r\Delta r)\frac{\partial T}{\partial x}|_{x+\Delta x}
\end{equation}
\]
-
\(q_x-q_{x+\Delta x}\)可得轴向净热传导。
-
径向热传导：
-
流入内表面的热量：
\[\begin{equation}q_r=-kA_r\frac{\partial T}{\partial r}|_r=-k(2\pi r\Delta x)\frac{\partial T}{\partial r}|_r
\end{equation}
\]
-
流出外表面的热量：
\[\begin{equation}q_{r+\Delta r}=-kA_{r+\Delta r}\frac{\partial T}{\partial r}|_{r+\Delta r}=-k(2\pi(r+\Delta r)\Delta x)\frac{\partial T}{\partial r}|_{r+\Delta r}
\end{equation}
\]
-
\(q_r-q_{r+\Delta r}\)可得净径向热传导。
-
对流：假设速度u只和r有关，而和x无关。
-
流入的热量：
\[\begin{equation}q_{conv}=\dot{m}C_{p}T|_{x}=(\rho A_{x}u)C_{p}T|_{x}=\rho(2\pi r\Delta r)uC_{p}T|_{x}
\end{equation}
\]
其中，\(\dot{m}=\rho A_x u\)是质量流量，\(C_[\)是比热容，\(T|_{x}\)表示x位置处的温度。
-
流出的热量：
\[\begin{equation}q_{conv+\Delta x}=\dot{m}C_pT|_{x+\Delta x}=(\rho A_{x+\Delta x}u)C_pT|_{x+\Delta x}=\rho(2\pi r\Delta r)uC_pT|_{x+\Delta x}
\end{equation}
\]
-
\(q_{conv}-q_{conv+\Delta x}\)可得对流带走的热量。
根据能量守恒，轴向轴向净热传导 + 径向净热传导 - 对流带走的热量 = 0：
\[\begin{equation}\begin{aligned}
&k(2\pi r\Delta r)[\frac{\partial T}{\partial x}|_{x+\Delta x}-\frac{\partial T}{\partial x}|_x]
-k(2\pi r\Delta x)\frac{\partial T}{\partial r}|_r
\\&+k(2\pi(r+\Delta r)\Delta x)\frac{\partial T}{\partial r}|_{r+\Delta r}
-\rho(2\pi r\Delta r)uC_p[T|_{x+\Delta x}-T|_x]=0
\end{aligned}
\end{equation}
\]
具体化简步骤略，主要是把上式除以\(2\pi \Delta x \Delta r\)，然后令\(\Delta x\)和\(\Delta r\)趋近于0，得到几个能求极限的分式的极限，然后回代，再把u的表达式带进去，除以r，最后改写一下。
回顾一下对物理系统进行数学建模建模的重点：
- 大多数物理系统的建模都会产生一阶或二阶偏微分方程。
- 建模涉及通过简化物理系统来开发数学解释。
这会产生具有边界/初始条件和明确参数值的合理确定性数学方程。
模型将根据选定的边界条件/初始条件和参数值计算任何给定时刻的变量。
- 相比之下，模拟则伴随着时间演化的混沌不确定性。
对微观子系统相互作用导致宏观系统行为的演化研究可以大致描述为模拟。
有限元方法
通用步骤和公式
Finite Element Method，这里就是把一个二维平面离散为一个个网格，如下图所示的小三角形网格。
有限元法的基本步骤为：
-
从强形式（原始数学模型、控制方程）开始构建公式
-
根据强形式推导出弱形式（乘以了测试函数，二阶导变一阶导）
-
单元近似
-
选择权重函数
-
构建公式：
\[\begin{equation}[\mathbb{K}]\{\mathrm{u}\}{=}\{\mathbb{F}\}
\end{equation}
\]
-
据公式和各类边界条件求解未知数。
考虑一个热传导问题，其强形式可能为：
- 控制方程：在区域\(\Omega\)内，有\(\frac{d}{dx}(kA\frac{dT}{dx})+Q=0\)；
- 边界条件：在x=0处，\(-\frac{d}{dx}kT=h\)；在\(x=L\)处，\(T=g\)。
强形式要求，场变量T必须连续，其导数也必须连续。
通过将上述控制方程乘以一个权重函数v，然后在整个区域\(\Omega\)上积分，可得：
\[\begin{equation}\int_\Omega v[\frac{d}{dx}(kA\frac{dT}{dx})+Q]dx=0
\end{equation}
\]
使用分部积分，易得：
\[\begin{equation}-\int_\Omega\frac{dv}{dx}kA\frac{dT}{dx}dx+[vAh]_{x=0}+\int_\Omega vQdx=0
\end{equation}
\]
使用弱形式的好处是，它降低了对解的光滑性/连续性要求，只需要解的不连续性是可积分的即可，因此更容易找到近似解。
那么，如何利用弱形式构建近似解呢？
设有如下二阶常微分方程：
\[\begin{equation}\frac{d^2u}{dx^2}+4u=8x^2
\end{equation}
\]
给定\(0<x<\frac{\pi}{4}\)，边界条件\(u(0)=u(\frac{\pi}{4})=0\)。
作为参考，其解析解为：\(u(x)=\cos(2x)+(1+\frac{\pi^2}{8})\sin(2x)+2x^2-1\)
假设U是该方程的解，我们可以假设一个近似解的形式，称为试函数（trial function）：
\[\begin{equation}U(x)=\varphi(x)=a_1\varphi_1+a_2\varphi_2+...+a_n\varphi_n=\sum_{i=0}^na_i\varphi_i
\end{equation}
\]
这里，\(a_i\)是待定的系数，而\(\varphi_i\)是已知的、我们自行定义的基函数，通常选择多项式函数。
这里，选择\(\varphi(x)=\sum_{i=1}^{n}a_{i}x^{i}(\frac{\pi}{4}-x)\)作为基函数。
下一步，使用伽辽金法（Galerkin's Method），选择一组已知的基函数\(\varphi_i (x)\)来构建试函数，并选择相同的基函数构建权函数，即令权函数（weight function）\(v(x)\)和试函数\(U(x)\)相同，使得\(v(x)=U(x)=\varphi(x)\)。
把试函数U代入原微分方程，得到残差：
\[\begin{equation}\varepsilon=|\frac{d^2U}{dx^2}+4U-8x^2|
\end{equation}
\]
这个残差实际上就是把原来函数的\(8x^2\)移动到左边了。如果U是精确解，那么代入U的残差计算结果就应当为0。
把上述残差乘以权重函数v，并在整个区域上积分，得到加权残差积分：
\[\begin{equation}R(x)=\int_0^{\pi/4}v(x)[\frac{d^2U}{dx^2}+4U-8x^2]dx
\end{equation}
\]
分部积分处理后的结果为：（因边界条件为0，忽略边界项）
\[\begin{equation}R(x)=\int_0^{\pi/4}(4vU-\frac{dv}{dx}\frac{dU}{dx}-8vx^2)dx
\end{equation}
\]
我们的目标是找到一组待定系数\(a_i\)，使得上述加权残差几分\(R(x)\)最小。由于我们使用了伽辽金法，令权函数v(x)=试函数U(x)，这意味着\(R(x)\)是\(a_i\)的函数，记作\(R(a_i)\)。
为最小化\(R_{ai}\)，需要求解\(\frac{dR(a_i)}{da_i}=0\)。
把试函数\(U(x)=\sum_{i=0}^na_i\varphi_i\)和权函数\(v(x)=U(x)\)代入上述\(R(x)\)的表达式，并求解\(\frac{dR(a_i)}{da_i}=0\)，得到一个关于待定系数\(a_i\)的线性方程组：
\[\begin{equation}[K]\{u\}=\{f\}
\end{equation}
\]
求解即可。
试函数的注意事项
- 完整性（completeness）：它必须能表示场变量的梯度，必须能表示场变量的任意常数值
- 兼容性（compatibility）：单元边界上的近似值必须是连续的。
收敛=完整性+兼容性
二维控制方程的正式推导
弱形式
设有如图所示的2D域，离散为许多个小的三角形网格。
原二维控制方程为：
\[\begin{equation}T\left(\frac{\partial^2w}{dx^2}+\frac{\partial^2w}{dy^2}\right)+f(x,y)=0
\end{equation}
\]
接下来，选择一个测试函数（或基函数）\(\phi_j (x,y)\)，把上面的强形式方程乘以该测试函数，然后对区域A做积分。
\[\begin{equation}\int_A\left[T\left(\frac{\partial^2w}{\partial x^2}+\frac{\partial^2w}{\partial y^2}\right)+f(x,y)\right]\phi_j(x,y)\mathrm{d}A=0
\end{equation}
\]
该式中，最关键的步骤是，对\(\int_AT\left(\frac{\partial^2w}{\partial x^2}+\frac{\partial^2w}{\partial y^2}\right)\phi_j\mathrm{d}A\)做积分。
为简洁，先令\(\Delta w=\frac{\partial^2w}{\partial x^2}+\frac{\partial^2w}{\partial y^2}\)，故该部分变成\(\int_AT\Delta w\phi_j\mathrm{d}A\)。
在1D问题中，可以直接使用分部积分转化，但二维需要用散度定理。
二维散度定理简写为：
\[\begin{equation}\int_A\nabla\cdot\mathbf{F}\mathrm{d}A=\int_\Gamma\mathbf{F}\cdot\mathbf{n}\mathrm{d}S
\end{equation}
\]
其中，\(\mathbf{n}=(n_x,n_y)\)是边界\(\Gamma\)的法外向。
注意，该公式和PPT里给的公式是等效的：
\[\begin{equation}\int_{A}\left(\frac{\partial G_{x}}{\partial x}+\frac{\partial G_{y}}{\partial y}\right)dA=\oint_{\Gamma}(G_{x}n_{x}+G_{y}n_{y})dS
\end{equation}
\]
这里，对二维散度定理进行详细分析：
通过闭合曲线\(\Gamma\)的流体的总通量（净流出量），等于该曲线所包围的区域A内所有源和汇的总和，即向量场\(\mathbf{F}\)的散度在区域A上的积分。
一项项分析：
-
\(\mathbf{F}\)是一个二维向量场，可写为\(\mathbf{F}(x,y)=\begin{pmatrix}F_x(x,y)\\F_y(x,y)\end{pmatrix}\)，其中\(F_x\)和\(F_y\)是标量函数。
-
\(\nabla\cdot\mathbf{F}\)是向量场\(\mathbf{F}\)的散度，衡量向量场在给定点处“源”或“汇”的强度，其计算公式为：
\[\begin{equation}\nabla\cdot\mathbf{F}=\frac{\partial F_x}{\partial x}+\frac{\partial F_y}{\partial y}
\end{equation}
\]
-
\(\int_A\nabla\cdot\mathbf{F}\mathrm{d}A\)表示在区域A上，对向量场\(\mathbf{F}\)的散度做面积分，即闭合曲线\(\Gamma\)所包围的区域A内，所有源和汇的总和。
-
\(\mathbf{F}\cdot\mathbf{n}\)：向量场\(\mathbf{F}\)和单位法向量\(\mathbf{n}\)的点积，表示\(\mathbf{F}\)在法向量上的分量，或垂直于曲线\(\Gamma\)的分量，可写作\(F_xn_x+F_yn_y\)。
-
\(\int_\Gamma\mathbf{F}\cdot\mathbf{n}\mathrm{d}S\)表示沿着闭合曲线\(\Gamma\)做线积分，即通过该闭合曲线\(\Gamma\)的流体的净流出量。
下一步，我们把\(\Delta w\phi_j\)用散度形式表示：
\[\begin{equation}\Delta w\phi_j=\nabla\cdot(\phi_j\nabla w)-\nabla\phi_j\cdot\nabla w
\end{equation}
\]
具体推导过程如下：
把上述等式代入\(\int_AT\Delta w\phi_j\mathrm{d}A\)，可得：
\[\begin{equation}\int_AT\Delta w\phi_jdA=\int_AT[\nabla\cdot(\phi_j\nabla w)-\nabla\phi_j\cdot\nabla w]dA
\end{equation}
\]
把该积分拆成两部分，然后对第一部分使用二维散度定理，把向量场\(\mathbf{F}\)替换为\(\phi_j \nabla w\)：
\[\begin{equation}T\int_A\nabla\cdot(\phi_j\nabla w)dA=T\oint_\Gamma(\phi_j\nabla w)\cdot\mathbf{n}dS
\end{equation}
\]
把两部分结合起来，得到表达式：
\[\begin{equation}\int_AT\Delta w\phi_j\mathrm{d}A=T\int_\Gamma\left(\nabla w\cdot\mathbf{n}\right)\phi_j\mathrm{d}S-T\int_A\left(\nabla w\cdot\nabla\phi_j\right)\mathrm{d}A
\end{equation}
\]
其中，\(\frac{\partial w}{\partial n}=\nabla w\cdot\mathbf{n}=\frac{\partial w}{\partial x}n_x+\frac{\partial w}{\partial y}n_y\)是w沿着边界\(\mathbf{n}\)（法外向）的方向导数。
PPT当中的公式是把式38拆开：
\[\begin{equation}
\begin{aligned}&\int_{A}T(\frac{\partial^{2}w}{dx^{2}}+\frac{\partial^{2}w}{dy^{2}})\phi_{j}(x,y)dA
\\&=T\oint_{\Gamma}\left(\frac{\partial w}{\partial x}n_{x}+\frac{\partial w}{\partial y}n_{y}\right)\phi_{j}dS
-T\int_{A}\left(\frac{\partial w}{\partial x}\frac{\partial\phi_{j}}{\partial x}+\frac{\partial w}{\partial y}\frac{\partial\phi_{j}}{\partial y}\right)dA
\end{aligned}
\end{equation}
\]
这里出现了一个边界项：\(\int_\Gamma\left(\nabla w\cdot\mathbf{n}\right)\phi_j\mathrm{d}S=\oint_\Gamma\left(\frac{\partial w}{\partial x}n_x+\frac{\partial w}{\partial y}n_y\right)\phi_jdS\)
应用Dirichlet或Neumann边界条件，这里不再讨论，直接令这一项=0。
这样，我们就能得到最终的弱形式表达式：
\[\begin{equation}\int_A\left(T\left(\frac{\partial w}{\partial x}\frac{\partial\phi_j}{\partial x}+\frac{\partial w}{\partial y}\frac{\partial\phi_j}{\partial y}\right)-f(x,y)\phi_j(x,y)\right)dA=0
\end{equation}
\]
域的离散和形函数
三角形网格由其节点函数值决定。
构建位移函数：
\[\begin{equation}w(x,y)=ax+by+c
\end{equation}
\]
代入形函数，根据网格节点函数值线性插值：
\[\begin{equation}w(x,y)=N_1(x,y)w_1+N_2(x,y)w_2+N_3w_3(x,y)
\end{equation}
\]
代入上图所示的各节点坐标值，构建矩阵：
\[\begin{equation}\begin{bmatrix}x_1&y_1&1\\x_2&y_2&1\\x_3&y_3&1\end{bmatrix}\begin{Bmatrix}
a \\b\\
c
\end{Bmatrix}=\begin{Bmatrix}
w_1 \\w_2\\
w_3
\end{Bmatrix}
\end{equation}
\]
下一步，我们把三角形网格的面积表示为顶点坐标：
\[\begin{equation}
\begin{aligned}A^e&=\frac{1}{2}\begin{bmatrix}x_1&y_1&1\\x_2&y_2&1\\x_3&y_3&1\end{bmatrix}
\\&=\frac{1}{2}\left((x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)\right)\neq0
\end{aligned}
\end{equation}
\]
这个面积公式是从把三角形的两条边表示为向量，然后求叉积，再除以2得到的。
接下来，应用克莱姆法则（Cramer's Rule）：
\[\begin{equation}x_i=\frac{\det(\mathbf{A}_i)}{\det(\mathbf{A})}
\end{equation}
\]
其中det(A)是系数矩阵A的行列式，\(A_i\)是通过把A的第i列替换为常数向量b得到的矩阵。
这样，获得三个参数：
\[\begin{equation}a=\frac{1}{2A^e}\begin{vmatrix}w_1&y_1&1\\w_2&y_2&1\\w_3&y_3&1\end{vmatrix},~b=\frac{1}{2A^e}\begin{vmatrix}x_1&w_1&1\\x_2&w_2&1\\x_3&w_3&1\end{vmatrix},~c=\frac{1}{2A^e}\begin{vmatrix}x_1&y_1&w_1\\x_2&y_2&w_2\\x_3&y_3&w_3\end{vmatrix}
\end{equation}
\]
把上述行列式展开（乘出来），回代到场变量的原始方程（式30），然后改写成形函数N的形式：
\[\begin{equation}\begin{gathered}N_{1}(x,y)=\frac{1}{2A^e}[(y_2-y_3)x+(x_3-x_2)y+(x_2y_3-x_3y_2)]\\N_{2}(x,y)=\frac{1}{2A^e}[(y_3-y_1)x+(x_1-x_3)y+(x_3y_1-x_1y_3)]\\N_{3}(x,y)=\frac{1}{2A^{e}}[(y_{1}-y_{2})x+(x_{2}-x_{1})y+(x_{1}y_{2}-x_{2}y_{1})]\end{gathered}
\end{equation}
\]
这个图表示，形函数\(N_i\)在对应节点i处为1，减小直到在非对应节点处为0。
下一步，对位移分别关于x和y求导，可得：
\[\begin{equation}\begin{Bmatrix}\frac{\partial w}{\partial x}\\\frac{\partial w}{\partial y}\end{Bmatrix}=\begin{bmatrix}\frac{\partial N_1}{\partial x}&\frac{\partial N_2}{\partial x}&\frac{\partial N_3}{\partial x}\\\frac{\partial N_1}{\partial y}&\frac{\partial N_2}{\partial y}&\frac{\partial N_3}{\partial y}\end{bmatrix}\begin{Bmatrix}w_1\\w_2\\w_3\end{Bmatrix}=[B]\{\Delta\}^e
\end{equation}
\]
这里，运动方程（kinematic matrix）[B]：
\[\begin{equation}[B]=\begin{bmatrix}\frac{\partial N_1}{\partial x}&\frac{\partial N_2}{\partial x}&\frac{\partial N_3}{\partial x}\\\frac{\partial N_1}{\partial y}&\frac{\partial N_2}{\partial y}&\frac{\partial N_3}{\partial y}\end{bmatrix}=\frac{1}{2A^e}\begin{bmatrix}(y_2-y_3)&(y_3-y_1)&(y_1-y_2)\\(x_3-x_2)&(x_1-x_3)&(x_2-x_1)\end{bmatrix}
\end{equation}
\]
刚度矩阵
回顾弱形式表达式，我们把它变成网格化累加，并采用矩阵形式：
\[\begin{equation}
\begin{aligned}&\int_A\left(T\left(\frac{\partial w}{\partial x}\frac{\partial\phi_j}{\partial x}+\frac{\partial w}{\partial y}\frac{\partial\phi_j}{\partial y}\right)-f(x,y)\phi_j(x,y)\right)dA
\\&=\sum_{e=1}^n\int_{A^e}\left(T\left(\frac{\partial w}{\partial x}\frac{\partial\phi_j}{\partial x}+\frac{\partial w}{\partial y}\frac{\partial\phi_j}{\partial y}\right)-f(x,y)\phi_j(x,y)\right)dA
\\&=\left.\int_A\left(\begin{bmatrix}\frac{\partial\phi_j}{\partial x}\frac{\partial\phi_j}{\partial y}\end{bmatrix}T\begin{Bmatrix}\frac{\partial w}{\partial x}\\\frac{\partial w}{\partial y}\end{Bmatrix} \right)-f(x,y)\phi_j(x,y)\right)dA=0
\end{aligned}
\end{equation}
\]
老规矩，令形函数=测试函数，然后令\(\begin{Bmatrix}\frac{\partial w}{\partial x}\\\frac{\partial w}{\partial y}\end{Bmatrix}=[B]\{\Delta\}^e\)（对插值公式\(w(x,y)=\sum_{i=1}^nN_i(x,y)w_i\)求偏导数），则上式等于：
\[\begin{equation}\int_{A}\left(\left[\frac{\partial\phi_{j}}{\partial x}\quad\frac{\partial\phi_{j}}{\partial y}\right]T[B]\{\Delta\}^{e}-f\phi_{j}\right)dA
\end{equation}
\]
考虑\(j=1,2,3\)三个节点，使用上述式子，叠加一下得：
\[\begin{equation}\int_A\left(\begin{bmatrix}\frac{\partial N_1}{\partial x}&\frac{\partial N_1}{\partial y}\\\frac{\partial N_2}{\partial x}&\frac{\partial N_2}{\partial y}\\\frac{\partial N_3}{\partial x}&\frac{\partial N_3}{\partial y}\end{bmatrix}T[B]\{\Delta\}^e-\begin{Bmatrix}N_1\\N_2\\N_3\end{Bmatrix}f(x,y)\right)dA
\end{equation}
\]
代入一下运动矩阵[B]和形函数矩阵[S]：
\[\begin{equation}\int_{A^e}([B]^TT[B]dA)\{\Delta\}^e-\int_{A^e}([S]^Tf(x,y)dA)
\end{equation}
\]
这样，我们就推导出了：
\[\begin{equation}[K]^e\{\Delta\}^e-\{F\}^e
\end{equation}
\]
刚度矩阵*位移=力，其中：
\[\begin{equation}[K]^{e}=\int_{A^{e}}([B]^{T}T[B]dA)\quad\{F\}^{e}=\int_{A^{e}}[S]^{T}f(x,y)dA
\end{equation}
\]
首先讨论刚度矩阵。
把上式简化一下：
\[\begin{equation}[K]^e=\int_{A^e}([B]^TT[B]dA)=[B]^TT[B]\int_{A^e}(dA)=A^eT[B]^T[B]
\end{equation}
\]
把运动矩阵[B]带进去，易得：
\[\begin{equation}[K]^e=\frac{T}{4A^e}\begin{bmatrix}(x_3-x_2)^2+(y_2-y_3)^2&(x_3-x_2)(x_1-x_3)+(y_2-y_3)(y_3-y_1)&(x_3-x_2)(x_2-x_1)+(y_2-y_3)(y_1-y_2)\\(x_3-x_2)(x_1-x_3)+(y_2-y_3)(y_3-y_1)&(x_1-x_3)^2+(y_3-y_1)^2&(x_1-x_3)(x_2-x_1)+(y_3-y_1)(y_1-y_2)\\(x_3-x_2)(x_2-x_1)+(y_2-y_3)(y_1-y_2)&(x_1-x_3)(x_2-x_1)+(y_3-y_1)(y_1-y_2)&(x_2-x_1)^2+(y_1-y_2)^2\end{bmatrix}
\end{equation}
\]
节点力
\[\{F\}^{e}=\int_{A^{e}}[S]^{T}f(x,y)dA
\]
类似地，用形函数插值顶点值，来表示网格面内任一点所受的力：
\[\begin{equation}f(x,y)\approx[N_1\quad N_2\quad N_3]\begin{pmatrix}f_1\\f_2\\f_3\end{pmatrix}=[S]\{f\}^e
\end{equation}
\]
\[\begin{equation}\{f\}^e=\begin{Bmatrix}f_1\\f_2\\f_3\end{Bmatrix}=\begin{Bmatrix}f_1(x_1,y_1)\\f_2(x_2,y_2)\\f_3(x_3,y_3)\end{Bmatrix}
\end{equation}
\]
回代入原式：
\[\begin{equation}\{F\}^e=\left(\int_{A^e}[S]^T[S]dA\right)\{f\}^e=[H]\{f\}^e
\end{equation}
\]
其中：
\[\begin{equation}
\begin{aligned}[H]&=\int_{A^e}[S]^T[S]dA=\int_{A^e}\begin{Bmatrix}N_1\\N_2\\N_3\end{Bmatrix}[N_1\quad N_2\quad N_3]dA
\\&\approx\int_{A^{e}}\begin{Bmatrix}{N_{1}}^{c}\\{N_{2}}^{c}\\{N_{3}}^{c}\end{Bmatrix}[{N_{1}}^{c}\quad{N_{2}}^{c}\quad{N_{3}}^{c}]dA
\\&=A^e\begin{Bmatrix}N_1^c\\N_2^c\\N_3^c\end{Bmatrix}[N_1^c\quad N_2^c\quad N_3^c]=\frac{1}{A^e}\begin{Bmatrix}\overline{N_1}^c\\\overline{N_2}^c\\\overline{N_3}^c\end{Bmatrix}\begin{bmatrix}\overline{N_1}^c&\overline{N_2}^c&\overline{N3}^c\end{bmatrix}
\end{aligned}
\end{equation}
\]
这一步注意角标的变化，这是如何做到的呢？
首先，注意到网格中心点的坐标可以这么表示：
\[\begin{equation}\begin{aligned}
x_{c} & = \frac{1}{3}\left(x_{1}+x_{2}+x_{3}\right)\\y_{c}&= \frac{1}{3}\left(y_{1}+y_{2}+y_{3}\right)
\end{aligned}
\end{equation}
\]
然后，我们把\(N_i (x,y)\)在单元中心点\((x_c, y_c)\)处取值，再令\(\overline{N_i}^c=A^e N_i^c\)，这样做是为了后续在面积\(A^e\)上积分方便，一乘就把系数干掉了。
可变积分函数（矩阵）用单元中心处的值来近似。如果单元尺寸足够小，这种近似是相当合理的。更严格的数学证明可在数值积分中找到，但超出了本模块的范围。
组装
节点处的连续性意味着由于单元内的线性插值，单元之间的连续性。
\[\begin{equation}\sum_{e=1}^n\int_A\left(T\left(\frac{\partial w}{\partial x}\frac{\partial\phi_j}{\partial x}+\frac{\partial w}{\partial y}\frac{\partial\phi_j}{\partial y}\right)-f(x,y)\phi_j(x,y)\right)dA
\end{equation}
\]
\[\begin{equation}\sum_{e=1}^{n}\left([K]^{e}\{\Delta\}^{e}-\{F\}^{e}\right)=[K]\{\Delta\}-\{F\}=0
\end{equation}
\]
线性形函数可以保证单元之间的连续性，即在相邻单元的公共边上，函数值是连续的。
这是因为共享边的两个三角形单元，在这条边上的值，均由这两个顶点唯一确定，而形函数又是线性的，所以两个单元共享的边上的函数值分布是相同的。
吐槽
真多啊，写的累，凌晨十二点半了。
