▲图1
如图1所示为某大桥的主桁架。主桁架的上弦杆与两个斜腹杆的中心线不交于一点,而是相隔一定距离,如果忽略它将导致误差。由于结点1和2都连在刚性很大的结点块上,如图2所示,因此可假定它们的线位移和转角都相等。
▲图2
如图3所示的张弦梁结构,钢索和撑杆都是铰接于主梁。在有限元模型中,梁、杆、索属于不同的单元类型,虽然这些结点具有相同的节点线位移,但截面转角不相同,此时我们可以在该处定义两个坐标一样的结点,然后指定这两个结点的线位移相等。
▲图3
由上述结构案例可知,两个自由度之间有某种约束,不失一般性,我们可以把它写成
\beta_1Q_1 + \beta_2Q_2 = \alpha \quad\cdots (1)
考虑约束条件
(1)下的势能泛函极值问题
\begin{split}
min.\quad &\Pi = \frac{1}{2} \mathbf Q^T \mathbf K \mathbf Q - \mathbf Q^T \mathbf F\\
s.t.\quad & \beta_1Q_1 + \beta_2Q_2 - \alpha=0\\
\end{split} \quad \cdots (2)
用罚函数将有约束问题转化为无约束问题。引入一个很大的正参数
C,构造新的泛函
\Pi_1 = \frac{1}{2} \mathbf Q^T \mathbf K \mathbf Q + \frac{1}{2} C(\beta_1Q_1 + \beta_2Q_2 - \alpha)^2 - \mathbf Q^T \mathbf F \quad\cdots (3)
令
\frac {\partial \Pi_1 }{\partial Q_i}=0,i=1,2,\cdots,n
得到新的平衡方程
\begin{bmatrix}
K_{11}+C\beta_1^2 & K_{12}+C\beta_1\beta_2 & \cdots & K_{1n} \\
K_{21}+C\beta_1\beta_2 & K_{22}C\beta_2^2 & \cdots & K_{2n} \\
\cdots & \cdots &\cdots &\cdots \\
K_{n1} & K_{n2} & \cdots & K_{nn} \\
\end{bmatrix}
\begin{Bmatrix}
Q_1 \\
Q_2 \\
\cdots \\
Q_n \\
\end{Bmatrix} =
\begin{Bmatrix}
F_1+C\alpha\beta_1 \\
F_2+C\alpha\beta_2 \\
\cdots \\
F_n \\
\end{Bmatrix} \quad\cdots (4)
例1
如图4所示,有一根忽略质量的刚性杆,它的一端铰接,其上还连接有根钢质杆和一根铝质杆,其右端作用有外力
P=30\times10^3N。采用两个单元进行建模,求1、2两点的位移。
▲图4
用两个单元对该问题进行建模,单元节点信息见下表
节点3、4处的边界条件为
Q_3=0和
Q_4=0。因为刚性杆保持直线,
Q_1和
Q_2的相对关系如图5所示,根据几何关系,得到相关节点的约束如下
▲图5
Q_1 -0.4Q_2=0 \quad \cdots (5)
单元①刚度矩阵为
\begin{split}
\mathbf k^1 &= \frac {200\times10^3 \times 1200}{4500}
\begin{bmatrix}
1 & -1 \\
-1 & 1 \\
\end{bmatrix} \\
&= 10^3 \begin{bmatrix}
53.33 & -53.33 \\
-53.33 & 53.33 \\
\end{bmatrix} \\
\end{split}
单元②刚度矩阵为
\begin{split}
\mathbf k^2 &= \frac {70\times10^3 \times 900}{3000}
\begin{bmatrix}
1 & -1 \\
-1 & 1 \\
\end{bmatrix} \\
&= 10^3 \begin{bmatrix}
21 & -21 \\
-21 & 21 \\
\end{bmatrix} \\
\end{split}
整体刚度矩阵为
\begin{split}
\mathbf K &= 10^3
\begin{bmatrix}
53.33& 0& -53.33& 0\\
0 & 21& 0& -21\\
-53.33& 0& 53.33& 0\\
0 & -21& 0& 21\\
\end{bmatrix}
\\
\end{split}
接下来修正矩阵
\mathbf K。选取一个远大于刚度系数的正数
C=53.33\times10^7,由于
Q_3=Q_4=0,将
C加到
\mathbf K中
K_{33}和
K_{44}的位置上。
\begin{split}
\mathbf K &= 10^3
\begin{bmatrix}
53.33& 0& -53.33& 0\\
0 & 21& 0& -21\\
-53.33& 0& 53.33+C& 0\\
0 & -21& 0& 21+C\\
\end{bmatrix}
\\
\end{split}
然后考虑(5)中给出的多点约束方程,注意到
\alpha=0,\beta_1=1,\beta_2=-0.4则由式(4)得到修正后的刚度矩阵为
\begin{split}
\mathbf K &= 10^3
\begin{bmatrix}
53.33+C& -0.4C & -53.33& 0\\
-0.4C & 21+0.16C & 0& -21\\
-53.33& 0& 53.33+C & 0\\
0 & -21& 0& 21+C\\
\end{bmatrix}
\\
\end{split}
修正后的整体等效节点荷载矩阵
\mathbf F =
\begin{Bmatrix}
0+0C \\
30\times 10^3+0C\\
0 \\
0 \\
\end{Bmatrix}
\\
例2
▲图6
如图6所示的结构,中间的铰接点不能看作拥有两个自由度的一个节点。因为连续梁的挠度函数在铰接点这里虽然连续但不可导,即在节点两边,不同单元的转角是不一样的。
▲图7
所以铰接点要建立两个节点,如图7所示。这样一来自由度1和自由度3对应的线位移必须相等,就需要建立约束关系
u_1-u_3=0选取一个远大于刚度系数的正数
C,修正矩阵
K.由(1)知,
\alpha=0,\beta_1=1,\beta_2=-1。修正后的刚度矩阵为
\mathbf K =
\begin{bmatrix}
\frac{12EI}{l^3}+C & -\frac{6EI}{l^2}& -C& 0\\
-\frac{6EI}{l^2} & \frac{4EI}{l}& 0& 0\\
-C& 0& \frac{12EI}{l^3}+C & \frac{6EI}{l^2}\\
0 & 0& \frac{6EI}{l^2}& \frac{4EI}{l}\\
\end{bmatrix}
数值验证参考:罚单元