【科研重器】从MATLAB到Python:Abaqus三维随机粗糙表面建模全攻略-致敬经典 Rand_surf,手把手教你实现高效六面体网格自动化生成

# -*- coding: utf-8 -*-"""Abaqus 3D 随机粗糙表面建模脚本参考来源：MATLAB 脚本 Rand_surf (Youngbin LIM)功能：基于 Weibull 分布和样条插值生成具有随机粗糙顶面的六面体实体模型"""from abaqus import *from abaqusConstants import *import numpy as npfrom scipy.special import gammafrom scipy.interpolate import RectBivariateSplinedef generate_bottom_up_hex_mesh(Lx, Ly, sn_x, sn_y, WB_a, WB_b, N, z_offset, num_layers):    """    核心计算函数：生成三维网格的节点坐标和单元拓扑关系    参数说明：    Lx, Ly: 模型在 X 和 Y 方向的长度    sn_x, sn_y: 原始随机采样的点数    WB_a: Weibull 分布的尺度参数 (Scale parameter)，影响起伏高度    WB_b: Weibull 分布的形状参数 (Shape parameter)，影响分布形态    N: 插值细分倍数，用于控制最终网格密度    z_offset: 基础厚度，确保粗糙面下方有实心基座    num_layers: Z 方向（厚度方向）的等分数/层数    """    # 1. 随机高度场生成    # U: 0-1 之间的均匀随机矩阵    U = np.random.rand(sn_y, sn_x)    # 计算均值偏置：Weibull 分布的期望值为 $a \cdot \Gamma(1 + 1/b)$    mean_offset = -WB_a * gamma(1 + 1 / WB_b)    # 利用逆变换采样生成符合 Weibull 分布的高度值 $Z = a \cdot [-\ln(1-U)]^{1/b}$    Z_raw = mean_offset + WB_a * (-np.log(1 - U)) ** (1 / WB_b)    # 2. 边界处理与插值    # 在四周填充一圈 0，确保模型边缘对齐（类似于周期性边界的简化处理）    Z_padded = np.pad(Z_raw, pad_width=1, mode='constant', constant_values=0)    # 定义插值的原始坐标点    x_p = np.linspace(0, Lx, sn_x + 2)    y_p = np.linspace(0, Ly, sn_y + 2)    # 使用双三次样条插值 (Bivariate Spline) 平滑表面    interp_func = RectBivariateSpline(x_p, y_p, Z_padded.T)    # 3. 生成最终网格坐标    nx, ny = (sn_x + 2) * N, (sn_y + 2) * N    x_new = np.linspace(0, Lx, nx)    y_new = np.linspace(0, Ly, ny)    # 计算顶层粗糙面的实际 Z 坐标    Z_top = interp_func(x_new, y_new).T + z_offset    # 4. 节点生成逻辑 (自底向上)    nodes = []    node_id_map = np.zeros((num_layers + 1, ny, nx), dtype=int)    node_count = 1    for k in range(num_layers + 1):        fraction = k / float(num_layers) # 当前层的高度比例        for j in range(ny):            for i in range(nx):                # Z 坐标从 0 按比例增长到顶层的随机高度                z = Z_top[j, i] * fraction                nodes.append([node_count, x_new[i], y_new[j], z])                # 记录节点 ID，方便后续构建单元                node_id_map[k, j, i] = node_count                node_count += 1    nodes = np.array(nodes)    # 5. 单元生成逻辑 (8节点六面体 C3D8)    elements = []    elem_count = 1    for k in range(num_layers): # 遍历每一层        for j in range(ny - 1): # 遍历 Y 方向            for i in range(nx - 1): # 遍历 X 方向                # 按照 Abaqus 六面体单元的节点编号顺序定义连接关系 (底面4点+顶面4点)                n1 = node_id_map[k, j, i]                n2 = node_id_map[k, j, i + 1]                n3 = node_id_map[k, j + 1, i + 1]                n4 = node_id_map[k, j + 1, i]                n5 = node_id_map[k + 1, j, i]                n6 = node_id_map[k + 1, j, i + 1]                n7 = node_id_map[k + 1, j + 1, i + 1]                n8 = node_id_map[k + 1, j + 1, i]                elements.append([elem_count, n1, n2, n3, n4, n5, n6, n7, n8])                elem_count += 1    elements = np.array(elements)    return nodes, elementsdef create_abaqus_model(nodes, elements, model_name='RoughSurfaceModel'):    """    Abaqus API 交互函数：将计算数据转化为 CAE 模型    """    # 清理已存在的同名模型    if model_name in mdb.models.keys():        del mdb.models[model_name]    myModel = mdb.Model(name=model_name)    # 创建 Part (3D, 可变形体)    myPart = myModel.Part(name='HexMeshPart', dimensionality=THREE_D, type=DEFORMABLE_BODY)    print('Adding nodes...')    # 逐一添加节点，必须指定 label 以保持与计算出的 ID 一致    for node in nodes:        node_id = int(node[0])        coords = (float(node[1]), float(node[2]), float(node[3]))        myPart.Node(coordinates=coords, label=node_id)    print('Adding elements...')    # 关键步骤：Abaqus 创建单元需要 MeshNode 对象，而不是整数 ID    # 这里通过字典建立 Label -> Node 对象的快速映射    node_dict = {}    for n in myPart.nodes:        node_dict[n.label] = n    # 逐一创建 HEX8 单元    for elem in elements:        elem_id = int(elem[0])        node_ids = [int(n) for n in elem[1:]]        # 将 ID 列表转换为 Node 对象元组        node_objects = tuple([node_dict[nid] for nid in node_ids])        myPart.Element(nodes=node_objects, elemShape=HEX8, label=elem_id)    # 设置单元集合 (Set)    all_elements = myPart.elements    myPart.Set(elements=all_elements, name='AllElements')    # 定义材料属性 (钢材示例)    myMaterial = myModel.Material(name='Steel')    myMaterial.Elastic(table=((210000.0, 0.3),)) # 杨氏模量, 泊松比    # 定义并指派截面属性 (Section)    mySection = myModel.HomogeneousSolidSection(        name='SolidSection',        material='Steel',        thickness=None    )    myPart.SectionAssignment(        region=myPart.sets['AllElements'],        sectionName='SolidSection'    )    # 创建装配体 (Assembly)    myAssembly = myModel.rootAssembly    myAssembly.DatumCsysByDefault(CARTESIAN)    myInstance = myAssembly.Instance(name='HexMesh-1', part=myPart, dependent=ON)    # 打印成功信息    print('=' * 60)    print('Mesh created successfully!')    print('Model name: %s' % model_name)    print('Total nodes: %d' % len(nodes))    print('Total elements: %d' % len(elements))    print('=' * 60)    return myModel, myPartif __name__ == '__main__':    # 建模参数配置    params = {        'Lx': 10.0,          # X 长度        'Ly': 10.0,          # Y 长度        'sn_x': 20,          # X 方向初始随机点数        'sn_y': 20,          # Y 方向初始随机点数        'WB_a': 0.3,         # Weibull 幅值系数        'WB_b': 2.0,         # Weibull 形状参数        'N': 2,              # 网格细分倍数 (nx = (sn_x+2)*N)        'z_offset': 5.0,     # 基座高度        'num_layers': 10     # 厚度方向分10层    }    # 执行生成流程    nodes, elements = generate_bottom_up_hex_mesh(**params)    model, part = create_abaqus_model(nodes, elements)