计算力学(四)弹性力学问题的变分解法(二)
加权残差法
现考虑一个微分方程\(L(u)+f=0\),其中\(L\)是微分算子,\(u\)是微分方程的解,\(g\)是已知函数.
构造一个试探解\(u_{app}\),想直接由试探解得到精确解通常很困难,所以我们想尽可能的减小试探解和精确解的差别,即减小残差(余量)\(R=L(u_{app})-L(u)=L(u_{app})+f\).\(u_{app}\)中有多个未知数,为了构造足够的方程,故而采用对残差加权的方法来判断残差什么情况下最小,选择不同的加权函数\(v\)就得到了不同的方法,包括
- 配置点法
- 子域法
- 最小二乘法
- 矩法
- 伽辽金法
其中伽辽金法的权函数就是试函数,具体的内容可以看这篇文章,后面只介绍伽辽金法.
强形式和弱形式
直接对残差进行积分,得到的方程形式 \[ \displaystyle\int_a^b v(L(u)+f)\mathrm dx=0\\ \] 就被称作强形式.对上式进行分部积分后就转化为了弱形式,本质上强形式和弱形式是一样的.强形式的主要缺点在于,所设的近似解函数需要有非常高阶的可微性.对应的,弱形式的优点在于可以通过分部积分将要求的可微次数减半,并且近似解函数只需满足Essential boundary condition,不必满足Natural boundary condition.(这里是直接引用了这篇文章,大佬写的太好了)下面结合例子来进行说明.
伽辽金法
现考虑悬臂梁上受到均布载荷问题.
列出其挠度方程 \[ EI\frac{\mathrm d\omega^4}{\mathrm d^4x}-q=0\\ \] 为了书写方便将其无量纲化 \[ \begin{align} \frac{\mathrm d(\frac{EI\omega}{ql^4})^4}{\mathrm d^4\frac{x}{l}}-1&=0\\ \frac{\mathrm du^4}{\mathrm d^4x}-1&=0\qquad 0\leq x\leq1\\ \end{align}\\ \] 相应的边界条件(Essential Boundary Conditions)为 \[ u\bigg|_{x=0}=0\quad,\quad u^{'}\bigg|_{x=0}=0\\ \] 构造近似解\(u=\sum a\phi\)其中\(a\)是系数,\(\phi\)即为选取的基函数,满足上面的边界条件(这里有问题,后面再说), 所以得到方程组 \[ \displaystyle\int_0^1\phi(u^{''''}-1)\mathrm dx=0\\ \] 这就是强形式.下面推导其弱形式 \[ \begin{align} LHS&=\displaystyle\int_0^1\phi\mathrm du^{'''}+\displaystyle\int_0^1\phi\mathrm dx\\ &=\phi u^{'''}\bigg|_0^1-\displaystyle\int_0^1\phi^{'}u^{'''}\mathrm dx+\displaystyle\int_0^1\phi\mathrm dx\\ &=\phi u^{'''}\bigg|_0^1-\phi^{'} u^{''}\bigg|_0^1+\displaystyle\int_0^1(\phi^{''}u^{''}+\phi)\mathrm dx\\ &=\phi u^{'''}\bigg|_{x=1}-\phi^{'} u^{''}\bigg|_{x=1}+\displaystyle\int_0^1(\phi^{''}u^{''}+\phi)\mathrm dx\\ \end{align}\\ \] 这里得到了近似解需要满足的另一类边界条件(essential boundary condition) \[ u^{'''}\bigg|_{x=1}=0\quad,\quad u^{''}\bigg|_{x=1}=0\\ \] 这里就回到了上面说到的问题,基函数同时满足这两类边界条件是很难的,强行构造则会出现其他的问题,例如构造\(x^5,x^6,x^7\cdots\)会损失低阶的部分.但让基函数满足Essential Boundary Conditions,再让近似解整体满足essential boundary condition是比较容易做到的.上面的近似解其实还需带入essential boundary condition找到系数的关系,使之整体满足边界条件.
现选取\(\phi_i=(1+x-e^{x})^i\)代入到弱形式中(呃呃,这里就不去找系数关系了,后面再解释),在matlab上计算得到2-6个基函数的结果
再选取\(\phi_i=(1-\cos x)^i\)同样代入到弱形式中,在matlab上计算得到2-5个基函数的结果
可以看到基函数选取的不同对结果的影响不小.
现在解释为啥不带入边界条件,找系数关系,呃呃主要因为太麻烦了,悬臂梁的均布载荷的解是多项式函数,我就不太想用\(x^i\)作基函数,选了\(\phi_i=(1+x-e^{x})^i\)然后惊奇的发现弱形式下自动满足边界条件了,然后强形式下得到的结果就不满足了,上面也就只给出了弱形式下结果.严谨的做法看这个回答.
有限元法解微分方程
从上面的结果可以看出,想要基函数在全区间上满足各种边界条件是非常困难的,而基函数如果是分段的就可以较好的满足这些条件.不过后面要写的有限元法是从最小势能原理推出来的,而不是从求解微分方程的角度出发,这里就不细说了,推荐一个网课数值分析 中科大讲的非常清楚.
简单来说,就是将待求解区间分为n段,对应n个基函数,每个基函数在对应区间上如下图所示(不一定是三角形的)
在其他区间上为零,这样就轻松的满足各类边界条件.同时这样处理让积分的过程也变的更加简单,只存在3种积分情况,并且都能容易的积出.系数的含义也更加明确了.
这个部分写的不是太严谨,借鉴了几位大佬的智慧,也有的地方是我自己的一些理解,希望看到的朋友不吝批评指正.(收藏的朋友顺便点点赞)
#! https://zhuanlan.zhihu.com/p/693721727