有限元入门 Part 1.

有限元入门 Part 1.

我终于又更新了!看了下上一篇文章已是两年前,当时看文章的同学现在都已经毕业了吧。翻了一下草稿箱,发现了这篇一年多前的未完成的文章哈哈哈。当时是我博士毕业答辩一周后写的,结果还没写完...现在认真补完,顺便提前祝大家新年快乐!

“一周前终于顺利答辩结束了,坑了一年的知乎专栏也准备继续更新。” <-- 这是我之前草稿的开场白lol

好了,开始正题。先总结一下之前的一系列文章:

首先,第一部分是变分法简介,介绍了极值化泛函数的方法。类似于函数求导,对泛函数求极值是对泛函数进行一阶变分并令其等于0,由此便可以转化为求解一个对应的偏微分方程(PDE)的问题。在大多数情况下,很难直接求得PDE的解析解,因此发展出了很多近似求解PDE的方法,即本专栏的第二部分。其中,在上一篇文章(weak form Galerkin)中提到了应用广泛的有限元方法是基于弱形式的。为了进一步举例说明,在这篇文章中,我们用小形变的线弹性模型举例。

首先,我们可以构造一个用来极值化的泛函数,即系统的总势能:总势能由系统内部应变能U和外力对系统做功所形成的势能V所构成:

\Pi = U+V

基于我们之前的线弹性模型的假设, 应变能密度 w_0=\frac{1}{2}\sigma_{ij}\epsilon_{ij} :

U=\int w_0d\Omega=\int\frac{1}{2}\sigma_{ij}\epsilon_{ij}d\Omega

V=-(\int b_i u_i d\Omega+\int t_i u_idA)

外力对系统做的功分为两类,第一种是体力( b_i )做功,所以是体积分,第二种是边界张力( t_i )做功,所以是边界面积分。外力如果做正功则系统势能降低,所以V加了一个负号。现在我们使用最小势能原理,对 \Pi 进行一阶变分(这里等效于对泛函数 \Pi 进行极值化):

\delta\Pi=\int \frac{1}{2}\sigma_{ij}\delta\epsilon_{ij}+\frac{1}{2}\epsilon_{ij}\delta\sigma_{ij}d\Omega-\int b_i\delta u_{i}d\Omega-\int t_i\delta u_i dA=0

然后带入应力-应变本构关系: \sigma_{ij}=C_{ijkl}\epsilon_{kl} , 由于 C_{ijkl} 的对称性,上式可以简化为:

\delta\Pi=\int \sigma_{ij}\delta\epsilon_{ij}d\Omega-\int b_i\delta u_{i}d\Omega-\int t_i\delta u_i dA=0

再带入应变-位移关系方程,对于小应变, \epsilon_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i}) 。利用矩阵的性质,一个对称矩阵 A_{ij} 和一个反对称矩阵 B_{ij} 的内积等于0( A_{ij}B_{ij}=0 ),则可得 \sigma_{ij}(u_{i,j}-u_{j,i})=0 。由此,上式可以简化为:

\delta\Pi=\int \sigma_{ij}\frac{1}{2}\delta(u_{i,j}+u_{j,i}+u_{i,j}-u_{j,i})d\Omega-\int b_i\delta u_{i}d\Omega-\int t_i\delta u_i dA=0

\delta\Pi=\int \sigma_{ij}\delta u_{i,j}d\Omega-\int b_i\delta u_{i}d\Omega-\int t_i\delta u_i dA=0

现在,其实我们已经得到弱形式的式子了,我们先把这个式子放在一边。这个式子是从极值化泛函数出发的。从之前的文章中我们知道,我们还可以从对应的PDE出发。对于固体力学,所解的PDE即平衡方程为:

-\sigma_{ij,j}-b_i=0

我们用弱形式的思想,将上式乘以试函数( \tilde{u}_i ),并在计算域上积分:

\int-(\sigma_{ij,j}+b_i)\tilde{u}_i d\Omega=0

利用分部积分,可以将上式转化为:

\int\sigma_{ij}\tilde{u}_{i,j}-b_i\tilde{u}_id\Omega-\int\sigma_{ij}n_j\tilde{u}_idA=0

对比之前变分法得到的 \delta\Pi 的式子,可以发现是一致的(在弱形式那篇文章中有论证一阶变分 \delta u_i 和试函数 \tilde{u}_i 的等效性), t_i=\sigma_{ij}n_j 是Neumann boundary condition, 而Dirichlet boundary condition对应于 \delta u_i=0 (或者换句话说, u_i=\hat{u}_i )。

到目前为止,我们用到的都是之前文章里介绍的知识,还没有提到有限元。在继续阅读前,推荐一下之前自己的一个回答,会对有限元的宏观上的理解有帮助:

有限元方法的核心思想是什么?www.zhihu.com图标

有限元的其中一个重要思想是用形函数作为检验函数。有限元中形函数是一个用来插值的函数,它使用单元节点值在单元内进行插值从而得到整个单元上的近似解。拿最简单的1D的单元进行举例,一个1D单元由两个节点构成,如果两个节点值用 u_1u_2 表示,两个节点上的形函数用 N_1(x)N_2(x) 表示,则这个单元上的任意一点的值可表示为: u(x)=u_1N_1(x)+u_2N_2(x) 。如果这个1D单元是从0到1的,那么带入两个节点的边界条件不难求得, N_1(x)=1-xN_2(x)=x

现在回到之前线弹性模型的例子,假设我们求解的是一个2D问题,并选用一阶三角形单元,每个单元有三个节点,每个节点有两个自由度(x和y方向的位移: u_x , u_y )。那么我们可以将位移场用节点位移和有限元形函数表示:

u_i=N_{ik}a_k

其中

u_i=[u_1(x,y)  \ \ \ \ u_2(x,y)]^T ,代表单元内的位移场,

N_{ik}= \left[   \begin{matrix}    N_1(x,y) & 0 & N_2(x,y) & 0 & N_3(x,y) & 0 \\    0 & N_1(x,y) & 0 & N_2(x,y) & 0 & N_3(x,y)   \end{matrix} \right]N_1, N_2, N_3为三个节点上的形函数,

a_k=\left[   \begin{matrix}    u_{1x} & u_{1y} & u_{2x} & u_{2y} & u_{3x} & u_{3y}   \end{matrix} \right]^Tu_{1x} 为节点1上的x方向的位移,其他分量以此类推。向量 a_k 为节点位移值向量,也是待求解的向量。

进而可得位移场的梯度: u_{i,j}=N_{ik,j}a_k 。带入之前的弱形式的式子可得(式子里的上标e代表该式子作用在某一个单元):

\begin{align} \delta\Pi&=\int \sigma_{ij}\delta u_{i,j}d\Omega^e-\int b_i\delta u_{i}d\Omega^e-\int t_i\delta u_i dA^e=0\\ &=\int C_{ijmn}\epsilon_{mn}\delta u_{i,j}d\Omega^e-\int b_i N_{ik}\delta a_k d\Omega^e -\int t_i N_{ik}\delta a_k dA^e\\ &=\int C_{ijmn}u_{m,n}\delta u_{i,j}d\Omega^e-\int b_i N_{ik}\delta a_k d\Omega^e -\int t_i N_{ik}\delta a_k dA^e\\ &=\int C_{ijmn}N_{mr,n} a_r N_{ik,j}\delta a_{k}d\Omega^e-\int b_i N_{ik}\delta a_k d\Omega^e -\int t_i N_{ik}\delta a_k dA^e \end{align}

之后可以把 \delta a_k 提出来并约去,另外由于 a_k 代表节点上的值,是固定的值,所以可以将其提到积分外,最后可以将上式简化为:

\int C_{ijmn}N_{mr,n} N_{ik,j}d\Omega^e a_r=\int b_i N_{ik} d\Omega^e +\int t_i N_{ik}dA^e

可以看出,积分 \int C_{ijmn}N_{mr,n} N_{ik,j}d\Omega^e 是一个二阶矩阵(由于有两个free indices --- k和r),可以将其记为 K^e_{kr}a_r 是节点值向量,也是待求的向量;等式右边只有一个free index,因此是一个向量,并且如果已经假定了形函数,是可以求解出来的,我们将整个等式右边记为 F^e_k 。这样,我们就可以将上式写成:

K^e_{kr}a_r=F^e_k

这就是大家经常看到的[k][U]=[F]的形式,我们想要求解的是节点值向量[U]。当然这只是一个单元内的,要得到全局解,还需要进行刚度矩阵的assembly。这个在下一篇文章中结合例子具体说明。这里总结一下,有限元是解节点上的值,然后再用形函数在单元内进行插值从而得到单元内任意一点的值。

本篇文章以线弹性模型为例,从之前文章中介绍的Weak form Galerkin出发,引入了有限元的重要概念形函数,并最终推导出了求解有限元的一般形式[K][U]=[F],从原理上解释了有限元方法,但还缺失很多细节。准备在下篇文章中,介绍如何用MATLAB写一个简单的有限元小程序,这样大家就会对有限元方法的细节有直观的认识。

现在毕业工作了,空余时间比以前多一些了,应该不会再鸽那么久了,尽量至少每个月更一篇吧哈哈。

发布于 2019-08-22

文章被以下专栏收录

    这是一个中二少年记录有关有限元,计算力学,固体力学的基础知识和自己的一些科研进展的专栏。