板壳结构matlab有限元编程（二）：Mindlin厚板单元程序代码讲解
- 作者优秀
- 优秀教师/意见领袖/博士学历/特邀专家/独家讲师
- 平台推荐
- 内容稀缺
导读：板壳结构是工程上常见的结构形式，应用也非常广泛，因此有限元仿真计算中板壳单元也是非常重要的单元。这里的板壳单元严格意义上可分为板单元和壳单元，板单元自由度包含挠度、法线转动，而壳单元自由度要比板单元多出两个面内的平动自由度。因此商业软件中并不单独设置板单元。但是从有限元理论发展的角度讲，有必要将板单元与壳单元区分开来，这样也有助于理论的理解。本系列博文将系统介绍板壳单元的有限元理论基础及Matlab编程实现，具体包括Kirchhoff薄板理论与Mindlin/Reissner 板壳理论，这是进行有限元离散的基础，基于上述理论建立的板壳结构的三大类方程，几何方程、物理方程和平衡方程，在三大类方程基础上离散得到的有限元相关方程，此外针对的物理问题包括了板壳结构的静力求解和模态分析的内容。
图3
图4
一、有限元方程的建立
基于Mindlin板理论的剪切变形板单元，与上一篇博文《板壳结构matlab有限元编程（一）：板单元》中所讲的基于Kirchhoff经典薄板弯曲理论的薄板单元最大区别就是，垂直于板中面的平面在板变形后将不再垂直于中面, 所以需要同时考虑板的剪切能量和弯曲能量。可以类比欧拉梁和铁木辛柯梁，将其关系拓展至二维，就是Kirchhoff经典板元和Mindin剪切板元的关系。
固体力学有限元方程的建立过程通常基于变分原理（如最小位能原理、虚功原理）或加权残值法（如伽辽金方法）。以下以最小位能原理为例，详细描述Mindlin板单元的有限元方程的推导过程：
对于线弹性固体，总势能泛函
由应变能U和外力功 W 组成：
因为考虑了剪切能量和弯曲能量，Mindlin板的应变能表达式为：
其中弯曲应力和应变为
剪切应力和应变为
是考虑实际的剪切应变沿厚度方向非均匀分布而引入的校正系数,为 5/6.将弯曲和剪切应力应变关系带入式(7-30), 可得
其中
板内的位移可以表示为
式中,
和
为板中面绕
轴和
轴的转角. 假设板中面没有板内位移,则对剪切变形板, 有
式中,
为由横向剪切变形造成的角度.
由于位能表达式中的
和
是各自独立的,所以对其进行的插值也可以各自独立进行. 这样其插值函数只要求
的连续性. 可采用适用于四边形单元的双线性等参形函数,
和可以分别插值表示为
式中,
为单元节点个数.其中
为双线性等参形函数。
下面以双线性等参形函数为例来说明进一步的计算过程. 使用高阶形函数也可以用类似的推导过程完成计算. 弯曲和剪切应变可以从位移中计算出
其中
将式(7-46)和式(7-47)带入式(7-45)可得
根据变分原理和弱形式改写（可参考任意一本有限元教材），从上述能量方程即可得弯曲板的单元刚度矩阵
式中, h 为板的厚度.
当板的厚度相对于其边长很小的时候, 剪切能量将远大于弯曲能量, 这种现象称为剪切锁死. 为避免此问题, 一般可以选用选择积分或者缩减积分的方法. 研究结果还表明, 在避免剪切锁死方面, Lagrange 单元通常具有较好的性能.
二、Mindlin剪切板单元代码实现
单元开发的核心是单元刚度矩阵的建立，根据上一节有限元公式的推导过程，得到Mindlin剪切板单元的Matlab代码如下：
function[K]=StiffnessMatrix(material,mesh,quadrature,ShapeOption)K=zeros(mesh.nn*3);D_f=material.E/(1-material.v^2).*[1 material.v 0material.v 1 00 0 (1-material.v)/2];%bendingD_t=[material.G 00 material.G];%shearalpha=5/6;% shear stress non-uniform modification factorbeta=material.t^3/12; %integral(z^2,-t/2,t/2)shape_order=4;for iel=1:mesh.neKElem=zeros(shape_order*3);for ig=1:size(quadrature.points,1)shapeFunction=ShapeFunction(quadrature.points(ig,1),quadrature.points(ig,2),ShapeOption);%型函数在积分点处的值elemCoordinates=[mesh.Nid(mesh.Eid(iel,2:end),2), mesh.Nid(mesh.Eid(iel,2:end),3)];%四个节点的坐标jacobian=Jacobian(shapeFunction,elemCoordinates);%2*2矩阵Nd=jacobian.globalDerivatives;%型函数的物理坐标导数B_f=zeros(3,shape_order*3);%bendingB_t=zeros(2,shape_order*3);%shearfor ib=1:shape_ordertemp=[0 Nd(ib,1) 0 %0 0 Nd(ib,2)0 Nd(ib,2) Nd(ib,1)];B_f(:,ib*3-2:ib*3)=temp;%bendingtemp=[Nd(ib,1) shapeFunction.fun(ib) 0Nd(ib,2) 0 shapeFunction.fun(ib)];B_t(:,ib*3-2:ib*3)=temp; %shearendKElem=KElem+quadrature.weights(ig)*(beta*B_f'*D_f*B_f+alpha*material.t*B_t'*D_t*B_t)*det(jacobian.matrix);% 弯曲刚度+剪切刚度endntot=mesh.Eid(iel,2:end)*3-2;%index of DOF%assemble to golobal stiffness matrixfor j=1:numel(ntot)for k=1:numel(ntot)K(ntot(j):ntot(j)+2,ntot(k):ntot(k)+2) = K(ntot(j):ntot(j)+2,ntot(k):ntot(k)+2) + KElem(j*3-2:j*3,k*3-2:k*3);endendend
shapeFunction.fun=1/4*[(1-xi)*(1-eta)(1+xi)*(1-eta)(1+xi)*(1+eta)(1-xi)*(1+eta)];shapeFunction.dfun=1/4*[-(1-eta),-(1-xi)1-eta,-(1+xi)1+eta,1+xi-(1+eta),1-xi];
三、我的Matlab有限元编程精品课
推荐大家关注我的原创视频课程里面《Matlab有元编程从入门到精通》目前加餐到第35期。我还会持续更新，强烈推荐学习者订阅。
可开电子发票，赠送答疑专栏
本课程为matlab有限元编程专题课，课程主要以案例的形式进行讲解，中间会穿插案例中所涉及到的有限元基本理论，案例不局限于力学问题的有限元求解，还会涉及传热学、电学等问题的有限元求解。
因为固体力学领域我最熟悉，所以我们从固体力学开始，所涉及的单元有杆单元，梁单元，平面三角形单元，薄板单元，厚板单元，四面体实体单元等等，力学问题有静力学问题，也有动力学问题，后期还会涉及材料非线性、几何非线性、接触非线性等非线性问题，内容丰富，不断更新完善。
此外，笔者为所有订阅用户提供知识圈答疑服务和VIP用户交流群。并附赠课程相关资料等（平台支持自行开具电子发票）。
希望大家关注笔者仿真秀专栏，点击文章末尾阅读原文即可关注，大家学习Matlab编程过程中，强烈推荐大家学习以下资料（在文章末尾点赞或再看，截图发到本公 众 号，回复 SimPC ，我会24小时发给你以下资料）。
