知识库
借助API进行结构优化

本文将为您介绍如何借助德儒巴软件第二代API实现结构优化。

在结构工程领域，寻找兼顾安全与经济的最优结构是一项极具挑战的工作。传统基于软件图形界面的人工反复调整、试错的方法相对耗时，并且容易陷入局部最优。

借助参数化建模和API，可以更高效的进行这项工作。

示例一压力容器壁厚优化

示例一为矩形压力容器，顶部设置了圆锥形开口。底部为压力容器的转角边设置了固定线支座。该模型共考虑两个荷载工况：自重和活荷载。结构优化的目的是在满足等效米塞斯应力不超过屈服强度的前提下，获得最轻的结构重量。

我们将面单元的厚度设置为全局参数控制。实现了参数化建模以后，我们可以借助API、编写python代码对面单元厚度进行优化。

示例一的python代码分为以下几个部分：

导入外部库并设置基本变量

在代码的第一部分，我们需要导入所需程序库，然后定义各类变量。此处设定屈服强度用于可视化展示，并需要给出板材厚度的待调整取值范围。


          
          
          # --------------------- Imports ---------------------
import itertools
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from dlubal.api import rfem


# --------------------- Configuration ---------------------
parameter_ranges = {
    't_1': {'min': 8, 'max': 12, 'step': 2},
    't_2': {'min': 4, 'max': 8, 'step': 2},
} # thicknesses in mm
f_y = 235  # yield strength in MPa

定义主要功能函数

完成了外部函数库的导入、基本变量的定义以后，我们可以定义主要功能对应的函数。首先是build_parameter_grid函数。该函数的功能是遍历输入的param_ranges字典，根据之前定义的变量最小值、最大值、步长生成所有可能的组合，并封装为结构化pandas.DataFrame输出，以便最终结果的导出、数据分析和图表绘制。


          
          
          def build_parameter_grid(param_ranges):
    """Build a complete mutation dataset from min/max/step values in mm."""
    columns = []
    value_lists = []
    for name, bounds in param_ranges.items():
        min_val = bounds['min']
        max_val = bounds['max']
        step = bounds['step']
        if step <= 0:
            raise ValueError(f"Step width for '{name}' must be positive.")
        if max_val < min_val:
            raise ValueError(f"Max value for '{name}' must be >= min value.")
        values = np.arange(min_val, max_val + step * 0.5, step)
        columns.append(name)
        value_lists.append(values)
    combinations = list(itertools.product(*value_lists))
    return pd.DataFrame(combinations, columns=columns)

另外我们还可以定义一个set_glpa函数，用来建立python脚本和RFEM6程序的通信，动态覆写RFEM中的全局参数。


          
          
          def set_glpa(dl_app, p_name, p_value,):
    """Sets the global parameter p_name to the specified value p_value"""
    params = dl_app.get_object_list([rfem.global_parameters.GlobalParameter()])
    for p in params:
        if p.name == p_name:
            p.value = p_value
            dl_app.update_object(p)
            check=True
            return True
    if not check:
        raise ValueError(f"  Error: Parameter '{p_name}' not found!")

解空间的清洗

为了避免最终设计出的结构出现“头重脚轻”的不合理搭配，我们可以定义一些约束，提取剔除这部分“上部壁厚大于下部壁厚”不合理的解。

另外，我们还需要进行单位换算。Dlubal API 底层严格遵循国际单位制 (SI，即米、牛顿、帕斯卡)。图形界面中为了方便输入的毫米 (mm)，在传入 API 前必须强制归一化为米 (m) / 1000。


          
          
          # Build parameter set
para_dfo = build_parameter_grid(parameter_ranges)
para_dfo = para_dfo[para_dfo['t_1'] >= para_dfo['t_2']]
para_df = para_dfo / 1000  # mm -> m
print(f"Parameter set:\n {para_dfo}")

RFEM循环计算

此后我们可以将数据传输给RFEM，进行有限元分析，获取当前变量组合对应的最大等效米塞斯应力和重量。在使用python操纵RFEM时，建议注意以下内容：

1. 使用上下文管理with语法。这样可以确保无论计算是否中断或崩溃，通信隧道都会在退出时被安全释放，防止后台残留进程、RFEM6被始终锁定。
2. 在全局参数被覆写并改变模型几何拓扑后，旧的有限元网格可能并没有被自动删除，最好手动delete_mesh，保证每次求解时模型的干净。
3. 在进行计算时，建议设定skip_warnings为True，这样可以允许程序在遇到非致命警告（如某些局部网格质量不佳）时强行越过，保证了无人值守情况下的“无监测运行”不被弹窗阻塞。
4. 设置容错机制。如果某个变体因为参数极端导致网格畸变或计算不收敛，循环过程不应该直接退出，而是将结果标记为 NaN，保证整个遍历矩阵能够完整跑完。


          
          
          # Connect to RFEM with current model and run parameter study
with rfem.Application() as rf_app:
    # Check connection and print model info
    app_info = rf_app.get_application_info()
    print("Application Info:", app_info)
    res_paras=pd.DataFrame(columns=['sigvm','sfz'])
    for i in para_df.index:
            # set global parameters
            print(f"Set mutation {i} with parameters {[x for x in para_df.loc[i]]}")
            for p in para_df.columns:
                # print(f"Set mutation {i} with parameter {p} to {para_df.loc[i, p]}")
                if not set_glpa(rf_app, p, para_df.loc[i, p]):
                    break
            # Remesh
            rf_app.delete_mesh()
            rf_app.generate_mesh(skip_warnings=True)
            # Run calculation
            calculation = rf_app.calculate_all(skip_warnings=True)
            # Results
            if calculation.succeeded:
                # Von Mises stresses
                res_sigvm = rf_app.get_results(
                    results_type=rfem.results.STATIC_ANALYSIS_SURFACES_EQUIVALENT_STRESSES_MISES_MESH_NODES,
                    filters=[
                        rfem.results.ResultsFilter(column_id='loading', filter_expression='DS1'),
                        ]      
                    ).data
                sigvm = res_sigvm['sigma_eqv_mises'].max()/10**6 # N/m^2 -> MPa                
                # self-weight
                sfz = rf_app.get_result_table(
                    table=rfem.results.ResultTable.STATIC_ANALYSIS_SUMMARY_TABLE,
                    loading=rfem.ObjectId(no=1,object_type=rfem.OBJECT_TYPE_LOAD_CASE)
                ).data
                sfz = float(sfz.loc[6].value)  / 10 # N -> kg (g=10 m/s^2 for gravity acceleration)

                res_paras.loc[i] = pd.Series({'sigvm':sigvm,'sfz':sfz})
            else:
                print(f"Calculation failed for mutation {i} with parameters {[x for x in para_df.loc[i]]}")
                res_paras.loc[i] = pd.Series({'sigvm':np.nan,'sfz':np.nan})

整理与输出结果

计算完成后，我们可以将包含所有参数组合（t1, t2）、质量（sfz）和最大应力（sigvm）的Pandas DataFrame写入 Excel 文件。确保即使脚本关闭，数据也可被二次读取或审查。

另外，我们还可以借助python代码，绘制帕累托最优解的决策图。寻找红线下方、最靠近图表左侧的那个带有标注的散点，该点对应的数据（如）即为当前参数空间内满足结构安全的最轻设计方案。

示例二网格无关性验证

示例二是一个网格无关性、收敛性验证研究。核心逻辑是通过不断加密有限元网格，观察结构临界荷载系数的变化趋势，从而在“计算精度”与“计算成本”之间寻找最优平衡点。

在示例二中，我们将网格尺寸设置为全局参数控制

示例二的python代码分为以下几个部分：

导入外部库并设置基本变量

在代码的第一部分，我们需要导入所需程序库，然后定义实验的边界条件与评价标尺。我们设置了一个单调递减数组，来作为有限元网格尺寸的可选值；此外，我们还在这里定义了网格收敛的准则，当相邻两次网格加密结果差异小于1%时，认为网格收敛。


          
          
          # --------------------- Imports ---------------------
# Standard library
import time

# Third-party
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from dlubal.api import rfem

# ------------------- Configuration -------------------
FE_SIZES_MM = [15, 12, 10, 9, 8, 7, 6, 5, 4, 3] # List of finite element sizes to sweep in millimetres
PARAM_NAME = "l_fe"                             # Name of the global parameter in the model
CONVERGENCE_THRESHOLD_PCT = 1.0                 # Convergence threshold in percentage
FCR_ANA = 1065                                  # Analytical critical load factor (for reference in the plot)

循环计算

这部分代码实现了参数注入、网格重生成、计算和提取数据的功能。在网格重生成后，我们可以通过get_mesh_statistics()获取有限元网格总数。之后通过time.perf_counter获取每次计算的耗时c_time。这两个指标将用来体现计算成本。

计算时，仍然将skip_warnings设置为True，保证计算可以正常运行。

每次计算完成后，通过结果表筛选器获取LC1的临界荷载系数，并与上一次循环进行对比。


          
          
          # ------------------- Functions -------------------
def set_glpa(dl_app, p_name, p_value,):
    """Sets the global parameter p_name to the specified value p_value"""
    params = dl_app.get_object_list([rfem.global_parameters.GlobalParameter()])
    for p in params:
        if p.name == p_name:
            p.value = p_value
            dl_app.update_object(p)
            check=True
            return True
    if not check:
        raise ValueError(f"  Error: Parameter '{p_name}' not found!")
    
# ------------------- Main execution -------------------
with rfem.Application() as rf_app:
    # Check connection and print model info
    app_info = rf_app.get_application_info()
    print("Application Info:", app_info)

    results = pd.DataFrame(columns=['l_fe_mm', 'n_elements', 'c_time_s', 'f_cr', 'delta_pct'])
    f_prev = None
    i = 0
    for l_fe in FE_SIZES_MM:
        i += 1
        l_fe_m = l_fe / 1000.0  # Convert mm to m for use in RFEM
        if not set_glpa(dl_app=rf_app, p_name=PARAM_NAME, p_value=l_fe_m):
            break
        # Remesh
        rf_app.delete_mesh()
        rf_app.generate_mesh(skip_warnings=True)
        n_elem = rf_app.get_mesh_statistics().surface_2D_finite_elements
        # Run calculation
        t0 = time.perf_counter()
        calculation = rf_app.calculate_all(skip_warnings=True)
        t1 = time.perf_counter()
        c_time = t1 - t0
        # Results
        if calculation.succeeded:
            f_cr = rf_app.get_results(
            results_type=rfem.results.STABILITY_ANALYSIS_CRITICAL_LOAD_FACTORS,
            filters=[
            rfem.results.ResultsFilter(column_id='loading', filter_expression='LC1'),
            ]      
            ).data.loc[0].f

            if f_prev is not None and f_prev != 0 and f_cr is not None:
                delta_pct = abs(f_cr - f_prev) / abs(f_prev) * 100.0
                delta_str = f"{delta_pct:>4.2f}"
            else:
                delta_pct = None
                delta_str = f"{'---':>4}"

            if f_cr is not None:
                print(f"Mutation {i}: l_fe = {l_fe:.1f} mm | n_elements = {n_elem} | f_cr = {f_cr:.2f} | delta = {delta_str} %")
            results.loc[i] = pd.Series({
                'l_fe_mm': l_fe, 'n_elements': n_elem, 'c_time_s': c_time, 
                'f_cr': f_cr, 'delta_pct': delta_pct
                })
            f_prev = f_cr
        else:
            print(f"Calculation failed for mutation {i} with mesh size {l_fe:.1f} mm")
            results.loc[i] = pd.Series({
                'l_fe_mm': l_fe, 'n_elements': n_elem, 'c_time_s': np.nan, 
                'f_cr': np.nan, 'delta_pct': np.nan
                })

# Calculate 1 / n_elements for the second plot
results['inv_n_elements'] = 1.0 / results['n_elements']
# Calculate relative deviation from minimum critical load factor
results['rel_dev_f_cr_min'] = (results['f_cr'] - results['f_cr'].min()) / results['f_cr'].min() * 100

数据可视化

最后，我们可以将网格数量、网格尺寸、相对误差、计算时间等用图表绘制出来。


          
          
          # ---------------------- Export Results ---------------------
results.to_excel('./ParaS-LBA_MeshConvergence-min_results.xlsx')

# ------------------- Plotting -------------------
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot 1: f_cr vs. n_elements
ax1 = axes[0, 0]
ax1.plot(results['n_elements'], results['f_cr'], 'o-', color='tab:blue', markersize=6)
ax1.axhline(y=FCR_ANA, color='red', linestyle='--', linewidth=1, label=f'$F_{{cr,ana}}$ = {FCR_ANA} kN')
ax1.set_xlabel('Number of Elements')
ax1.set_ylabel('Critical Load Factor [kN]')
ax1.set_title('Critical Load Factor vs. Number of Elements')
ax1.grid(True, alpha=0.3)
ax1.legend()

# Plot 2: f_cr vs. 1/n_elements
ax2 = axes[0, 1]
ax2.plot(results['inv_n_elements'], results['f_cr'], 'o-', color='tab:blue', markersize=6)
ax2.axhline(y=FCR_ANA, color='red', linestyle='--', linewidth=1, label=f'$F_{{cr,ana}}$ = {FCR_ANA} kN')
ax2.set_xlabel('1 / Number of Elements')
ax2.set_ylabel('Critical Load Factor [kN]')
ax2.set_title('Critical Load Factor vs. 1 / Number of Elements')
ax2.grid(True, alpha=0.3)
ax2.legend()

# Plot 3: Relative Deviation vs. Calculation Time
ax3 = axes[1, 0]
ax3.plot(results['c_time_s'], results['rel_dev_f_cr_min'], 'o-', color='tab:blue', markersize=6)
ax3.set_xlabel('Calculation Time [s]')
ax3.set_ylabel('Relative Deviation [%]')
ax3.set_title('Relative Deviation vs. Calculation Time')
ax3.grid(True, alpha=0.3)

# Plot 4: Relative Deviation vs. l_fe_mm
ax4 = axes[1, 1]
ax4.plot(results['l_fe_mm'], results['rel_dev_f_cr_min'], 'o-', color='tab:blue', markersize=6)
ax4.set_xlabel('FE Mesh Size [mm]')
ax4.set_ylabel('Relative Deviation [%]')
ax4.set_title('Relative Deviation vs. FE Mesh Size')
ax4.grid(True, alpha=0.3)

fig.suptitle('Mesh Convergence Analysis', fontsize=14, fontweight='bold')
fig.tight_layout()
fig.savefig('./ParaS-LBA_MeshConvergence-min_results.png', dpi=200)
plt.show()