有限元分析薄板挠度(附C程序)

1问题描述某周边简支非均匀的矩形（或圆形）板在均布载荷作用下挠度过大。

结合实际，提出集中改进设计方案，并进行对比分析。

2.问题分析不均匀板有两种主要的情况，结构不均匀和材料不均匀，结构不均匀是指板的厚度不是常量，材料不均匀体现在板的弹性模量和泊松比是变化的。

另外，有的板可以是以上两种情况的混合情形。

不均匀板与均匀板的有限元问题有哪些差别呢？下面从均匀板问题推导出非均匀板有限元问题的解决方法。

2.1应力应变先以结构不均匀板为例来讨论。

假设一矩形板长为2，宽为2，厚度沿x ，y 不均匀，由一函数()h ,h x y =描述，但仍然符合薄板假设。

对于均匀板，显然h 是一个常数。

设挠度为()=x,y ωω，则板内应变向量可以表示为{}2222211==z 12x x y y xy xy x z y x y ρεεεωεγγ⎧⎫⎧⎫∂⎪⎪⎪⎪∂⎪⎪⎪⎪⎧⎫⎪⎪⎪⎪∂⎪⎪⎪⎪⎪⎪=-⎨⎬⎨⎬⎨⎬∂⎪⎪⎪⎪⎪⎪⎩⎭⎪⎪⎪⎪∂⎪⎪⎪⎪∂∂⎪⎪⎪⎪⎩⎭⎩⎭应力应变关系为{}1p z D σρ⎧⎫⎡⎤=⎨⎬⎣⎦⎩⎭弯矩扭矩矩阵{}{}()()h ,2h ,2x y x y M zdz σ-=⎰这里就体现出不均匀板和均匀板的区别了。

积分完毕后，可以得到{}[]1M D ρ⎧⎫=⎨⎬⎩⎭其中薄板的弯曲系数矩阵[]()()()321,1012101/2Eh x y D μμμμ⎡⎤⎢⎥=⎢⎥-⎢⎥-⎣⎦是关于薄板总体坐标的函数，所以对各个分单元都是不同的。

各单元的弯曲系数矩阵可以采用单元中心处的代替。

那么就可以得出一系列的弯曲系数矩阵[]D ei 。

如果单元划分得足够细，是可以代替真实解的。

2.2单元分析可以将板分为边长为0.25的矩形小单元，每一个单元都是一样的。

对于任何一个单元的节点，都有3项独立的位移{}i i i xi i yi i w w w y w x δθθ⎧⎫⎪⎪⎪⎪⎧⎫⎪⎪⎛⎫∂⎪⎪⎪⎪==⎨⎬⎨⎬⎪∂⎝⎭⎪⎪⎪⎪⎩⎭⎪⎪∂⎛⎫⎪⎪- ⎪∂⎪⎪⎝⎭⎩⎭位移模式()22312345672233389101112,w x y x y x xy y x x y xy y x y xy αααααααααααα=+++++++++++形状函数矩阵是一个112⨯的行向量()[],kl mn N x y N N N N =⎡⎤⎣⎦其中222222222222222211128111111i i i i i i i i i i i i i x x y y x x y y x y N a b a b a b x x y y y y x x y y x x y x a b b a b a ⎛⎫⎡⎛⎫⎛⎫=++++--⎡⎤ ⎪⎪⎪⎣⎦⎢⎝⎭⎝⎭⎣⎝⎭⎤⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫++--++-⎥ ⎪⎪ ⎪ ⎪ ⎪⎪⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭⎥⎦(),,,i k l m n =单元刚度矩阵[][][][]1212ee TS k B D B dxdy ⨯=⎰很明显，积分式中包含了弹性系数矩阵，而不同单元的弹性系数矩阵是不同的，所以，即便单元划分相同，得到的单元刚度矩阵也不同。

对于均匀板，相同形式的单元，刚度矩阵是相同的。

均匀与非均匀的差别，完全体现在弹性系数矩阵上。

但是非均匀板的一些结果可以间接地运用。

矩形单元四节点单元刚度矩阵是一个规律性很强的对称矩阵。

矩阵中待求的独立元素只有21个。

1425637101111084231195632121516172021115132018421614211956317202112151671011201815131082119161400(,)360(1)00000000k k k k k k k k k k k k k k k k k k k Eh x y k k k k k k k ab k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k μ-----------------对称1421195630k k k k k k k ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥-⎣⎦如果以单元中心点的参数代替单元的参数，那么非均匀的薄板单元刚度矩阵与均匀板差别不大，即各个单元的刚度矩阵都需要利用单元中心的参数来代替。

这个结论对材料不均匀板也是实用的。

下面运用这种近似来求解非均匀薄板问题。

3问题求解求解的模型仍然是边长为2的正方形薄板，材料的弹性模量E 25000000000Pa =，泊松比0.3μ=。

图1 薄板单元3.1 均匀板将板分为均匀的正方形单元44⨯，88⨯等，可以得到精度不等的数值解。

厚度为0.035的薄板，在载荷100000000P q a =作用下可以计算出板的挠度。

可以看出挠度是完全对称的，并且最大挠度出现在板的几何中心上。

图2 均匀板均布载荷下变形（整体变形和中线变形）3.2 结构非均匀板通过一个均匀板的计算，可以看出对称均匀板受到均布载荷后，挠度也是完全对称的。

但是对于非均匀板，受到均布载荷挠度如何分布，最大挠度在什么地方，是最关键的问题。

只有找到最大挠度的位置，以及最大挠度是多少才能采取相应的解决办法。

下面解决一个结构非均匀薄板。

设板在x 方向上厚度从0.02线性变化到0.05（平均厚度为0.035，与上一均匀板厚度相同），厚度关于中面对称，其余参数和前面的均匀板相同。

结构如图。

图3 不均匀板模型同样，将板分为88⨯（64个）正方形单元，在均布载荷100000000P q a =作用下，板的形变已经完全不一样了。

最大挠度不再出现在几何中心处。

很显然，挠度最大处向板较薄的一侧推移了。

值得注意的是，这种不均匀薄板在薄弱环节挠度往往增大得很快，平均厚度相同的薄板受到相同载荷，不均匀板的挠度往往大得多。

图4 不均匀板在均布载荷下变形但是新的问题出现了。

实际情况确实是挠度最大处向板较薄的方向推移了，可是目前的网格还是比较粗，所有的位移量都是针对节点计算的，并不能精确反映出现的位置。

所以，要想较精确地得到具体位置，网格需要加密。

但针对现有的精度，我们可以提出几种解决大挠度的方案。

1）在最大挠度处添加支撑（固支或简支）如果在最大挠度处添加支撑，那么该处的挠度就已固定。

简支时，转角自由；固支时，转角为0或者给定。

假如限制最大挠度处挠度为0.1。

下面分别就固支和简支情况给出数值解。

图5 薄弱环节采用简支图6 薄弱环节采用固支2）在最大挠度出添加固支肋条在最大挠度出使挠度为0，可以得到下面的情形。

可以通过不断改变添加的位置找到最合适的方法。

显然，图中不是最好的方法。

图7 薄弱处添加肋条3）施加外力例如给薄板最大挠度出施加一集中载荷500000N ，方向与均布载荷相反。

可以看到最大挠度（0.12）也有明显减小。

图8 薄弱处施加外力4）改变材质（材质分布）不均匀薄板本来就可以是材质的不均匀。

当改变材料性质可行时，那么这也是一种好的方法。

下面通过改变材料的性质改变薄板的形变。

首先应该清楚，在薄板较薄的位置，结构的刚度较小。

如果结构不能改变，那么需要适当地调整材料的性质。

采用不同的材料，或者对材料进行不同的处理。

那么，下面我们让薄板较薄的半边弹性模量为275Gpa，看看会有什么发生。

图9 改变材质（分布）通过对比发现，薄板的最大挠度随着较薄区域的刚度的增加，最大挠度明显减小；如果全域都改变材料，使用弹性模量较大的材料，挠度会更小，但没有通过直接增加薄区域弹性模量的办法有效。

4总结1. 不均匀板主要有两种情况，结果不均匀和材质不均匀；2. 文中展示了如何求解不均匀薄板的有限元问题。

主要是通过有限单元的平均参数代替单元参数。

理论上，当单元足够小的时候，是可以逼近精确解的。

不论是结构不均匀还是材质不均匀，都是适用；3.解决挠度过大的方案，仍然是通过改变结构、改变外部条件和改变材质（材质分布）等方法来解决。

求解的方法就是使用不均匀板的有限元方法；4.这种方法不能给出一个优化的方法，只能是校核；5.求解过程是基于C语言完成的。

单元数目不多时,可以直接计算.当单元数目很多，矩阵的维数很大，通过手算是不可能的。

所有的工作都必须借助计算机完成。

事实上，用有限元解决问题，单元数目都非常多。

最初的将整个过程用C语言完成，是耗费一定时间的。

当整个程序建立之后，使用分析就很方便了。

由于编程序和调程序花费了一定时间，所以分析上显得很不足。

但是通过完成作业，很清楚地认识到有限元理论和应用的区别，以及可以熟练地掌握有限元方法的求解过程、编程计算和问题解决。

尽管解决的方案简单，但对自己提高很多。

（后面附部分程序代码）五、附录1、单元节点记录EL[NE][4]={1,2,11,10,2,3,12,11,3,4,13,12,4,5,14,13,5,6,15,14,6,7, 16,15,7,8,17,16,8,9,18,17,10,11,20,19,11,12,21,20,12,13,22,21,13,14,23,22,14,15,24,23,15 ,16,25,24,16,17,26,25,17,18,27,26,19,20,29,28,20,21,30,29,21,22,31,30,22,23,32,31,23,24,33,32,24 ,25,34,33,25,26,35,34,26,27,36,35,28,29,38,37,29,30,39,38,30,31,40,39,31,32,41,40,32,33,42,41,33 ,34,43,42,34,35,44,43,35,36,45,44,37,38,47,46,38,39,48,47,39,40,49,48,40,41,50,49,41,42,51,50,42 ,43,52,51,43,44,53,52,44,45,54,53,46,47,56,55,47,48,57,56,48,49,58,57,49,50,59,58,50,51,60,59,51 ,52,61,60,52,53,62,61,53,54,63,62,55,56,65,64,56,57,66,65,57,58,67,66,58,59,68,67,59,60,69,68,60 ,61,70,69,61,62,71,70,62,63,72,71,64,65,74,73,65,66,75,74,66,67,76,75,67,68,77,76,68,69,78,77,69 ,70,79,78,70,71,80,79,71,72,81,80};2、刚度矩阵EK0=E*h*h*h/(360*a*b*(1-um*um));for(i=0;i<=11;i=i+3){EK[i][i]=21-6*um+30*a*a/b/b+30*b*b/a/a;EK[i+1][i+1]=8*b*b-8*um*b*b+40*a*a;EK[i+2][i+2]=8*a*a-8*um*a*a+40*b*b;//k1—k3}EK[1][0]=3*b+12*um*b+30*a*a/b;EK[4][3]=EK[1][0];EK[7][6]=-EK[4][3];EK[10][9]=EK[7][6];//k4EK[2][0]=-(3*a+12*um*a+30*b*b/a);EK[5][3]=-EK[2][0];EK[8][6]=EK[5][3];EK[11][9]=EK[2][0];//k5EK[2][1]=-30*um*a*b;EK[5][4]=-EK[2][1];EK[8][7]=EK[2][1];EK[11][10]=EK[5][4];//k6EK[3][0]=-21+6*um-30*b*b/a/a+15*a*a/b/b;EK[9][6]=EK[3][0];//k7EK[4][1]=-8*b*b+8*um*b*b+20*a*a;EK[10][7]=EK[4][1];//k8EK[5][2]=-2*a*a+2*um*a*a+20*b*b;EK[11][8]=EK[5][2];//k9EK[3][1]=-3*b-12*um*b+15*a*a/b;EK[4][0]=EK[3][1];EK[9][7]=-EK[3][1];EK[10][6]=EK[9][7];//k10EK[3][2]=3*a-3*um*a+30*b*b/a;EK[9][8]=-EK[3][2];EK[5][0]=EK[9][8];EK[11][6]=EK[3][2];//k11EK[6][0]=21-6*um-15*b*b/a/a-15*a*a/b/b; EK[9][3]=EK[6][0];//k12EK[7][1]=2*b*b-2*um*b*b+10*a*a;EK[10][4]=EK[7][1];//k13EK[8][2]=2*a*a-2*um*a*a+10*b*b;EK[11][5]=EK[8][2];//k14EK[7][0]=-3*b+3*um*b+15*a*a/b;EK[10][3]=EK[7][0];EK[6][1]=-EK[10][3];EK[9][4]=EK[6][1];//k15EK[6][2]=-3*a+3*um*a+15*b*b/a;EK[9][5]=-EK[6][2];EK[8][0]=EK[9][5];EK[11][3]=EK[6][2];//k16EK[6][3]=-21+6*um+15+b*b/a/a-30*a*a/b/b; EK[9][0]=EK[6][3];//k17EK[7][4]=-2*b*b+2*um*b*b+20*a*a;EK[10][1]=EK[7][4];//k18EK[11][2]=-8*a*a+8*um*a*a+20*b*b;EK[8][5]=EK[11][2];//k19EK[10][0]=3*b-3*um*b+30*a*a/b;EK[9][1]=-EK[10][0];EK[7][3]=EK[10][0];EK[6][4]=EK[9][1];//k20EK[6][5]=-3*a-12*um*a+15*b*b/a;EK[8][3]=EK[6][5];EK[9][2]=-EK[8][3];EK[11][0]=EK[9][2];//k21for(i=0;i<=11;i++){for(j=i;j<=11;j++){EK[i][j]=EK[j][i];}}for(i=0;i<=11;i++){for(j=0;j<=11;j++){EK[i][j]=EK[i][j]*EK0;}}3、刚度矩阵叠加double SK[NN][NN];//整体刚度矩阵243*243for(i=0;i<=NN-1;i++){for(j=0;j<=NN-1;j++){SK[i][j]=0;}}for(i=0;i<=NE-1;i++){for(j=0;j<=3;j++){for(k=3*EL[i][j],m=3*j;k<=3*EL[i][j]+2,m<=3*j+2;k++,m++) {for(l=3*EL[i][0],n=0;l<=3*EL[i][0]+2,n<=2;l++,n++){SK[k][l]=SK[k][l]+EK[m][n];}for(l=3*EL[i][1],n=3;l<=3*EL[i][1]+2,n<=5;l++,n++){SK[k][l]=SK[k][l]+EK[m][n];}for(l=3*EL[i][2],n=6;l<=3*EL[i][2]+2,n<=8;l++,n++){SK[k][l]=SK[k][l]+EK[m][n];}for(l=3*EL[i][3],n=9;l<=3*EL[i][3]+2,n<=11;l++,n++){SK[k][l]=SK[k][l]+EK[m][n];}}}}4、高斯-赛德尔迭代double A=0;do{R=L;L=0;for(j=1;j<=NN-1;j++){A=A+SK[0][j]*u[j];}u[0]=Q[0]-A;u[0]=u[0]/SK[0][0];A=0;for(i=1;i<=NN-2;i++){for(j=0;j<i;j++){A=A+SK[i][j]*u[j];}for(j=i+1;j<=NN-1;j++){A=A+SK[i][j]*u[j];}u[i]=Q[i]-A;u[i]=u[i]/SK[i][i];A=0;}for(j=0;j<=NN-2;j++){A=A+SK[NN-1][j]*u[j];}u[NN-1]=Q[NN-1]-A;u[NN-1]=u[NN-1]/SK[NN-1][NN-1];A=0;for(i=0;i<=NN-1;i++){L=L+u[i]*u[i];}L=sqrt(L);}while (abs(L-R)>=0.0000000001);。