冰溜子，也称为冰锥，通常是在雪融化成水，顺屋檐或其他高处滴下时，由于气温骤降至冰点以下，这些雪水会迅速冻结，形成上粗下细的锥形冰柱。
本文按照"模型简化-划分单元-组装整体刚度矩阵和整体节点力矩阵"的常规思路来建立冰溜子的有限元模型。
模型简化
对于一根悬挂的冰溜子，将其简化为一维线性变截面杆模型，荷载就是其自重。如图1所示
▲图1
单元划分
利用最小势能原理来推导线性变截面杆的刚度矩阵和考虑自重时的等效节点荷载。
▲图2
如图2所示，线性变截面杆件两个节点的位移分别为
u_1,u_2，两端截面积分别为
A_1,A_2.冰的密度为ρ，由线性插值得位移场为
u(x)=(1-\frac{x}{l})u_1 + \frac{x}{l}u_2 = \mathbf N \mathbf q^e \quad \cdots (1)
其中
\mathbf N叫做形状函数矩阵,
\mathbf q^e叫做节点位移列阵。
\mathbf N =[1-\frac{x}{l}, \frac{x}{l}]
\quad \cdots (2)
\mathbf q^e =
\begin{Bmatrix}
u_1 \\
u_2 \\
\end{Bmatrix} \quad \cdots (3)
同理，任意一个截面得面积为
A(x)=(1-\frac{x}{l})A_1 + \frac{x}{l}A_2 = \mathbf N \mathbf A \quad \cdots (4)
其中
\mathbf A =
\begin{Bmatrix}
A_1 \\
A_2 \\
\end{Bmatrix}
应变
\begin{split}
\epsilon(x) &= \frac{du(x)}{dx} \\
&= \frac{d}{dx} \mathbf N(x) \mathbf q^e \\
&= \mathbf B \mathbf q^e
\end{split}\quad \cdots (5)
其中，
\mathbf B是应变矩阵
\mathbf B =[-\frac{1}{l}, \frac{1}{l}]
\quad \cdots (6)
应力
\sigma(x) = E \epsilon(x) = E \mathbf B \mathbf q^e
\quad \cdots (7)
单元势能的表达式
\begin{split}
\Pi &= U - W_e \\
&= \frac{1}{2} \int_V \sigma(x)\epsilon(x) dV - \int_V \rho g u(x)dV \\
&= \frac{1}{2} \mathbf q^{eT}\int_0^l E\mathbf B^T \mathbf B A(x)dx \mathbf q^{e} - \mathbf q^{eT}\int_0^l \rho g\mathbf N^T \mathbf N \mathbf A dx \\
\end{split} \quad \cdots (8)
(8)作变分运算
\delta \Pi = \int_0^l E\mathbf B^T \mathbf B A(x)dx \mathbf q^{e} - \int_0^l\rho g \mathbf N^T \mathbf N \mathbf A dx \quad \cdots (9)
记
\mathbf k^e = \int_0^l E \mathbf B^T \mathbf B A(x) dx
\mathbf f^e = \int_0^l \rho g\mathbf N^T \mathbf N \mathbf A dx
令
\delta \Pi = 0，得
\mathbf k^e \mathbf q^{e}= \mathbf f^e \quad \cdots (10)
经计算得
\begin{split}
\mathbf k^e &= E\int_0^l
\begin{Bmatrix}-\frac{1}{l} \\ \frac{1}{l} \\ \end{Bmatrix}
\begin{bmatrix}-\frac{1}{l} & \frac{1}{l} \\ \end{bmatrix}
[(1-\frac{x}{l})A_1 + \frac{x}{l}A_2 ]dx \\
&=\frac{E}{l}\frac{A_1+A_2}{2}
\begin{bmatrix}1 & -1 \\ -1 & 1 \\\end{bmatrix} \\
\end{split} \quad \cdots (11)
\begin{split}
\mathbf f^e &= \rho g\int_0^l
\begin{Bmatrix}1-\frac{x}{l} \\ \frac{x}{l}\\ \end{Bmatrix}
\begin{bmatrix}1-\frac{x}{l} & \frac{x}{l}\\ \end{bmatrix}
\begin{Bmatrix}A_1 \\ A_2\\ \end{Bmatrix} dx \\
&=\frac{\rho g l}{6}
\begin{Bmatrix}2A_1 + A_2 \\ A_1 + 2A_2 \end{Bmatrix} \\
\end{split} \quad \cdots (12)
对于线性变截面杆件，第一步是将杆件转化为一个阶梯状模型。这个模型由若干个离散的单元组成，每个单元都简化成等截面，由(11)可知，截面面积取单元两端截面积的平均值。比如，我们用四个单元来建立这根杆的模型，如图3a所示
▲图3
图3b即为由此得到的4个单元5个节点的有限元模型。此时单元
i(i=1,2,3,4)刚度矩阵为
\mathbf k^{i} = \frac{EA_i}{l_i}
\begin{Bmatrix}
1 & -1\\
-1 & 1\\
\end{Bmatrix}
若单元已近似为等截面，则单元
i(i=1,2,3,4)等效节点力为
(\mathbf f^e)_i = \frac{\rho gA_il_i}{2}
\begin{bmatrix}
1\\
1\\
\end{bmatrix}
组装刚度矩阵
在一维问题中，每个节点只有一个自由度，图2b中的五个节点的有限元模型就有五个自由度。其中
Q_1=0，此即为边界条件。由于每个单元有两个节点，因此，可以简单地对单元节点连接的信息进行表述，如图4所示。这样，根据单元节点连接信息就可以建立起局部-整体的节点对应关系。
▲图4
例如，单元1刚度矩阵在整体刚度矩阵得位置为
\mathbf K = \frac{EA_1}{l_1}
\begin{bmatrix}
1& -1& 0& 0& 0\\
-1 & 1& 0& 0& 0\\
0 & 0& 0& 0& 0\\
0 & 0& 0& 0& 0\\
0 & 0& 0& 0& 0\\
\end{bmatrix}
单元2刚度矩阵在整体刚度矩阵得位置为
\mathbf K = \frac{EA_2}{l_2}
\begin{bmatrix}
0 & 0& 0& 0& 0\\
0 & 1& -1& 0& 0\\
0 & -1& 1& 0& 0\\
0 & 0& 0& 0& 0\\
0 & 0& 0& 0& 0\\
\end{bmatrix}
单元3刚度矩阵在整体刚度矩阵得位置为
\mathbf K = \frac{EA_3}{l_3}
\begin{bmatrix}
0 & 0& 0& 0& 0\\
0 & 0& 0& 0& 0\\
0 & 0& 1& -1& 0\\
0 & 0& -1& 1& 0\\
0 & 0& 0& 0& 0\\
\end{bmatrix}
单元4刚度矩阵在整体刚度矩阵得位置为
\mathbf K = \frac{EA_4}{l_4}
\begin{bmatrix}
0 & 0& 0& 0& 0\\
0 & 0& 0& 0& 0\\
0 & 0& 0& 0& 0\\
0 & 0& 0& 1& -1\\
0 & 0& 0& -1& 1\\
\end{bmatrix}
单元1等效节点力在整体等效节点力矩阵得位置为
\mathbf F = \frac{\rho gA_1l_1}{2}
\begin{Bmatrix}
1\\
1\\
0\\
0\\
0\\
\end{Bmatrix}
单元2等效节点力在整体等效节点力矩阵得位置为
\mathbf F = \frac{\rho gA_2l_2}{2}
\begin{Bmatrix}
0\\
1\\
1\\
0\\
0\\
\end{Bmatrix}
单元3等效节点力在整体等效节点力矩阵得位置为
\mathbf F = \frac{\rho gA_3l_3}{2}
\begin{Bmatrix}
0\\
0\\
1\\
1\\
0\\
\end{Bmatrix}
单元4等效节点力在整体等效节点力矩阵得位置为
\mathbf F = \frac{\rho gA_4l_4}{2}
\begin{Bmatrix}
0\\
0\\
0\\
1\\
1\\
\end{Bmatrix}
整体刚度矩阵
K\mathbf K = E
\begin{bmatrix}
\frac{A_1}{l_1} & -\frac{A_1}{l_1}& 0& 0& 0\\
-\frac{A_1}{l_1} & (\frac{A_1}{l_1}+\frac{A_2}{l_2}) & -\frac{A_2}{l_2} &0& 0\\
0 & -\frac{A_2}{l_2}& (\frac{A_2}{l_2}+\frac{A_3}{l_3})& 0& 0\\
0 & 0& -\frac{A_3}{l_3}& (\frac{A_3}{l_3}+\frac{A_4}{l_4})& \frac{A_4}{l_4}\\
0 & 0& 0& -\frac{A_4}{l_4}& \frac{A_4}{l_4}\\
\end{bmatrix}
整体等效节点力矩阵
\mathbf F = \frac{\rho g}{2}
\begin{Bmatrix}
A_1l_1\\
A_1l_1+A_2l_2\\
A_2l_2+A_3l_3\\
A_3l_3+A_4l_4\\
A_4l_4\\
\end{Bmatrix}
用"划行划列法"处理边界条件之后，得到的有限元平衡方程为
\begin{bmatrix}
(\frac{A_1}{l_1}+\frac{A_2}{l_2}) & -\frac{A_2}{l_2} &0& 0\\
-\frac{A_2}{l_2}& (\frac{A_2}{l_2}+\frac{A_3}{l_3})& 0& 0\\
0& -\frac{A_3}{l_3}& (\frac{A_3}{l_3}+\frac{A_4}{l_4})& \frac{A_4}{l_4}\\
0& 0& -\frac{A_4}{l_4}& \frac{A_4}{l_4}\\
\end{bmatrix}
\begin{Bmatrix}
q_2\\
q_3\\
q_4\\
q_5\\
\end{Bmatrix}=
\begin{Bmatrix}
A_1l_1+A_2l_2\\
A_2l_2+A_3l_3\\
A_3l_3+A_4l_4\\
A_4l_4\\
\end{Bmatrix}
