计算不收敛？接触问题有限元Matlab编程求解攻略(下)

详细信息

作者优秀: 优秀教师/意见领袖/博士学历/特邀专家/独家讲师
平台推荐: 内容稀缺

11月前浏览1061

导读：本次课程主要围绕接触问题的有限元Matlab编程求解进行介绍，主要分为两个部分，一是接触理论有限元求解的基本原理和算法，二是利用Matlab对接触问题有限元求解算法进行实现，并将计算结果与商业有限元计算结果进行对比，验证了程序算法的可靠性。

一、接触问题有限元求解原理及算法（Part1）

本章节内容笔者在《计算不收敛？接触问题有限元Matlab编程求解攻略(上)》已经介绍，今天给大家继续讲第二部分内容。见下文：

二、Matlab编程实现接触问题的有限元求解（Part2）

根据第一节介绍的接触问题有限元算法，本节利用Matlab对接触问题有限元求解算法进行实现。本文以图1所示的两物体接触的平面应力问题为对象介绍接触问题的有限元Matlab编程求解过程。图1中上部块体上表面受均布荷载作用，左图为Matlab计算结果，有图为abaqus模拟结果，经比较，二者应力分布形式基本一致，但因采用算法的不同计算结果略有差异。本文的接触问题的有限元算法在计算精度方面还有待提高，但作为接触问题有限元编程入门来讲还是有一定借鉴意义，姑且抛砖引玉吧，欢迎大家批评指正！

图1 两物体接触的平面应力问题Matlab有限元编程求解

图2 两物体接触的平面应力问题Abaqus求解

本程序基于四节点平面应力单元的代码进行开发，四节点平面应力单元相关的理论介绍和代码实现可参考文章《四节点八节点四边形单元悬臂梁的Matlab有限元编程》或仿真秀视频教程《Matlab有限元编程从入门到精通20讲》中对应的四节点四边形单元的课程。

《Matlab有限元编程从入门到精通20讲》

1、主程序main.m

Main.m是主程序，部分代码展示如下。由于接触问题是两个独立的物体相互接触的过程，所以主程序对于两物体的节点单元信息、边界条件以及各自的弹性刚度矩阵都是分开定义。其中Geometry()返回节点单元信息(P,E)、荷载施加节点(q0_node)和边界约束施加节点(Boundary)，通过Patch=1或2分别定义上下物体的上述信息；然后主程序通过while循环，分100个荷载步逐步施加荷载，并计算每个荷载步对应的节点位移d，然后通过d_final=d_final+d与上一步计算的节点位移d叠加。节点位移的计算也是通过Kd=F求得，而这里的K即本篇part1中提到的切线刚度矩阵，包含了两物体各自的弹性刚度和接触刚度两部分。两物体各自的弹性刚度矩阵通过Linear_Elasticity()函数得到，接触刚度矩阵通过Contact_Formualtion()。Linear_Elasticity()函数可参考文章《四节点八节点四边形单元悬臂梁的Matlab有限元编程》或仿真秀视频教程《Matlab有限元编程从入门到精通20讲》中对应的四节点四边形单元的课程。本篇重点讲一下，接触刚度矩阵的编程求解方法。

Patch=1;%底部块体[P1,E1,q0_node1,Boundary1] = Geometry(Patch, 0, 0);Patch=2;%上部块体[P2,E2,q0_node2,Boundary2] = Geometry(Patch, 0, 0);%定义节点和单元Node=[P1;P2];Element=[E1;E2+size(P1,1)];time=1;d_final=zeros(2*size(Node,1),1);while(time <100)%建立下部块体刚度阵Patch=1;[K_Patch_1, F_Patch_1] = Linear_Elasticity(P1, E1, q0_node1,Boundary1,Patch);%建立上部块体刚度阵Patch=2;[K_Patch_2, F_Patch_2] = Linear_Elasticity(P2, E2,q0_node2,Boundary2, Patch);%将两个全局patch的刚度矩阵和荷载矩阵组装为一体K_Global=[K_Patch_1,zeros(size(K_Patch_1,1),size(K_Patch_2,2)) ;zeros(size(K_Patch_2, 1),size(K_Patch_1,2)) ,K_Patch_2];F_Global=[F_Patch_1; F_Patch_2];%考虑接触约束Penalty=1e7;%惩罚因子[K_Contact RC_force] = Contact_Formualtion(P1, P2,E1,E2, Penalty, dis_1, dis_2);K_Global=K_Global+K_Contact;d=linsolve(K_Global, F_Global);d_final=d_final+d;time=time+1;end

2、接触刚度矩阵求解程序Contact_Formualtion()

接触刚度矩阵计算是求解接触问题的核心，具体算法如下:

1）搜索接触：在预定义的主面和从面之间遍历所有可能的node-to-segement组合；

2）判断是否施加接触约束，需要同时满足两个条件：条件一，node与segement的距离gN<0 ；条件二， node在segement上投影点的参数坐标满足；

3）根据公式19计算刚度矩阵。

在实际问题中，并非所有这些接触面上的点都发生接触力学行为，只是定义一个可能发生接触的节点集合，在计算中只在此集合中进行搜索判断是否发生接触。接触算法的核心就是搜所可能发生接触的物体确定实际被激活的接触面。一旦激活接触，就必须建立基于穿透函数 (5) 的局部接触约束条件，进而产生接触刚度。接触刚度矩阵求解函数Contact_Formualtion()首先通过接触面几何定义函数Contact_Geometry()返回主面和从面的节点信息和对应的自由度的全局索引值，索引信息方便被激活的接触对的刚度矩阵组装至整体接触刚度矩阵中。然后程序对所有从节点和主面线段（segment）的组合进行搜索，通过cntelm2d()函数来判断是否形成接触对，同时判断接触对是否被激活从而贡献刚度，并返回对应的接触对刚度矩阵和接触力。接着将接触对刚度矩阵组装至全局接触刚度矩阵中。因此，cntelm2d()函数对于整体刚度矩阵的计算至关重要，接下来一小节就介绍接触对刚度矩阵的计算程序。

function [K_Contact RC_force] = Contact_Formualtion(P1, P2, E1,E2,Penalty, dis_1, dis_2) K_Contact=zeros(2*(size(P1,1)+size(P2,1))); %initialize the Output;    SD=2;%空间维度%slave接触面的节点及对应的全局索引[ P_s, patch_s, Global_DOF_ID_s]=  Contact_Geometry(0, dis_1, dis_2);  [ P_m, patch_m, Global_DOF_ID_m]=  Contact_Geometry( 1, dis_1, dis_2);  RC_force=zeros(size(P_s,1)+size(P_m,1),2);nnd_s=size(P_s,1);%节点数ne_s=nnd_s-1;%单元数nnd_m=size(P_m,1);%节点数ne_m=nnd_m-1;%单元数%对接触面进行搜索得到node-to-segement接触对for n_i=1:nnd_s %从面节点进行循环nodeID_s=Global_DOF_ID_s(2*n_i-1:2*n_i);  for e_i=1:ne_m %主面segement进行循环ID_m1=  Global_DOF_ID_m(2*e_i-1:2*e_i);ID_m2=  Global_DOF_ID_m(2*(e_i+1)-1:2*(e_i+1));iDof=[ID_s,ID_m1,ID_m2];xs=P_s(n_i,1);ys=P_s(n_i,2);x1=P_m(e_i,1);x2=P_m(e_i+1,1);y1=P_m(e_i,2);y2=P_m(e_i+1,2);ELXY=[xs,x1,x2;ys,y1,y2];            [FORCE, STIFF] = cntelm2d(Penalty, ELXY);%返回接触刚度矩阵和接触力if isempty(STIFF)continue;end            for row=1:size(STIFF,1)for col=1:size(STIFF,2)r=iDof(row);c=iDof(col);                 K_Contact(r,c)=K_Contact(r,c)+STIFF(row,col);endendRC_force(n_i,:)=FORCE(1:2);RC_force(5+e_i,:)=FORCE(3:4);RC_force(5+e_i+1,:)=FORCE(5:6);endend   end

3、接触对刚度矩阵求解程序cntelm2d()

接触对刚度矩阵可根据part1中的公式19进行计算，具体的实现程序如下所示。程序首先求得segement线段的单位切向量和法向量，根据公式3计算出从节点x1到主面segement的距离gNs，程序中用GAPN表示；接下来就根据gNs的正负来判断是否发生穿透，如果为负值，说明发生node to segement接触，根据公式11计算从节点x1对应的主面上投影接触点的参数坐标ξ，程序中用alpha表示。然后通过判断alpha是否在[0,1]之间来判断从节点是否落在主线段segement上，如果落在主线段上，则根据公式1计算Ns，程序中用NN表示，然后根据公式19计算接触对刚度矩阵。这里需要补充介绍接触力的求解方法。对于从节点在接触点ξ处的接触力pN根据公式13给出，该接触力分配到主线段segement上的两节点力为

（20）

将接触对的接触力卸载一起可以得到整体接触对的接触力向量为

(21)

function [FORCE, STIFF] = cntelm2d(OMEGAN, ELXY)%***********************************************************************% SEARCH CONTACT POINT AND RETURN STIFFNESS AND RESIDUAL FORCE % IF CONTACTED FOR NORMAL CONTACT%***********************************************************************%ZERO = 0.D0; ONE = 1.D0; EPS = 1.E-6; P05 = 0.05; FORCE=[]; STIFF=[];XT = ELXY(:,3)-ELXY(:,2); %切向量XLEN = norm(XT);%x1 x2的距离if XLEN < EPS, return; end%% UNIT NORMAL AND TANGENTIAL VECTORXT = XT/XLEN;XN = [-XT(2); XT(1)];%% NORMAL GAP FUNCTION Gn = (X_s - X_1).NGAPN = (ELXY(:,1)-ELXY(:,2))'*XN;%公式3%% CHECK IMPENETRATION CONDITIONif (GAPN >= ZERO) || (GAPN <= -XLEN), return; end%% NATURAL COORDINATE AT CONTACT POINTALPHA = (ELXY(:,1) - ELXY(:,2))'*XT/XLEN;%公式11%% OUT OF SEGMENTif (ALPHA > ONE+P05) || (ALPHA < -P05), return; end%% CONTACT OCCURS IN THIS SEGMENTXLAMBN = -OMEGAN*GAPN; %公式%% DEFINE VECTORSNN = [XN; -(ONE-ALPHA)*XN; -ALPHA*XN];%2*3 公式16% TT = [XT; -(ONE-ALPHA)*XT; -ALPHA*XT];% PP = [ZERO; ZERO; -XN; XN];% QQ = [ZERO; ZERO; -XT; XT];CN = NN;%去掉-之后的差别不大% CN = NN - GAPN*QQ/XLEN;% CT = TT + GAPN*PP/XLEN;%% CONTACT FORCEFORCE = XLAMBN*CN;% FORM STIFFNESSSTIFF = OMEGAN*(CN*CN');end