有限元方法

绪论

有限单元法

有限元法：基于近代计算机快算发展而发展出的求解偏微分方程边值问题近似解的数值技术

  • 核心思想：离散化+数值近似

固体力学控制微分方程

  • 梁弯曲问题, 挠度ω(x)\omega(x)

    • EJd4ω(x)dx4−q=0EJ\frac{\mathrm{d}^4\omega(x)}{\mathrm{d}x^4}-q=0

    • 悬臂梁边界条件, 固定端挠度转角为零，自由端弯矩剪力为零

    • {x=0: ω=0,ω′(x)=0x=l: ω′′(x)=0,ω′′′(x)=0\begin{cases} x=0: \ \omega=0,\omega\prime(x)=0 \\ x=l: \ \omega\prime\prime(x)=0,\omega\prime\prime\prime(x)=0 \end{cases}

  • 薄板弯曲问题, 挠度

    • D∇4ω(x,y)=q(x,y)D \nabla^4 \omega(x,y)=q(x,y)

    • 其中 D=Eh3/12(1−u2)D=E h^3 / 12(1-u^2)
    • ∇4=∂4∂x4+∂4∂x2∂y2+∂4∂y4\nabla^4=\frac{\partial^4}{\partial x^4}+\frac{\partial^4}{\partial x^2 \partial y^2}+\frac{\partial^4}{\partial y^4}
    • 薄板四边简支边界条件，挠度弯矩为零

      • {x=0, a: w=0,wx′′(x,y)=0y=0, b: w=0,wy′′(x,y)=0\begin{cases} x=0,~a:~w=0,w_x^{''}(x,y)=0 \\ y=0,~b:~w=0,w_y^{''}(x,y)=0 \end{cases}

弹性力学三维问题

几何方程: 位移应变关系

正应变-uvw/xyz, 各回各家各找各妈

切应变-yz, xz, xy, uvw交换身体再相加

εx=∂u∂x\varepsilon_x=\frac{\partial u}{\partial x} εy=∂v∂y\varepsilon_y=\frac{\partial v}{\partial y} εz=∂w∂z\varepsilon_z=\frac{\partial w}{\partial z}
γyz=∂w∂y+∂v∂z\gamma_{yz} = \frac{\partial w}{\partial y}+\frac{\partial v}{\partial z} γxz=∂u∂z+∂w∂x\gamma_{xz} = \frac{\partial u}{\partial z}+\frac{\partial w}{\partial x} γxy=∂u∂y+∂v∂x\gamma_{xy} = \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}

物理方程: σ=Dε\sigma = D\varepsilon, 或ε=D−1σ\varepsilon=D^{-1}\sigma

D−1=1E[1−u−u000−u1−u000−u−u10000002(1+μ)0000002(1+μ)0000002(1+μ)]D^{-1}=\frac{1}{E}\begin{bmatrix} &1 &-u &-u &0 &0 &0 \\ &-u &1 &-u &0 &0 &0 \\ &-u &-u &1 &0 &0 &0 \\ &0 &0 &0 &2(1+\mu) &0 &0 \\ &0 &0 &0 &0 &2(1+\mu) &0 \\ &0 &0 &0 &0 &0 &2(1+\mu) \\ \end{bmatrix}

平衡方程: σij,j+Fi=0\sigma_{ij,j}+F_i = 0

∂σx∂x+∂τyx∂y+∂τzx∂z+Fbx=0∂σy∂y+∂τxy∂x+∂τzy∂z+Fby=0∂σz∂z+∂τxz∂x+∂τyz∂y+Fbz=0\begin{aligned} \frac{\partial \sigma_x}{\partial x}+\frac{\partial \tau_{yx}}{\partial y}+\frac{\partial \tau_{zx}}{\partial z}+F_{bx}&=0 \\ \frac{\partial \sigma_y}{\partial y}+\frac{\partial \tau_{xy}}{\partial x}+\frac{\partial \tau_{zy}}{\partial z}+F_{by}&=0 \\ \frac{\partial \sigma_z}{\partial z}+\frac{\partial \tau_{xz}}{\partial x}+\frac{\partial \tau_{yz}}{\partial y}+F_{bz}&=0 \end{aligned}

偏微分方程的数值解法

  1. 有限元方法（FEM - finite Element Method）
  2. 边界元方法（BEM - Boundary Element Method）
  3. 加权残值方法（MWR - Method of Weighted Residuals）
  4. 有限差分法（FDM - Finite Differential Method）
  5. 无网格法（Mesh Free Method）

有限元基础

应变分析

  • 位移矢量v0=vv=li+mj+nk\bold v^0=\frac{\bold v}{v}=l\bold i+m\bold j+n\bold k

    • l=cos(v,i) , m=cos(v,j) , n=cos(v,k)

  • 在外力下，P点 - > P’点， P点的位移矢

    • u={u v w }T(x,y,z)u=\{ u~v~w~\}^T(x,y,z)

  • 应变六分量

    • ε={εx , εy , εz , γyz , γzx , γxy }T\varepsilon=\{ \varepsilon_x \ , \ \varepsilon_y \ , \ \varepsilon_z \ , \ \gamma_{yz} \ , \ \gamma_{zx} \ , \ \gamma_{xy}\ \}^T

  • 几何方程

    εx=∂u∂x\varepsilon_x=\frac{\partial u}{\partial x} εy=∂v∂y\varepsilon_y=\frac{\partial v}{\partial y} εz=∂w∂z\varepsilon_z=\frac{\partial w}{\partial z}
    γyz=∂w∂y+∂v∂z\gamma_{yz} = \frac{\partial w}{\partial y}+\frac{\partial v}{\partial z} γxz=∂u∂z+∂w∂x\gamma_{xz} = \frac{\partial u}{\partial z}+\frac{\partial w}{\partial x} γxy=∂u∂y+∂v∂x\gamma_{xy} = \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}

应变

引入算子

E(∇)=[∂∂x000∂∂z∂∂y0∂∂y0∂∂z0∂∂x00∂∂z∂∂y∂∂x0]E(\nabla) = \begin{bmatrix} &\frac{\partial}{\partial x} &0 &0 &0 &\frac{\partial}{\partial z} &\frac{\partial}{\partial y} \\ &0 &\frac{\partial}{\partial y} &0 &\frac{\partial}{\partial z} &0 &\frac{\partial}{\partial x} \\ &0 &0 &\frac{\partial}{\partial z} &\frac{\partial}{\partial y} &\frac{\partial}{\partial x} &0 \\ \end{bmatrix}

梯度矢量 ∇=i∂∂x+j∂∂y+k∂∂z=(∂∂x ∂∂y ∂∂z)T\nabla = \mathbf i \frac{\partial}{\partial x}+\mathbf j \frac{\partial}{\partial y}+\mathbf k \frac{\partial}{\partial z} = (\frac{\partial}{\partial x} ~ \frac{\partial}{\partial y}~ \frac{\partial}{\partial z})^T

位移列阵u={u v w}Tu=\{ u \ v \ w \} ^T

应变=算子的转置矩阵 x 位移列阵 ， 即ε=ET(∇)u\varepsilon = E^T(\nabla)u

应力

σ={σx σy σz τyz τzx τxy}T\sigma = \{\sigma_x ~\sigma_y ~ \sigma_z ~\tau_{yz} ~\tau_{zx} ~\tau_{xy}\}^T
体力(外力): f={fx fy fz}T\mathbf f = \{f_x ~ f_y ~ f_z\}^T

平衡方程

E(∇)σ+f=0E(\nabla)\sigma+f=0

{∂σx∂x+∂τyx∂y+∂τzx∂z+fx=0∂σy∂y+∂τxy∂x+∂τzy∂z+fy=0∂σz∂z+∂τxz∂x+∂τyz∂y+fz=0\begin{cases} \frac{\partial \sigma_x}{\partial x}+\frac{\partial \tau_{yx}}{\partial y}+\frac{\partial \tau_{zx}}{\partial z}+f_x&=0 \\ \frac{\partial \sigma_y}{\partial y}+\frac{\partial \tau_{xy}}{\partial x}+\frac{\partial \tau_{zy}}{\partial z}+f_y&=0 \\ \frac{\partial \sigma_z}{\partial z}+\frac{\partial \tau_{xz}}{\partial x}+\frac{\partial \tau_{yz}}{\partial y}+f_z&=0 \end{cases}

应力边界条件

Pˉ=E(v)σ\bar P = E(\bold v)\sigma

一个面上应力可分解为一个正应力, 两个剪应力分量

  • 物体表面的面力矢量: Pˉ={Pˉx Pˉy Pˉz}\bar P = \{ \bar P_x ~ \bar P_y ~ \bar P_z \}
  • 此表面处外法线方向: v=(l , m , n)\mathbf v = (\mathbf l ~,~\mathbf m ~,~ \mathbf n)

平衡方程

{Pˉx=σxl+τxym+τxznPˉx=τyzl+σym+τyznPˉx=τzxl+τzym+σzn\begin{cases} \bar P_x = \sigma_x l + \tau_{xy}m+\tau_{xz}n \\ \bar P_x = \tau_{yz} l + \sigma_y m+\tau_{yz}n \\ \bar P_x = \tau_{zx} l + \tau_{zy}m+\sigma_{z}n \end{cases}

本构方程

σ=Dε\sigma = D \varepsilon

  • σ={σx σy σz τyz τzx τxy}T\sigma = \{\sigma_x ~\sigma_y ~ \sigma_z ~\tau_{yz} ~\tau_{zx} ~\tau_{xy}\}^T
  • ε={εx εy εz γyz γzx γxy}T\varepsilon=\{\varepsilon_x~\varepsilon_y~ \varepsilon_z ~\gamma_{yz} ~ \gamma_{zx} ~ \gamma_{xy}\}^T

弹性矩阵D

D=E(1−μ)(1+u)(1−2μ)[1μ1−μμ1−μ000μ1−μ1μ1−μ000μ1−μμ1−μ10000001−2μ2(1−μ)0000001−2μ2(1−μ)0000001−2μ2(1−μ)]D = \frac{E(1-\mu)}{(1+u)(1-2\mu)} \begin{bmatrix} &1 &\frac{\mu}{1-\mu} &\frac{\mu}{1-\mu} &0 &0 &0 \\ &\frac{\mu}{1-\mu} &1 &\frac{\mu}{1-\mu} &0 &0 &0 \\ &\frac{\mu}{1-\mu} &\frac{\mu}{1-\mu} &1 &0 &0 &0 \\ &0 &0 &0 &\frac{1-2\mu}{2({1-\mu)}} &0 &0 \\ &0 &0 &0 &0 &\frac{1-2\mu}{2({1-\mu)}} &0 \\ &0 &0 &0 &0 &0 &\frac{1-2\mu}{2({1-\mu)}} \\ \end{bmatrix}

或

D−1=1E[1−u−u000−u1−u000−u−u10000002(1+μ)0000002(1+μ)0000002(1+μ)]D^{-1}=\frac{1}{E} \begin{bmatrix} &1 &-u &-u &0 &0 &0 \\ &-u &1 &-u &0 &0 &0 \\ &-u &-u &1 &0 &0 &0 \\ &0 &0 &0 &2(1+\mu) &0 &0 \\ &0 &0 &0 &0 &2(1+\mu) &0 \\ &0 &0 &0 &0 &0 &2(1+\mu) \\ \end{bmatrix}