考虑一段一端由线性弹簧支撑的细长梁，如图1所示，弹簧刚度为
r，则弹簧的应变能为
▲图1
U_k = \frac 12 r(\omega_l)^2
其中，
\omega_l是梁在
x=l处的挠度。
梁的势能为
\Pi = \frac 12 \int_V \sigma \epsilon dV + \frac 12 r(\omega_l)^2 \quad \cdots (1)
现在来求
\omega_l▲图2
如图2所示的两节点梁单元，
\omega_1,\theta_1,\omega_2,\theta_2分别为四个节点自由度。物理坐标系
(xOy)和自然坐标系
(\xi)的线性映射关系
x = (1-\xi)x_1 +\xi x_2
在节点1的位置时
\xi=0,在节点2的位置时
\xi=1。挠度由三次多项式表示为
\omega = b_1 + b_2\xi + b_3\xi^2 + b_4\xi^3 \quad \cdots (2)
梁的转角
\begin{split}
\theta &= -\frac {\partial \omega }{\partial x} \\
&= -\frac {\partial \omega }{\partial \xi}\frac {\partial \xi }{\partial x} \\
&= -\frac 1l(b_2+2b_3\xi + 3b_4\xi^2) \\
\end{split}\quad \cdots (3)
其中，
l= x_2-x_1为单元长度。在节点1的位置时
\xi=0,由(2)(3)可得
\omega_1 = b_1\quad \cdots (4.1)
\theta_1 = -\frac 1l b_2 \quad \cdots (4.2)
在节点2的位置时
\xi=1,由(2)(3)可得
\omega_2 = b_1 + b_2 + b_3 + b_4 \quad \cdots (4.3)
\theta_2 = -\frac 1l (b_2 + 2b_3 + 3b_4 )\quad \cdots (4.4)
(4.1)~(4.4)联立解得
\begin{split}
b_1 &= \omega_1 \\
b_2 &= -l \theta_1 \\
b_3 &= -3\omega_1 + 3\omega_2 + 2l\theta_1 + l\theta_2 \\
b_4 &= 2 \omega_1 -2\omega_2 -l\theta_1 - l\theta_2
\end{split} \quad \cdots (5)
(5)代入(2)得
\omega = N_1\omega_1 + N_2l \theta_1 + N_3\omega_2 + N_4l \theta_2 \quad \cdots (6)
其中
\begin{split}
N_1 &= 1-3\xi + 2\xi^3 \\
N_2 &= -\xi + 2\xi^2 -\xi^3 \\
N_3 &= 3\xi^2 -2\xi^3 \\
N_4 &= \xi^2 -\xi^3
\end{split} \quad \cdots (7)
记
\mathbf N=[N_1, N_2l, N_3, N_4l],(6)写成矩阵形式
\omega = \mathbf N \mathbf q^e \quad \cdots (8)
由(3)可知，梁的广义应变(曲率)
\begin{split}
\epsilon &= \frac {\partial^2 \omega }{\partial x^2} \\
&= \frac {\partial }{\partial \xi }(\frac {\partial \omega }{\partial\xi }
\frac {\partial \xi }{\partial x }) \\
&= \frac {1}{l^2}\frac {\partial^2 \omega }{\partial \xi^2} \\
&= \frac {1}{l^2}(\frac {\partial^2 N_1 }{\partial \xi^2}\omega_1+ \frac {\partial^2 N_2 }{\partial \xi^2}l \theta_1 + \frac {\partial^2 N_3 }{\partial \xi^2}\omega_2 + \frac {\partial^2 N_4 }{\partial \xi^2}l \theta_2) \\
&= \frac {1}{l^2}[-6+12\xi, (4-6\xi)l, 6-12\xi, (2-6\xi)l]\begin{Bmatrix}
v_1 \\
\theta_1 \\
v_2 \\
\theta_2 \\
\end{Bmatrix} \\
\end{split}
记
\mathbf B = \frac {1}{l^2}[-6+12\xi, (4-6\xi)l, 6-12\xi, (2-6\xi) ]
则
\epsilon = \mathbf B \mathbf q^e
其中
\mathbf B是应变矩阵。应力
\sigma = E \epsilon= E \mathbf B \mathbf q^e
由(7)(8)可知，(1)中的
\omega_l就是
\xi=1时
\omega的值。
\begin{split}
N_1(1) &= 1-3 + 2 = 0 \\
N_2(2) &= -1 + 2 -1 =0 \\
N_3(3) &= 3-2=1 \\
N_4(4) &= 1-1=0
\end{split}
\omega_l = [0,0,1,0] \mathbf q^e \quad \cdots (10)
梁单元势能的表达式
\begin{split}
\Pi &= U + U_r \\
&= \frac{1}{2} \mathbf q^{eT} \mathbf k_1^e \mathbf q^{e} +
\frac{1}{2} \mathbf q^{eT} \mathbf k_2^e \mathbf q^{e}
\end{split} \quad \cdots (11)
其中 记
\mathbf k_1^e = EI \int_0^1\mathbf B^T \mathbf B d\xi
\mathbf k_2^e = r \mathbf N(\xi=1)^T \mathbf N(\xi=1)
\mathbf k^e = \mathbf k_1^e + \mathbf k_2^e
经计算得
\mathbf k^e =
\frac {EI}{l^3} \begin{bmatrix}12 & 6l & -12 & 6l \\ 6l & 4l^2 & -6l & 2l^2 \\
-12 & 6l & 12 & -6l \\
6l & 2l^2 & -6l & 4l^2 \\
\end{bmatrix}+
\begin{bmatrix}0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\
0 & 0 & r & 0 \\
0 & 0 & 0 & 0 \\
\end{bmatrix}
对于具有弹性支承的单元，只需将
\mathbf k_2^e 加入到常规单元刚度矩阵中。
▲图3
对于图3所示得梁，弹簧刚度为
r= \frac {EI}{l^3}。只划分一个单元时，其有限元平衡方程为
\frac {EI}{l^3} \begin{bmatrix}12 & 6l & -12 & 6l \\ 6l & 4l^2 & -6l & 2l^2 \\
-12 & 6l & 12+1 & -6l \\
6l & 2l^2 & -6l & 4l^2 \\
\end{bmatrix}
\begin{Bmatrix}\omega_1 \\ \theta_1 \\ \omega_2 \\ \theta_2 \\
\end{Bmatrix}=
\begin{Bmatrix}R_1 \\ R_2 \\ -P \\ 0 \\
\end{Bmatrix}
考虑边界条件之后
\frac {EI}{l^3} \begin{bmatrix}
13 & -6l \\
-6l & 4l^2 \\
\end{bmatrix}
\begin{Bmatrix} \omega_2 \\ \theta_2 \\
\end{Bmatrix}=
\begin{Bmatrix} -P \\ 0 \\
\end{Bmatrix}
粤公网安备44030502008569号