一次C3D4单元的debbug过程
20230326使用python编程:
完成C3D4的类实现
完成了C3D4尽力分析的计算模板Script_for_C3D4_Static_Analysis.py
完成pyvista+tkinter的小型后处理
使用悬臂梁模型和abaqus对比结果,发现对不上,自己的后处理显示很小. abaqus节点力1000N
自编结果:
后来提取其中一个单元的刚度矩阵发现也对不上.
# -------------可视化结果-
# ic(element_list[267].label)
# ic(element_list[267].Ke)
# # 读取abaqus生成的刚度矩阵文件
from EF2D.tools import *
# e268_stiff=read_mtx("L:\\EF2D\\tests\\data\\c3d4_abaqus\\e268stiff_STIF2.mtx")
# ic(e268_stiff)
# diff=e268_stiff-element_list[267].Ke
# imShowMat(diff,"diff")使用RAO的书中的四面体单元例子:
然后在验证自己的C3D4单元:
import numpy as np
from EF2D.Elements import C3D4
import scipy.io as sio
import matplotlib.pyplot as plt
from typing import List,Dict
from icecream import ic
import sys
import os
from EF2D.tools import *
from EF2D import *
n1=Node(label=1,x=0,y=0,z=0)
n2=Node(label=2,x=0.10,y=0,z=0)
n3=Node(label=3,x=0,y=0.15,z=0)
n4=Node(label=4,x=0,y=0,z=0.2)
initNumpyPrintOptions()
e=C3D4(label=1,node1=n1,node2=n2,node3=n3,node4=n4,E=207e9,nu=0.3)
ic(e.B)
# ic| e.B: array([[-10. , 0. , 0. , 10. , 0. , 0. , -0. , 0. , 0. , -0. , 0. , 0. ],
# [ 0. , -6.6667, 0. , 0. , -0. , 0. , 0. , 6.6667, 0. , 0. , 0. , 0. ],
# [ 0. , 0. , -5. , 0. , 0. , 0. , 0. , 0. , -0. , 0. , 0. , 5. ],
# [ -6.6667, -10. , 0. , -0. , 10. , 0. , 6.6667, -0. , 0. , 0. , -0. , 0. ],
# [ 0. , -5. , -6.6667, 0. , 0. , -0. , 0. , -0. , 6.6667, 0. , 5. , 0. ],
# [ -5. , 0. , -10. , 0. , 0. , 10. , -0. , 0. , -0. , 5. , 0. , -0. ]])
ic(e.Ve)
# ic| e.Ve: 0.0005
ic(e.D)
# 1e11*array([[2.7865, 1.1942, 1.1942, 0. , 0. , 0. ],
# [1.1942, 0. , 1.1942, 0. , 0. , 0. ],
# [1.1942, 1.1942, 2.7865, 0. , 0. , 0. ],
# [0. , 0. , 0. , 0.7962, 0. , 0. ],
# [0. , 0. , 0. , 0. , 0.7962, 0. ],
# [0. , 0. , 0. , 0. , 0. , 0.7962]])
ic(e.Ke*1e-10)
ic(e.Ke[1,7])
# ic| e.Ke*1e-10: array([[ 1.6697, 0.6635, 0.4976, -1.3933, -0.2654, -0.199 , -0.1769, -0.3981, 0. , -0.0995, 0. , -0.2986],
# [ 0.6635, 0.4976, 0.3317, -0.3981, -0.3981, 0. , -0.2654, 0. , -0.1327, 0. , -0.0995, -0.199 ],
# [ 0.4976, 0.3317, 0.9233, -0.2986, 0. , -0.3981, 0. , -0.199 , -0.1769, -0.199 , -0.1327, -0.3483],
# [-1.3933, -0.3981, -0.2986, 1.3933, 0. , 0. , 0. , 0.3981, 0. , 0. , 0. , 0.2986],
# [-0.2654, -0.3981, 0. , 0. , 0.3981, 0. , 0.2654, 0. , 0. , 0. , 0. , 0. ],
# [-0.199 , 0. , -0.3981, 0. , 0. , 0.3981, 0. , 0. , 0. , 0.199 , 0. , 0. ],
# [-0.1769, -0.2654, 0. , 0. , 0.2654, 0. , 0.1769, 0. , 0. , 0. , 0. , 0. ],
# [-0.3981, 0. , -0.199 , 0.3981, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.199 ],
# [ 0. , -0.1327, -0.1769, 0. , 0. , 0. , 0. , 0. , 0.1769, 0. , 0.1327, 0. ],
# [-0.0995, 0. , -0.199 , 0. , 0. , 0.199 , 0. , 0. , 0. , 0.0995, 0. , 0. ],
# [ 0. , -0.0995, -0.1327, 0. , 0. , 0. , 0. , 0. , 0.1327, 0. , 0.0995, 0. ],
# [-0.2986, -0.199 , -0.3483, 0.2986, 0. , 0. , 0. , 0.199 , 0. , 0. , 0. , 0.3483]])
# ic| e.Ke[1,7]: 0.0
# 20250328,悬臂梁结构个abaqus验证,结果对比不上
## array([[2.7865, 1.1942, 1.1942, 0. , 0. , 0. ],
# [1.1942, 0. , 1.1942, 0. , 0. , 0. ],
# [1.1942, 1.1942, 2.7865, 0. , 0. , 0. ],
# [0. , 0. , 0. , 0.7962, 0. , 0. ],
# [0. , 0. , 0. , 0. , 0.7962, 0. ],
# [0. , 0. , 0. , 0. , 0. , 0.7962]])最后破案了,弹性矩阵(2,2)忘记赋值了!!!!!!!!!!!!
后续:修改了之后,还是错误,最后通过提取abaqus刚度矩阵进行对比,发现是刚度矩阵不对.
解决方案:重新推导单元刚度矩阵,发现之前的推导有误,现在已经修正了,结果和abaqus对上了.