平面应力问题的Matlab有限元编程问题

shenxiong

已知条件如图所示，上下边界受均不拉应力载荷1Mpa；
有限元编程：以单元1为例，节点号1 2 4 3，
我们知道，弹性力学里面有个应力边界条件问题，现在底部边界法线的方向余弦，l=0(与x轴正半轴夹角的余弦值)，m=-1（与y轴正半轴夹角的余弦值）
表达式：l*sigmax+m*taoyx=0,m*sigmay+l*taoxy=-1;解得下边界节点的应力（理论值）：水平方向切应力taoyx=0，竖直方向应力sigmay=1Mpa,

问题1,有限元matlab编程计算结果:下面边界的节点竖直方向应力结果是不是接近1Mpa，水平方向切应力接近零才近乎正确吧。。。
问题2，我们知道，有限元的刚度矩阵是奇异矩阵，也即是说刚度矩阵的行列式是接近或为零的，
           一般非齐次线性方程组K.u=f,施加位移边界条件的话，对刚度矩阵进行了相关变化，使得奇异变成了非奇异矩阵，可以求解位移u中未知的值，依据是：线性代数的克拉默法则，系数矩阵的行列式值不为0时，方程组有唯一解，位移u有唯一解；（逆否命题：解不唯一的话，行列式值为0；）
         现在这个平面应力问题，上下边界受拉，刚度矩阵奇异，解方程K.u=f,方程组的解还是唯一的吗？不唯一的话，但是我的Matlab求解的位移u值却给出了一组，数量级可接受的。

tonnyw

问题1,有限元matlab编程计算结果:下面边界的节点竖直方向应力结果是不是接近1Mpa，水平方向切应力接近零才近乎正确吧。。。
>>I think you are right on this.
问题2，我们知道，有限元的刚度矩阵是奇异矩阵，也即是说刚度矩阵的行列式是接近或为零的，
           一般非齐次线性方程组K.u=f,施加位移边界条件的话，对刚度矩阵进行了相关变化，使得奇异变成了非奇异矩阵，可以求解位移u中未知的值，依据是：线性代数的克拉默法则，系数矩阵的行列式值不为0时，方程组有唯一解，位移u有唯一解；（逆否命题：解不唯一的话，行列式值为0；）
         现在这个平面应力问题，上下边界受拉，刚度矩阵奇异，解方程K.u=f,方程组的解还是唯一的吗？不唯一的话，但是我的Matlab求解的位移u值却给出了一组，数量级可接受的。
>>This type of problem is called pure Neumann problem which doesn't have displacement boundary condition. The solution is not unique, which means adding an arbitrary constant c to any solution u will produce another solution as well. The reason that Matlab doesn't give your error message even though the stiffness matrix is singular is because of round-off error. In this case if you use Matlab to calculate the determinant, you still get a very small number instead of 0. So Matlab still considers the stiffness matrix is not singular. You can also check the singular behavior of the stiffness matrix by checking its condition number. Usually in handling the pure Neumann problem, you can add small perturbation on the diagonal terms, which means you attach soft springs to each nodes.

shenxiong

tonnyw 发表于 2013-4-15 21:27
问题1,有限元matlab编程计算结果:下面边界的节点竖直方向应力结果是不是接近1Mpa，水平方向切应力接近零才 ...

这么牛x，谢谢你的帮助，我回去试一试

pasuka

shenxiong 发表于 2013-4-16 09:02
这么牛x，都是英文有点不是很懂

斑竹在美国，办公室电脑没有中文输入法，只能打英文
又，如果是全日制非在职的研究生，这点简单的英语也看不懂的话，真有点让人失望

shenxiong

pasuka 发表于 2013-4-16 14:16
斑竹在美国，办公室电脑没有中文输入法，只能打英文
又，如果是全日制非在职的研究生，这点简单的英语也 ...

我只是感觉很奇怪，为啥都用英文回复

shenxiong

pasuka 发表于 2013-4-16 14:16
斑竹在美国，办公室电脑没有中文输入法，只能打英文
又，如果是全日制非在职的研究生，这点简单的英语也 ...

我只是感觉很奇怪，为啥都用英文回复

shenxiong

pasuka 发表于 2013-4-16 14:16
斑竹在美国，办公室电脑没有中文输入法，只能打英文
又，如果是全日制非在职的研究生，这点简单的英语也 ...

我只是感觉很奇怪，为啥都用英文回复

shenxiong

pasuka 发表于 2013-4-16 14:16
斑竹在美国，办公室电脑没有中文输入法，只能打英文
又，如果是全日制非在职的研究生，这点简单的英语也 ...

我只是感觉很奇怪，为啥都用英文回复

shenxiong

pasuka 发表于 2013-4-16 14:16
斑竹在美国，办公室电脑没有中文输入法，只能打英文
又，如果是全日制非在职的研究生，这点简单的英语也 ...

我只是感觉很奇怪，为啥都用英文回复

snowwave02

tonnyw 发表于 2013-4-15 21:27
问题1,有限元matlab编程计算结果:下面边界的节点竖直方向应力结果是不是接近1Mpa，水平方向切应力接近零才 ...

请问版主，这个微小的扰动刚度应该加多少？对平面应力单元有没有一个确定的公式？最近也遇到这个问题了，在abaqus的减缩积分中，abaqus加入了沙漏修正，我的程序没加，所以导致了刚度矩阵奇异。

xiaofei170

可以看看，了解下。