有限元基础编程——耦合边界条件处理

引言

本期内容主要讨论有限元编程过程中，出现多点之间约束耦合的关系（如斜支座），就会出现带约束方程的表达式，这时使用直接消去法、置"1"法、乘大数法时会有所限制，读者可根据本篇的例子进行计算，更能明白限制之处。本文先从通用的含有约束线性方程组求解入手，带着大家了解拉格朗日乘子法和罚函数的基础概念，然后使用有限元方法开始解决斜支座问题。理论部分详看曾攀老师的《有限元分析基础教程》P179, 本文基于Matlab语言进行数值实现，参照知乎Mathematica风格^[1]。

线性方程组求解

  在Matlab中进行线性方程组求解时，可使用反斜杠“\”进行求解，解析解为： 

位移约束耦合

约束关系：    

出现类似这样的耦合关系的时候，直接用高斯消去法进行求解时，会和真实解出现误差。这时候采用拉格朗日乘子法或者罚函数方法进行求解可得到满足约束关系的解析解或者近似解。

拉格朗日乘子法

约束关系中的系数矩阵C、d可写为：

然后将原系数矩阵进行重新组装：

然后进行矩阵求解：

上式的解满足例题中的约束关系，同时也可算出拉格朗日乘子    ，工程中即为约束内力。

罚函数法

该方法将引入一个较大的系数来形成一个罚函数(Penalty Function approach)，从而取得近似解。这里仍要用到上述约束方程的系数矩阵    ，罚函数的系数可设为原系数矩阵最大元素最大值的1e6倍（注意不是1e6这个数值哦），例子中的罚系数可取    。罚函数法的刚度矩阵可写为：

约束内力可由近似解的误差乘以大数    获得，即：

所得结果与拉格朗日乘子值一致，罚函数求解过程中没有改变原系数矩阵的大小，维数不变，变化后的K还是一个对称矩阵，为求解方程带来极大的便利性。

平面应力状态斜支座边界约束问题

例：

书中给出的解析解为：

拉格朗日乘子法得出的结果为：

罚函数法得出的结果为：

支反力为：

Matlab源码

clear
%---------约束条件----------------------
% u4-v4=0
%---------------------------------------
%-------拉格朗日乘子法------------------
Q = [0 0 1 -1]';
f = 0;
K_Lambda = [[K Q];[Q' 0]];
F_Lambda = [F;f];
U_Lambda = K_Lambda\F_Lambda
%-------罚函数法-----------------------
Penalty = 1.0e6*E;
K_Penalty = K+Penalty*Q*Q';
F_Penalty = F+Penalty*f*Q;
U_Penalty = K_Penalty\F_Penalty
RF = (U_Penalty(3)-U_Penalty(4)-f)*Penalty

引用链接

[1] : https://zhuanlan.zhihu.com/p/114555696