Mindlin板单元 单元刚度矩阵程序计算问题

zzjo

参考王勖成的《有限单元法》第10.4节的Mindlin板单元的推导，编程计算了板的弯曲问题，为什么计算出来的节点位移（挠度w和两个转角rotx，roty）与用ANSYS（采用shell143）计算出来的差很多呢？
我采用的是四节点四边形单元，形函数为：Ni=(1+ks(i)*ks)(1+it(i)*it)/4，Kb和Ks都采用2*2高斯积分。
算例是一个1m长的方板，厚度是h=0.05m，这算不算中厚板呢，是不是因为还是属于薄板，导致所说的“剪切锁死”。我算出的单元刚度矩阵里，Ks的值比Kb大，所以它占得权重很大。
希望能得到指教。

tonnyw

1# zzjo

Is it possible that in Ansys you are using reduced integration scheme?

You might just take one element and make a comparison of element stiffness matrix between yours and Ansys.

pasuka

参考王勖成的《有限单元法》第10.4节的Mindlin板单元的推导，编程计算了板的弯曲问题，为什么计算出来的节点位移（挠度w和两个转角rotx，roty）与用ANSYS（采用shell143）计算出来的差很多呢？
我采用的是四节点四边 ...
zzjo 发表于 2010-12-13 09:28

去找找A.J.M. Ferreira的《MATLAB Codes for Finite Element Analysis》的电子版，里面附带四节点中厚板单元的Matlab源代码，再做一下比较
ANSYS的壳单元，每节点6个自由度，板元部分也许是MITC4或RDKT

zzjo

1# zzjo

Is it possible that in Ansys you are using reduced integration scheme?

You might just take one element and make a comparison of element stiffness matrix between yours and Ansys.
tonnyw 发表于 2010-12-15 13:45

I don't use educed integration scheme,but I don't know whether it is used in Ansys.When I  use Ansys,all the selection is default.
I have made a comparison of element stiffness matrix, a part of them is the same,such as  row w<-to-> line w, the other is different,such as row Rotx<-to-> line Rotx.
I have a question, whether the element stiffness matrix ( mine and Ansys)must be all the same.

zzjo

3# pasuka
这本书不错，可是我不懂MatlAB，但我看了他的理论推导，最后那个B矩阵，我有跟他差个负号。但是我改用他的B阵后，算出来的结果还是不对。我就不明白了，单元刚度矩阵只跟这个B阵有关吧，为什么就是算的不对呢？
单刚跟ANSYS输出的单刚比较，只有一部分是不相同的，就是在整个单元刚度矩阵（12*12）中，每个节点的两个转角（sitax，sitay）所在行与列对应的那个值。其他的都是一样的。

zzjo

这个问题还没解决，实在不知道问题出在哪里，现在我把我的推导过程添加上来，希望得到大家解答，谢谢

pasuka

3# pasuka
这本书不错，可是我不懂MatlAB，但我看了他的理论推导，最后那个B矩阵，我有跟他差个负号。但是我改用他的B阵后，算出来的结果还是不对。我就不明白了，单元刚度矩阵只跟这个B阵有关吧，为什么就是算的 ...
zzjo 发表于 2010-12-16 18:50

二个转角对应的刚度值互换就正确了？
Mindlin板转角对应的坐标轴不要搞混了
绕X轴的转角应该对应y轴
绕y轴的转角应该对应x轴
缩减积分和完全积分的问题，ANSYS里面也许是非协调板单元

tonnyw

6# zzjo

As you can see from the attached figure, the formulation in Ansys is based on assumed strain method. I think this type of plate element is called MITC4 plate element which is locking free. I don't think you need to use reduced integration scheme.

In forming the element stiffness matrix, the strain components should be the ones as shown in the attached figure which is the linear interpolation of the strains directly calculated at the tying points A, B, C, and D.

zzjo

谢谢tonnyw和pasuka的解答，现在刚度阵的问题已经解决了。但是又遇到一个奇怪的现象。我计算了《MATLAB Codes for Finite Element Analysis》中给出的算例，，结果跟它和ANSYS比较都很好，然后我把算例改了一下，施加不同的荷载，考虑不同的约束，得到的结果跟ANSYS比较，有好有坏，不知道是什么原因。我把计算结果拿出来，希望大家能帮忙分析分析，特别是做板弯曲问题的大侠，到底是什么原因呢？

zzjo

最奇怪的是我计算的结点转角位移跟ANSYS得出的吻合的非常好，但是结点的挠度却有一半跟ANSYS的吻合很好，而另一半却不好，而且是交替出现的。不知道是不是求解方程的问题，我用了两种解方程的方法，一种是全选主元高斯消去，一种是LDLT分解解对称正定方程，两种解法的结果只是稍有不同而已，但这种奇怪现象没有消除。

pasuka

谢谢tonnyw和pasuka的解答，现在刚度阵的问题已经解决了。但是又遇到一个奇怪的现象。我计算了《MATLAB Codes for Finite Element Analysis》中给出的算例，，结果跟它和ANSYS比较都很好，然后我把算例改了一下，施加 ...
zzjo 发表于 2010-12-29 19:25

那本书里面的算例的结果lz和自己的程序比较过吗？那本书上的板单元是最基本的教学用代码，改边界条件和载荷条件，算出来的结果不理想当然是可能的，收敛速度慢，加密网格之后的结果呢？
四边简支、固支方板受集中力或均布载荷都是有解析解和级数解的，和自己的程序对比过了吗？
Matlab是非常简单的，会Fortran 的话，1个礼拜肯定可以上手，输出单刚、总刚和右端向量比较一下，又不是很困难的事情。
btw，请多动动脑子，多上网搜索资料，特别是英文资料，类似于Mindlin板的有限元matlab或Fortran代码网上肯定有的
又，老板的课题组就没有指定小老板或者博士师兄什么的带你吗？就只有lz一个人瞎捣鼓吗？

tonnyw

谢谢tonnyw和pasuka的解答，现在刚度阵的问题已经解决了。但是又遇到一个奇怪的现象。我计算了《MATLAB Codes for Finite Element Analysis》中给出的算例，，结果跟它和ANSYS比较都很好，然后我把算例改了一下，施加 ...
zzjo 发表于 2010-12-29 19:25

What do you mean you already solved the problem of element stiffness matrix? Did you get the same stiffness matrix as the one from Ansys?

zzjo

恩，我自己当然思考过，现在也一直在试验，只是结果确实比较怪。我只是想把原因找到，为什么这样。
现在的情况是，我的程序没有问题，单刚、总刚和右端向量比较过，用自己的程序算的书里给出的算例，结果和他给出的是一样的，而且跟ＳＮＳＹＳ算的也一样。只是改变约束条件后，我算的结果就不好了。
在王勖成那本书里有个算例，角点用柱子支撑的方形薄板承受均布荷载ｑ。他采用的是矩形板元，我用自己的程序（Ｍｉｎｄｌｉｎ板元＋线性等参元）计算，结点的挠度结果是不正确的。
他里面有句话：“上述矩形板元不能推广到一般的四边形板元，因为经过坐标变换的一般四边形单元不能满足常应变准则，收敛性很差不能用于实际计算”。不是这个原因吧

zzjo

12# tonnyw
单元刚度阵跟ＡＮＳＹＳ的不一样，但是跟《MATLAB Codes for Finite Element Analysis》中算的刚度阵是一样的。正如您们说是，ＡＮＳＹＳ里采用的是协调元，我计算的刚度阵当然不跟他一样了。
我现在处理的方法跟《MATLAB Codes for Finite Element Analysis》是一样的。
The sti&#64256;ness integral is solved by considering for the Q4 element, 2×2　Gauss points for the bending contribution and 1 point for the shear contribution.
This selective integration proved to be one of the simplest remedies for avoidingshear locking [4, 12].

tonnyw

12# tonnyw
单元刚度阵跟ＡＮＳＹＳ的不一样，但是跟《MATLAB Codes for Finite Element Analysis》中算的刚度阵是一样的。正如您们说是，ＡＮＳＹＳ里采用的是协调元，我计算的刚度阵当然不跟他一样了。
我现在 ...
zzjo 发表于 2010-12-30 12:47

According to what you mentioned before, the element in Ansys is based on assumed strain method which is not compatible. The element formulation in Matlab codes for finite element analysis is compatible.

tonnyw

谢谢tonnyw和pasuka的解答，现在刚度阵的问题已经解决了。但是又遇到一个奇怪的现象。我计算了《MATLAB Codes for Finite Element Analysis》中给出的算例，，结果跟它和ANSYS比较都很好，然后我把算例改了一下，施加 ...
zzjo 发表于 2010-12-29 19:25

I always like playing with one element first. How about choosing one element and fixing two corner nodes to see how the results look like?

By the way, how did you enforce the displacement boundary condition?

zzjo

16# tonnyw
我这样做过，取一个单元比较单刚，取四个单元比较单刚和合成的总刚。现在不是刚度阵地问题了，是我采用的Mindlin板元适合什么样的求解问题。厚板（h/L>1/5）的情况可以计算，但当板的厚度变薄，我通过“The sti&#64256;ness integral is solved by considering for the Q4 element, 2×2　Gauss points for the bending contribution and 1 point for the shear contribution.”这样处理，也可以解决薄板四边简支或固支的问题，但是却不能解决只在薄板边界上约束几个点的问题（如固支四个角点）。这是困惑的地方。也许是这种“This selective integration proved to be one of the simplest remedies for avoidingshear locking ”简单的处理方法不能正确处理在薄板边界上约束几个点的问题。如图示的情况（或者受均布荷载），约束1、2、52三点。

tonnyw

16# tonnyw
我这样做过，取一个单元比较单刚，取四个单元比较单刚和合成的总刚。现在不是刚度阵地问题了，是我采用的Mindlin板元适合什么样的求解问题。厚板（h/L>1/5）的情况可以计算，但当板的厚度变薄，我通过 ...
zzjo 发表于 2010-12-31 11:00

What  I am saying is that you take one element and fix two corner nodes. I didn't ask you to compare the stiffness matrix.

zzjo

Ok,我照你说的是了下，如果固定同一个边上的两点，可以计算出另两个结点的位移，但是固定对角线的两个点，就不能求解了，解方程错误。分别算了下两种情况刚度阵的秩，前一种情况是满秩的12，后一种情况的秩是11，行列式为零，固定一个角点也这种情况。为什么会这样呢？

tonnyw

Ok,我照你说的是了下，如果固定同一个边上的两点，可以计算出另两个结点的位移，但是固定对角线的两个点，就不能求解了，解方程错误。分别算了下两种情况刚度阵的秩，前一种情况是满秩的12，后一种情况的秩是11，行列 ...
zzjo 发表于 2010-12-31 15:05

Each node has three degrees of freedom and the total dofs are 12. The stiffness matrix is always singular before applying boundary conditions. So tell us how you enforce the boundary conditions.