静力学FEM12.30

1.静力学方程

考虑图所示变截面弹性杆的静态响应。这是线性应力分析或线弹性问题的一个例子，我们需要求杆内的应力分布σ(x)。

应力由物体的变形产生，而变形由物体内各点的位移u(x)表征。位移导致用ε(x)表示的应变；应变是一个无量纲变量。杆受到分布力b(x)或集中力作用。这个力可能由重力（如果杆垂直放置而非水平放置）、磁力或热应力引起；在一维情况下，我们将考虑单位长度上的体力，所以b(x)的单位是力/长度。此外，载荷可以规定在杆端，在这些位置位移未被规定；这些载荷被称为牵引力，用t表示。这些载荷的单位是力/面积，当乘以面积时，就得到所施加的力。

2.加权余量法

静力学方程

\[\frac{\mathrm{d}}{\mathrm{d}x}\left(AE\frac{\mathrm{d}u}{\mathrm{d}x}\right)+b(x)=0,\quad0<x<l \]

\[\sigma\left(0\right)=\left(E\frac{\mathrm{d}u}{\mathrm{d}x}\right)=-\bar{t}\\ \]

\[u\left(l\right)=\bar{u} \]

积分弱形式

\[\int_0^lw\left[\frac{\mathrm{d}}{\mathrm{d}x}\left(AE\frac{\mathrm{d}u}{\mathrm{d}x}\right)+b\right]\mathrm{d}x=0, \quad\forall w\quad \mathrm{(a)} \]

\[\left[wA\left(E\frac{\mathrm{d}u}{\mathrm{d}x}+\bar{t}\right)\right]_{x=0}=0,\quad\forall w\mathrm{(b)} \]

（4）整理成

\[\int_0^lw\frac{\mathrm{d}}{\mathrm{d}x}\left(AE\frac{\mathrm{d}u}{\mathrm{d}x}\right)\mathrm{d}x+\int_0^lwb\mathrm{d}x=0,\quad\forall w \]

分布积分

\[\int_{0}^{l} w \frac{d}{dx} \left( AE \frac{du}{dx} \right) dx = \left[ w \left( AE \frac{du}{dx} \right) \right]_{0}^{l} - \int_{0}^{l} \left( AE \frac{du}{dx} \right) \frac{dw}{dx} dx \]

（7 ）代入（6）

\[\left.\left(wAE\frac{\mathrm{d}u}{\mathrm{d}x}\right)\right|_0^l-\int_0^l\frac{\mathrm{d}w}{\mathrm{d}x}AE\frac{\mathrm{d}u}{\mathrm{d}x}\mathrm{d}x+\int_0^lwb\mathrm{d}x=0,\quad\forall w\text{且}w(l)=0 \]

把 \(\sigma=E\frac{\mathrm{d}u}{\mathrm{d}x}\)代入（8）

\[(wA\sigma)_{x=l}-(wA\sigma)_{x=0}-\int_0^l\frac{\mathrm{d}w}{\mathrm{d}x}AE\frac{\mathrm{d}u}{\mathrm{d}x}\mathrm{d}x+\int_0^lwb\mathrm{d}x=0,\quad\forall w\text{且}w(l)=0 \]

因为 $ w(l) = 0 $，且由 (2）式得到 $ (wA\sigma){x=0} = -(wAt) $(应力边界条件)，将这些条件代入 (8) 式，得

\[\int_0^l\frac{\mathrm{d}w}{\mathrm{d}x}AE\frac{\mathrm{d}u}{\mathrm{d}x}\mathrm{d}x=\int_0^lwb\mathrm{d}x+\left(wA\bar{t}\right)_{x=0},\quad\forall w\text{且}w(l)=0 \]

应力边界条件自然满足， 自然边界条件。

3.伽辽金法

伽辽金

将图 中的杆件划分为 \(n_e\) 个杆单元。那么，在整个域内权函数插值为

\[w(x) = N(x)w \]

其中

\[N(x) = [N_1(x) \quad N_2(x) \quad \cdots \quad N_{n_e}(x)] \]

\[w = [w_1 \quad w_2 \quad \cdots \quad w_{n_e}]^T \]

场函数的域内插值为

\[u(x) = N(x)u \]

其中

\[u = [u_1 \quad u_2 \quad \cdots \quad u_n]^T \]

场函数必须满足强制边界条件，即令\(x = l\)的节点位移

\[u_1=\bar{u}_1 \]

其余节点的位移是未知的，将由接下来的伽辽金弱形式求解。

在强制边界条件上，要求\(w(l)=0\)，所以须令

\[w_1 = 0 \]

其余节点的权函数是任意的。

单元和总体的矩阵是由集成矩阵关联的，如(2.18)式所示，可得到单元的节点信息，即

\[w^e = L^e w, \quad u^e = L^e u \]

所以，单元的权函数与场函数插值表达式为

\[w^e(x) = N^e(x) w^e \]

\[u^e(x) = N^e(x) u^e \]

那么，在整个杆件域\([0, l]\)上的积分，转换为在每个单元域积分的叠加。进而，将(5.9)式代入(3.20)式，得

\[\sum_{e = 1}^{n_e} \int_0^l \left( \frac{d w^e}{d x} \right)^T A E \left( \frac{d u^e}{d x} \right) d x - \int_0^l \left( w^e \right)^T b d x - \left( w^e \right)^T A F \left. \frac{d u^e}{d x} \right|_{x = 0} ^{x = l} = 0 \]

其中

\[\frac{d w^e}{d x} = B^e w^e \]

\[\frac{d u^e}{d x} = B^e u^e \]

将(5.11)式代入(5.10)式，得到

\[\sum_{e = 1}^{n_e} \left( w^e \right)^T \left( \int_0^l \left( B^e \right)^T A E B^e u^e d x - \int_0^l \left( N^e \right)^T b d x - \left( \left( N^e \right)^T A F \left. \frac{d u^e}{d x} \right|_{x = 0} ^{x = l} \right) \right) = 0 \]

其中，第一项\(u^e\)不含积分变量，可移到积分符号外，即

\[\int_0^l \left( B^e \right)^T A E B^e u^e d x = \int_0^l \left( B^e \right)^T A E B^e d x u^e \]

那么(5.12)式可写为

\[\left( w^e \right)^T \left( \int_0^l \left( B^e \right)^T A E B^e d x u^e - \int_0^l \left( N^e \right)^T b d x - \left( \left( N^e \right)^T A F \left. \frac{d u^e}{d x} \right|_{x = 0} ^{x = l} \right) \right) = 0 \]

对(5.14)式，定义两个非常重要的矩阵：

(1) 单元刚度矩阵

\[K^e = \int_0^{l^e} (B^e)^T A E B^e dx = \int_0^{l^e} (B^e)^T A E B^e dx \]

(2) 单元等效节点力列阵

\[f^e = \int_0^{l^e} (N^e)^T b dx + \left. [(N^e)^T A F] \right|_{x = 0} ^{x = l^e} = \int_0^{l^e} (N^e)^T b dx + \left. [(N^e)^T A F] \right|_{x = 0} ^{x = l^e} \]

将(5.15)式与(5.16)式代入(5.14)式，并使用(5.8)式，得

\[w^T \left( \sum_{e = 1}^{n_e} (L^e)^T K^e L^e u - \sum_{e = 1}^{n_e} (L^e)^T f^e \right) = 0 \]

在(5.17)式的推导中，\(w\)与\(x\)无关，所以可将其提到求和符号的外边。于是，得到了杆件的总体刚度矩阵与载荷列阵

\[K = \sum_{e = 1}^{n_e} (L^e)^T K^e L^e \]

\[f = \sum_{e = 1}^{n_e} (L^e)^T f^e \]

将(5.18)式代入(5.17)式，得到

\[w^T (K u - f) = 0, \quad \forall w \text{除了} w_1 = 0 \]

为了求解(5.19)式，令

\[r = K u - f \]

其中，\(r\)称为不平衡力列阵。(5.19)式可写为

\[w^T r = 0, \quad \forall w \text{除了} w_1 = 0 \]

针对图5 - 1所示的结构，(5.21)式可以展开为

\[w_{2} r_{2} + w_{3} r_{3} + \cdots + w_{n} r_{n} = 0 \]

由于\(w\)的任意性，可知\(r_{2} = r_{3} = \cdots = r_{n} = 0\)。此时\(r_1\)未知，事实上，\(r_1\)是1节点处的不平衡力，也就是所谓的约束反力。于是，(5.20)式可写为

\[r = \begin{bmatrix} r_1 \\ 0 \\ \vdots \\ 0 \end{bmatrix} = \begin{bmatrix} k_{11} & k_{12} & \cdots & k_{1n} \\ k_{21} & k_{22} & \cdots & k_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ k_{n1} & k_{n2} & \cdots & k_{nn} \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ u_n \end{bmatrix} - \begin{bmatrix} f_1 \\ f_2 \\ \vdots \\ f_n \end{bmatrix} \]

重新整理，得

\[\begin{bmatrix} k_{11} & k_{12} & \cdots & k_{1n} \\ k_{21} & k_{22} & \cdots & k_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ k_{n1} & k_{n2} & \cdots & k_{nn} \end{bmatrix} \begin{bmatrix} \bar{u}_1 \\ u_2 \\ \vdots \\ u_n \end{bmatrix}= \begin{bmatrix} f_1 + r_1 \\ f_2 \\ \vdots \\ f_n \end{bmatrix} \]

(5.24)式未知数为\(\bar{u}_1, u_2, \cdots, u_n\)与\(r_1\)，共\(n\)个，因此(5.24)式可解。至此，强制边界条件对场函数和权函数的附加要求即(5.6)式与(5.7)式，同时在(5.24)式中得到满足。

4.形函数构造

节点线性单元

对于杆单元，为了满足完备性要求，至少要选取一次线性多项式，即

\[\phi^e(x)=\alpha_0+\alpha_1 x \]

在1和2节点分别满足

\[\phi^e(x_1^e)=\phi_1^e x \quad \text{与} \quad \phi^e(x_2^e)=\phi_2^e x \]

此时得到两个线性方程，恰好对应(4.2)式中同样数量的待定系数\(\alpha_0\)与\(\alpha_1\)，因此这些待定系数可由(4.3)式解得。

写成矩阵形式

\[\phi^e(x)=\left[1 \quad x\right]\left[\begin{array}{l}\alpha_0 \\ \alpha_1\end{array}\right]=p(x) \alpha^e \]

采用节点值表示待定系数

\[\begin{array}{l} \phi^e(x_1^e)=\phi_1^e x=\alpha_0+\alpha_1 x_1^e \rightarrow\left[\begin{array}{l} 1 \\ x_1^e \end{array}\right]\left[\begin{array}{l} \alpha_0 \\ \alpha_1 \end{array}\right] \\ \phi^e(x_2^e)=\phi_2^e x=\alpha_0+\alpha_1 x_2^e \rightarrow\left[\begin{array}{l} 1 \\ x_2^e \end{array}\right]\left[\begin{array}{l} \alpha_0 \\ \alpha_1 \end{array}\right] \end{array} \]

由(4.5)式解出

\[\alpha^e=(A^e)^{-1} \phi^e \]

其中，\((A^e)^{-1}\)是\(A^e\)的逆矩阵，它可按下式计算得到

\[(A^e)^{-1}=\frac{1}{\left|A^e\right|}\left(A^e\right)^* \]

\(\left|A^e\right|\)是\(A^e\)的行列式。

\((A^e)^*\)是\(A^e\)的伴随矩阵，它的元素\((A_{ij}^e)^*\)是\(A^e\)的元素\(A_{ij}^e\)的代数余子式。那么，对于(4.5)式中的\(A^e\)，其逆矩阵为

\[(A^e)^{-1}=\frac{1}{l^e}\left[\begin{array}{cc} x_2^e & -x_1^e \\ -1 & 1 \end{array}\right]=\frac{1}{l^e}\left[\begin{array}{cc} l^e & -x_1^e \\ -1 & 1 \end{array}\right] \]

其中，\(l^e\)为杆单元长度。

将(4.6)式代入(4.4)式，得

\[\phi^e(x)=N^e(x) \phi^e \]

文中首先定义了单元形函数矩阵\(N^e(x)=[N_1^e(x) \quad N_2^e(x)] = p(x)(A^e)^{-1}\)，并通过将\((4.8)\)式代入\((4.9)\)式，得到：

\[N^e = [N_1^e(x) \quad N_2^e(x)] = p(x)(A^e)^{-1} = \left[1 \quad x\right]\left[\begin{array}{cc} \frac{x_2^e}{l^e} & -\frac{x_1^e}{l^e} \\ -\frac{1}{l^e} & \frac{1}{l^e} \end{array}\right] = \left[\frac{l^e - x}{l^e} \quad \frac{x}{l^e}\right] \]

在\((4.10)\)式中，\(N_1^e(x)\)与\(N_2^e(x)\)分别为单元在1节点与2节点的形函数，且这些形函数是线性函数，具有以下属性：

\[N_1^e(x_1^e) = 1, \quad N_1^e(x_2^e) = 0 \]

\[N_2^e(x_1^e) = 0, \quad N_2^e(x_2^e) = 1 \]

简记为：

\[N_i^e(x_j^e) = \delta_{ij} \]

其中，\(\delta_{ij}\)为Kronecker符号，即：

\[\delta_{ij} = \begin{cases} 1, & i = j \\ 0, & i \neq j \end{cases} \]

在第3章中的“弱”形式中，需要求导函数与场函数的导数，如下：

\[\frac{d\phi^e}{dx} = \frac{d(N^e \phi^e)}{dx} = \frac{dN^e}{dx} \phi^e = \frac{dN_1^e}{dx} \phi_1^e + \frac{dN_2^e}{dx} \phi_2^e \]

写成矩阵形式：

\[\frac{d\phi^e}{dx} = \left[\begin{array}{ll} \frac{dN_1^e}{dx} & \frac{dN_2^e}{dx} \end{array}\right]\left[\begin{array}{l} \phi_1^e \\ \phi_2^e \end{array}\right] = B^e \phi^e \]

其中，\(B^e\)称为单元应变矩阵，即：

\[B^e = \left[\begin{array}{ll} \frac{dN_1^e}{dx} & \frac{dN_2^e}{dx} \end{array}\right] = \left[\begin{array}{ll} -\frac{1}{l^e} & \frac{1}{l^e} \end{array}\right] \]

5.等参变换

(4.4)式与(4.17)式的位移插值函数构造于杆单元的全局坐标系中。有限元积分格式的推导需要在单元的域内进行积分运算，但是对于三角形、四边形或者实体单元，由于被积函数的项数太多，不宜整理成显式，往往采用数值积分。数值积分的方法有Gauss积分、Irons积分与Hammer积分等。为了数值积分公式表达的方便，这些数值积分方法一般都在规则的积分域中给出积分点和积分权系数。因此，单元积分格式的推导也应该在规则的积分域中进行，这个新引入的规则积分域被称为母单元。那么，除了母单元域内位移与节点位移的位移插值关系以外，还需要引入杆单元局部坐标系与规则母单元自然坐标系的插值关系，也即坐标变换关系。

在杆单元的母单元中，坐标的变换范围为节点1处的 -1到节点2处的 +1，如图7 - 1所示。坐标变换的插值表达式可待定

\[x = \alpha_0+\alpha_1 \xi \]

将\((\xi=-1, x = x_1)\)与\((\xi=+1, x = x_2)\)两点代入(7.38)式，可得到插值函数

\[N_1^{\xi} = \frac{1}{2}(1-\xi), \quad N_2^{\xi} = \frac{1}{2}(1+\xi) \]

那么，(7.38)式可写为

\[x(\xi) = [N_1^{\xi}(\xi) \quad N_2^{\xi}(\xi)] \left[\begin{array}{c}x_1 \\ x_2\end{array}\right] = N^{\xi}(\xi) x^e \]

如果将(4.2)式的单元位移插值函数也在母单元内插值，即

\[u^e(\xi) = \alpha_0+\alpha_1 \xi \]

类似的有\((\xi=-1, u = u_1)\)和\((\xi=-1, u = u_2)\)。于是，(7.41)式可写为

\[u^e = [N_1^{\xi}(\xi) \quad N_2^{\xi}(\xi)] \left[\begin{array}{c}u_1 \\ u_2\end{array}\right] = N^{\xi}(\xi) u^e \]

可以看到，坐标变换关系式(7.40)与位移的插值表达式(7.42)在形式上是相同的，即\(N^{\xi}=N^e\)。如果坐标变换和函数插值采用相同的节点，并且采用相同的插值函数，则称这种变换为等参变换，这种单元称为等参单元，如图7 - 1所示。

在母单元中，场函数是用自然坐标表述的，又因为在自然坐标中积分限是规则化的-1到+1，因此希望能够在自然坐标内按照规则化的数值积分进行单元的积分运算。为此，还需要建立自然坐标系与全局坐标系之间导数、微元的变换关系。

\[\frac{d N_i^e(x(\xi))}{d x} = \frac{d N_i^e(x)}{d x} = J(\xi) \frac{d N_i^e(\xi)}{d \xi} \]

其中，\(J(\xi)\)称为雅可比(Jacobian)系数，且得到微元间的变换关系\(d x = J(\xi) d \xi\)，将(7.40)式两端对\(\xi\)求导数，\(J(\xi)\)可以表达为自然坐标的函数

\[J(\xi) = \frac{d x}{d \xi} = \frac{d \bar{N}^e(\xi) x^e}{d \xi} \]

进一步求解(7.44)式，得到

\[J(\xi) = \frac{d x}{d \xi} = \frac{x_2^e - x_1^e}{2} = \frac{l^e}{2} \]

所以两节点杆单元的\(J(\xi)\)是常数。于是，整理(7.43)式，得到自然坐标系下\(N_i^e(x)\)对全局坐标\(x\)导数的显式表达式

\[\frac{d N_i^e(x)}{d x} = \frac{1}{J(\xi)} \frac{d N_i^e(\xi)}{d \xi} \]

于是，(5.15)式在全局坐标系下的单元积分可以转化到自然坐标系下，将(7.46)式代入(5.15)式，得

\[K^e = A E \int_{x_1^e}^{x_2^e} [B^e(x)]^T B^e(x) d x \]

\[= A E \int_{-1}^{1} \left[ \frac{d N_1^e(x)}{d x} \frac{d N_2^e(x)}{d x} \right]^T \left[ \frac{d N_1^e(x)}{d x} \frac{d N_2^e(x)}{d x} \right] J(\xi) d \xi \]

于是，得到了自然坐标系下的单元刚度矩阵

\[K^e = AE \int_{-1}^{1} [B^e(\xi)]^T B^e(\xi) \frac{1}{J(\xi)} d\xi \]

由于(7.47)式的\(B^e(\xi)\)与\(J(\xi)\)是常数，所以积分号内的矩阵是常数矩阵，不必采用数值积分，直接解析积分即可。将(7.39)式与(7.45)式代入(7.47)式，得

\[K^e = \frac{2AE}{l^e} \int_{-1}^{1} \left[ \begin{array}{cc} -\frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & -\frac{1}{2} \end{array} \right]^T \left[ \begin{array}{cc} -\frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & -\frac{1}{2} \end{array} \right] d\xi \\= \frac{2AE}{l^e} \int_{-1}^{1} \left[ \begin{array}{cc} \frac{1}{4} & -\frac{1}{4} \\ -\frac{1}{4} & \frac{1}{4} \end{array} \right] d\xi \\= \frac{AE}{l^e} \left[ \begin{array}{cc} 1 & -1 \\ -1 & 1 \end{array} \right] \]

如果是具有3个节点的等参杆单元，形函数取自然坐标系内的一维3节点拉格朗日插值多项式，即

\[N_1^e = \frac{1}{2}\xi(\xi - 1), \quad N_2^e = (1 - \xi^2), \quad N_3^e = \frac{1}{2}\xi(\xi + 1) \]

那么，单元应变矩阵\(B^e\)为

\[B^e = \left[ \begin{array}{ccc} \frac{dN_1^e}{d\xi} & \frac{dN_2^e}{d\xi} & \frac{dN_3^e}{d\xi} \end{array} \right] = \left[ \begin{array}{ccc} \xi - 1 & -2\xi & \xi + 1 \end{array} \right] \]

Jacobi系数为

\[J(\xi) = \frac{dx}{d\xi} = \frac{dN^e(\xi) x^e}{d\xi} = \left( \frac{\xi - 1}{2} \right) x_1^e + (1 - \xi^2) x_2^e + \left( \frac{\xi + 1}{2} \right) x_3^e \]

于是，得到了自然坐标系下的单元刚度矩阵

\[K^e = AE \int_{-1}^{1} [B^e(\xi)]^T B^e(\xi) \frac{1}{J(\xi)} d\xi \]

由于(7.47)式的\(B^e(\xi)\)与\(J(\xi)\)是常数，所以积分号内的矩阵是常数矩阵，不必采用数值积分，直接解析积分即可。将(7.39)式与(7.45)式代入(7.47)式，得

\[K^e = \frac{2AE}{l^e} \int_{-1}^{1} \left[ \begin{array}{cc} -\frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & -\frac{1}{2} \end{array} \right]^T \left[ \begin{array}{cc} -\frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & -\frac{1}{2} \end{array} \right] d\xi \]

\[= \frac{2AE}{l^e} \int_{-1}^{1} \left[ \begin{array}{cc} \frac{1}{4} & -\frac{1}{4} \\ -\frac{1}{4} & \frac{1}{4} \end{array} \right] d\xi = \frac{AE}{l^e} \left[ \begin{array}{cc} 1 & -1 \\ -1 & 1 \end{array} \right] \]

如果是具有3个节点的等参杆单元，形函数取自然坐标系内的一维3节点拉格朗日插值多项式，即

\[N_1^e = \frac{1}{2}\xi(\xi - 1), \quad N_2^e = (1 - \xi^2), \quad N_3^e = \frac{1}{2}\xi(\xi + 1) \]

那么，单元应变矩阵\(B^e\)为

\[B^e = \left[ \begin{array}{ccc} \frac{dN_1^e}{d\xi} & \frac{dN_2^e}{d\xi} & \frac{dN_3^e}{d\xi} \end{array} \right] = \left[ \begin{array}{ccc} \xi - 1 & -2\xi & \xi + 1 \end{array} \right] \]

Jacobi系数为

\[J(\xi) = \frac{dx}{d\xi} = \frac{dN^e(\xi) x^e}{d\xi} = \left( \frac{\xi - 1}{2} \right) x_1^e + (1 - \xi^2) x_2^e + \left( \frac{\xi + 1}{2} \right) x_3^e \]

6.高斯积分

线性静态有限元分析是有限元分析的最基本内容，通过求解结构的平衡方程，可以求得结构的位移，进而得到结构的应变、应力与约束反力。第2、3、5、6章通过不同方法得到了由单元平衡方程组装而成的结构平衡方程，在静态有限元分析之前需要先求解积分形式的单元刚度矩阵与载荷列阵。通常被积函数的表达式项数很多，进行精确的解析积分十分困难，只能用数值积分方法近似求积。

6.1 数值积分方法

数值积分的基本思想是把积分化为求和，如下式所示：

\[\int_{\xi_1}^{\xi_2} F(\xi) d\xi = \sum_{i} \alpha_i F(\xi_i) \]

\[\int_{\eta_1}^{\eta_2} \int_{\xi_1}^{\xi_2} F(\xi, \eta) d\xi d\eta = \sum_{i} \sum_{j} \alpha_{ij} F(\xi_i, \eta_j) \]

一维问题标准形式可写成

\[\int_{-1}^{+1} F(\xi) d\xi \]

注意积分上、下限分别为 +1，-1。

6.2 高斯积分法

如果事先不确定取样点的位置，对每个取样点有两个待定的量，即积分点的位置\(\xi_i\)和权系数\(H_i\)。如果有\(n\)个取样点就有\(2n\)个条件，可以确定一个\(2n - 1\)阶多项式，用这个多项式近似被积函数并精确积分，积分可写成

\[\int_{-1}^{+1} F(\xi) d\xi=\sum_{i = 1}^{n} H_i F(\xi_i) \]

表7 - 1列出了前三阶高斯积分的积分点坐标和对应的权系数。显然，选取相同数目取样点，高斯给出更高的积分精度。在有限元法中几乎毫无例外地都用高斯法进行数值积分。

表7 - 1 高斯积分的积分点坐标和权系数

积分点数\(n\)	积分点坐标\(\xi_i\)	积分权系数\(H_i\)
1	\(\xi_1 = 0.00000000\)	\(H_1 = 2.00000000\)
2	\(\xi_1, \xi_2=\pm 0.57735027\)	\(H_1, H_2 = 1.00000000\)
3	\(\xi_1, \xi_3=\pm 0.77459667\)	\(H_1, H_3 = 0.55555556\)
	\(\xi_2 = 0.00000000\)	\(H_2 = 0.88888889\)

正因为构造的被积函数是\(2n - 1\)次多项式，因此\(n\)个积分点的高斯积分可达\(2n - 1\)阶的精度。也就是说如果\(F(\xi)\)是不超过\(2n - 1\)次的多项式，积分结果将是精确的。

6.3 单元刚度矩阵与载荷列阵的解析积分和数值积分

对于两节点杆单元，\(K^e\)可以解析积分出具体数值，而不必采用数值积分。将

将(4.16)式代入(6.21)式，得

\[K^e = EA \int_{x_1^e}^{x_2^e} \left( \frac{1}{l^e} \left[ \begin{array}{cc} -1 & 1 \\ 1 & -1 \end{array} \right] \right)^T \left( \frac{1}{l^e} \left[ \begin{array}{cc} -1 & 1 \\ 1 & -1 \end{array} \right] \right) dx \]

\[= \frac{EA}{(l^e)^2} \int_{x_1^e}^{x_2^e} \left[ \begin{array}{cc} -1 & 1 \\ 1 & -1 \end{array} \right] \left[ \begin{array}{cc} -1 & 1 \\ 1 & -1 \end{array} \right] dx = \frac{EA}{l^e} \left[ \begin{array}{cc} 1 & -1 \\ -1 & 1 \end{array} \right] \]

对于3节点2次杆单元（见4.3节），单元应变矩阵\(B^e\)为

\[B^e = \frac{dN^e(x)}{dx} = \left[ \begin{array}{ccc} \frac{dN_1^e(x)}{dx} & \frac{dN_2^e(x)}{dx} & \frac{dN_3^e(x)}{dx} \end{array} \right] = \frac{1}{(l^e)^2} \left[ \begin{array}{ccc} 2x - (x_2^e + x_3^e) & -4x + 2(x_2^e + x_3^e) & 2x - (x_2^e + x_3^e) \end{array} \right] \]

令

\[x_c = \frac{x_2^e + x_3^e}{2} \]

由于\(x_2^e - x_1^e = l^e\)，\(x_3^e - x_1^e = x_3^e - x_2^e = l^e/2\)，那么(7.10)式可进一步表达为

\[B^e = \frac{1}{(l^e)^2} \left[ \begin{array}{ccc} (x - x_c) - \frac{l^e}{2} & -2(x - x_c) & (x - x_c) + \frac{l^e}{2} \end{array} \right] \]

将(7.12)式代入下式的单元刚度矩阵：

\[K^e = \int_{x_1^e}^{x_3^e} (B^e)^T A E B^e dx \]

其中

\[(B^e)^T B^e = \left[ \begin{array}{ccc} \frac{4(x - x_c)^2}{(l^e)^4} & \frac{8(x - x_c)}{(l^e)^3} & \frac{32(x - x_c)^2}{(l^e)^4} \\ \frac{8(x - x_c)}{(l^e)^3} & \frac{32(x - x_c)^2}{(l^e)^4} & \frac{64(x - x_c)^3}{(l^e)^5} \\ \frac{16(x - x_c)^2}{(l^e)^4} & \frac{32(x - x_c)^3}{(l^e)^5} & \frac{1}{(l^e)^4} \end{array} \right] \]

\[+ \left[ \begin{array}{ccc} \frac{16(x - x_c)^2}{(l^e)^4} & \frac{8(x - x_c)}{(l^e)^3} & \frac{32(x - x_c)^2}{(l^e)^4} \\ \frac{8(x - x_c)}{(l^e)^3} & \frac{32(x - x_c)^2}{(l^e)^4} & \frac{64(x - x_c)^3}{(l^e)^5} \\ \frac{4(x - x_c)^2}{(l^e)^4} & \frac{8(x - x_c)}{(l^e)^3} & \frac{1}{(l^e)^2} \end{array} \right] \]

将(7.14)式代入(7.13)式并积分，得

\[K^e = AE \left[ \begin{array}{ccc} \frac{(l^e)^2}{12} \left( \frac{4(x - x_c)}{(l^e)^2} - \frac{1}{l^e} \right) & \frac{4x(x - 2x_c)}{(l^e)^3} - \frac{32(x - x_c)^3}{3(l^e)^4} \\ \frac{4x(x - 2x_c)}{(l^e)^3} - \frac{32(x - x_c)^3}{3(l^e)^4} & \frac{32(x - x_c)^3}{(l^e)^4} - \frac{64(x - x_c)^3}{(l^e)^5} \\ \frac{16(x - x_c)^3}{(l^e)^4} - \frac{4x(x - 2x_c)}{(l^e)^3} & \frac{32(x - x_c)^3}{(l^e)^4} - \frac{32(x - x_c)^3}{(l^e)^4} \\ \frac{16(x - x_c)^3}{(l^e)^4} - \frac{4x(x - 2x_c)}{(l^e)^3} & \frac{32(x - x_c)^3}{(l^e)^4} - \frac{32(x - x_c)^3}{(l^e)^4} \\ \frac{(l^e)^2}{12} \left( \frac{4(x - x_c)}{(l^e)^2} + \frac{1}{l^e} \right) \end{array} \right]_{x_1^e}^{x_3^e} \]

(7.15)式可进一步整理为

\[K^e = \frac{AE}{3l^e} \left[ \begin{array}{ccc} 7 & -8 & 1 \\ -8 & 16 & -8 \\ 1 & -8 & 7 \end{array} \right] \]

显然，以上的解析积分是比较繁琐的，但是仍然可以手动推导，只是工作量较大。事实上，绝大多数的有限元单元都必须通过数值积分来计算。(7.14)式仅仅是3个自由度的杆单元，尚且如此复杂，而平面三角形单元、四边形单元等，其刚度矩阵\(K^e\)只能通过数值积分得到，而推导显式表达式的计算量将不可估计。

接下来，采用高斯数值积分求解(7.14)式，令

\[f(x) = EA \left[ \begin{array}{ccc} \frac{4}{(l^e)^2} (x - x_c) - \frac{1}{l^e} & \frac{8}{(l^e)^3} (x - x_c) & \frac{4}{(l^e)^2} (x - x_c) + \frac{1}{l^e} \\ \frac{8}{(l^e)^3} (x - x_c) & \frac{32}{(l^e)^4} (x - x_c)^2 & \frac{8}{(l^e)^3} (x - x_c) \\ \frac{4}{(l^e)^2} (x - x_c) + \frac{1}{l^e} & \frac{8}{(l^e)^3} (x - x_c) & \frac{4}{(l^e)^2} (x - x_c) + \frac{1}{l^e} \end{array} \right] \]

为了将积分区间\([x_1^e, x_3^e]\)换到标准积分区间\([-1,1]\)，可经变换

\[x = \frac{x_3^e - x_1^e}{2} \xi + \frac{x_3^e + x_1^e}{2} = \frac{l^e}{2} \xi + x_c, \quad x \in [x_1^e, x_3^e] \]

且

\[dx = \frac{l^e}{2} d\xi \]

此时

\[K^e = \int_{x_1^e}^{x_3^e} f(x) dx = \int_{-1}^{1} f\left( \frac{l^e}{2} \xi + x_c \right) \frac{l^e}{2} d\xi = \frac{l^e}{2} \int_{-1}^{1} g(\xi) d\xi \]

其中，\(g(\xi) = f\left( \frac{l^e}{2} \xi + x_c \right) \frac{l^e}{2}\)。

首先采用单点高斯积分

\[K^e = \frac{l^e}{2} \int_{-1}^{1} g(\xi) d\xi = \frac{l^e}{2} H_1 g(\xi_1) \]

其中

\[\xi_1 = 0, \quad H_1 = 2.0 \]

(7.22)

当\(\xi_1 = 0\)时

\[g(0) = \frac{EA}{(l^e)^2} \left[ \begin{array}{ccc} 1 & 0 & -1 \\ 0 & 0 & 0 \\ -1 & 0 & 1 \end{array} \right] \]

将(7.22)式与(7.23)式代入(7.21)式，得到单点积分的刚度矩阵结果为

\[K^e = \frac{EA}{l^e} \left[ \begin{array}{ccc} 1 & 0 & -1 \\ 0 & 0 & 0 \\ -1 & 0 & 1 \end{array} \right] \]

与精确解(7.16)相比较，可知采用单点积分误差较大。

(7.20)式的2点高斯积分公式为

\[K^e = \frac{l^e}{2} \int_{-1}^{1} g(\xi) d\xi = \frac{l^e}{2} [H_1 g(\xi_1) + H_2 g(\xi_2)] \]

其中

\[\xi_1 = -0.57735027, \quad H_1 = 1 \]

\[\xi_2 = +0.57735027, \quad H_2 = 1 \]

对于\(\xi_1=-0.57735027\)

\[g(\xi_1)=\frac{EA}{(l^e)^2} \left[ \begin{array}{ccc} 4.6427344 & -4.9760678 & 0.3333333 \\ -4.9760678 & 5.3333333 & -0.3572655 \\ 0.3333333 & -0.3572655 & 0.0239226 \end{array} \right] \]

对于\(\xi_2=+0.57735027\)

\[g(\xi_2)=\frac{EA}{(l^e)^2} \left[ \begin{array}{ccc} 0.0239226 & -0.3572655 & 0.3333333 \\ -0.3572655 & 5.3333333 & -4.9760678 \\ 0.3333333 & -4.9760678 & 4.6427344 \end{array} \right] \]

将(7.27)式、(7.28)式代入(7.25)式，得到2点高斯积分的刚度矩阵为

\[K^e=\frac{EA}{l^e} \left[ \begin{array}{ccc} 2.3333333 & -2.6666667 & 0.3333333 \\ -2.6666667 & 5.3333333 & -2.6666667 \\ 0.3333333 & -2.6666667 & 2.3333333 \end{array} \right] \]

由此可见，采用2点高斯积分已经可以得到精确结果。

(7.20)式的3点高斯积分公式为

\[K^e=\frac{l^e}{2} \int_{-1}^{1} g(\xi) d\xi=\frac{l^e}{2} [H_1 g(\xi_1)+H_2 g(\xi_2)+H_3 g(\xi_3)] \]

其中

\[\xi_1=-0.77459667, \quad H_1=0.55555556 \]

\[\xi_2=0, \quad H_2=0.88888889 \]

\[\xi_3=+0.77459667, \quad H_3=0.55555556 \]

对于\(\xi_1=-0.77459667\)

\[g(\xi_1)=\frac{EA}{(l^e)^2} \left[ \begin{array}{ccc} 6.49838658 & -7.89838653 & 1.39999994 \\ -7.89838653 & 9.59999977 & -1.70161325 \\ 1.39999994 & -1.70161325 & 0.30161330 \end{array} \right] \]

对于\(\xi_2=0\)

\[g(\xi_2)=\frac{EA}{(l^e)^2} \left[ \begin{array}{ccc} 1 & 0 & -1 \\ 0 & 0 & 0 \\ -1 & 0 & 1 \end{array} \right] \]

对于\(\xi_3 = +0.77459667\)，

\[g(\xi_3)=\frac{EA}{(l^e)^2} \left[ \begin{array}{ccc} 0.30161330 & -1.70161341 & 1.40000007 \\ -1.70161341 & 9.60000027 & -7.89838653 \\ 1.40000007 & -7.89838653 & 6.49838658 \end{array} \right] \]

将(7.32)式~(7.34)式代入(7.30)式，得到

\[K^e=\frac{EA}{l^e} \left[ \begin{array}{ccc} 2.33333333 & -2.66666667 & 0.33333333 \\ -2.66666667 & 5.33333333 & -2.66666667 \\ 0.33333333 & -2.66666667 & 2.33333333 \end{array} \right] \]

便得到\(K^e\)的数值表达式，其与2点高斯积分结果一致。

由7.1节可知：当积分点的个数为\(n\)时，高斯积分具有\(2n - 1\)次的代数精度，即对于一切不高于\(2n - 1\)次的多项式，高斯积分都精确成立；若多项式次数高于\(2n - 1\)次，则不成立。那么，当积分点个数为1时，高斯积分具有1次代数精度；当积分点个数为2时，高斯积分具有3次代数精度。所以，对于3节点杆单元的2次多项式被积函数(7.17)，需要2个积分点即可精确积分，这个结论也得到了上文的验证。

对于有3个及3个以上节点的单元，内部节点自由度可以在单元层次凝聚掉，而只保持外部节点自由度参加系统方程的集成，以提高计算效率。

之所以采用数值积分，是因为对于自由度较多的单元(如板壳结构单元、平面或实体单元)解析推导积分的过程太过繁琐，而并不是因为被积函数不可解析积分。事实上，有限元单元常用的Lagrange插值函数或Hermite插值函数都是解析可积的。

以上是单元刚度矩阵的积分，对于(6.22)式载荷列阵的积分相对容易，可以采用解析积分或者数值积分。此处，将2节点杆单元插值形函数代入(6.22)式的\(f_{x}^{e}\)项，同时假设\(b\)为常数，得到

\[f_{x}^{e}=\int_{x_1^{e}}^{x_2^{e}} (N^e)^T b dx = \frac{b l^e}{2} \left[ \begin{array}{c} 1 \\ 1 \end{array} \right] \]

可以看出，(7.36)式的积分相当于将体力集中于杆单元的两个节点自由度上。另外，对于(6.22)式的\(f_{r}^{e}\)项，可由插值函数的基本属性\(N^e(x_r^{e})=\delta_{ij}\)，得

\[f_{r}^{e}=(N^e)^T A \left. \frac{du^e}{dx} \right|_{x = x_r^{e}} = \left[ \begin{array}{c} 0 \\ A \left. \frac{du^e}{dx} \right|_{x = x_r^{e}} \end{array} \right] \]

上式说明，作用于单元节点上的载荷，直接将其放置到单元载荷列阵相应自由度上即可。

7.方程求解