MATLAB和Python网格桁架框架构件刚度载荷位移和受力微分方程

53 篇文章 0 订阅
39 篇文章 0 订阅

🎯要点

  1. 数学​方法​:🎯一维线性边界值问题:🖊高斯求积法则 | 🖊洛巴托求积法则 | 🖊矩阵插值和微分计算 | 🖊在细化网格上生成值。🎯二维边界值问题:构建二维网格:🖊生成几何边界 | 🖊生成扩展转置连接矩阵 | 🖊生成互补网格 | 🖊绘制标记网格图和互补标记网格 | 🖊从定义几何尺寸生成网格求解:使用逐元素策略、部分(或全部)矢量化方法或稀疏矩阵的坐标形式解刚度矩阵,使用全部矢量化方法解全局力矢量,计算方程系数,计算狄利克雷边界条件,解有限元矩阵和偏微分方程。
  2. 构件分析方法​:🎯矩阵分析:🖊平面桁架位移和形变状态:提取结构的刚度矩阵,绘制变形形状,节点位移后,计算杆件的内力和桁架 | 🖊 空间桁架位移和形变状态:提取结构的刚度矩阵,计算节点载荷下 | 🖊平面框架位移和形变状态:提取局部刚度矩阵和旋转矩阵,计算出内力,计算钢框架承载下 | 🖊空间框架位移和形变状态:提取结构刚度矩阵,绘制变形形状,计算钢框架承载下,计算承受横向载荷和力矩的空间钢框架 | 🖊网格结构位移和形变:提取结构刚度矩阵,绘制形变形状,计算承受集中力矩和点载荷,计算钢具有相同的圆形横截面网架结构承受两个集中力矩 | 🖊构件载荷位移和形变:承受集中和均匀荷载的钢网架结构 | 🖊支撑沉降平面框架位移:外力影响下形变 | 🖊温度因素产生的构件位移:构件的拉伸或压缩的应力 | 🖊非标件框架位移:提取结构的刚度矩阵,计算给定载荷下 | 🖊释放构件位移:节点载荷情况下 | 🖊弹性支撑构件位移:均匀分布载荷的节点。🎯线性弹性分析:🖊双节点和三节点桁架:计算刚度矩阵 | 🖊弹簧和不同横截面的一维结构:计算外力条件下的位移,计算节点载荷下位移和法向应力 | 🖊平面桁架:计算外力下的形变 | 🖊网格二维结构:计算平面应力、平面应变和轴对称单元的刚度矩阵,计算节点力矢量,绘制形变 | 🖊平面应力计算情景:受拉伸表面载荷作用下,均匀弯曲和表面牵引载荷作用下,平面应变结构,轴对称结构。🎯弹塑性分析 | 🎯有限变形和超弹性 | 🎯有限应变。

🍇Python平面应力下变形位移求解器

考虑长度为 50 毫米、宽度为 10 毫米、厚度为 1 毫米的条带。使用 200 N 的力,该条带在垂直的平面应变张力下变形。条带的底部保持固定。带材的弹性特性为E=100MPa,ν=0.48。

在此,我使用线性弹性方程:
σ i j = 2 μ ε i j + λ ε k k δ i j \sigma_{i j}=2 \mu \varepsilon_{i j}+\lambda \varepsilon_{k k} \delta_{i j} σij=2μεij+λεkkδij
拉梅常数可以根据已知的弹性性质来确定。通过分配平面应变并假设 1 方向的自由位移,我可以执行快速计算,得出:最大顶部位移 = 8 毫米。

import numpy as np
import math
from matplotlib import pyplot as plt
def shape(xi):
	x,y = tuple(xi)
	N = [(1.0-x)*(1.0-y), (1.0+x)*(1.0-y), (1.0+x)*(1.0+y), (1.0-x)*(1.0+y)]
	return 0.25*np.array(N)

print('create mesh')

mesh_ex = 9
mesh_ey = 49
mesh_lx = 10.0
mesh_ly = 50.0
# derived
mesh_nx      = mesh_ex + 1
mesh_ny      = mesh_ey + 1
num_nodes    = mesh_nx * mesh_ny
num_elements = mesh_ex * mesh_ey
mesh_hx      = mesh_lx / mesh_ex
mesh_hy      = mesh_ly / mesh_ey
nodes = []
for y in np.linspace(0.0, mesh_ly, mesh_ny):
	for x in np.linspace(0.0, mesh_lx, mesh_nx):
		nodes.append([x,y])
nodes = np.array(nodes)
conn = []
for j in range(mesh_ey):
	for i in range(mesh_ex):
		n0 = i + j*mesh_nx
		conn.append([n0, n0 + 1, n0 + 1 + mesh_nx, n0 + mesh_nx])

print ('material model - plane strain')
E = 100.0
v = 0.48
C = E/(1.0+v)/(1.0-2.0*v) * np.array([[1.0-v,     v,     0.0],
								      [    v, 1.0-v,     0.0],
								      [  0.0,   0.0,   0.5-v]])

print('create global stiffness matrix')
K = np.zeros((2*num_nodes, 2*num_nodes))
q4 = [[x/math.sqrt(3.0),y/math.sqrt(3.0)] for y in [-1.0,1.0] for x in [-1.0,1.0]]
B = np.zeros((3,8))
for c in conn:
	xIe = nodes[c,:]
	Ke = np.zeros((8,8))
	for q in q4:
		dN = gradshape(q)
		J  = np.dot(dN, xIe).T
		dN = np.dot(np.linalg.inv(J), dN)
		B[0,0::2] = dN[0,:]
		B[1,1::2] = dN[1,:]
		B[2,0::2] = dN[1,:]
		B[2,1::2] = dN[0,:]
		Ke += np.dot(np.dot(B.T,C),B) * np.linalg.det(J)
	for i,I in enumerate(c):
		for j,J in enumerate(c):
			K[2*I,2*J]     += Ke[2*i,2*j]
			K[2*I+1,2*J]   += Ke[2*i+1,2*j]
			K[2*I+1,2*J+1] += Ke[2*i+1,2*j+1]
			K[2*I,2*J+1]   += Ke[2*i,2*j+1]

print('assign nodal forces and boundary conditions')
f = np.zeros((2*num_nodes))
for i in range(num_nodes):
	if nodes[i,1] == 0.0:
		K[2*i,:]     = 0.0
		K[2*i+1,:]   = 0.0
		K[2*i,2*i]   = 1.0
		K[2*i+1,2*i+1] = 1.0
	if nodes[i,1] == mesh_ly:
		x = nodes[i,0]
		f[2*i+1] = 20.0
		if x == 0.0 or x == mesh_lx:
			f[2*i+1] *= 0.5

print('solving linear system')
u = np.linalg.solve(K, f)
print('max u=', max(u))
print('plotting displacement')
ux = np.reshape(u[0::2], (mesh_ny,mesh_nx))
uy = np.reshape(u[1::2], (mesh_ny,mesh_nx))
xvec = []
yvec = []
res  = []
for i in range(mesh_nx):
    for j in range(mesh_ny):
        xvec.append(i*mesh_hx + ux[j,i])
        yvec.append(j*mesh_hy + uy[j,i])
        res.append(uy[j,i])
t = plt.tricontourf(xvec, yvec, res, levels=14, cmap=plt.cm.jet)
plt.scatter(xvec, yvec, marker='o', c='b', s=2)
plt.grid()
plt.colorbar(t)
plt.axis('equal')
plt.show()

下图显示了运行 Python 代码的结果,最大位移为6.75 mm,与上面的近似计算结果一致。

参阅一:计算思维

参阅二:亚图跨际

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值