利用拓扑优化ESO算法进行悬臂梁的优化代码示例 使用Python语言实现

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

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

# 定义悬臂梁问题的参数
nelx, nely = 80, 40  # 模型的尺寸
volfrac = 0.4  # 材料的体积分数
penal = 3.0  # 惩罚参数p
rmin = 3.0  # 最小半径r

# 初始化设计变量和有限元分析参数
x = volfrac * np.ones(nely * nelx, dtype=float)
xmin = np.zeros(nely * nelx, dtype=float)
xmax = np.ones(nely * nelx, dtype=float)
dc = np.zeros(nely * nelx, dtype=float)
ke = np.array([[1, 0, -1, 0], [0, 0, 0, 0], [-1, 0, 1, 0], [0, 0, 0, 0]])  # 单元刚度矩阵
edofMat = np.zeros((nelx * nely, 4), dtype=int)
iK = np.zeros((2 * (nelx + 1) * (nely + 1),), dtype=int)
jK = np.zeros((8 * nelx * nely,), dtype=int)
K = np.zeros((8 * nelx * nely,), dtype=float)
F = np.zeros((2 * (nely + 1) * (nelx + 1),), dtype=float)
U = np.zeros((2 * (nely + 1) * (nelx + 1),), dtype=float)
fixeddofs = np.union1d(np.where(np.abs(x - 1) < 1e-5)[0],
                        np.array([2 * (nely + 1) * (nelx + 1) - 2, 2 * (nely + 1) * (nelx + 1) - 1]))

# 定义有限元分析的函数
def FE(x):
    # 初始化刚度矩阵和载荷向量
    K *= 0
    F *= 0

    for elx in range(nelx):
        for ely in range(nely):
            n1 = (nely + 1) * elx + ely
            n2 = (nely + 1) * (elx + 1) + ely
            edof = np.array([2 * n1, 2 * n1 + 1, 2 * n2, 2 * n2 + 1], dtype=int)
            edofMat[nely * elx + ely, :] = edof
            K00 = x[n1] ** penal
            K01 = x[n1] * x[n2] ** penal
            K10 = K01
            K11 = x[n2] ** penal
            K[np.ix_(edof, edof)] += ke * np.array([[K00, K01], [K10, K11]])

    # 设置载荷向量
    F[1] = -1

    # 移除固定自由度
    K, F = K[np.ix_(np.setdiff1d(range(2 * (nely + 1) * (nelx + 1)), fixeddofs),),
            np.ix_(np.setdiff1d(range(2 * (nely + 1) * (nelx + 1)), fixeddofs),)] \
        , F[np.setdiff1d(range(2 * (nely + 1) * (nelx + 1)), fixeddofs)]

    # 解算有限元方程
    U[np.setdiff1d(range(2 * (nely + 1) * (nelx + 1)), fixeddofs),] = spsolve(K, F)

    # 计算材料敏感度分布
    for elx in range(nelx):
        for ely in range(nely):
            n1 = (nely + 1) * elx + ely
            n2 = (nely + 1) * (elx + 1) + ely
            edof = np.array([2 * n1, 2 * n1 + 1, 2 * n2, 2 * n2 + 1], dtype=int)
            u_e = U[edof]
            dc[n1] += (x[n1] ** penal) * np.dot(u_e.T, np.dot(ke, u_e))
            dc[n2] += (x[n2] ** penal) * np.dot(u_e.T, np.dot(ke, u_e))

    return np.sum((F * U)[0])

# 定义拓扑优化的函数
def topopt(nelx, nely, volfrac, penal, rmin):
    xold = np.ones(nely * nelx, dtype=float)
    x = np.ones(nely * nelx, dtype=float)
    loop = 0
    change = 1
    while change > 0.01 and loop < 50:
        loop += 1
        xold[:] = x

        # 计算材料敏感度分布
        dc *= 0
        obj = FE(x)

        # 计算最小密度和最大密度
        xmin = np.zeros(nely * nelx, dtype=float)
        xmax = np.ones(nely * nelx, dtype=float)
        xPhys = x.copy()
        for i in range(50):
            # 计算移动平均滤波器
            xPhys = np.minimum(np.maximum(xPhys, 0.001), 0.999)
            if rmin > 0:
                nfilter = int(np.ceil(2 * rmin / 2))
                H = np.zeros((nfilter * 2 + 1, nfilter * 2 + 1))
                for i1 in range(-nfilter, nfilter + 1):
                    for j1 in range(-nfilter, nfilter + 1):
                        r = np.sqrt(i1 ** 2 + j1 ** 2)
                        if r <= rmin:
                            H[i1 + nfilter, j1 + nfilter] = 0.5 * (rmin + np.cos(np.pi * r / rmin) * rmin)
                        else:
                            H[i1 + nfilter, j1 + nfilter] = 0.5 * (rmin + 1.5 * (r - rmin))
                H /= np.sum(H)

                # 应用移动平均滤波器
                xPhys = np.asarray([np.sum(H * xPhys[il: il + 2 * nfilter + 1, jl: jl + 2 * nfilter + 1])
                                    for il in range(nelx) for jl in range(nely)])

            # 更新设计变量
            xnew = np.zeros(nely * nelx, dtype=float)
            for i in range(nely * nelx):
                l = np.maximum(i - nely, 0)
                r = np.minimum(i + nely, nely * nelx - 1)
                t = np.maximum(i % nely - 1, 0)
                b = np.minimum(i % nely + 1, nely - 1)
                nbrs = [i, l, r, t, b]
                dcnbrs = dc[nbrs]
                xnbrs = xPhys[nbrs]
                if x[i] < 0.001:
                    xnew[i] = 0
                elif x[i] > 0.999:
                    xnew[i] = 1
                else:
                    # 计算敏感度加权平均密度
                    xnew[i] = np.sum(xnbrs * dcnbrs) / np.sum(dcnbrs)

            # 通过投影算法确保密度在[0,1]之间
            xnew = np.minimum(np.maximum(xnew, 0), 1)
            change = np.max(np.abs(x - xnew))
            x[:] = xnew

        # 打印迭代信息
        print("Iteration: {0:3d}  Obj: {1:.4f}  Vol: {2:.4f}  ch: {3:.4f}".format(
            loop, obj, np.mean(x), change))

    return x

# 运行拓扑优化
xOpt = topopt(nelx, nely, volfrac, penal, rmin)