abaqus中形函数问题

zbh

求助各位大神：

我计算了一个简单的问题，就是一个1mm*1mm的二维平板，下边固定，右上角有一个45角度方向的0.1mm的位移。计算模型见下图


abaqus的计算输入文件job-1.inp和计算结果输出文件job-1.dat见附件。

我做这么简单计算的目的就是要验证abaqus形函数是否正确，如果正确，用形函数和单元4个结点的应力值（未应力磨平处理的）做内插，应该就得到应用内高斯积分点处的应力值。

但是通过编制python脚本，由单元4个结点处应力值通过形函数矩阵N内插得到的单元内4个高斯积分点处的应力值，与abaqus计算得到的高斯积分点略有不同，即第三个和第四个积分点处的应力值对调了，我查了很长时间，应该是我的形函数有问题，请版主给看看，问题出在哪里，谢谢您的关注。

zbh

test.py脚本内容如下：

# 计算过程
# 第一步： abaqus job=job-1
# 第二步： 打开job-1.dat，取单元6的4个未磨平的结点应力stress_node和单元内的gauss积分点应力stress_gauss
# 第三步： 按照二维等参元的积分点位置排序1./sqrt(3.)*[(-1,-1),(1,-1),(1,1),(-1,1)],计算等参元的形函数
# 第四步： 由单元4个结点处应力通过形函数内插得到单元内的积分点处应力值

import numpy as np

# 第二步
# -----------------------------------------------------------------------------
# 由abaqus软件计算得到的stress_node应力值，由job-1.dat文件拷贝得到，所以数据是准确的
# 以单元6为例，单元4个结点号分别为7，8，12，11，未经磨平处理的应力如下(平面应力)
stress_node=np.matrix([ [   1705.,      2.0837E+04,   5551.], \
                        [   1065.,      1.8822E+04,   4857.], \
                        [  -985.1,      1.8305E+04,   4142.], \
                        [  -332.8,      2.0331E+04,   4844.]])


# 由abaqus软件计算得到的stress_gauss应力值，由job-1.dat文件拷贝得到，所以数据是准确的
# gauss积分点处应力值(平面应力)
stress_gauss=np.matrix([ [  1138.,      2.0303E+04,   5255. ],\
                         [  767.5,      1.9139E+04,   4853. ],\
                         [ -39.49,      2.0010E+04,   4845. ],\
                         [ -414.6,      1.8842E+04,   4441. ]])

# 第三步：
# -----------------------------------------------------------------------------
g3 = 1.0/np.sqrt(3.0)
ksieta=np.matrix([[-g3,g3,g3,-g3],[-g3,-g3,g3,g3]]) # integral point's position

for ip in range(0,4): # loop over 4 integral point
        ksi = ksieta[0,ip]
        eta = ksieta[1,ip]
# 二维4结点等参元的形函数
        N = 0.25*np.matrix([(1-ksi)*(1-eta),(1+ksi)*(1-eta),(1+ksi)*(1+eta),(1-ksi)*(1+eta)])
        print N


# 形函数矩阵

N=np.matrix([[ 0.62200847,  0.16666667,  0.0446582,   0.16666667],\
             [ 0.16666667,  0.62200847,  0.16666667,  0.0446582 ],\
             [ 0.0446582,   0.16666667,  0.62200847,  0.16666667],\
             [ 0.16666667,  0.0446582,   0.16666667,  0.62200847]])

# 第四步
# -----------------------------------------------------------------------------
# 为了上面形函数矩阵的正确性，将单元4个结点处未经磨平的应力值，通过形函数矩阵内插，应该得到单元内gauss积分点处的应力值
# 由于 形函数矩阵N * 单元4个结点处应力值 => 单元4个Gauss积分点处应力值
stress_gauss=np.dot(N,stress_node)

stress_gauss = N * stress_node

print stress_gauss

# 根据上面的stress_gauss，发现与下面的第3行和第4行正好相反
# [[  1138.5649843   20303.7589709    5254.57664704]
# [   767.56010732  19139.05608368   4852.9194919 ]
# [  -414.56497702  18841.90808726   4441.09011661]
# [   -39.46010008  20010.27764111   4845.41393839]]

# 由此推断上面的形函数矩阵N有问题，请版主给看看，问题出在哪里，谢谢您的关注。

cartographer

可能是应为这个吧

sunx

cartographer 发表于 2017-2-24 13:11
可能是应为这个吧

你好，请问你这个是什么资料

cartographer

sunx 发表于 2019-4-14 16:01
你好，请问你这个是什么资料

帮助文档啊。