【游戏流体力学基础及Unity代码（十四）】舌尖上的有限元Galerkin法

伊索谗言 · 发表于 2021-3-18 14:17

本文主要介绍有限元方法中的Galerkin方法及如何用其解微分方程。
今晚吃啥，一直是一个困扰人类的千古难题。各路仁人志士为此耗费了无数的下午的时光，也没能得出统一的结论。而本文为计算出今晚究竟会吃到啥，提出了一种切实可行的方法，取名为“舌尖上的有限元法”。所有代码在https://github.com/clatterrr/FluidSimulationTutorialsUnity/tree/main/python/FiniteElement
今晚吃啥

假设有一条笔直的一维街道，从x = 0处开始，每隔dx = 1米处就有一个卖包子的晚餐店，每份晚餐包含的包子数量u就是自身所在坐标轴数值。如在x = 0处的就叫“零号店”，每顿晚餐有0个包子，在x = 1处的就叫“一号店“，每顿晚餐1个包子，以此类推。
而你的下午的学习生活工作主要集中在x0处。在度过这个美好下午后，你今晚会去哪家店吃，受到与这家店距离的影响。如果x0刚好在整数处，那么你只会去这里的晚餐店吃。如果不在整数处，那么你的左右各有一家店，你晚餐去这两家店的频率之比就是离这两家店距离的反比。比如你在x0 = 2.7处，那么离二号店距离0.7，离三号店距离0.3。那么你去某一晚去二号店的概率就是0.3，去三号店的概率就是0.7。那么平均每顿晚餐你可以吃上的包子数量U，就是0.3x2+ 0.7x3 = 2.7个包子。
写成公式如下。那个圆中间插一竖那个符号叫phi，代表你去某家店的概率，phi1就是去左边那家店的概率。那个反着的3叫做xi，代表你离你左边那家店的距离。xi的取值范围为[0,1]，因为距离与你超过1的晚餐店你都不会去。第三个式子，x是生活区的坐标，x1是左边店的坐标。

$\phi_1 = 1 - \xi \qquad \phi_2 = \xi \tag{1} \qquad \xi = x - x_1$
剧透一下，这里的phi就是有限元中插值函数，或者说是Galerkin方法中的权函数，或者基底函数等，或者形函数Shape Function。有的教程会把phi写成N，不过它们是同一个东西。
平均每顿晚餐吃到的包子数U = 去左边包子店的概率乘以左边店的每顿包子数 + 去右边包子店的概率乘以右边店每顿包子数，这句话翻译成数学公式如下。

$U = \begin{bmatrix} 1 - \xi & \xi \end{bmatrix}\begin{bmatrix} u_1 \\u_2 \end{bmatrix} = \phi^T u \tag{2}$
但是，如果你开始不知道每个晚餐店提供多少包子的话，就得自己去计算了。接下来都假设你并不知道每个晚餐的店的包子数量。但是你发现，如果你生活的地点x0是2.7，那么可以吃2.7个包子，而如果是2.8的话，那么你每晚可以吃到2.8个包子。换句话，这意味着生活地点只要前进0.1，那么每晚包子数量也能0.1。再换句话，每晚包子数量增加量除以生活地点前进的距离再减去一等于零。再再换句话：

$\frac{d U}{dx} - 1 = 0 \tag{3}$
上式是个伟大的定律，在这定律的背后，必定蕴藏着每个包子店晚餐到底提供多少包子的真相。你的任务就是去探索真相。而浮出幕后的第二个线索是而且x = 1处的1号包子店提供1个包子，这也叫初始条件。
而且你发现，这个定律不受与你去哪家店的概率影响，只要前进0.1，那么包子数量就能增强0.1。我们可以利用这个方法来探索另一个伟大的定律。也就是

$\int_{x1}^{x2} \phi(\frac{U}{dx} - 1)dx = 0 \tag{4}$
x1是左边店的坐标，x2是右边店的坐标。
这个方法也叫做Galerkin法。Galerkin法是有限元中的一种加权残差法。除了Galerkin法还有点装配法和子域法等。好了，你可以把这段话给忘了，继续往下看。
上面这个式子4也可以换成如下的式子。记住dx是生活地点前进的距离，dU是多吃到的包子数量，公式3已经说明它们比值为1。

$\int_{x1}^{x2} \phi(dU) = \int_{x1}^{x2} \phi dx \tag{5}$
再来看一下公式2并求导。注意不带上标的u是指晚餐包子铺提供的包子数量。

$dU = d(\phi^T u) = d(\phi^T)u = \begin{bmatrix} \phi_{1x} & \phi_{2x} \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \end{bmatrix}\tag{6}$
对phi求导得到如下。

$\phi_{1x} = -1 \qquad \phi_{2x} = 1 \tag{7}$
上面这个式子上实际上在说，生活地点每前进0.1，那么去左边店的概率就减少0.1，所以phi{1x} = -1。而去右边店的概率就增加0.1，所以phi{2x} = 1。不带上标的u指每个包子店卖的包子。公式5和6翻译成人话就是
距离左边店的距离变化量乘 phi1x 乘左边店的包子数量 + 距离右边店的距离变化量乘 phi2x 乘右边店的包子数量等于包子的总变化量，也就是 0.1x(-1)x2 + 0.1x1x3 = 0.1。
将公式5带入公式4，得到如下式子。注意我已将积分范围设置在左右两家店之间，让左边的店坐标为0，右边店坐标为1。那么生活区的坐标dx其实就是距左边店的距离dxi。这也是使用Galerkin方法解一阶偏微分方程的一般形式。

$\int_{0}^{1}\begin{bmatrix} \phi_{1} \\ \phi_{2} \end{bmatrix}\begin{bmatrix} \phi_{1x} & \phi_{2x} \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \end{bmatrix}d \xi = \int_{0}^{1}\begin{bmatrix} \phi_{1} \\ \phi_{2} \end{bmatrix}d\xi \tag{8}$
上式的未知数就是左边店与右边店的包子数量u1和u2。对于我们这个问题来说，首先将phi{1x}和phi{2x}代入得

$\int_{0}^{1}\begin{bmatrix} -\phi_{1} & \phi_{1}\\ -\phi_{2} & \phi_{2}\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \end{bmatrix}d \xi = \int_{0}^{1}\begin{bmatrix} \phi_{1} \\ \phi_{2} \end{bmatrix}d\xi \tag{9}$
然后积分

$\int_0^1 \phi_1 d\xi = \int_0^1(1-\xi)d\xi = \frac{1}{2} \qquad \int_0^1 \phi_2 d\xi = \int_0^1\xi d\xi = \frac{1}{2} \tag{10}$
将10代入9得

$\begin{bmatrix} -1/2 & 1/2\\ -1/2 & 1/2\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \end{bmatrix}= \begin{bmatrix} 1/2 \\ 1/2 \end{bmatrix} \tag{11}$
上面虽然只有两个相同的线性方程组，但别忘了第二个线索即初始条件u1 = 1。由此解得u2 = 2。如果也想顺手知道u3和u4的值，那么上式也可写成下面的样子

$\begin{bmatrix} -1/2 & 1/2 & 0 & 0\\ -1/2 & 1/2 & 0 & 0 \\0 & 0 & 0 & 0 \\0 & 0 & 0 & 0\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\u_3 \\u_4 \end{bmatrix}= \begin{bmatrix} 1/2 \\ 1/2 \\ 0 \\ 0 \end{bmatrix} \tag{12}$
上面的式子看起来很沙雕也没什么用，实际上确实很沙雕也没什么用。但是如果将其与下面两个式子相加

$\begin{bmatrix} 0 & 0 & 0 & 0 \\0 &-1/2 & 1/2 & 0\\ 0 & -1/2 & 1/2 & 0 \\0 & 0 & 0 & 0\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\u_3 \\u_4 \end{bmatrix}= \begin{bmatrix} 0 \\ 1/2 \\ 1/2 \\ 0 \end{bmatrix} \\ \begin{bmatrix} 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 \\0 & 0 & -1/2 & 1/2 \\ 0 & 0 & -1/2 & 1/2 \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\u_3 \\u_4 \end{bmatrix}= \begin{bmatrix} 0 \\ 0 \\ 1/2 \\ 1/2 \\ \end{bmatrix}$
就得到了

$\begin{bmatrix} -1/2 & 1/2 & 0 & 0 \\-1/2 & 0 & 1/2 & 0\\ 0 & -1/2 & 0 & 1/2 \\0 & 0 & -1/2 & 1/2\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\u_3 \\u_4 \end{bmatrix}= \begin{bmatrix} 1/2 \\ 1 \\ 1 \\ 1/2 \end{bmatrix} \tag{13}$
这就是平衡矩阵，也可以写成下面的形式

$Ku = b \tag{14}$
在固体力学中，上面的方程也叫刚度方程。K叫做刚度矩阵(stiffness matrix)，刚度矩阵相加也叫装配。如果在处理与弹簧有关的有限元问题中，u就是节点位移，b就是节点力。
13式虽然是4x4的矩阵，但只有三个不同的线性方程组。因此还要一个初始条件。为了求出结果，我们假设u1 = 5，那么上式改写成如下

$\begin{bmatrix} 1 & 0 & 0 & 0 \\-1/2 & 0 & 1/2 & 0\\ 0 & -1/2 & 0 & 1/2 \\0 & 0 & -1/2 & 1/2\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\u_3 \\u_4 \end{bmatrix}= \begin{bmatrix} 5 \\ 1 \\ 1 \\ 1/2 \end{bmatrix} \tag{15}$
初始条件可以是某个点的值，但也可是某处的梯度。如果u1 = u2，那么矩阵应该如下设置

$\begin{bmatrix} 1 & -1 & 0 & 0 \\-1/2 & 0 & 1/2 & 0\\ 0 & -1/2 & 0 & 1/2 \\0 & 0 & -1/2 & 1/2\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\u_3 \\u_4 \end{bmatrix}= \begin{bmatrix} 0 \\ 1 \\ 1 \\ 1/2 \end{bmatrix}\tag{16}$
求K的逆矩阵时，暂时用高斯消元法即可。由于主对角线上大部分为零，所以将矩阵k和矩阵b整体上移一位就可用普通高斯消元法。由于上面矩阵K大部分元素都是零，因此又被称为稀疏矩阵。如何利用稀疏矩阵的特性快速求它的逆矩阵，则是另外的话题了。完整代码文件为fem1d.py。其中组装矩阵的部分代码如下。注意我的代码文件中u是从0开始编号的。
for i in range(0,nmax-1):
kmat[i:i+2,i:i+2] += [[-0.5,0.5],[-0.5,0.5]]
bmat[i:i+2] += [a*0.5,a*0.5]今晚吃啥二维函数季

你好不容易算完了，你的朋友们也都已经出去吃晚餐了，你也正打算出去解决晚餐的时候，你突然想到城市现在一般都是二维的，这样城市被化为一块块的正方形，你应该居住在某个正方形之间。正方形的顶点上各有一个晚餐店，因此你周围有四个晚餐店。晚餐提供包子的数量随x轴和y轴变化。
这样，并且假设店的分布位置如下。那么你离一号店和四号店的x轴距离为x，离1号店和二号店的y轴距离为y。下图括号外的数是晚餐店的编号，而括号内的数是实际的配给的包子数量。
那么你去四个店的概率分别为

$\phi_1 = (1-x)(1-y) \qquad \phi_2 = x(1-y) \qquad\phi_3 = (1-x)y \qquad \phi_4 = xy \tag{17}$
上面的概率，或者说是权函数加起来刚好为一。假设每个店提供的包子数量为u，、那么你平均每晚吃到的包子数为

$U = \phi_1u_1+ \phi_2u_2 + \phi_3u_3 + \phi_4u_4 = \phi^T u\tag{18}$
初始条件可能如下，U沿着x轴或y轴都是一次函数，且ab为常数

$\qquad \frac{\partial U}{\partial x} + a\frac{\partial U}{\partial y} + b = 0 \tag{19}$
然后Galerkin法

$\int_{0}^{1}\int_{0}^{1} \phi(\frac{\partial U}{\partial x} + a\frac{\partial U}{\partial y} + b)dxdy = 0 \tag{20}$

$\int_{0}^{1}\int_{0}^{1} \phi\frac{\partial U}{\partial x} dxdy + a\int_{0}^{1}\int_{0}^{1} \phi\frac{\partial U}{\partial y} dxdy = -b\int_{0}^{1}\int_{0}^{1} \phi dxdy \tag{21}$
对于上式左边，继续代入并展开

$\int_0^1\int_{0}^{1}\begin{bmatrix} \phi_{1} \\ \phi_{2} \\\phi_{3} \\ \phi_{4} \end{bmatrix}\begin{bmatrix} \phi_{1x} & \phi_{2x} &\phi_{3x} & \phi_{4x} \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\u_3 \\u_4 \end{bmatrix}dxdy \\+a\int_0^1\int_{0}^{1}\begin{bmatrix} \phi_{1} \\ \phi_{2} \\\phi_{3} \\ \phi_{4} \end{bmatrix}\begin{bmatrix} \phi_{1y} & \phi_{2y} &\phi_{3y} & \phi_{4y} \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\u_3 \\u_4 \end{bmatrix}dxdy = -b\int_0^1\int_{0}^{1}\begin{bmatrix} \phi_{1} \\ \phi_{2} \\\phi_{3} \\ \phi_{4} \end{bmatrix}dxdy \tag{22}$
算phi的积分与微分。如果你对二重积分不熟悉，可以看看我的一个回答https://www.zhihu.com/question/44875342/answer/821701089，主题仍然是舌尖上的二重积分。

$\phi_1 = (1-x)(1-y) \qquad \phi_2 = x(1-y) \qquad\phi_3 = (1-x)y \qquad \phi_4 = xy\\ \frac{\partial \phi_1}{\partial x} = y - 1 \qquad \frac{\partial \phi_2}{\partial x} = 1-y\qquad \frac{\partial \phi_3}{\partial x} = -y\qquad \frac{\partial \phi_4}{\partial x} = y \\ \frac{\partial \phi_1}{\partial y} = x - 1 \qquad \frac{\partial \phi_2}{\partial y} = -x\qquad \frac{\partial \phi_3}{\partial y} = 1-x\qquad \frac{\partial \phi_4}{\partial y} = x\\ \iint \phi_1 = \frac{1}{4} \qquad \iint \phi_2 = \frac{1}{4} \qquad\iint \phi_3 = \frac{1}{4} \qquad\iint \phi_4 = \frac{1}{4}$
接下来就是爆肝的算参数时刻了。

$\begin{bmatrix} -1/6 & 1/6 & -1/12 & 1/12 \\-1/6 & 1/6 & -1/12 & 1/12 \\ -1/12 & 1/12 & -1/6 & 1/6 \\ -1/12 & 1/12 & -1/6 & 1/6 \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\u_3 \\u_4 \end{bmatrix} + a\begin{bmatrix} -1/6 & -1/12 & 1/6 & 1/12 \\-1/12 & -1/6 & 1/12 & 1/6 \\ -1/6 & -1/12 & 1/6 & 1/12 \\ -1/12 & -1/6 & 1/12 & 1/6\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\u_3 \\u_4 \end{bmatrix} =- b\begin{bmatrix} 1/4 \\1/4 \\ 1/4 \\ 1/4\end{bmatrix}\tag{23}$
稍微解释一下，左边矩阵最最左上角的元素是这样计算的

$\int_0^1 \int_0^1\phi_1 \phi_{1x} = \int_0^1 \int_0^1(1-x)(1-y)(y-1) = -\frac{1}{6} \tag{24}$
同样要解上面的矩阵，必须给出一些点的确切值。完整代码文件为fem2d.py。据我测试，似乎是需要三个初始条件才能正确求解出逆矩阵，但我暂时想不通为什么是三个。
如果把式子19稍微改一下，把对y微分改为对时间t微分，那么就成了线性对流方程

$\frac{\partial U}{\partial t} + a\frac{\partial U}{\partial x} = 0\tag{25}$
源代码稍经修改，结果如下。从左往右为x轴，从上往下为t轴。U为1的地方，在第零秒还在第一个网格，在第三秒就到了第2个网格。完整代码文件为fem2dadvection.py
今晚吃啥二元函数季

你的朋友们都吃完晚餐回来了，他们告诉你晚餐不仅提供包子，也开始提供馒头了。本着严谨认真的态度，你决定继续研究下去。为了化简这个问题，你继续把城市拍扁成一维的。晚餐店提供包子馒头关系式如下

$\frac{\partial U}{\partial x} + a\frac{\partial U}{\partial y} + b\frac{\partial V}{\partial x} = 0 \qquad \frac{\partial V}{\partial y} = 0\tag{26}$
其中U是平均每晚能吃到的包子数量，V是平均每晚吃到的馒头数量。

$U = (1-x)(1-y)u_1 + x(1-y) u_2 + (1-x)yu_3 + xyu_4 = \phi^Tu\\ V = (1-x)(1-y)v_1 + x(1-y) v_2 + (1-x)yv_3 + xyv_4 = \phi^Tv\\$ 虽然包子和馒头食材不同，但这和你晚上去哪个晚餐店没什么关系。

$\int_{0}^{1}\int_{0}^{1} \phi(\frac{\partial U}{\partial x} + a\frac{\partial U}{\partial y} + b\frac{\partial V}{\partial x})dxdy = 0 \tag{27}$
替换

$\int_{0}^{1}\int_{0}^{1} \phi\frac{d\phi_T}{dx}udxdy + a\int_{0}^{1}\int_{0}^{1} \phi\frac{d\phi_T}{dy}udxdy +\int_{0}^{1}\int_{0}^{1} b\phi\frac{d\phi_T}{dx}vdxdy = 0 \tag{28}$
与之上一节相比，并不需要额外计算什么东西，最后算得如下式子。

$\begin{bmatrix} -1/6 & 1/6 & -1/12 & 1/12 \\-1/6 & 1/6 & -1/12 & 1/12 \\ -1/12 & 1/12 & -1/6 & 1/6 \\ -1/12 & 1/12 & -1/6 & 1/6 \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\u_3 \\u_4 \end{bmatrix} + a\begin{bmatrix} -1/6 & -1/12 & 1/6 & 1/12 \\-1/12 & -1/6 & 1/12 & 1/6 \\ -1/6 & -1/12 & 1/6 & 1/12 \\ -1/12 & -1/6 & 1/12 & 1/6\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\u_3 \\u_4 \end{bmatrix} \\+ b\begin{bmatrix} -1/6 & 1/6 & -1/12 & 1/12 \\-1/6 & 1/6 & -1/12 & 1/12 \\ -1/12 & 1/12 & -1/6 & 1/6 \\ -1/12 & 1/12 & -1/6 & 1/6 \end{bmatrix}\begin{bmatrix} v_1 \\ v_2 \\v_3 \\v_4 \end{bmatrix} = 0 \tag{29}$
如果剧透原函数为U = x + y，且V = x。那么我们知道的初始条件式子a = 1，b = -2。并且由于dV/dy = 0，那么v1 = v3而且v2 = v4。再加上初始条件后最后算得，u2 = u3 = v2 = v4 = 1，u1 = v1 = v3 = 0，u4 = 2。完美契合原函数。你实在不信，拿那个4x4矩阵第一行检验一下也行，如下

$(-\frac{1}{6}u_1 + \frac{1}{6}u_2 - \frac{1}{12}u_3 + \frac{1}{12}u_4) + a(-\frac{1}{6}u_1 - \frac{1}{12}u_2 + \frac{1}{6}u_3 + \frac{1}{12}u_4)\\ + b(-\frac{1}{6}v_1 + \frac{1}{6}v_2 - \frac{1}{12}v_3 + \frac{1}{12}v_4) = 0$
如果你无法在一秒钟之内想出这样的验证方法的话，你应该补一下线性代数与矩阵的知识。这部分编程留给读者了。
今晚吃啥二次函数季

终于算出来了每个晚餐店的馒头包子配给，这使你充满了决心。但你突然发现你其实对馒头没什么兴趣，也不愿意去对面街道吃包子。于是继续把包子店拍扁成一维的，而晚餐点仍然只搭配包子这一种小吃。但是吧，包子种类越来越多，比如喷香鲜肉包，多汁榨菜包等，每个晚餐店提供的包子随坐标轴变化呈二次函数关系。为了找出这个二次函数，你需要三个点来确定这个二次函数。也就是说这次你把每三个店首尾相连的分为一组，你每晚可能去这三店中的一个。也就是如果你住在x = 0到x = 2之间，那么就会受到零一二号店的影响。而如果住在x = 2到x = 4之间，就会受到二三四号店的影响，依次类推。
不过二维函数的权函数有两种，第一种类似于二阶精度前向后向差分，计算点为三点中的最左边或最右边。这种权函数的推导方法为，假设每晚吃到包子数量为

$U = c_1x^2 + c_2x + c_3 \tag{30}$
那么

$U(0) = c_3 \qquad U(1) = c_1 + c_2 + c_3 \qquad U(2) = 4c_1 + 2c_2 + c_3 \tag{31}$
令c = Au，那么

$A = \begin{bmatrix} 0 & 0 & 1 \\ 1 & 1 &1 \\ 4 & 2 & 1\end{bmatrix}^{-1} = \begin{bmatrix} 1/2 & -1 & 1/2 \\ -2/3 & 2 &-1/2 \\ 1 & 0 & 0\end{bmatrix} \tag{32}$
那么

$\phi_1 = \frac{1}{2}\xi^2 - \frac{3}{2}\xi + 1 \qquad \phi_2 = - \xi^2 + 2\xi \qquad \phi_3 = \frac{1}{2}\xi^2-\frac{1}{2}\xi \tag{33}$
那么U就应该用如下式子表达

$U = \phi_1u_1 + \phi_2u_2 + \phi_3 u_3 = \phi^Tu \tag{34}$
验证一下式子的正确性。假设U = x^2 + x + 1，那么U0 = 1，U2 = 3，U3为7。也就是

$\xi^2 + \xi + 1 = (\frac{1}{2}\xi^2 - \frac{3}{2}\xi + 1) + 3(- \xi^2 + 2\xi) + 7( \frac{1}{2}\xi^2-\frac{1}{2}\xi) \tag{35}$
如果选择U1 = 3，U2 = 7，U3 = 13，最后算出来的就是(x+1)^2 + (x+1) + 1，也是正确的。毕竟这次最左边的店需要往左移动一个。但是你也可以将31式子改成算U(0)和U(0.5)和U(1)的，翻译成人话就是现在每隔0.5个网格就有一个包子店，因此每个网格包括首尾有三个包子店，这时候你的活动范围就是这一个网格里的三个包子店。最后算出的权函数也有所不同。这种方法的x的积分域是0到2。不过有的教程喜欢用xi来代替x，xi的积分域是0到1，这样x = 2xi，dx = 2dxi。
另一种权函数类似于中心差分的计算方法，计算的位置在中间。针对一维二次函数的权函数为

$\phi_1 = \frac{1}{2}\xi(\xi - 1) \qquad \phi_2 = 1 -\xi^2 \qquad \phi_3 = \frac{1}{2}\xi(\xi+1) \tag{36}$
再算如果U[0] = 1，U[1] = 3，U[2] = 7，那么上式计算出(x + 1)^2 + (x + 1) + 1，代入x = 0算出的是第二个值3。这种方法的积分域是-1到1。实际算参数的时候，这两种权函数没什么差别。
不过现在我们暂时用第一种形式的权函数。把假设忘了，我们现在只知道下面的公式，U是晚餐店搭配包子的数量

$\frac{\partial^2 U}{\partial x^2} - 2 = 0 \tag{37}$
且零号店搭配一个包子即U(0) = 1，这时候如何反求出U呢？我们仍然可以给上面积分一个权函数

$\int_{0}^{2} \phi\frac{dU^2}{dx} = 2\int_{0}^{2} \phi dx \tag{38}$
然后分部积分

$[\phi \frac{dU}{d\xi}]_0^2 - \int_0^2 \frac{d\phi}{d\xi}\frac{dU}{d\xi} d\xi = 2\int_0^2 \phi d \xi \tag{39}$
为了验证上面式子正确性，继续手算。先对权函数求导和积分

$\frac{\partial \phi_1}{\partial \xi} = \xi - \frac{3}{2} \qquad \frac{\partial \phi_2}{\partial \xi} = -2\xi + 2 \qquad \frac{\partial \phi_2}{\partial \xi} = \xi - \frac{1}{2}\\ \int_0^2 \phi_1d\xi = \frac{1}{3} \qquad \int_0^2 \phi_2d\xi = \frac{4}{3} \qquad \int_0^2 \phi_3d\xi = \frac{1}{3}$
再次假设我们提前拿到剧本，知道U = x^2 + x + 1，那么dU = 2x + 1，然后将phi1代入

$((\frac{1}{2}\xi^2 - \frac{3}{2}\xi + 1)(2\xi +1))|_0^2 - \int_0^2(\xi - \frac{3}{2})(2\xi+1)d\xi = 2\int_0^2 (\frac{1}{2}\xi^2 - \frac{3}{2}\xi + 1)d\xi \tag{40}$
计算可得

$(-1) - (-\frac{5}{3}) = 2(\frac{1}{3}) \tag{41}$
确实是正确的。代入phi2或phi3一样也不会出错。推导公式时每一步都要保证正确是很重要的。
继续把假设的解析解忘了，现在我们知道式子分部积分后39式肯定是对的，然后再将34式代入39式，得到

$[\phi \frac{d\phi^T}{d\xi}u]_0^2 - \int_0^2 \frac{d\phi}{d\xi}\frac{d\phi^T}{d\xi}u d\xi = 2\int_0^2 \phi d \xi \tag{42}$
然后展开

$(\begin{bmatrix} \phi_{1} \\ \phi_{2} \\ \phi_{3} \end{bmatrix}\begin{bmatrix} \phi_{1x} & \phi_{2x} & \phi_{3x} \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix})_0^2- \\\int_0^2 \begin{bmatrix} \phi_{1x} \\ \phi_{2x} \\ \phi_{3x} \end{bmatrix}\begin{bmatrix} \phi_{1x} & \phi_{2x} & \phi_{3x} \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix}dx = 2\int_0^2 \begin{bmatrix} \phi_{1} \\ \phi_{2} \\ \phi_{3} \end{bmatrix} dx \tag{43}$
计算可得。计算这一步不会出现最终的有限元程序中，这些值在编写有限元程序前就需要手动计算好，然后作为参数填入有限元程序。虽然这是用上述第一种权函数算出来的，但如果用式子35的第二种权函数，也可以算得完全一样的值。注意第二种权函数积分范围为-1到1。

$\begin{bmatrix} 3/2 & -2 & 1/2 \\ 0 & 0 & 0\\ 1/2 & -2 & 3/2 \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix}- \begin{bmatrix} 7/6 & -4/3 & 1/6 \\ -4/3 & 8/3 & -4/3\\ 1/6 &-4/3 & 7/6 \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix} = 2\begin{bmatrix} 1/3 \\ 4/3 \\ 1/3 \end{bmatrix} \tag{44}$
合并一下，得到

$\begin{bmatrix} 1/3 & -2/3 & 1/3 \\ 4/3 & -8/3 & 4/3\\ 1/3 &-2/3 & 1/3 \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix} = a\begin{bmatrix} 1/3 \\ 4/3 \\ 1/3 \end{bmatrix} \tag{45}$
很好，饿着肚子算半天算出来个这么个穷酸矩阵，表面上看有三个方程组，其实它们都是一样的。一个方程组没法求出三个变量，因此还需要额外两个初始条件。干脆把它们化成如下形式

$\begin{bmatrix} 1/3 & -2/3 & 1/3 \\ 0 & 0 & 0 \\ 0 & 0 & 0\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix} = a\begin{bmatrix} 1/3\\ 0 \\ 0 \end{bmatrix} \tag{46}$
这个方程可以用于解下面的二阶偏导数，其中a是常量，对于我们的问题而言a = 2。

$\frac{\partial^2 U}{\partial x^2} - a = 0 \tag{47}$
最终用于解一维二次函数的矩阵如下。在解之前还要把k矩阵最后两行替换上初始条件

$\begin{bmatrix} 1/3 & -2/3 & 1/3 & 0 & 0\\ 0 &1/3 & -2/3 & 1/3 & 0\\ 0& 0 &1/3 & -2/3 & 1/3 \\0& 0& 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ u_3 \\u_4 \\ u_5\end{bmatrix} = a\begin{bmatrix} 1/3 \\ 1/3 \\ 1/3 \\ 0 \\ 0 \end{bmatrix} \tag{48}$
上面就是用有限元法Galerkin法解一维二次函数的一般形式了。完整代码文件为fem1dx2.py。二次函数稍微变化一下就成了流体力学中的扩散方程。
突然你意识到你自己晚餐还没吃上，很想抓起身边某应急食品来尽享美味。但用Galerkin法解NavierStokes方程的准备都做好了，怎么能不立刻开始呢？
这篇所有代码都在https://github.com/clatterrr/FluidSimulationTutorialsUnity/tree/main/python/FiniteElement。下一篇用有限元Galerkin法解NavierStokes方程，敬请期待！可压缩流体那几篇预计再跳票几个月。
上一篇：光影帽子：【游戏流体力学基础及Unity代码（十三）】压力泊松方程,python写SIMPLE算法
参考

《计算流体力学有限元方法及其编程详解》作者：毕超。很棒的解释得很详细的一本书。也是从解微分方程开始解释有限元，而不是从杆弹簧模型之类开始。书中附有matlab代码。
https://zhuanlan.zhihu.com/p/264901969
https://www.bilibili.com/video/BV1d4411i7Wr?from=search&seid=7135824331022518280 有关固体力学中的有限元方法慕课，讲的非常详细很棒
《Finite Element Procedures》
剩下的倒没什么推荐的，因为有限元的书里有50%是使用商业软件的，40%是那种满篇公式初学者很难理解的。剩下的不错的基本上都是在说固体力学。我倒是觉得与其一开始就说杆单元弹簧单元什么的，不如从一阶二阶偏微分方程开始说起，多举一些手算的例子，然后把刚度矩阵形函数这些吓人的名词换成比较亲切的比喻，这样对初学者更容易理解一些。

		自动登录	找回密码
密码			立即注册

[笔记] 【游戏流体力学基础及Unity代码（十四）】舌尖上的有限元Galerkin法

本帖子中包含更多资源

浏览过的版块