考虑一个有限元模型的势能泛函
\Pi = \frac{1}{2} \mathbf Q^T \mathbf K \mathbf Q - \mathbf Q^T \mathbf F \quad\cdots (1)
其中，
\mathbf Q是整体节点位移矩阵，
\mathbf K是整体刚度矩阵。已知沿支座的自由度1的方向产生位移为
a_1，即
Q_1=a_1。那么求势能泛函(1)的极值变成了在约束条件
Q_1=a_1下的极值问题。
\begin{split}
min.\quad &\Pi = \frac{1}{2} \mathbf Q^T \mathbf K \mathbf Q - \mathbf Q^T \mathbf F\\
s.t.\quad & Q_1-a_1 = 0 \\
\end{split} \quad \cdots (2)
用罚函数将有约束问题转化为无约束问题。引入一个很大的正参数
C，构造新的泛函
\Pi_1 = \frac{1}{2} \mathbf Q^T \mathbf K \mathbf Q + \frac{1}{2} C(Q_1-a_1)^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 & K_{12} & \cdots & K_{1n} \\
K_{21} & K_{22} & \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+Ca_1 \\
F_2 \\
\cdots \\
F_n \\
\end{Bmatrix} \quad\cdots (4)
这里可以看到，为处理
Q_1=a_1，需要对以上方程进行修正，即将一个大数
C加到
\mathbf K的第1个对角元上以及将
Ca_1加到
F_1上，式(4)的解即为节点位移列阵
\mathbf Q。节点1处的支反力为
R_1 = -C(Q_1-a_1) \quad\cdots (5)
需要说明的是:这里所说的罚函数法只是一种近似的方法，最终求解的精度，特别是支反力的求解精度，取决于
C的选取。考察方程组(4)中的第一个方程
(K_{11}+C)Q_1 +K_{12}Q_2 +\cdots + K_{1n}Q_n = F_1+Ca_1 \quad\cdots (6)
(6)两边除以C
(\frac {K_{11}}{C}+1)Q_1 +\frac {K_{12}}{C}Q_2 +\cdots + \frac {K_{1n}}{C}Q_n = \frac {F_1 }{C}+a_1 \quad\cdots (7)
从(7)式中我们可以发现，如果
C比刚度系数
K_{11},K_{12},…,K_{1n}大得多时，那么有
Q_1\approx a_1。
F是施加在支座处的外力(如果存在的话)，而且
F/C通常是一个很小的值。
例1
图1中，有一个外力
P=60\times10^3N在作用在杆中点，求杆的位移、应力和支座支反力(取
E=20\times10^3N/mm)。
▲图1
在这个问题中，我们首先应该确定杆与墙之间是否发生接触。假设墙不存在，那么点B的位移为
Q_B=1.8mm.从这个结果可以看出接触是存在的，因为边界条件发生了变化，即点B的位移是给定的
1.2mm，所以需要重新求解。建立两个单元的有限元模型如图1b所示，边界条件为
Q_1=0和
Q_3=1.2mm，结构刚度矩阵为
\mathbf K =\frac {20\times 10^3 \times 250 }{150}
\begin{bmatrix}
1 & -1 & 0 \\
-1 & 2 & -1 \\
0 & -1 & 1 \\
\end{bmatrix}
整体节点力列阵为
\mathbf F =
\begin{Bmatrix}
0 \\
60\times10^3\\
0 \\
\end{Bmatrix}
按照上述方法，取
C=\frac {2}{3}\times10^9,则修正后的方程如下
\frac {10^5 }{3}
\begin{bmatrix}
20001 & -1 & 0 \\
-1 & 2 & -1 \\
0 & -1 & 20001 \\
\end{bmatrix}
\begin{Bmatrix}
Q_1\\
Q_2\\
Q_3\\
\end{Bmatrix} =
\begin{Bmatrix}
0 \\
60\times10^3\\
80\times10^8\\
\end{Bmatrix}
解得
\mathbf Q =
\begin{Bmatrix}
7.49985\times10^-5\\
1.500045\\
1.200015\\
\end{Bmatrix}
支反力为
\begin{split}
R_1 &= -C\times 7.49985\times10^-5 = 49.999\times10^-3 \\
R_3 &= -C(1.200015-1.2) = -10.001\times10^3\\
\end{split}
精确解为
R_1 =-50\times10^-3N,
R_3 =-10\times10^3N。由罚函数法得到的结果有一个小的近似误差。
例2
▲图2
如图2所示的四杆桁架结构。对于每个单元，给定
E=29.5\times10^6,
A=1.0。节点2竖直下移
0.12。用先处理法建立如图3所示的有限元模型，得到的平衡方程组为
▲图3
\frac {29.5\times10^6}{600}
\begin{bmatrix}
15.0 & 0& 0& 0\\
0 & 20.0& 0& -20 \\
0 & 0& 22.68& 5.76 \\
0 & -20.0& 5.76& 24.32 \\
\end{bmatrix}
\begin{Bmatrix}
Q_1\\
Q_2\\
Q_3\\
Q_4\\
\end{Bmatrix} =
\begin{Bmatrix}
20000 \\
0\\
0\\
-25000\\
\end{Bmatrix}
考虑支座移动的影响，需要修正平衡方程组。在与指定位移值的自由度相对应的结构刚度矩阵对角线元素上加上一个大常数
C。一般来讲，
C可以取为修正前刚度矩阵的最大对角线元素的
10^4倍。此外，载荷列阵的对应位置还要加上力
Ca，其中
a为指定的位移值。在本例当中，对于自由度2来说，
a=-0.12，因而载荷列阵的第2行上要加上一个大小为
-0.12C的力。修正后的有限元方程组为
\frac {29.5\times10^6}{600}
\begin{bmatrix}
15.0 & 0& 0& 0\\
0 & 20.0+C& 0& -20 \\
0 & 0& 22.68& 5.76 \\
0 & -20.0& 5.76& 24.32 \\
\end{bmatrix}
\begin{Bmatrix}
Q_1\\
Q_2\\
Q_3\\
Q_4\\
\end{Bmatrix} =
\begin{Bmatrix}
20000 \\
-0.12C\\
0\\
-25000\\
\end{Bmatrix}
