离散化1.2.2.2

现在市场上的有限元分析软件普遍都是基于位移法进行求解的其数学原理就是上节的最小势能原理变分提法。

回顾上节推导过程首先需要根据位移条件确定可能位移的范围其次根据假设的可能位移代入几何方程和本构方程得到基于位移的应力应变的表达式最后代入

\delta

δII0的方程求解相关待定参数。

线性有限元求解的一般格式基本上也是一样的除此之外线性有限元会多一些其他步骤比如结构的离散化、位移模式假定等。

下图为线性有限元求解的一般格式。

1.2.2.1

将结构进行离散化建模是有限元分析工作中的第一步某典型的结构如下图所示其中结构为一个圆盘将其离散成若干四边形网格单元如下图所示。

结构离散化合适与否一定程度上决定使用有限元分析方法解决工程问题的精度。

1.2.2.2

在最小势能原理中首先需要根据位移条件确定可能位移的范围但是在通常的分析中我们为了方便处理往往将可能位移用多项式来逼近这个过程就是单元的位移模式选择。

我们以一阶二维三角形单元为例二维三角形单元有6各位移自由度。

用多项式拟合可能位移如下式所示

(1-54)

u(x,y)a0​a1​xa2​yv(x,y)b0​b1​xb2​y​(1-54)

b形函数

多项式位移模式有六个待定参数二维三角形单元有六个节点位移变量所以六个待定参数可以表示成六个节点位移变量的形式如下式待定系数满足下式。

(1-55)

\end{bmatrix}^{-1}\frac{[M_{ij}\cdot

M_{ij}

(−1)ij称为代数余子式矩阵的逆是其代数余子式组成的伴随矩阵除以该矩阵的行列式。

(1-58)

​a1​b1​c1​​a2​b2​c2​​a3​b3​c3​​

​−1⋅

a_0\frac{1}{A}(a_1u_1a_2u_2a_3u_3)\\

\overline

a_1\frac{1}{A}(b_1u_1b_2u_2b_3u_3)\\

\overline

a_2\frac{1}{A}(c_1u_1c_2u_2c_3u_3)

a0​A1​(a1​u1​a2​u2​a3​u3​)a1​A1​(b1​u1​b2​u2​b3​u3​)a2​A1​(c1​u1​c2​u2​c3​u3​)​(1-59)

(1-60)

\frac{1}{A}(a_1u_1a_2u_2a_3u_3)\frac{1}{A}(b_1u_1b_2u_2b_3u_3)x\frac{1}{A}(c_1u_1c_2u_2c_3u_3)y\\

\frac{1}{A}(a_1b_1xc_1y)u_1\frac{1}{A}(a_2b_2xc_2y)u_2\frac{1}{A}(a_3b_3xc_3y)u_3\\

N_1(x,y)u_1N_2(x,y)u_2N_3(x,y)u_3

u(x,y)​a0​a1​xa2​yA1​(a1​u1​a2​u2​a3​u3​)A1​(b1​u1​b2​u2​b3​u3​)xA1​(c1​u1​c2​u2​c3​u3​)yA1​(a1​b1​xc1​y)u1​A1​(a2​b2​xc2​y)u2​A1​(a3​b3​xc3​y)u3​N1​(x,y)u1​N2​(x,y)u2​N3​(x,y)u3​​(1-60)

(1-61)

\frac{1}{A}(a_1v_1a_2v_2a_3v_3)\frac{1}{A}(b_1v_1b_2v_2b_3v_3)x\frac{1}{A}(c_1v_1c_2v_2c_3v_3)y\\

\frac{1}{A}(a_1b_1xc_1y)v_1\frac{1}{A}(a_2b_2xc_2y)v_2\frac{1}{A}(a_3b_3xc_3y)v_3\\

N_1(x,y)v_1N_2(x,y)v_2N_3(x,y)v_3

v(x,y)​a0​a1​xa2​yA1​(a1​v1​a2​v2​a3​v3​)A1​(b1​v1​b2​v2​b3​v3​)xA1​(c1​v1​c2​v2​c3​v3​)yA1​(a1​b1​xc1​y)v1​A1​(a2​b2​xc2​y)v2​A1​(a3​b3​xc3​y)v3​N1​(x,y)v1​N2​(x,y)v2​N3​(x,y)v3​​(1-61)

(1-62)

u(x,y)[u(x,y)v(x,y)​]​[N1​(x,y)0​0N1​(x,y)​N2​(x,y)0​0N2​(x,y)​N3​(x,y)0​0N3​(x,y)​]

​u1​v1​u2​v2​u3​v3​​

应变和位移的关系由几何方程确定在本文二维的例子中应变仅有三个分量那么应变与位移的关系可以如下式所示。

(1-63)

应力应变关系又成本构方程在线性范围内应力与应变的关系写成矩阵如下(本例中为二维三角形单元应力分量只有三个)

(1-64)

\frac{E}{1-\mu^2}\begin{bmatrix}

\mu

\boldsymbol{\varepsilon}(x,y)\\

\mathbf

​D(x,y)⋅ε(x,y)D(x,y)⋅B(x,y)⋅qeS(x,y)⋅qe​(1-64)

1.2.2.5

II\frac{1}{2}\int_{\Omega}\sigma_{ij}\varepsilon_{ij}d\Omega-\int_{\Omega}

b_{i}u_{i}d\Omega

\frac{1}{2}\int_{\Omega}\boldsymbol{\sigma}^T\boldsymbol{\varepsilon}\mathbf

\mathbf

II​21​∫Ω​σij​εij​dΩ−∫Ω​bi​ui​dΩ−∫A​pi​ui​dA21​∫Ω​σTεdΩ−∫Ω​bTudΩ−∫A​pTudA​(1-65)

将前述单元的应力、应变矩阵等代入上式建立单元的总势能表达式如下所示

(1-66)

\frac{1}{2}\int_{\Omega}\boldsymbol{\sigma}^T\boldsymbol{\varepsilon}\mathbf

\mathbf

\frac{1}{2}\int_{\Omega}\boldsymbol{q^e}^T\boldsymbol{B^TD^TBq^e}\mathbf

\mathbf

\frac{1}{2}\boldsymbol{q^e}^T(\int_{\Omega}\boldsymbol{B^TD^TB

d}\Omega)\boldsymbol{q^e}-(\int_{\Omega}

\mathbf

\frac{1}{2}\boldsymbol{q^e}^T\boldsymbol{K^e}\boldsymbol{q^e}-\boldsymbol{P^e}^T\boldsymbol{q^e}

II​21​∫Ω​σTεdΩ−∫Ω​bTudΩ−∫A​pTudA21​∫Ω​qeTBTDTBqedΩ−∫Ω​bTNqedΩ−∫A​pTNqedA21​qeT(∫Ω​BTDTBdΩ)qe−(∫Ω​bTNdΩ∫A​pTNdA)qe21​qeTKeqe−PeTqe​(1-66)

回顾之前最小势能原理真实位移是使得势能取得最小值上式是单元的总势能那么真实的位移使得势能取到最小值也就是说总势能一阶变分为零即

(1-67)

q^e(\boldsymbol{K^e}\boldsymbol{q^e}-\boldsymbol{P^e})\delta

q^e0

δII∂qe∂II​δqe(Keqe−Pe)δqe0(1-67)

\delta

\boldsymbol{K^e}\boldsymbol{q^e}-\boldsymbol{P^e}0\tag{1-68}

Keqe−Pe0(1-68)

对于物体来说一般有多个单元那么就可以建立多个单元平衡方程如下所示

(1-69)

\boldsymbol{K^e_{(i)}}\boldsymbol{q^e_{(i)}}-\boldsymbol{P^e_{(i)}}0\tag{1-69}

按照顺序将所有的单元的节点位移组成总体节点位移列阵相应的节点等效载荷和单元刚度矩阵组成总体节点等效载荷列阵和整体刚度矩阵以下图为例下图总共有两个单元得到两个单元平衡方程

(1-70)

\end{bmatrix}\end{aligned}\tag{1-70}

​k111​k211​k311​k411​k511​k611​​k121​k221​k321​k421​k521​k621​​k131​k231​k331​k431​k531​k631​​k141​k241​k341​k441​k541​k641​​k151​k251​k351​k451​k551​k651​​k161​k261​k361​k461​k561​k661​​

​u1​v1​u2​v2​u3​v3​​

​px11​py11​px21​py21​px31​py31​​

​k112​k212​k312​k412​k512​k612​​k122​k222​k322​k422​k522​k622​​k132​k232​k332​k432​k532​k632​​k142​k242​k342​k442​k542​k642​​k152​k252​k352​k452​k552​k652​​k162​k262​k362​k462​k562​k662​​

​u2​v2​u3​v3​u4​v4​​

​px22​py22​px32​py32​px42​py42​​

​​(1-70)

p^1_{y2}p^2_{y2}\\p^1_{x3}p^2_{x3}\\

p^2_{y4}\\

​k111​k211​k311​k411​k511​k611​00​k121​k221​k321​k421​k521​k621​00​k131​k231​k331​k112​k431​k212​k531​k312​k631​k412​k512​k612​​k141​k241​k341​k122​k441​k222​k541​k322​k641​k422​k522​k622​​k151​k251​k351​k132​k451​k232​k551​k332​k651​k432​k532​k632​​k161​k261​k361​k142​k461​k242​k561​k342​k661​k442​k542​k642​​00k152​k252​k352​k452​k552​k652​​00k162​k262​k362​k462​k562​k662​​

​px11​py11​px21​px22​py21​py22​px31​px32​py31​py32​px42​py42​​

​(1-71)

p^1_{y2}p^2_{y2}\\p^1_{x3}p^2_{x3}\\

p^2_{y4}\\

​k111​k211​k311​k411​k511​k611​00​k121​k221​k321​k421​k521​k621​00​k131​k231​k331​k112​k431​k212​k531​k312​k631​k412​k512​k612​​k141​k241​k341​k122​k441​k222​k541​k322​k641​k422​k522​k622​​k151​k251​k351​k132​k451​k232​k551​k332​k651​k432​k532​k632​​k161​k261​k361​k142​k461​k242​k561​k342​k661​k442​k542​k642​​00k152​k252​k352​k452​k552​k652​​00k162​k262​k362​k462​k562​k662​​

​Rx1​Ry1​px21​px22​py21​py22​px31​px32​py31​py32​px42​py42​​

​(1-72)

p^1_{y2}p^2_{y2}\\p^1_{x3}p^2_{x3}\\

p^2_{y4}\\

​k331​k112​k431​k212​k531​k312​k631​k412​k512​k612​​k341​k122​k441​k222​k541​k322​k641​k422​k522​k622​​k351​k132​k451​k232​k551​k332​k651​k432​k532​k632​​k361​k142​k461​k242​k561​k342​k661​k442​k542​k642​​k152​k252​k352​k452​k552​k652​​k162​k262​k362​k462​k562​k662​​

​u2​v2​u3​v3​u4​v4​​

​px21​px22​py21​py22​px31​px32​py31​py32​px42​py42​​

​P​(1-73)

[k131​k231​​k141​k241​​k151​k251​​k161​k261​​]

​u2​v2​u3​v3​​

qe后通过一系列的回代得到单元内任意位置的位移、应变、应力矩阵而式1-74用来求解约束支反力。

(1-75)

u(x,y)N(x,y)⋅qeε(x,y)B(x,y)⋅qeσ(x,y)S(x,y)⋅qe​(1-75)