计算力学(六)坐标变换和高斯积分

发表于 分类于 ， ， 阅读次数： 本文字数： 8.1k 阅读时长 ≈ 7 分钟

坐标变换

在实际应用中,单元可能处于整体坐标系中的任意位置,具有任意大小,将整体坐标系中的运算转化到局部坐标系中,不同的单元有公共的基准,方便程序的设计和实现.

一维坐标变换

对于变换\(x=\displaystyle\sum_{i=1}^3N_i(\xi)x_i\),其逆变换\(\xi=\xi(x)\)的存在条件是 \[ \frac{\mathrm dx}{\mathrm d\xi}\neq0\\ \] 一维坐标变换示意图

上图中 \[ x=\displaystyle\sum_{i=1}^3N_i(\xi)x_i=\frac{1}{2}(\xi^2-\xi)x_1+\frac{1}{2}(\xi^2+\xi)x_2+(1-\xi^2)x_3\\ \frac{\mathrm dx}{\mathrm d\xi}=\sum_{i=1}^3\frac{\partial N_i}{\partial\xi}x_1=(x_1+x_2-2x_3)\xi+\frac{1}{2}(x_2-x_1)\\ 若2x_3=x_1+x_2\quad\frac{\mathrm dx}{\mathrm d\xi}=\frac{1}{2}(x_2-x_1)\neq0\quad逆变换始终存在\\ 若2x_3\neq x_1+x_2\quad由\frac{\mathrm dx}{\mathrm d\xi}=0得到\quad\xi=\frac{0.5(x_2-x_1)}{2x_3-x_1-x_2}时逆变换不存在\\ \] 插值函数在两个坐标系中的导数关系 \[ \frac{\mathrm du}{\mathrm dx}=\frac{\mathrm du}{\mathrm d\xi}\cdot\frac{1}{\dfrac{\mathrm dx}{\mathrm d\xi}}=\frac{\mathrm du}{\mathrm d\xi}\cdot\frac{1}{\displaystyle\sum_{i=1}^3\dfrac{\mathrm dN_i}{\mathrm d\xi}x_i}\\ \] 积分微元关系 \[ \mathrm dx=\frac{\mathrm dx}{\mathrm d\xi}\mathrm d\xi=(\sum_{i=1}^3\frac{\mathrm dN_i}{\mathrm d\xi}x_i\mathrm)d\xi\\ \]

二维坐标变换

正方形单元

二维中函数偏导数的关系 \[ \frac{\partial u}{\partial \xi}=\frac{\partial u}{\partial x}\frac{\partial x}{\partial \xi}+\frac{\partial u}{\partial y}\frac{\partial y}{\partial \xi}\\ \frac{\partial u}{\partial \eta}=\frac{\partial u}{\partial x}\frac{\partial x}{\partial \eta}+\frac{\partial u}{\partial y}\frac{\partial y}{\partial \eta}\\ \] 偏导数的变换关系为 \[ \frac{\partial}{\partial \xi}=\frac{\partial}{\partial x}\frac{\partial x}{\partial \xi}+\frac{\partial}{\partial y}\frac{\partial y}{\partial \xi}\\ \frac{\partial}{\partial \eta}=\frac{\partial}{\partial x}\frac{\partial x}{\partial \eta}+\frac{\partial}{\partial y}\frac{\partial y}{\partial \eta}\\ \] 写成矩阵形式 \[ \begin{bmatrix} \dfrac{\partial}{\partial\xi}\\ \dfrac{\partial}{\partial\eta} \end{bmatrix}=\textbf{J} \begin{bmatrix} \dfrac{\partial}{\partial x}\\ \dfrac{\partial}{\partial y} \end{bmatrix}\\ \] 其中\(\textbf{J}\)为Jacobi矩阵 \[ \textbf{J}= \begin{bmatrix} \dfrac{\partial x}{\partial \xi}&\dfrac{\partial y}{\partial \xi}\\ \dfrac{\partial x}{\partial \eta}&\dfrac{\partial y}{\partial \eta} \end{bmatrix}\\ \] 偏导数的逆变换形式 \[ \begin{bmatrix} \dfrac{\partial}{\partial x}\\ \dfrac{\partial}{\partial y} \end{bmatrix}=\textbf{J}^{-1} \begin{bmatrix} \dfrac{\partial}{\partial \xi}\\ \dfrac{\partial}{\partial \eta} \end{bmatrix}=\frac{1}{|\textbf{J}|} \begin{bmatrix} \dfrac{\partial y}{\partial \eta}&-\dfrac{\partial y}{\partial \xi}\\ -\dfrac{\partial x}{\partial \eta}&\dfrac{\partial x}{\partial \xi} \end{bmatrix} \begin{bmatrix} \dfrac{\partial}{\partial \xi}\\ \dfrac{\partial}{\partial \eta} \end{bmatrix}\\ \] 其中\(|\textbf{J}|=\dfrac{\partial x}{\partial\xi}\dfrac{\partial y}{\partial\eta}-\dfrac{\partial y}{\partial\xi}\dfrac{\partial x}{\partial\eta}\).

当函数是插值基函数\(N_i\)​时 \[ \begin{bmatrix} \dfrac{\partial N_i}{\partial x}\\ \dfrac{\partial N_i}{\partial y} \end{bmatrix}=\textbf{J}^{-1} \begin{bmatrix} \dfrac{\partial N_i}{\partial \xi}\\ \dfrac{\partial N_i}{\partial \eta} \end{bmatrix}\\ \] 函数对整体坐标的偏导数 \[ \dfrac{\partial u}{\partial x}=\sum\dfrac{\partial N_i}{\partial x}u_i,\quad\dfrac{\partial u}{\partial y}=\sum\dfrac{\partial N_i}{\partial y}u_i\\ \] 积分微元关系 \[ \mathrm dx\mathrm dy=\det\textbf{J}\cdot\mathrm d\xi\mathrm d\eta\\ \]

三角形单元

由于三个面积坐标\(L_1,L_2,L_3\)中只有两个独立,任取其中两个作为坐标变换中独立的面积坐标.令 \[ \xi=L_1\qquad\eta=L_2\quad由此得到\quad L_3=1-\xi-\eta\\ \] 令函数\(f\)为面积坐标的函数 \[ f=f(L_1,L_2,L_3)\\ \] 得到函数\(f\)对独立面积坐标和原始面积坐标的偏导数之间的关系 \[ \frac{\partial f}{\partial\xi}=\frac{\partial f}{\partial L_1}\frac{\partial L_1}{\partial\xi}+\frac{\partial f}{\partial L_2}\frac{\partial L_2}{\partial\xi}+\frac{\partial f}{\partial L_3}\frac{\partial L_3}{\partial\xi}=\frac{\partial f}{\partial L_1}-\frac{\partial f}{\partial L_3}\\ \frac{\partial f}{\partial\eta}=\frac{\partial f}{\partial L_1}\frac{\partial L_1}{\partial\eta}+\frac{\partial f}{\partial L_2}\frac{\partial L_2}{\partial\eta}+\frac{\partial f}{\partial L_3}\frac{\partial L_3}{\partial\eta}=\frac{\partial f}{\partial L_2}-\frac{\partial f}{\partial L_3}\\ \] 得到Jacobi矩阵 \[ \textbf{J}=\begin{bmatrix} \dfrac{\partial x}{\partial \xi}&\dfrac{\partial y}{\partial \xi}\\ \dfrac{\partial x}{\partial \eta}&\dfrac{\partial y}{\partial \eta} \end{bmatrix}= \begin{bmatrix} \displaystyle\sum(\dfrac{\partial N_i}{\partial L_1}-\dfrac{\partial N_i}{\partial L_3})x_i&\displaystyle\sum(\dfrac{\partial N_i}{\partial L_1}-\dfrac{\partial N_i}{\partial L_3})y_i\\ \displaystyle\sum(\dfrac{\partial N_i}{\partial L_2}-\dfrac{\partial N_i}{\partial L_3})x_i&\displaystyle\sum(\dfrac{\partial N_i}{\partial L_2}-\dfrac{\partial N_i}{\partial L_3})y_i \end{bmatrix}\\ \] 插值基函数对整体坐标的偏导数 \[ \begin{bmatrix} \dfrac{\partial N_i}{\partial \xi}\\ \dfrac{\partial N_i}{\partial \eta} \end{bmatrix}=\textbf{J}^{-1} \begin{bmatrix} \dfrac{\partial N_i}{\partial L_1}-\dfrac{\partial N_i}{\partial L_3}\\ \dfrac{\partial N_i}{\partial L_2}-\dfrac{\partial N_i}{\partial L_3} \end{bmatrix}\\ \] 积分微元关系 \[ \mathrm dx\mathrm dy=\det\textbf{J}\cdot\mathrm d\xi\mathrm d\eta\\ \]

三维情况

Jacobi矩阵 \[ \textbf{J}=\begin{bmatrix} \dfrac{\partial x}{\partial \xi}&\dfrac{\partial y}{\partial \xi}&\dfrac{\partial z}{\partial \xi}\\ \dfrac{\partial x}{\partial \eta}&\dfrac{\partial y}{\partial \eta}&\dfrac{\partial z}{\partial \eta}\\ \dfrac{\partial x}{\partial \zeta}&\dfrac{\partial y}{\partial \zeta}&\dfrac{\partial z}{\partial \zeta}\\ \end{bmatrix}= \begin{bmatrix} \displaystyle\sum\dfrac{\partial N_i}{\partial \xi_1}x_i&\displaystyle\sum\dfrac{\partial N_i}{\partial \xi_1}y_i&\displaystyle\sum\dfrac{\partial N_i}{\partial \xi_1}z_i\\ \displaystyle\sum\dfrac{\partial N_i}{\partial \eta_1}x_i&\displaystyle\sum\dfrac{\partial N_i}{\partial \eta_1}y_i&\displaystyle\sum\dfrac{\partial N_i}{\partial \eta_1}z_i\\ \displaystyle\sum\dfrac{\partial N_i}{\partial \zeta_1}x_i&\displaystyle\sum\dfrac{\partial N_i}{\partial \zeta_1}y_i&\displaystyle\sum\dfrac{\partial N_i}{\partial \zeta_1}z_i\\ \end{bmatrix}\\ \] 插值基函数对整体坐标的偏导数 \[ \begin{bmatrix} \dfrac{\partial N_i}{\partial x}\\ \dfrac{\partial N_i}{\partial y}\\ \dfrac{\partial N_i}{\partial z}\\ \end{bmatrix}=\textbf{J}^{-1} \begin{bmatrix} \dfrac{\partial N_i}{\partial \xi}\\ \dfrac{\partial N_i}{\partial \eta}\\ \dfrac{\partial N_i}{\partial \zeta}\\ \end{bmatrix}\\ \]

高斯积分

有限元方法经过坐标变换后,常出现如下形式的积分 \[ \displaystyle\int_{-1}^{1}f(\xi)\mathrm d\xi\qquad\displaystyle\int_{-1}^{1}\displaystyle\int_{-1}^{1}f(\xi,\eta)\mathrm d\xi\mathrm d\eta \qquad\displaystyle\int_{-1}^{1}\displaystyle\int_{-1}^{1}\displaystyle\int_{-1}^{1}f(\xi,\eta,\zeta)\mathrm d\xi\mathrm d\eta\mathrm d\zeta\\ \] 一般通过数值积分的方法算出积分的结果,最常用的是高斯积分,用少量的积分点函数值加权相加得到积分的结果,它具有较高的精度和稳定性

  • \(n\)点高斯积分具有\(2n-1\)阶代数精度
  • 高斯积分的权系数均为正,因而是稳定的
高斯点和权重

可以用勒让德多项式来构造高斯点和高斯权重系数,具体内容看数值分析的教材.设高斯权重为\(H\)

在整体坐标中的积分,经坐标变换转化为高斯积分 \[ \displaystyle\int_{x_1}^{x_2}f(x)\mathrm dx=\displaystyle\int_{-1}^{1}f(\xi)\frac{\mathrm dx}{\mathrm d\xi}\mathrm d\xi=\sum_{i=1}^nH_if(\xi_i)\\ \displaystyle\iint_D f(x,y)\mathrm dx\mathrm dy=\displaystyle\int_{-1}^{1}\displaystyle\int_{-1}^{1}f(\xi,\eta)\det{\textbf{J}}\mathrm d\xi\mathrm d\eta=\sum_{i=1}^n\sum_{j=1}^nH_iH_jf(\xi_i,\eta_j)\\ \displaystyle\iiint_\Omega f(x,y,z)\mathrm dx\mathrm dy\mathrm dz=\displaystyle\int_{-1}^{1}\displaystyle\int_{-1}^{1}\displaystyle\int_{-1}^{1}f(\xi,\eta,\zeta)\det{\textbf{J}}\mathrm d\xi\mathrm d\eta\mathrm d\zeta=\sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^nH_iH_jH_kf(\xi_i,\eta_j,\zeta_k)\\ \] 其中\(D\)是平面四边形区域,\(\Omega\)是空间六面体区域.

三角形上数值积分

设函数\(f(x,y)\)定义在三角形区域上,其在该区域上的积分通过坐标变换,转换成以面积坐标为积分变元的积分 \[ I=\displaystyle\iint_\Delta f(x,y)\mathrm dx\mathrm dy=\frac{1}{2}\displaystyle\int_0^1\displaystyle\int_0^{1-L_1}f(L_1,L_2,L_3)\det{\textbf{J}}\mathrm dL_1\mathrm dL_2\\ \] 这类积分用Hammer积分公式,与高斯积分公式类似,有积分点和权重系数

Hammer积分 \[ I=\frac{1}{2}\sum_{i=1}^{n}H_if(L_1^{(i)},L_3^{(i)},L_3^{(i)})\det\textbf{J}\\ \]

十四点积分

对于三点高斯积分在二维情况中一共需要计算9个点的函数值与权重的乘积,到了三维就需要计算27个点的函数值与权重的乘积,而两点高斯积分则只需要计算8个点,为了计算更高效同时精度不变,在三维情况中可以挑选14个点进行积分,下面给处这14个点与相应的权重. \[ \displaystyle\int_{-1}^{1}\displaystyle\int_{-1}^{1}\displaystyle\int_{-1}^{1}f(\xi,\eta,\zeta)\det{\textbf{J}}\mathrm d\xi\mathrm d\eta\mathrm d\zeta=B[f(-b,0,0)+f(b,0,0)+f(0,-b,0)+f(0,b,0)+f(0,0,-b)+f(0,0,b)]\\+C[f(-c,-c,-c)+f(c,c,c)+f(c,-c,-c)+f(-c,c-c)+f(-c,-c,c)+f(c,c,-c)+f(-c,c,c)+f(c,-c,c)]\\ 其中\quad B=\frac{320}{361},\quad C=\frac{121}{361},\quad b=\sqrt\frac{19}{30},\quad c=\sqrt\frac{19}{33}\\ \] 可以看出,这14个点中,6个点在每个面中心,8个点靠近各个顶点.(ps.全网都没找到这是怎么推出来的,有大佬知道的话麻烦告诉我一下,找了好久)

下一篇应该能到刚度矩阵,代码明天再补.

#! https://zhuanlan.zhihu.com/p/694799344