8节点矩形单元一致质量矩阵推导？

shipsj

       最近要做个结构动力学方面的东西，遇到点问题，现描述如下
       将一致质量矩阵通过缩放因子的方法转化为对角矩阵，推导的时候不知道出了什么问题和书上的结果不一样，现将matlab程序附上，等待高手指点。具体的可以看王勖成的书上474页
       注：此程序换成4节点9节点矩形单元的形函数时，结果和书上的一样，就是这个8节点的算出来不一样
-----------------------------------------------------------------------
clear
syms x y cofx cofe;
%corner node 1,2,3,4 shape function
Ni=1/4*(1+cofx*x)*(1+cofe*y)*(cofx*x+cofe*y-1);
N1=subs(Ni,{cofx,cofe},{-1,-1});
N2=subs(Ni,{cofx,cofe},{1,-1});
N3=subs(Ni,{cofx,cofe},{1,1});
N4=subs(Ni,{cofx,cofe},{-1,1});
%middle node 5,6,7,8 shape function
N5=1/2*(1-x^2)*(1-y);
N6=1/2*(1-y^2)*(1+x);
N7=1/2*(1-x^2)*(1+y);
N8=1/2*(1-y^2)*(1-x);
%Nsh is composed by shape functions
%Nsh=[0 N1 0 N2 0 N3 0 N4 0 N5 0 N6 0 N7 0 N8;N1 0 N2 0 N3 0 N4 0 N5 0 N6 0 N7 0 N8 0]
%We can simply vertify through u or v degree of freedom
Nsh=[N1 N2 N3 N4 N5 N6 N7 N8];
Masstr=Nsh.' * Nsh;
Massco=int(int(Masstr,x,-1,1),y,-1,1)
%vector scalar equal consistent mass matrix of the sum of each row
vm=sym(zeros(1,length(Massco)));
for i=1:length(Massco), vm(i) = sum(Massco(i,); end
vm
vs=sym(zeros(1,length(Massco)));
for i=1:length(Massco), vs(i) = Massco(i,i); end
sum(vs) .\ vs
终于有了结果啊！感谢caoer对本版的支持，tonnyw一致跟我保持探讨，pasuka最后的指点，感谢大家的参与！

yiqing

相应背景是什么啊，结构力学吗

tonnyw

Then the book has something wrong.

shipsj

这个是动力学有限元解法的基础理论。蓝科维奇和王勖成的书上结果都一样啊，结果应该不会有问题吧，节点处等效质量为总质量1/36，中点2/9

tonnyw

Could you put here the equations you are using? I think your mass matrix is okay. The problem is how you get the factor.

shipsj

vs[  2/15,  2/15,  2/15,  2/15, 32/45, 32/45, 32/45, 32/45]这个是我求出来的质量矩阵的对角线的上的量，如果质量矩阵没问题，那么说明缩放因子不是一个常量？
因为书上给定[  1/36, 1/36,  1/36, 1/36, 2/9,  2/9, 2/9,  2/9]，中节点上的等效质量是角节点的8倍，而上面的却是16/3倍，如果将vs乘以一个常量是得不出这个结果的。不知tonnyw兄有何见解
好说说的我方程，对于动力学问题，工程上往往采用达朗贝尔原理，变成一个“静力”问题，平衡方程中多了惯性力和阻尼力，如下
sigma[ij,j]+f-rho*u[i,tt]-mu*u[i,t] = 0
sigma为应力，f为体力，u为位移函数，rho质量密度，mu阻尼系数[ij,j]表述ij分量对j的偏导，[i,tt]为i分量对时间的二次偏导（不知道论坛的公式编辑语法，写的好是繁杂啊）
转化成有限元形式为
M*diff(a,t,t)+C*diff(a,t)+K*a=Q
u=Na, N为形函数，a为位移关于时间的变量
M有单元质量矩阵组装而成，C，K同理
单元质量矩阵由int(transpose(N)*rho*N,in element body），由此得出的单元质量矩阵为协调矩阵或称一致质量矩阵。现在我想做的就是把质量集中在节点上，把它变成集中质量矩阵，即矩阵对角化处理。方法有很多，常用的有两种：1）每行主元素m_ii=sum(m_ij),(i=1..n)，主元素的值为所在行所有元素的数值和，位置在对角线上；2）每行主元素等于该行主元素乘以缩放因子factor，非主元素为零。factor根据质量守恒确定。
我上面的程序，vm为第一种方法，vs为第二种方法

shipsj

怎么都没人顶啊！自己先顶一个，占个位

caoer

This is a good question....coz rare people care about the assembly of macroelement.
I would highlight this entry and let the mess talk about this.

My suggestion is that forgetting the 8-node element, which is not used in the modern
FEM package any more.

shipsj

谢谢斑竹好心提醒，我现在主要是想了解一下基本理论，为后续的软件计算垫点低，心里也好有个谱，质量集中的方法，现在用的是也不多，但在前期估算时还是很节省时间的，期待大家能一起讨论一下，把这个问题给解决了

tonnyw

8# caoer

I think 8-noded element is still used in FEM package. For instance, Adina still uses it.

caoer

8# caoer

I think 8-noded element is still used in FEM package. For instance, Adina still uses it.
tonnyw 发表于 2010-7-2 23:22

A book mentioned this 8-node 2d element is not good as 9-node, so several mainstream packages have abandon it. unfortunately I can't recall the name of
textbook. maybe Adina keeps it for some kinds of purposes or just make it work
as a tradeoff.

BTW, about the Adina, I am curious what is the development language used? fortran
or c?

tonnyw

Well. First of all, I think Adina is a very good software. I don't know what kind of coding language Adina is using. As you know, Adina is originally from NONSAP written in Fortran.

As far as I know, eight-noded element has good performance and that's why it is called serendipity element.

tonnyw

6# shipsj

I cannot find anything wrong with what you did. I cannot get the results as Zienkiewicz. This baffles me. Maybe someone else could give a try.

shipsj

13#tonnyw
这个也让我纠结了好几天了，和几个同学也讨论了一下，也无终而果

tonnyw

14# shipsj

If you can get the lumped mass matrix for nine-noded element, I see no reason you cannot get the one for eight-noded element.

shipsj

是啊，这也是我纳闷的地方，如果方法不对的话，其它模型结果也不对才对啊，这是我9-noded element的结果[ 1/36, 1/36, 1/36, 1/36, 1/9, 1/9, 1/9, 1/9, 4/9]，和书上的一样啊

shipsj

clear
syms x y cofx cofe;
%corner node 1,2,3,4 shape function
Ni=1/4*cofx*x*(1+cofx*x)*cofe*y*(1+cofe*y);
N1=subs(Ni,{cofx,cofe},{-1,-1});
N2=subs(Ni,{cofx,cofe},{1,-1});
N3=subs(Ni,{cofx,cofe},{1,1});
N4=subs(Ni,{cofx,cofe},{-1,1});
%middle node 5,6,7,8 shape function
N5=-1/2*(1-x^2)*y*(1-y);
N6=1/2*x*(1-y^2)*(1+x);
N7=1/2*(1-x^2)*y*(1+y);
N8=-1/2*x*(1-y^2)*(1-x);
N9=(1-x^2)*(1-y^2);
%Nsh is composed by shape functions
%Nsh=[0 N1 0 N2 0 N3 0 N4 0 N5 0 N6 0 N7 0 N8;N1 0 N2 0 N3 0 N4 0 N5 0 N6 0 N7 0 N8 0]
%We can simply vertify through u or v degree of freedom
Nsh=[N1 N2 N3 N4 N5 N6 N7 N8 N9];
Masstr=Nsh.' * Nsh;
Massco=int(int(Masstr,x,-1,1),y,-1,1)
%vector scalar equal consistent mass matrix of the sum of each row
vm=sym(zeros(1,length(Massco)));
for i=1:length(Massco), vm(i) = sum(Massco(i,); end
vm;
sum(vm) .\ vm
%%%%%%%%%%%%%%%%%%%
vs=sym(zeros(1,length(Massco)));
for i=1:length(Massco), vs(i) = Massco(i,i); end
sum(vs) .\ vs
这是九节点的程序，和上面没什么两样啊。这个文章中的质量矩阵我也能由这个方法得到，当然三角形的程序会麻烦一点，但思路也是这样，如果按书上说的思路我觉得没什么问题，一开始我以为形函数给错了，仔细检查了也没发现问题，很疑惑

pasuka

是啊，这也是我纳闷的地方，如果方法不对的话，其它模型结果也不对才对啊，这是我9-noded element的结果[ 1/36, 1/36, 1/36, 1/36, 1/9, 1/9, 1/9, 1/9, 4/9]，和书上的一样啊
shipsj 发表于 2010-7-3 16:09

这个结果和高斯积分点有关
对于9节点矩形平面单元，采用2*2和3*3高斯积分，角点都是1/36，而边中点和中心点是4/36
对于8节点矩形平面单元，采用2*2和3*3高斯积分，角点分别是1/36、3/76，边中点是8/36、16/76
具体参看：库克，《有限元分析的概念和应用》

pasuka

所以lz说“中节点上的等效质量是角节点的8倍，而上面的却是16/3倍”正好是精确积分的结果，如果使用2*2点高斯积分就应该和王勖成书上的结果一致

shipsj

18# pasuka
[quote][/quote]这个结果和高斯积分点有关
我提的问题还没有用高斯积分的方法去实现，不过你的资料和数据好像和Zinkiewicz的书上不一样了，我去找下这边书