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

关键词

MATLAB | Python | 数学 | 物理 | 受力 | 载荷 | 构件 | 钢结构 | 非标件 | 节点 | 位移 | 内力 | 元素 | 几何 | 积分 | 算法 | 矩阵 | 网格 | 二维 | 一维 | 边界值 | 绘图 | 矢量 | 稀疏矩阵 | 刚度 | 平面 | 空间 | 桁架 | 框架 | 承载 | 形变 | 力矩 | 横截面 | 沉降 | 支撑 | 释放 | 弹性 | 线性 | 外力 | 法向应力 | 应变 | 轴对称 | 拉伸 | 弯曲 | 牵引 | 弹塑性 | 微分方程 | 数值

🏈指点迷津 | Brief

🎯要点

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

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

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

在此，我使用线性弹性方程：

\sigma_{i j}=2 \mu \varepsilon_{i j}+\lambda \varepsilon_{k k} \delta_{i j}

拉梅常数可以根据已知的弹性性质来确定。通过分配平面应变并假设 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，与上面的近似计算结果一致。

PreviousMATLAB雨刮通风空调模糊器和发电厂电力聚变器卷积神经 NextPython氮氧甲烷乙烷乙烯丙烯气体和固体热力学模型计算

Last updated 6 months ago

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()