提示:文章写完后,目录可以自动生成,如何生成可参考右边的帮助文档
解决bug:IndexError: index 1 is out of bounds for axis 0 with size 1
一、有限元算法
def Bernoulli3DFEA(Node, Element, Material, force, constrain, distributedLoad):
"""
计算三维Bernoulli梁有限元分析中每个节点的位移:
# 详细说明:
初始化全局刚度矩阵和载荷向量。
遍历每个单元,计算局部刚度矩阵,并将其添加到全局刚度矩阵。
应用载荷和约束条件。
解线性方程组以获得位移向量。
# 参数:
Node: 节点坐标的NumPy数组。
Element: 单元节点索引的NumPy数组。
Material: 材料属性的NumPy数组。
force: 集中载荷信息的NumPy数组。
constrain: 约束条件的NumPy数组。
distributedLoad: 均布载荷信息的NumPy数组。
# 返回:
U: 每个节点的位移的NumPy数组。
"""
numDofPerNode = 6
numDofTotal = Node.shape[0] * numDofPerNode
# 初始化全局刚度矩阵和力矢量
KKG = sp.lil_matrix((numDofTotal, numDofTotal), dtype=float)
FFG = np.zeros((numDofTotal, 1), dtype=float)
# 遍历每个元素以计算局部刚度矩阵
for iEle in range(Element.shape[0]):
Node1 = int(Element[iEle, 0]) - 1
Node2 = int(Element[iEle, 1]) - 1
# 确保正确访问材质数组
if Material.shape[1] != 6:
raise ValueError("Material array should have exactly 6 columns corresponding to E, G, A, Iy, Iz, and J.")
E = Material[iEle, 0]
G = Material[iEle, 1]
A = Material[iEle, 2]
Iy = Material[iEle, 3]
Iz = Material[iEle, 4]
J = Material[iEle, 5]
# 计算节点坐标和元素长度
x1, y1, z1 = Node[Node1]
x2, y2, z2 = Node[Node2]
L = np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2 + (z2 - z1) ** 2)
# 局部刚度矩阵
KLocal = np.zeros((12, 12), dtype=float)
KLocal[0, 0] = KLocal[6, 6] = A * E / L
KLocal[0, 6] = KLocal[6, 0] = -A * E / L
KLocal[1, 1] = KLocal[7, 7] = 12 * E * Iz / L ** 3
KLocal[1, 7] = KLocal[7, 1] = -12 * E * Iz / L ** 3
KLocal[1, 5] = KLocal[5, 1] = 6 * E * Iz / L ** 2
KLocal[1, 11] = KLocal[11, 1] = 6 * E * Iz / L ** 2
KLocal[2, 2] = KLocal[8, 8] = 12 * E * Iy / L ** 3
KLocal[2, 8] = KLocal[8, 2] = -12 * E * Iy / L ** 3
KLocal[2, 4] = KLocal[4, 2] = -6 * E * Iy / L ** 2
KLocal[2, 10] = KLocal[10, 2] = -6 * E * Iy / L ** 2
KLocal[3, 3] = KLocal[9, 9] = G * J / L
KLocal[3, 9] = KLocal[9, 3] = -G * J / L
KLocal[4, 4] = KLocal[10, 10] = 4 * E * Iy / L
KLocal[4, 10] = KLocal[10, 4] = 2 * E * Iy / L
KLocal[5, 5] = KLocal[11, 11] = 4 * E * Iz / L
KLocal[5, 11] = KLocal[11, 5] = 2 * E * Iz / L
KLocal[4, 8] = KLocal[8, 4] = 6 * E * Iy / L ** 2
KLocal[5, 7] = KLocal[7, 5] = -6 * E * Iz / L ** 2
# 其他刚度矩阵项
KLocal[6, 7] = 12 * E * Iz / L ** 3
KLocal[6, 11] = -6 * E * Iz / L ** 2
KLocal[7, 11] = -12 * E * Iz / L ** 3
KLocal[8, 10] = -6 * E * Iy / L ** 2
KLocal[4, 2] = KLocal[2, 4] = -6 * E * Iy / L ** 2
KLocal[5, 11] = KLocal[11, 5] = 2 * E * Iz / L
# 将上半部分镜像到下半部分以获得对称的刚度矩阵
for i in range(6, 12):
for j in range(6, 12):
KLocal[i, j] = KLocal[i - 6, j - 6]
# 变换矩阵
l = (x2 - x1) / L
m = (y2 - y1) / L
n = (z2 - z1) / L
T = np.array([[l, m, n, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-m, l, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-n, 0, l, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, l, m, n, 0, 0, 0, 0, 0, 0],
[0, 0, 0, -m, l, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, -n, 0, l, 0, 0, 0, 0, 0, 0]])
# 全局刚度矩阵
T = np.block([[T, np.zeros((6, 6))], [np.zeros((6, 6)), T]])
KGlobal = T.T @ KLocal @ T
# 全局自由度
index = np.r_[np.arange(6 * Node1, 6 * Node1 + 6), np.arange(6 * Node2, 6 * Node2 + 6)]
# 装配到全局刚度矩阵
for i in range(12):
for j in range(12):
KKG[index[i], index[j]] += KGlobal[i, j]
# 应用分布式负载
for distLoad in distributedLoad:
elemIdx, dir, mag = distLoad
elemIdx = int(elemIdx) - 1 # 转换为从 0 开始的索引
Node1 = int(Element[elemIdx, 0]) - 1
Node2 = int(Element[elemIdx, 1]) - 1
L = np.linalg.norm(Node[Node2] - Node[Node1])
# 在指定方向上分配负载
if dir == 1: # X-direction
loadVector = mag * L / 2 * np.array([1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0])
elif dir == 2: # Y-direction
loadVector = mag * L / 2 * np.array([0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0])
elif dir == 3: # Z-direction
loadVector = mag * L / 2 * np.array([0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0])
else:
continue
# 将分布荷载应用于全局力矢量
FFG[index] += loadVector.reshape(-1, 1)
# 施加集中负载
for iLoad in range(force.shape[0]):
NodeDof = (int(force[iLoad, 0]) - 1) * 6 + int(force[iLoad, 1]) - 1
FFG[NodeDof] = force[iLoad, 2]
# 应用约束
KKG, FFG = apply_constraints(KKG, FFG, constrain)
# 求解位移
U = sla.spsolve(KKG.tocsr(), FFG)
return U
二、bug 解析
为了解决 IndexError: index 1 is out of bounds for axis 0 with size 1 的问题,我们需要确保 Material、Element 或其他数组的大小与它们在代码中的使用保持一致。这个错误通常是因为数组的访问超出了其定义的范围。
三、可能的原因
1.元素数量与材料属性不匹配:如果 Element 数组中的单元数量多于 Material 数组中的行数,就会发生索引超界的问题。
2.访问未定义的索引:例如,尝试访问 Material[iEle] 时,iEle 的值可能比 Material 的行数多。
四、解决方案
1.确保每个 Element 都有相应的 Material 数据。
2.检查循环中的索引是否在数组的范围内。
import numpy as np
# 节点坐标 (x, y, z)
Node = np.array([
[0, 0, 0], # Node 1
[1, 0, 0], # Node 2
[1, 1, 0], # Node 3
[0, 1, 0] # Node 4
])
# 单元定义 (连接节点索引)
Element = np.array([
[1, 2, 3, 4] # Element 1 connects Node 1, 2, 3, and 4
])
# 材料属性 (E, ν)
Material = np.array([
[210000, 0.3] # Material for Element 1
])
# 节点载荷 (节点索引, Fx, Fy, Fz)
Load = np.array([
[1, 1000, 0, 0], # Node 1: Fx = 1000N
[2, 0, 0, 0], # Node 2: No force
[3, 0, 0, 0], # Node 3: No force
[4, 0, 0, 0] # Node 4: No force
])
# 节点约束 (节点索引, U, V, W)
Constr = np.array([
[1, 1, 1, 1], # Node 1: Fully constrained (U, V, W)
[2, 0, 1, 1], # Node 2: Constrained in V, W
[3, 1, 1, 0], # Node 3: Constrained in U, V
[4, 1, 0, 1] # Node 4: Constrained in U, W
])
# 分布载荷 (单元索引, 方向, qx, qy, qz)
DistributedLoad = np.array([
[1, 1, 100, 0, 0] # Element 1, q = 100N/m in x-direction
])
def Bernoulli3DFEA(Node, Element, Material, Load, Constr, DistributedLoad):
# 输出 Material 的形状
print(f"Material shape: {Material.shape}")
# 确保 Element 和 Material 的数量匹配
if len(Element) != len(Material):
print("Error: The number of elements does not match the number of material entries.")
return
# 遍历元素,访问其材料属性
for iEle in range(len(Element)):
try:
E = Material[iEle, 0]
ν = Material[iEle, 1]
print(f"Element {iEle + 1} Material: E={E}, ν={ν}")
except IndexError:
print(f"Error: Material index {iEle} out of range.")
return
# 返回某种形式的结果
U = np.zeros(Node.shape)
return U
# 运行函数
U = Bernoulli3DFEA(Node, Element, Material, Load, Constr, DistributedLoad)
print(f"Displacements:\n{U}")
1. 解释
1.匹配 Element 和 Material 的数量:通过在 Bernoulli3DFEA 函数中检查 Element 和 Material 的数量是否匹配,确保不会访问超出范围的材料属性。
2.使用 try 和 except 捕获 IndexError:如果发生索引超界错误,将打印错误信息并返回。这样可以帮助识别输入中的问题。
3.错误检查:在访问数组之前,确保索引在数组的范围内。
2. 测试输出
Material shape: (1, 2)
Element 1 Material: E=210000, ν=0.3
Displacements:
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
3. 如何添加更多单元
# 新的单元定义
Element = np.array([
[1, 2, 3, 4], # Element 1
[2, 3, 4, 1] # Element 2
])
# 对应的材料属性
Material = np.array([
[210000, 0.3], # Material for Element 1
[210000, 0.3] # Material for Element 2
])
重点!: 注意每个索引都有对应的匹配项
在这个扩展示例中,Element 数组有两个单元,每个单元都有相应的 Material 属性。这种匹配确保了不会发生索引超界错误。