Note_Fem_C3D8单元公式
单元几何
单元使用八节点，每个节点有三个平移自由度; 局部坐标使用(r,s,t)坐标系，r,s,t取值范围-1到+1；使用等参数单元推导刚度矩阵
- Hexahedron (brick) element faces
Face 1 1 – 2 – 3 – 4 face Face 2 5 – 8 – 7 – 6 face Face 3 1 – 5 – 6 – 2 face Face 4 2 – 6 – 7 – 3 face Face 5 3 – 7 – 8 – 4 face Face 6 4 – 8 – 5 – 1 face
- 积分点布置（abaqus）
单元刚度矩阵Ke
节点形函数为，ri,si,ti是节点i的(r,s,t)坐标值：
几何插值：
位移场插值：
参考abaqus理论手册，一阶六面体也是用这个插值函数：
应变分量推导：
节点的应变微分算子为：
使用链式法则推掉形函数对局部坐标的微分：
J是雅可比矩阵，形式为：
J进一步分解：
形函数对r,s,t的偏导：
至此B矩阵推导完毕
采用各向同性材料，则弹性矩阵为：
单元刚度矩阵Ke的公式为：
完全积分方案：每个方向上使用2个积分点进行积分，满足精确积分，共有8个积分点进行数值积分。
单元质量矩阵
推导公式参见《计算固体力学》page289
- 集中质量矩阵
集中质量矩阵的思路是将单元质量平均分配到每个节点自由度。因此：
单元总质量为：
因此，单元质量矩阵为：
- 一致质量矩阵
等效单元节点力
类似C3D4单元，针对热载荷，体积力(N/m^3)，作用在单元面1的面力N/m*m的等效节点力公式子为：
静力分析问题公式
由C3D8组成的系统的静力分析问题，可以归结为求解下列公式，可以得到节点位移\(\{Q\}_{N\times1}\)：
动力分析公式，特征值问题见C3D4单元
B-Bar修正
参考:
https://www.jishulink.com/post/1948589 https://zhuanlan.zhihu.com/p/27784056916?utm_psn=1880581810448221613 https://zhuanlan.zhihu.com/p/27663938155 https://www.jishulink.com/post/357900 https://www.fangzhenxiu.com/post/10174006/ https://mp.weixin.qq.com/s/3Akfa6T3vvomC7PiUVeX5Q 首先我这里使用的B矩阵并不是文中的标准矩阵，要有个T矩阵进行行变换。
Bmu：标准应变矩阵 Bmy：我的应变矩阵
>> T*Bmy==Bmu
ans =
1 1 1
1 1 1
1 1 1
1 1 1
1 1 1
1 1 1
>> T
T =
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 0 1 0
0 0 0 0 0 1
0 0 0 1 0 0
>> Bmy
Bmy =
[ b1, 0, 0]
[ 0, b2, 0]
[ 0, 0, b3]
[ b2, b1, 0]
[ 0, b3, b2]
[ b3, 0, b1]
>> Bmu
Bmu =
[ b1, 0, 0]
[ 0, b2, 0]
[ 0, 0, b3]
[ 0, b3, b2]
[ b3, 0, b1]
[ b2, b1, 0]
B-bar修正：
代码实现
def B_modify(self,r,s,t):
"""使用b-bar技术生成B=B_dil+B_dev"""
# step1计算B矩阵
b_dev,b_dil=[],[]
T=np.array([[1,0,0,0,0,0],
[0,1,0,0,0,0],
[0,0,1,0,0,0],
[0,0,0,0,1,0],
[0,0,0,0,0,1],
[0,0,0,1,0,0]])
for i in range(8):
nr=self.Ni_r(i)
ns=self.Ni_s(i)
nt=self.Ni_t(i)
b=np.array([nr(r,s,t),ns(r,s,t),nt(r,s,t)])
nx,ny,nz=sc.linalg.solve(self.J(r,s,t),b)
bi=np.array([[nx,0,0],
[0,ny,0],
[0,0,nz],
[ny,nx,0],
[0,nz,ny],
[nz,0,nx]],dtype=float)
# step2转换到标准B矩阵
bi_standard=T.dot(bi)
# step3 B矩阵变换
bi_dil=(1/3.0)*np.array([[nx,ny,nz],
[nx,ny,nz],
[nx,ny,nz],
[0,0,0],
[0,0,0],
[0,0,0]],dtype=float)
bi_dev=bi_standard-bi_dil
b_dev.append(bi_dev)
b_dil.append(bi_dil)
b8_dev=np.concatenate(b_dev,axis=1)
b8_dil=np.concatenate(b_dil,axis=1)
return b8_dev,b8_dil
def Ke_bbar(self):
"""返回b-bar修正后的单元刚度矩阵"""
ke=np.zeros((24,24))
# 对B_dev进行全积分
for i in range(self.gauss_w.shape[0]):
# 积分点坐标,权重
r,s,t=self.gauss_p[i]
w=self.gauss_w[i]
detJ=self.detJ(r,s,t)
B_dev,_=self.B_modify(r,s,t)
B_devT=B_dev.T
ke+=w*detJ*np.dot(B_devT,np.dot(self.D,B_dev))
# B_dil进行缩减积分
w,r,s,t=8.0,0.0,0.0,0.0
_,B_dil=self.B_modify(r,s,t)
B_dilT=B_dil.T
detJ=self.detJ(r,s,t)
ke+=w*detJ*np.dot(B_dilT,np.dot(self.D,B_dil))
return ke
使用abaqus建模单个C3D8单元，提取刚度矩阵，使用b-bar技术修正自编C3D8单元，以abaqus单元刚度为基准，前后单元刚度矩阵差别如下：
- 使用b-bar修正后：
- 修正前：
使用单个单元进行验证：
自编C3D8的节点位移结果：
abaqus节点位移结果：
基本一致