利用拓扑优化ESO算法进行悬臂梁的优化代码示例以及给出简化图

以下是一个基于拓扑优化ESO算法的悬臂梁优化的代码示例：

import numpy as np
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve

# 定义常数
E = 1  # 弹性模量
nu = 0.3  # 泊松比
rho = 1  # 密度
l = 1  # 梁长度
w = 0.1  # 梁宽度
h = 0.1  # 梁高度
p = -1  # 外载荷

# 定义初始布局
n_elx, n_ely = 60, 10  # 元素个数
nel = n_elx * n_ely  # 总元素数
nnx, nny = n_elx + 1, n_ely + 1  # 节点数
nn = nnx * nny  # 总节点数
ndof = 2  # 自由度个数
nnp = nn * ndof  # 总自由度个数
x = np.zeros(nnp)

# 定义初始密度分布
volfrac = 0.5  # 初始密度
x[:nnp//2] = volfrac

# 定义有限元模型
KE = np.array([[12, 3, -6, -3], [3, 12, 3, -6], [-6, 3, 12, -3], [-3, -6, -3, 12]]) / (12 * (1 - nu ** 2) * h * w)
KE1 = np.array([[-6, -3, 6, -3], [-3, 0, 3, 0], [6, 3, -6, 3], [-3, 0, 3, 0]]) / (2 * l)
KE2 = np.array([[12, -6, 3, -12], [-6, 4, -2, 6], [3, -2, 4, -3], [-12, 6, -3, 12]]) * E * h ** 3 / (12 * (1 - nu ** 2) * l ** 3)

edofMat = np.zeros((nel, 8), dtype=int)
for ely in range(n_ely):
    for elx in range(n_elx):
        el = ely * n_elx + elx
        n1 = (n_elx + 1) * ely + elx
        n2 = (n_elx + 1) * (ely + 1) + elx
        edofMat[el] = np.array([2 * n1, 2 * n1 + 1, 2 * n2, 2 * n2 + 1, 2 * (n1 + 1), 2 * (n1 + 1) + 1, 2 * (n2 + 1), 2 * (n2 + 1) + 1])

# 定义ESO算法
def ElementStiffness(x):
    # 计算元素刚度矩阵
    penal = 3
    KE_penalty = penal * x ** (penal - 1) * KE
    return KE_penalty

def ElementMass(x):
    # 计算元素质量矩阵
    rho_penalty = rho * x ** penal
    return rho_penalty * KE2

def ComputeK(x):
    # 计算总刚度矩阵
    K = lil_matrix((nnp, nnp))
    for el in range(nel):
        edof = edofMat[el]
        KE_penalty = ElementStiffness(x[edof[:4]])
        K[np.ix_(edof[:4], edof[:4])] += KE_penalty
        K[np.ix_(edof[:4], edof[4:])] -= KE1
        K[np.ix_(edof[4:], edof[:4])] -= KE1.T
        ME_penalty = ElementMass(x[edof[:4]])
        K[np.ix_(edof[:4], edof[:4])] += ME_penalty
    return K

def ESO(x, Emin=1e-9, Emax=1e9, rho_min=0.001, rho_max=1, beta=1):
    # 进行ESO优化
    loop = 0
    change = 1
    while change > 0.01 and loop < 100:
        loop += 1
        x_old = x.copy()
        K = ComputeK(x)
        U = spsolve(K, p)
        for el in range(nel):
            edof = edofMat[el]
            x[edof[:4]] = max(Emin, min(Emax, x[edof[:4]] + beta * (p @ U[edof[:4]]) / (KE @ U[edof[:4]])))
        x = np.minimum(np.maximum(x, rho_min), rho_max)
        change = np.linalg.norm(x - x_old, np.inf)

    return x

# 进行优化
x_opt = ESO(x)

# 绘制简化图
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

fig, ax = plt.subplots()

for ely in range(n_ely):
    for elx in range(n_elx):
        el = ely * n_elx + elx
        xc = (elx + 0.5) / n_elx
        yc = (ely + 0.5) / n_ely
        if x_opt[edofMat[el, :4]].sum() / 4 > 1e-4:
            ax.add_patch(Rectangle((xc - 0.5 / n_elx, yc - 0.5 / n_ely), 1 / n_elx, 1 / n_ely, facecolor='black'))
        else:
            ax.add_patch(Rectangle((xc - 0.5 / n_elx, yc - 0.5 / n_ely), 1 / n_elx, 1 / n_ely, facecolor='white'))
        
ax.set_xlim([0, 1])
ax.set_ylim([0, 1])
ax.set_aspect('equal')
plt.show()

下面是一个简化的悬臂梁示意图，其中黑色矩形表示材料存在的区域，白色矩形表示空洞区域。