基于拓扑优化ESO算法的悬臂梁优化的Python代码示例

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

import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve
from numba import jit

# 定义悬臂梁的几何参数和材料属性
L = 1.0  # 悬臂梁长度
H = 0.2  # 悬臂梁高度
W = 0.01  # 悬臂梁宽度
E = 1e7  # 弹性模量
nu = 0.3  # 泊松比
rho = 2600  # 密度
g = 9.8  # 重力加速度
q = 10000  # 集中载荷
nelx = 60  # 横向单元数
nely = 12  # 纵向单元数
volfrac = 0.5  # 材料体积分数
penal = 3.0  # 材料惩罚系数
rmin = 3.0  # 最小半径
ft = 2  # 滤波类型，1：Sensitivity，2：Density

# 计算单元尺寸和网格尺寸
dx = L / nelx
dy = H / nely
n = (nelx + 1) * (nely + 1)
m = nelx * nely

# 构建节点和单元编号
inds = np.arange(n).reshape((nely+1, nelx+1))
elems = np.zeros((m, 4), dtype=int)
for ely in range(nely):
    for elx in range(nelx):
        elem = ely * nelx + elx
        n1 = inds[ely, elx]
        n2 = inds[ely, elx+1]
        n3 = inds[ely+1, elx+1]
        n4 = inds[ely+1, elx]
        elems[elem] = [n1, n2, n3, n4]

# 初始化密度和位移场
rho0 = np.ones(m) * volfrac
x = np.ones(m) * volfrac

# 定义滤波器
def get_kernel(ft):
    if ft == 1:
        kernel = np.array([[0, 1, 0],
                           [1, 1, 1],
                           [0, 1, 0]])
    elif ft == 2:
        kernel = np.array([[0, 1, 0],
                           [1, 4, 1],
                           [0, 1, 0]])
    else:
        raise ValueError('Invalid filter type')
    return kernel

kernel = get_kernel(ft)

# 定义有限元分析函数
@jit(nopython=True)
def fem(x, rho0, penal):
    K = np.zeros((n, n))
    f = np.zeros(n)
    for elem in range(m):
        n1, n2, n3, n4 = elems[elem]
        x1, y1 = n1 % (nelx+1), n1 // (nelx+1)
        x2, y2 = n2 % (nelx+1), n2 // (nelx+1)
        x3, y3 = n3 % (nelx+1), n3 // (nelx+1)
        x4, y4 = n4 % (nelx+1), n4 // (nelx+1)
        X = np.array([[x1, x2, x3, x4],
                      [y1, y2, y3, y4]])
        B = np.array([[-1, 1, 0, 0],
                      [0, 0, 1, -1],
                      [1, -1, -1, 1],
                      [0, 0, -1, 1],
                      [-1, 0, 1, 0],
                      [0, -1, 0, 1]])
        D = E / (1 - nu**2) * np.array([[1, nu, 0],
                                        [nu, 1, 0],
                                        [0, 0, (1 - nu) / 2]])
        detX = np.linalg.det(X)
        invX = np.linalg.inv(X)
        Bhat = np.dot(invX.T, B)
        Ke = np.dot(np.dot(Bhat.T, D), Bhat) * detX * x[elem]**penal
        K[n1, n1] += Ke[0, 0]
        K[n1, n2] += Ke[0, 1]
        K[n1, n3] += Ke[0, 2]
        K[n1, n4] += Ke[0, 3]
        K[n2, n1] += Ke[1, 0]
        K[n2, n2] += Ke[1, 1]
        K[n2, n3] += Ke[1, 2]
        K[n2, n4] += Ke[1, 3]
        K[n3, n1] += Ke[2, 0]
        K[n3, n2] += Ke[2, 1]
        K[n3, n3] += Ke[2, 2]
        K[n3, n4] += Ke[2, 3]
        K[n4, n1] += Ke[3, 0]
        K[n4, n2] += Ke[3, 1]
        K[n4, n3] += Ke[3, 2]
        K[n4, n4] += Ke[3, 3]
        f[elem] = rho0[elem] * g * detX * x[elem]
    fixed_nodes = np.where((inds == 0) | (inds == nelx))[0]
    free_nodes = np.setdiff1d(np.arange(n), fixed_nodes)
    Kff = K[free_nodes][:, free_nodes]
    Kfc = K[free_nodes][:, fixed_nodes]
    u = np.zeros(n)
    uf = spsolve(Kff, -Kfc.dot(u[fixed_nodes]))
    u[free_nodes] = uf
    s = np.dot(K, u) - f
    return u, s

# 定义灵敏度分析函数
@jit(nopython=True)
def sensitivity(x, rho0, penal):
    dc = np.zeros(m)
    for elem in range(m):
        n1, n2, n3, n4 = elems[elem]
        x1, y1 = n1 % (nelx+1), n1 // (nelx+1)
        x2, y2 = n2 % (nelx+1), n2 // (nelx+1)
        x3, y3 = n3 % (nelx+1), n3 // (nelx+1)
        x4, y4 = n4 % (nelx+1), n4 // (nelx+1)
        X = np.array([[x1, x2, x3, x4],
                      [y1, y2, y3, y4]])
        B = np.array([[-1, 1, 0, 0],
                      [0, 0, 1, -1],
                      [1, -1, -1, 1],
                      [0, 0, -1, 1],
                      [-1, 0, 1, 0],
                      [0, -1, 0, 1]])
        D = E / (1 - nu**2) * np.array([[1, nu, 0],
                                        [nu, 1, 0],
                                        [0, 0, (1 - nu) / 2]])
        detX = np.linalg.det(X)
        invX = np.linalg.inv(X)
        Bhat = np.dot(invX.T, B)
        Ke = np.dot(np.dot(Bhat.T, D), Bhat) * detX * x[elem]**penal
        s = np.dot(Ke, u[[n1, n2, n3, n4]])
        dc[elem] = -rho0[elem] * x[elem]**(penal-1) * s
    return dc

# 定义拓扑优化函数
def topo_opt(rho0, x, volfrac, rmin, penal):
    loop = 0
    change = 1
    while change > 0.01 and loop < 200:
        # 有限元分析
        u, s = fem(x, rho0, penal)
        # 灵敏度分析
        dc = sensitivity(x, rho0, penal)
        # 计算滤波后的灵敏度
        if ft == 1:
            dc[:] = np.dot(kernel, dc) / np.sum(kernel)
        elif ft == 2:
            kernel *= np.nanmax(dc) / np.sum(dc)
            dc[:] = np.dot(kernel, dc)
        # 计算移动限制
        rho = np.maximum(0, np.minimum(rho0 + 0.05, np.maximum(0, x - 0.05)))
        # 计算更新步长
        l1 = 0
        l2 = 1e9
        while l2 - l1 > 1e-4:
            lmid = 0.5 * (l2 + l1)
            xnew = np.maximum(0, np.maximum(x - lmid * dc, rho))
            if np.sum(xnew) - volfrac * m > 0:
                l1 = lmid
            else:
                l2 = lmid
        change = np.max(np.abs(x - xnew))
        x = xnew
        loop += 1
    return x

# 进行拓扑优化
x = topo_opt(rho0, x, volfrac, rmin, penal)

# 输出结果
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
for elem in range(m):
    n1, n2, n3, n4 = elems[elem]
    x1, y1 = n1 % (nelx+1), n1 // (nelx+1)
    x2, y2 = n2 % (nelx+1), n2 // (nelx+1)
    x3, y3 = n3 % (nelx+1), n3 // (nelx+1)
    x4, y4 = n4 % (nelx+1), n4 // (nelx+1)
    if x[elem] > 0.5:
        ax.fill([x1*dx, x2*dx, x3*dx, x4*dx], [H-y1*dy, H-y2*dy, H-y3*dy, H-y4*dy], 'k')
    else:
        ax.fill([x1*dx, x2*dx, x3*dx, x4*dx], [H-y1*dy, H-y2*dy, H-y3*dy, H-y4*dy], 'w')
ax.set_xlim(0, L)
ax.set_ylim(0, H)
ax.set_aspect('equal', 'box')
plt.show()

该代码使用了NumPy、SciPy和Numba库，可以通过pip安装。在运行前需要安装这些库，并将代码保存为.py文件，然后在命令行中运行python file.py即可。代码中已经调整了一些参数，可以根据实际需求进行修改。