有限元分析中剪切锁定(shear locking)

8月前浏览177

平面问题的剪切锁定

剪切锁定的实质就是，在理论上没有剪切变形的单元中发生了剪切变形，这个剪切变形也被称为伴生剪切。该现象一般发生在受纯弯状态下，主要的现象是模型弯曲刚度变大，导致变形偏小。如果能直接得到控制微分方程的解析解，则不会出现剪切锁死的现象，但是如果采用有限元方法，由于网格划分则会出现剪切锁死现象，得到的结果将出现误差。

▲图1

如图1(b)所示，这个模型是材料力学中的平截面假定，截面受到弯矩作用发生纯弯变形，上表面受到水平压应力作用，发生收缩变形，下表面受到水平拉应力作用，发生拉伸变形，并且竖线直线长度不变，始终与上下表面垂直。此时没有发生剪切变形，剪应力为0，该截面变形为理论上的纯弯变形。如图1(c)，这个模型是有限元方法截面变形图，可以观察到该截面上表面也是压缩变形，下表面也是拉伸变形，但是竖向直线有了增大或减小，即发生了剪切变形，消耗了部分变形能量，导致单元过硬，弯曲变形减小，与理论计算有所偏差。

有限元方法可以通过细化网格来有效改善剪切锁死现象，通过增加模型的网格密度，特别是在潜在发生剪切锁死的区域进行网格细化，可以提高模型的准确性和稳定性，减少剪切锁死的发生。但是会导致单元节点增大,使计算的工作量大大增加。也可以通过采用缩减积分方法，减少单元的积分点，降低数值计算中的刚度，从而减轻剪切锁死现象的发生。但是如果网格的尺寸不到位，结果的误差还是会比较大，这与有限元的沙漏现象有关。采用非协调单元也是缓解剪切锁死现象的方法，但是在单元形状比较畸形时，计算结果将会出现很大误差。

中厚板问题的剪切锁定

由于Kirchhoff板理论忽略了剪切变形，但要求挠度的导数连续，给构造协调单元带来了不少麻烦。为此，采用考虑剪切变形的Mindlin板理论来克服，即节点转动和挠度相互独立。这种方法比较简单，精度较好，并且能利用等参变换，得到任意四边形甚至曲边四边形单元，因而实用价值较高。

▲图2

如图2所示的4结点四边形板单元。根据Mindlin板理论的假设，板内任意一点的位移由三个广义位移，和完全确定。为了与有限元的结点位移相对应，采用的位移数组为

式中和是跟挠度无关的独立的转角位移。其结点位移数组为

其中，

位移插值

写成矩阵形式

Mindlin板理论考虑了横向剪切变形，因此应变有5个分量，即

写成矩阵形式

其中

而其中的子矩阵为

而相应的应力分量也有5个，它们与应变之间的关系是

式中弹性矩阵为

应力与结点位移之间的关系是

单元刚度矩阵为

由此可知，单元剪切刚度矩阵是的函数，而单元弯曲刚度矩阵是的函数，当板的厚度与其边长相比变得非常小时，剪切能量占主导地位。这被称为剪切锁定。解决这个问题的一种方法是缩减积分。例如，如果使用8节点四边形单元时，计算弯曲刚度矩阵将使用个高斯点进行积分，而剪切刚度矩阵则仅使用进行积分。

来源：数值分析与有限元编程

高斯消元法的新思路

用Python做数据分析之后得到启发,对高斯消元法进行改进。大致思路仍然是将矩阵化为上三角矩阵，只是算法不同。大致思路仍然是将矩阵化为上三角矩阵，只是算法不同。假设线性方程组的增广矩阵为其中，系数矩阵可逆，确保方程组有唯一解。接下来通过行变换，将第一行主对角元素置为1。importnumpyasnp#Ab表示增广矩阵Ab=np.array([[11,12,13,14],[21,22,23,24],[31,32,33,34]],dtype=np.float32)#操作1，第一行主对角元素置为1Ab[0:1,0:4]=1/Ab[0,0]*Ab[0:1,0:4]print(Ab)经过操作1，下面引入一种特殊的矩阵，它是由相同维数的列向量和行向量通过矩阵乘法得到，其特殊之处就是矩阵的秩为1(各行成比例)。比如要将第一列除主对角之外的元素变为0，需要借助上述的矩阵。其中这个时候想必你也知道右边的列向量和行向量是从哪来的了。importnumpyasnp#Ab表示增广矩阵Ab=np.array([[11,12,13,14],[21,22,23,24],[31,32,33,34]],dtype=np.float32)#操作1,第一行主对角元素置为1Ab[0:1,0:4]=1/Ab[0,0]*Ab[0:1,0:4]#操作2#第一步构造一个秩为1的矩阵temp1=Ab[0:1,0:4]#提取第一行temp2=Ab[1:3,0:1]##提取第一列的第二行和第三行temp3=np.zeros((1,1),dtype=np.float32)temp4=np.vstack((temp3,temp2))#扩展r1=np.dot(temp4,temp1)#r1是秩为1的矩阵Ab=Ab-r1print(Ab)输出结果为：[[1.1.09090911.18181821.2727273][0.-0.90909195-1.8181839-2.727272][0.-1.8181839-3.636364-5.454544]]重复上述操作，将第二行主对角元素置为1，第二列主对角以下的元素化为0。第三行主对角元素置为1时结束，系数矩阵就成了上三角矩阵，用回代法就得到了方程组的解。以下是完整的代码importnumpyasnpdeftriang_solver(Ab,n):foriinrange(n):Ab[i:i+1,i:n+1]=1/Ab[i,i]*Ab[i:i+1,i:n+1]Ab[i+1:n,i:n+1]=Ab[i+1:n,i:n+1]-np.matmul(Ab[i+1:n,i:i+1],Ab[i:i+1,i:n+1])x=np.zeros(n)#回代foriinrange(n,0,-1):x[i-1]=Ab[i-1,n]-np.dot(Ab[i-1,i:n],x[i:n])returnxif__name__=="__main__":#Ab表示增广矩阵Ab=np.array([[4,2,-1,5],[1,4,1,12],[2,-1,4,12]],dtype=np.float32)x=triang_solver(Ab,3)print(x)来源：数值分析与有限元编程