以下是一个基于 Evolutionary Structural Optimization 算法的悬臂梁优化的 Python 代码示例:
```python
import numpy as np
import matplotlib.pyplot as plt
# 定义初始参数
L = 1.0 # 悬臂梁长度
H = 0.2 # 悬臂梁高度
W = 0.01 # 悬臂梁宽度
E = 1.0 # 杨氏模量
p = 1.0 # 外载荷
N = 50 # 网格数量
volfrac = 0.4 # 材料体积分数
penal = 3.0 # 惩罚因子
rmin = 1.5 # 最小半径
# 定义初始材料分布
x = np.linspace(0, L, N)
y = np.linspace(0, H, N)
X, Y = np.meshgrid(x, y)
rho = volfrac * np.ones((N, N))
# 定义有限元模型
nelx = N - 1
nely = N - 1
ndof = 2 * (N * (N + 1) // 2) # 节点自由度数
nele = nelx * nely # 单元数量
edof = 8 # 单元自由度数
dofs = np.arange(ndof)
# 定义单元节点编号
iK = np.array([
[0, 1, 2, 3],
[1, 4, 5, 2],
[2, 5, 6, 7],
[3, 2, 7, 8],
[2, 5, 8, 7],
[7, 8, 9, 10],
[8, 11, 12, 9],
[9, 12, 13, 14],
[10, 9, 14, 15],
[9, 12, 15, 14],
[14, 15, 16, 17],
[15, 18, 19, 16],
[16, 19, 20, 21],
[17, 16, 21, 22],
[16, 19, 22, 21]
])
# 定义单元刚度矩阵
KE = np.array([
[12, 3, -6, -12, 3, -6, -3, -3],
[3, 12, 12, -3, -6, 3, -3, -6],
[-6, 12, 24, -6, -12, 6, -6, -12],
[-12, -3, -6, 12, -3, -6, 3, 3],
[3, -6, -12, -3, 12, 3, -3, 6],
[-6, 3, 6, -6, 3, 12, -12, -3],
[-3, -3, -6, 3, -3, -12, 12, 6],
[-3, -6, -12, 3, 6, -3, 6, 12]
]) / (E * W * H ** 3 / 12)
# 定义加权函数
def get_w(x, y):
if x < 0.5 and y < 0.1:
return 10.0
else:
return 1.0
# 定义计算单元刚度矩阵的函数
def get_KE(x1, x2, x3, x4, y1, y2, y3, y4):
B = np.array([
[y2 - y4, 0, y4 - y3, 0, y3 - y2, 0, y1 - y3, 0, y2 - y1, 0, y4 - y2, 0],
[0, x4 - x2, 0, x3 - x4, 0, x2 - x3, 0, x1 - x3, 0, x4 - x1, 0, x2 - x4],
[x4 - x2, y2 - y4, x3 - x4, y4 - y3, x2 - x3, y3 - y2, x1 - x3, y3 - y1, x4 - x1, y1 - y4, x2 - x4, y4 - y2]
]) / ((x1 - x3) * (y2 - y4) - (x2 - x4) * (y1 - y3))
D = E * np.array([
[1, 0, 0],
[0, 1, 0],
[0, 0, 0.5]
])
return B.T @ D @ B * (x1 - x3) * (y2 - y4) / 2
# 定义计算刚度矩阵和载荷向量的函数
def get_K_and_f(rho):
K = np.zeros((ndof, ndof))
f = np.zeros(ndof)
for elx in range(nelx):
for ely in range(nely):
n1 = (N + 1) * ely + elx
n2 = (N + 1) * (ely + 1) + elx
n3 = (N + 1) * ely + elx + 1
n4 = (N + 1) * (ely + 1) + elx + 1
if rho[ely, elx] > 0.001:
edof_list = np.array([
2 * n1, 2 * n1 + 1, 2 * n2, 2 * n2 + 1,
2 * n3, 2 * n3 + 1, 2 * n4, 2 * n4 + 1
])
xe = np.array([x[elx], x[elx], x[elx + 1], x[elx + 1]])
ye = np.array([y[ely], y[ely + 1], y[ely], y[ely + 1]])
KE_el = get_KE(xe[iK], ye[iK])
K[np.ix_(edof_list, edof_list)] += KE_el * rho[ely, elx] ** penal
f[[2 * n2 + 1, 2 * n4 + 1]] += p * W * H / 2 * get_w(x[elx] + L / N / 2, y[ely] + H / N / 2) * rho[ely, elx] * (L / N) ** 2 / 2
return K, f
# 定义计算柔度函数的函数
def get_J(rho):
K, f = get_K_and_f(rho)
u = np.linalg.solve(K, f)
return np.dot(u, f)
# 定义计算灵敏度函数的函数
def get_sensitivity(rho):
K, f = get_K_and_f(rho)
u = np.linalg.solve(K, f)
dsdrho = -penal * rho ** (penal - 1)
KE_list = np.array([get_KE(xe[iK], ye[iK]) for xe, ye in zip(X.flatten()[dofs.reshape(-1, edof)[:, :-1]], Y.flatten()[dofs.reshape(-1, edof)[:, :-1]])])
B_list = np.array([np.array([
[y2 - y4, 0, y4 - y3, 0, y3 - y2, 0, y1 - y3, 0, y2 - y1, 0, y4 - y2, 0],
[0, x4 - x2, 0, x3 - x4, 0, x2 - x3, 0, x1 - x3, 0, x4 - x1, 0, x2 - x4],
[x4 - x2, y2 - y4, x3 - x4, y4 - y3, x2 - x3, y3 - y2, x1 - x3, y3 - y1, x4 - x1, y1 - y4, x2 - x4, y4 - y2]
]) / ((x1 - x3) * (y2 - y4) - (x2 - x4) * (y1 - y3))] for x1, x2, x3, x4, y1, y2, y3, y4 in zip(X.flatten()[dofs.reshape(-1, edof)[:, :-1]], X.flatten()[dofs.reshape(-1, edof)[:, 1:]], X.flatten()[dofs.reshape(-1, edof)[:, 2:]], X.flatten()[dofs.reshape(-1, edof)[:, 3:]], Y.flatten()[dofs.reshape(-1, edof)[:, :-1]], Y.flatten()[dofs.reshape(-1, edof)[:, 1:]], Y.flatten()[dofs.reshape(-1, edof)[:, 2:]], Y.flatten()[dofs.reshape(-1, edof)[:, 3:]])])
s = np.sum(B_list * (K @ u)[dofs.reshape(-1, edof)], axis=(1, 2)) * dsdrho * (L / N) ** 2
return s.reshape((N, N))
# 定义计算设计变量的函数
def get_design_variable(rho):
eta = 0.5
l1 = 0.0
l2 = 1e9
move = 0.2
while l2 - l1 > 1e-4:
lmid = 0.5 * (l1 + l2)
rho_new = np.maximum(0, np.maximum(rho - move, np.minimum(1, np.minimum(rho + move, rho * np.sqrt(-get_sensitivity(rho) / lmid / eta)))))
if np.sum(rho_new) - volfrac * N ** 2 > 0:
l1 = lmid
else:
l2 = lmid
return rho_new
# 开始迭代
for it in range(100):
rho = get_design_variable(rho)
J = get_J(rho)
print('Iteration:', it, ', Objective:', J)
plt.clf()
plt.imshow(-rho.T, cmap='gray', extent=[0, L, 0, H])
plt.colorbar()
plt.savefig('iteration_%03d.png' % it)
```
在这个示例中,我们定义了初始参数和初始材料分布,然后定义了有限元模型和单元刚度矩阵。接下来,我们定义了加权函数、计算单元刚度矩阵的函数、计算刚度矩阵和载荷向量的函数、计算柔度函数的函数和计算灵敏度函数的函数。最后,我们定义了计算设计变量的函数,并开始迭代。在每次迭代中,我们使用计算出的设计变量计算并更新柔度函数和灵敏度函数,然后使用设计变量更新材料分布。我们迭代了 100 次,并在每次迭代中保存材料分布的图像。