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

1问题描述某周边简支非均匀的矩形（或圆形）板在均布载荷作用下挠度过大。结合实际，提出集中改进设计方案，并进行对比分析。2.问题分析不均匀板有两种主要的情况，结构不均匀和材料不均匀，结构不均匀是指板的厚度不是常量，材料不均匀体现在板的弹性模量和泊松比是变化的。另外，有的板可以是以上两种情况的混合情形。不均匀板与均匀板的有限元问题有哪些差别呢？下面从均匀板问题推导出非均匀板有限元问题的解决方法。2.1应力应变先以结构不均匀板为例来讨论。假设一矩形板长为2，宽为2，厚度沿x，y不均匀，由一函数h,hxy描述，但仍然符合薄板假设。对于均匀板，显然h是一个常数。设挠度为=x,y，则板内应变向量可以表示为2222211==z12xxyyxyxyxzyxy应力应变关系为1pzD弯矩扭矩矩阵h,2h,2xyxyMzdz这里就体现出不均匀板和均匀板的区别了。积分完毕后，可以得到1MD其中薄板的弯曲系数矩阵3210,10121001/2EhxyD是关于薄板总体坐标的函数，所以对各个分单元都是不同的。各单元的弯曲系数矩阵可以采用单元中心处的代替。那么就可以得出一系列的弯曲系数矩阵Dei。如果单元划分得足够细，是可以代替真实解的。2.2单元分析可以将板分为边长为0.25的矩形小单元，每一个单元都是一样的。对于任何一个单元的节点，都有3项独立的位移iiixiiyii位移模式22312345672233389101112,wxyxyxxyyxxyxyyxyxy形状函数矩阵是一个112的行向量,klmnNxyNNNN其中222222222222222211128111111iiiiiiiiiiiiixxyyxxyyxyNabababxxyyyyxxyyxxyxabbaba,,,iklmn单元刚度矩阵1212eeTSkBDBdxdy很明显，积分式中包含了弹性系数矩阵，而不同单元的弹性系数矩阵是不同的，所以，即便单元划分相同，得到的单元刚度矩阵也不同。对于均匀板，相同形式的单元，刚度矩阵是相同的。均匀与非均匀的差别，完全体现在弹性系数矩阵上。但是非均匀板的一些结果可以间接地运用。矩形单元四节点单元刚度矩阵是一个规律性很强的对称矩阵。矩阵中待求的独立元素只有21个。1425637101111084231195632121516172021115132018421614211956317202112151671011201815131082119161400(,)360(1)000000000kkkkkkkkkkkkkkkkkkkEhxykkkkkkkabkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk对称1421195630kkkkkkk其中22221222222232425622227222822921021121630/30/8840884031230/31230/3021630/15/8820222031215/3330/kbaabkbbakaabkbbabkaabakabkbaabkbbakabbkbbabkaaba22221222213222142152162222172221822219220221k21615/15/221022103315/3315/21615/30/222088203330/31215/baabkbbakaabkbbabkaabakbaabkbbakaabkbbabkaaba如果以单元中心点的参数代替单元的参数，那么非均匀的薄板单元刚度矩阵与均匀板差别不大，即各个单元的刚度矩阵都需要利用单元中心的参数来代替。这个结论对材料不均匀板也是实用的。下面运用这种近似来求解非均匀薄板问题。3问题求解求解的模型仍然是边长为2的正方形薄板，材料的弹性模量E25000000000Pa，泊松比0.3。图1薄板单元3.1均匀板将板分为均匀的正方形单元44，88等，可以得到精度不等的数值解。厚度为0.035的薄板，在载荷100000000Pqa作用下可以计算出板的挠度。可以看出挠度是完全对称的，并且最大挠度出现在板的几何中心上。图2均匀板均布载荷下变形（整体变形和中线变形）3.2结构非均匀板通过一个均匀板的计算，可以看出对称均匀板受到均布载荷后，挠度也是完全对称的。但是对于非均匀板，受到均布载荷挠度如何分布，最大挠度在什么地方，是最关键的问题。只有找到最大挠度的位置，以及最大挠度是多少才能采取相应的解决办法。下面解决一个结构非均匀薄板。设板在x方向上厚度从0.02线性变化到0.05（平均厚度为0.035，与上一均匀板厚度相同），厚度关于中面对称，其余参数和前面的均匀板相同。结构如图。图3不均匀板模型同样，将板分为88（64个）正方形单元，在均布载荷100000000Pqa作用下，板的形变已经完全不一样了。最大挠度不再出现在几何中心处。很显然，挠度最大处向板较薄的一侧推移了。值得注意的是，这种不均匀薄板在薄弱环节挠度往往增大得很快，平均厚度相同的薄板受到相同载荷，不均匀板的挠度往往大得多。图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];//k