zbh 发表于 2017-2-23 23:38:48

abaqus中形函数问题

本帖最后由 zbh 于 2017-2-24 08:55 编辑

求助各位大神:

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


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

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

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



zbh 发表于 2017-2-24 07:22:08

本帖最后由 zbh 于 2017-2-24 07:26 编辑

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([ ,\
                         ,\
                         [ -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
      eta = ksieta
# 二维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行正好相反
# [
# [   767.5601073219139.05608368   4852.9194919 ]
# [-414.5649770218841.90808726   4441.09011661]
# [   -39.4601000820010.27764111   4845.41393839]]

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

cartographer 发表于 2017-2-24 13:11:25

可能是应为这个吧

sunx 发表于 2019-4-14 16:01:42

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

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

cartographer 发表于 2019-4-14 22:31:48

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

帮助文档啊。
页: [1]
查看完整版本: abaqus中形函数问题