第十篇使用8、14、20节点6面体的立方体弹性固体的三维分析

引言

本篇程序的整体路线基本与第七篇是类似的，所以大家可以以第七篇为参考。程序由8-14-20节点的六面体网格划分而成，其中的内置几何子程序hexahedron_xz产生以yz平面开始向y方向的节点和单元标号，具体可见下图。

具体的局部节点单元标号参照下图，以xz平面为顺时针编号，之后向y方向扩展

一个20节点六面体单元需要的积分点为27个，为了方便考虑”折减积分“，将积分点减少为8个

计算实例

节点值为20；
网格划分x方向为1，y方向为3，z方向为2，积分点为8，单元性质类型由两类；
单元性质分别有e=100.0、v=0.3(序号1)和e=50.0和v=0.3(序号2)两种；
单元对应的性质类型序号为1、2、1、2、1、2；
三个方向坐标点值为 x：0.0，0.5；
y：0.0，1.0，2.0，3.0；
z：0.0，-1.0，-2.0；
约束结点有46个（1代表自由；0代表约束）例如节点1 x，y方向约束，z方向自由；
荷载节点为8个，1、2、3、14、15、20、21、22的z方向荷载分别为(0.0417、-0.1667、0.0417、-0.1667、-0.1667、0.0417、-0.1667、0.0417)；
没有固定位移和转角值。

代码

import numpy as np
import A
type_2d='plane'
element='hexahedron'
dir='y'
nxe=1
nye=3
nze=2
nst=6
np_types=2
loaded_nodes=8
fixed_freedoms=0
nod=20
nprops=2
nr=46
nodof=3
ndim=3
nip=8
nels=np.array([0])
nn=np.array([0])
A.mesh_size(element,nod,nels,nn,nxe,nye,nze)
ndof=nod*nodof
if type_2d=='axisymmetric':nst=4
nels=int(nels)
nn=int(nn)
#初始化定义数组
nf=np.ones((nodof,nn),dtype=np.int64)
points=np.ones((nip,ndim))
g=np.ones((ndof,1),dtype=np.int64)
g_coord=np.ones((ndim,nn))
fun=np.ones((nod,1))
coord=np.ones((nod,ndim))
jac=np.ones((ndim,ndim))
g_num=np.ones((nod,nels))
der=np.ones((ndim,nod))
deriv=np.ones((ndim,nod))
bee=np.zeros((nst,ndof))
km=np.ones((ndof,ndof))
eld=np.ones((ndof,1))
weights=np.ones((nip,1))
g_g=np.ones((ndof,nels))
prop=np.ones((nprops,np_types))
num=np.ones((nod,1),dtype=np.int64)
#x_coords=np.ones((nxe+1,1))
#y_coords=np.ones((nxe+1,1))
etype=np.ones((nels,1),dtype=np.int64)
gc=np.ones((ndim,1))
dee=np.zeros((nst,nst))
sigma=np.ones((nst,1))
if np_types==1:etype[:,0]=1
else:etype[:,0]=(1,2,1,2,1,2)
if np_types==1:prop[:nprops,np_types-1]=(1.0e6,0.3)
else:prop[0,:]=(100,50)prop[1,:]=(0.3,0.3)
x_coords=np.array([0.0,0.5])
y_coords=np.array([0.0,1.0,2.0,3.0])
z_coords=np.array([0.0,-1.0,-2.0])
#读取nr
if nr!=0:Dim_1=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,16,18,19,20,23,25,28,30,31,32,33,35,37,38,39,42,44,47,49,50,51,52,54,56,57,58,61,63,66,68,69,70]nf_value=np.array([[0,1,1,0,1,0,1,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],\[0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,1,1,0,0,0,1,1,0,0,1,1,1,1,0,0,0,1,1,0,0,1,1,1,1,0,0,0],\[1,1,1,1,1,1,1,1,1,1,0,0,0,1,1,0,0,1,1,1,1,0,0,0,1,1,0,0,1,1,1,1,0,0,0,1,1,0,0,1,1,1,1,0,0,0]])m=0for i in Dim_1:for j in range(1,nodof+1):nf[j-1,i-1]=nf_value[j-1,m]m=m+1
#form nf
A.formnf(nf)
neq=int(max(nf.reshape(nf.shape[0]*nf.shape[1],1)))
kdiag=np.zeros((neq,1),dtype=np.int64)
loads=np.zeros((neq+1,1))
##累加单元去发现全局尺寸
for iel in range(1,nels+1):A.hexahedron_xz(iel,x_coords,y_coords,z_coords,coord,num)A.num_to_g(num,nf,g)g_num[:,iel-1]=num[:,0]g_coord[:,num[:,0]-1]=np.transpose(coord[:,:])g_g[:,iel-1]=g[:,0]A.fkdiag(kdiag,g)
for i in range(1,neq):kdiag[i]=kdiag[i]+kdiag[i-1]
kv=np.zeros((kdiag[neq-1,0],1),dtype=float)
print('一共有'+str(neq)+'等式和天际线存储个数为'+str(kdiag[neq-1]))
#单元刚度积分和安装
A.sample(element,points,weights)
for iel in range(1,nels+1):A.deemat(dee,prop[0,etype[iel-1]-1],prop[1,etype[iel-1]-1])num[:,0]=g_num[:,iel-1]coord[:,:]=np.transpose(g_coord[:,num[:,0]-1])g[:,0]=g_g[:,iel-1]km[:]=0for i in range(1,nip+1):#A.shape_fun(fun,points,i)A.shape_der(der,points,i)jac=np.dot(der,coord)det=np.linalg.det(jac)A.invert(jac)deriv=np.dot(jac,der)A.beemat(bee,deriv)km[:]=km[:]+np.dot(np.dot(np.transpose(bee),dee),bee)*det*weights[i-1]
#call fsparvA.fsparv(kv,km,g,kdiag)
###读荷载和位移
if loaded_nodes!=0:Dim_2 = [1,2,3,14,15,20,21,22]load_value=np.array([[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0.0417,-0.1667,0.0417,-0.1667,-0.1667,0.0417,-0.1667,0.0417]])m=0for i in Dim_2:for j in range(1,nodof+1):loads[nf[j-1,i-1]-1]=load_value[j-1,m]m=m+1
#####方程求解
if fixed_freedoms!=0:node=np.ones((fixed_freedoms,1),dtype=np.int64)sense=np.ones((fixed_freedoms,1),dtype=np.int64)value=np.ones((fixed_freedoms,1))no=np.ones((fixed_freedoms,1),dtype=np.int64)node[:,0]=(1,3)sense[:,0]=(2,1)value[:,0]=(-0.001,-0.005)
####位移赋值for i in range(1,fixed_freedoms+1):no[i-1]=nf[sense[i-1]-1,node[i-1]-1]kv[kdiag[no-1]-1]=kv[kdiag[no-1]-1]+1e20loads[no-1,0]=kv[kdiag[no-1,0]-1,0]*value
##荷载增量累加
A.sparin(kv,kdiag)
A.spabac(kv,loads,kdiag)
loads[neq]=0
print('节点 x-位移 y-位移 z-位移')
for k in range(1,nn+1):print(k,end=' ')for m in range(1,nodof+1):print('{:9.4e}'.format(loads[nf[m-1,k-1]-1,0]),end=' ')print()
#取得单元中心点的力矩
nip=1
points=np.ones((nip,ndim))
weights=np.ones((nip,1))
A.sample(element,points,weights)
print('积分点nip='+str(nip)+'应力为')
print('单元 x-坐标     y-坐标        z-坐标      sig_x       sig_y         sig_z       tau_xy\tau_yz      tau_zx')
for iel in range(1,nels+1):A.deemat(dee,prop[0,etype[iel-1]-1],prop[1,etype[iel-1]-1])num[:,0]=g_num[:,iel-1]coord[:,:]=np.transpose(g_coord[:,num[:,0]-1])g[:,0]=g_g[:,iel-1]eld=loads[g-1,0]for i in range(1,nip+1):A.shape_fun(fun,points,i)A.shape_der(der,points,i)for i in range(1,ndim+1):gc[i-1]=np.dot(np.transpose(fun[:,0]),coord[:,i-1])jac=np.dot(der,coord)A.invert(jac)deriv=np.dot(jac,der)A.beemat(bee,deriv)sigma[:]=np.dot(dee[:],np.dot(bee[:],eld[:]))print(iel,end='  ')for i in range(1,ndim+1):print('{:9.4e}'.format(gc[i-1,0]),end='  ')for i in range(1,nst+1):print('{:9.4e}'.format(sigma[i-1,0]),end='  ')print( )

终端输出结果

使用8、14、20节点6面体的立方体弹性固体的三维分析(python,有限元)相关推荐

8节点矩形弹性轴对称弹性固体的非轴对称分析(python,有限元)
第九篇 8节点矩形弹性轴对称弹性固体的非轴对称分析介绍本篇讨论一个轴对称或者反轴对称的固体遭受一个非对称的简谐波荷载的情况将引入三个新的参数,lth表示表示承受的对应简谐波荷载的振幅值,ifla ...
Delphi中预想不到的代码楼主zswang（伴水清清)(专家门诊清洁工）2002-05-16 14:20:38 在 Delphi / VCL组件开发及应用提问
Delphi中预想不到的代码楼主zswang(伴水清清)(专家门诊清洁工)2002-05-16 14:20:38 在 Delphi / VCL组件开发及应用提问 No.1 Delphi中的In ...
时间格式的转换例如：(2021-05-10 14:20:43) 转为( 2021年5月10日 14时20分43秒)
console.log(name('2021-02-10 14:20:43'));function name(date) {const arr = date.split(/[-: ]/)return ...
基于改进的 IEEE24 节点电力系统和比利时 20 节点天然气系统通过电转气和燃气轮机耦合
1 概述做综合能源经常需要用到算例,本文展示一个很棒的算例,还要数据. 2 算例如图所示的电-气互联综合能源系统中,IEEE24节点系统有10台发电机组;节点2.18和22为燃气轮机,分别与天然气 ...
Gerrit version 2.14.20 is now available
Gerrit version 2.14.20 is now available 今天发现2.14版本的又有个个更新了.现在2.14更新到2.14.17版本了. gerrit.war 历史版本下载各个 ...
【MATLAB深度学习工具箱】学习笔记--体脂估计算例再分析：拟合神经网络fitnet里面的函数】
介绍上一篇 [MATLAB深度学习工具箱]学习笔记--体脂估计算例再分析:拟合神经网络fitnet里面的数据结构]_bear_miao的博客-CSDN博客原文链接如下[MATLAB深度学习工具箱]学 ...
python有限元传热求解_Python进行有限元编程-平面应力问题（三节点三角形单元）...
参考书籍是:<有限元方法基础教程>(国际单位制)(第五版) 章节为:第6章建立平面应力和平面应变刚度方程重要的事情继续强调(有限元的基本计算流程): Step 1: 选择单元类型. S ...
从零开始刷Leetcode——字符串（13.14.20.28）
文章目录 13. 罗马数字转整数 14. 最长公共前缀 20. 有效的括号 28. 实现 strStr() 13. 罗马数字转整数给定一个罗马数字,将其转换成整数.输入确保在 1 到 3999 的范 ...
7-4 宿舍谁最高？（20 分) map+结构体的应用
7-4 宿舍谁最高? (20 分) 学校选拔篮球队员,每间宿舍最多有4个人.现给出宿舍列表,请找出每个宿舍最高的同学.定义一个学生类Student,有身高height,体重weight等. 输入格式: ...
谷哥学术2022年2月份资源分享下载列表14/20
资源名称下载地址关键词简约图文风家具装修PPT模板_PPT 16_9_2022-02-09.pptx https://download.csdn.net/download/tysonchiu/7 ...

使用8、14、20节点6面体的立方体弹性固体的三维分析(python,有限元)

第十篇使用8、14、20节点6面体的立方体弹性固体的三维分析

引言

计算实例

代码

终端输出结果

使用8、14、20节点6面体的立方体弹性固体的三维分析(python,有限元)相关推荐

最新文章

热门文章

使用8、14、20节点6面体的立方体弹性固体的三维分析(python,有限元)

第十篇 使用8、14、20节点6面体的立方体弹性固体的三维分析

引言

计算实例

代码

终端输出结果

使用8、14、20节点6面体的立方体弹性固体的三维分析(python,有限元)相关推荐

最新文章

热门文章

第十篇使用8、14、20节点6面体的立方体弹性固体的三维分析