案例四：DynamicCoupled 动态气动弹性入门¶

动态气动弹性关心时间历程。结构位移、速度、气动力和尾迹都会随时间推进。这个 Notebook 用一个很小的 DynamicCoupled 算例，看翼尖周期力激励下的响应；为了让页面能快点跑完，只取 32 个时间步。

动态问题比静态问题多了惯性、阻尼和时间离散。结构侧可以先写成

$$ \mathbf{M}\ddot{\mathbf{q}}+\mathbf{C}\dot{\mathbf{q}}+\mathbf{K}\mathbf{q} = \mathbf{f}_a(t)+\mathbf{f}_\mathrm{ext}(t), $$

$\mathbf{M}$、$\mathbf{C}$、$\mathbf{K}$ 分别是质量、阻尼和刚度，$\mathbf{f}_a(t)$ 是非定常气动力，$\mathbf{f}_\mathrm{ext}(t)$ 是外加载荷。DynamicCoupled 会在每个时间步里处理结构动力学、非定常 UVLM 和 FSI 子迭代。

0. 导入依赖¶

这一格只负责导入标准库、NumPy/Matplotlib、SHARPy 和 Notebook 显示工具。把导入单独放出来后，读者可以先确认环境是否能找到 sharpy，再继续往下看模型和求解器。

In [1]:

Copied!





import contextlib
import io
import json
import os
import shutil
import uuid
from pathlib import Path

os.environ.setdefault("MPLCONFIGDIR", str(Path("_sharpy_runs") / ".mplconfig"))

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from cycler import cycler

import sharpy.sharpy_main
import sharpy.utils.generate_cases as gc
from IPython.display import HTML, display
import contextlib
import io
import json
import os
import shutil
import uuid
from pathlib import Path

os.environ.setdefault("MPLCONFIGDIR", str(Path("_sharpy_runs") / ".mplconfig"))

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from cycler import cycler

import sharpy.sharpy_main
import sharpy.utils.generate_cases as gc
from IPython.display import HTML, display

0.1 图表样式¶

这里统一图表字体、字号、线宽和配色。后面的物理代码不需要再反复设置这些视觉参数；如果要改论文或讲义图的风格，也只需要改这一格。

In [2]:

Copied!





NATURE_COLORS = {
    "blue": "#3B6FB6",
    "teal": "#2A9D8F",
    "orange": "#D97904",
    "red": "#C44536",
    "purple": "#7B5EA7",
    "grey": "#6B7280",
    "light_grey": "#E5E7EB",
    "dark": "#111827",
}
NATURE_SEQUENCE = [
    NATURE_COLORS["blue"],
    NATURE_COLORS["orange"],
    NATURE_COLORS["red"],
    NATURE_COLORS["teal"],
    NATURE_COLORS["purple"],
]

mpl.rcParams.update({
    "figure.dpi": 150,
    "figure.facecolor": "white",
    "axes.facecolor": "white",
    "font.family": "serif",
    "font.serif": ["Times New Roman", "Times", "Nimbus Roman", "DejaVu Serif", "serif"],
    "mathtext.fontset": "stix",
    "font.size": 10,
    "axes.titlesize": 12,
    "axes.labelsize": 11,
    "xtick.labelsize": 10,
    "ytick.labelsize": 10,
    "legend.fontsize": 10,
    "axes.unicode_minus": False,
    "svg.fonttype": "none",
    "pdf.fonttype": 42,
    "axes.spines.top": False,
    "axes.spines.right": False,
    "axes.linewidth": 0.8,
    "xtick.direction": "out",
    "ytick.direction": "out",
    "xtick.major.width": 0.7,
    "ytick.major.width": 0.7,
    "lines.linewidth": 1.9,
    "lines.markersize": 4.8,
    "axes.prop_cycle": cycler(color=NATURE_SEQUENCE),
    "legend.frameon": False,
})
NATURE_COLORS = {
    "blue": "#3B6FB6",
    "teal": "#2A9D8F",
    "orange": "#D97904",
    "red": "#C44536",
    "purple": "#7B5EA7",
    "grey": "#6B7280",
    "light_grey": "#E5E7EB",
    "dark": "#111827",
}
NATURE_SEQUENCE = [
    NATURE_COLORS["blue"],
    NATURE_COLORS["orange"],
    NATURE_COLORS["red"],
    NATURE_COLORS["teal"],
    NATURE_COLORS["purple"],
]

mpl.rcParams.update({
    "figure.dpi": 150,
    "figure.facecolor": "white",
    "axes.facecolor": "white",
    "font.family": "serif",
    "font.serif": ["Times New Roman", "Times", "Nimbus Roman", "DejaVu Serif", "serif"],
    "mathtext.fontset": "stix",
    "font.size": 10,
    "axes.titlesize": 12,
    "axes.labelsize": 11,
    "xtick.labelsize": 10,
    "ytick.labelsize": 10,
    "legend.fontsize": 10,
    "axes.unicode_minus": False,
    "svg.fonttype": "none",
    "pdf.fonttype": 42,
    "axes.spines.top": False,
    "axes.spines.right": False,
    "axes.linewidth": 0.8,
    "xtick.direction": "out",
    "ytick.direction": "out",
    "xtick.major.width": 0.7,
    "ytick.major.width": 0.7,
    "lines.linewidth": 1.9,
    "lines.markersize": 4.8,
    "axes.prop_cycle": cycler(color=NATURE_SEQUENCE),
    "legend.frameon": False,
})

0.1.1 图表辅助函数¶

style_axes、label_panel 和 direct_label 只做重复排版工作。把它们和全局样式分开，后面读图代码时就不会被样式细节打断。

In [3]:

Copied!





def style_axes(ax, *, xlabel=None, ylabel=None, title=None, zero_line=False, grid_axis="y"):
    if xlabel:
        ax.set_xlabel(xlabel)
    if ylabel:
        ax.set_ylabel(ylabel)
    if title:
        ax.set_title(title, loc="left", pad=7)
    if zero_line:
        ax.axhline(0.0, color=NATURE_COLORS["grey"], linewidth=0.8, zorder=0)
    if grid_axis:
        ax.grid(True, axis=grid_axis, color=NATURE_COLORS["light_grey"], linewidth=0.6)
        ax.set_axisbelow(True)
    else:
        ax.grid(False)
    ax.tick_params(length=3)
    return ax

def label_panel(ax, label):
    ax.text(-0.14, 1.06, label, transform=ax.transAxes, fontsize=12, fontweight="bold", va="bottom", clip_on=False)

def direct_label(ax, x, y, text, color, *, dx=6, dy=0):
    ax.annotate(
        text,
        xy=(x, y),
        xytext=(dx, dy),
        textcoords="offset points",
        color=color,
        fontsize=10,
        va="center",
        clip_on=False,
    )

def finish_figure(fig):
    fig.tight_layout()
    plt.show()
def style_axes(ax, *, xlabel=None, ylabel=None, title=None, zero_line=False, grid_axis="y"):
    if xlabel:
        ax.set_xlabel(xlabel)
    if ylabel:
        ax.set_ylabel(ylabel)
    if title:
        ax.set_title(title, loc="left", pad=7)
    if zero_line:
        ax.axhline(0.0, color=NATURE_COLORS["grey"], linewidth=0.8, zorder=0)
    if grid_axis:
        ax.grid(True, axis=grid_axis, color=NATURE_COLORS["light_grey"], linewidth=0.6)
        ax.set_axisbelow(True)
    else:
        ax.grid(False)
    ax.tick_params(length=3)
    return ax

def label_panel(ax, label):
    ax.text(-0.14, 1.06, label, transform=ax.transAxes, fontsize=12, fontweight="bold", va="bottom", clip_on=False)

def direct_label(ax, x, y, text, color, *, dx=6, dy=0):
    ax.annotate(
        text,
        xy=(x, y),
        xytext=(dx, dy),
        textcoords="offset points",
        color=color,
        fontsize=10,
        va="center",
        clip_on=False,
    )

def finish_figure(fig):
    fig.tight_layout()
    plt.show()

0.2 运行目录与静默执行¶

SHARPy 会写 .fem.h5、.aero.h5、.sharpy 等案例文件。这里把它们统一放到 _sharpy_runs/，并屏蔽求解器内部的冗长终端输出，让 Notebook 页面只留下我们主动打印的关键信息。

In [4]:

Copied!





RUN_ROOT = Path("_sharpy_runs") / Path(__name__).stem
RUN_ROOT.mkdir(parents=True, exist_ok=True)

@contextlib.contextmanager
def silence_process_output():
    # Silence subprocess/C-level stdout and stderr during SHARPy internals.
    stdout_fd = os.dup(1)
    stderr_fd = os.dup(2)
    with open(os.devnull, "w") as devnull:
        try:
            os.dup2(devnull.fileno(), 1)
            os.dup2(devnull.fileno(), 2)
            yield
        finally:
            os.dup2(stdout_fd, 1)
            os.dup2(stderr_fd, 2)
            os.close(stdout_fd)
            os.close(stderr_fd)
RUN_ROOT = Path("_sharpy_runs") / Path(__name__).stem
RUN_ROOT.mkdir(parents=True, exist_ok=True)

@contextlib.contextmanager
def silence_process_output():
    # Silence subprocess/C-level stdout and stderr during SHARPy internals.
    stdout_fd = os.dup(1)
    stderr_fd = os.dup(2)
    with open(os.devnull, "w") as devnull:
        try:
            os.dup2(devnull.fileno(), 1)
            os.dup2(devnull.fileno(), 2)
            yield
        finally:
            os.dup2(stdout_fd, 1)
            os.dup2(stderr_fd, 2)
            os.close(stdout_fd)
            os.close(stderr_fd)

0.3 翼型中弧线¶

naca4_airfoil 只生成 NACA 四位数翼型的中弧线。这个教学算例不追求高精度翼型几何，重点是让 AerodynamicInformation 有一个清楚、可重复的二维截面定义。

In [5]:

Copied!





def naca4_airfoil(code="4415", n=80):
    # Return one SHARPy-compatible camber line array for a NACA 4-digit airfoil.
    m = int(code[0]) / 100
    p = int(code[1]) / 10
    x = np.linspace(0.0, 1.0, n)
    yc = np.zeros_like(x)
    if m > 0 and p > 0:
        left = x < p
        yc[left] = m / p**2 * (2 * p * x[left] - x[left] ** 2)
        yc[~left] = m / (1 - p) ** 2 * ((1 - 2 * p) + 2 * p * x[~left] - x[~left] ** 2)
    airfoil = np.zeros((1, n, 2))
    airfoil[0, :, 0] = x
    airfoil[0, :, 1] = yc
    return airfoil
def naca4_airfoil(code="4415", n=80):
    # Return one SHARPy-compatible camber line array for a NACA 4-digit airfoil.
    m = int(code[0]) / 100
    p = int(code[1]) / 10
    x = np.linspace(0.0, 1.0, n)
    yc = np.zeros_like(x)
    if m > 0 and p > 0:
        left = x < p
        yc[left] = m / p**2 * (2 * p * x[left] - x[left] ** 2)
        yc[~left] = m / (1 - p) ** 2 * ((1 - 2 * p) + 2 * p * x[~left] - x[~left] ** 2)
    airfoil = np.zeros((1, n, 2))
    airfoil[0, :, 0] = x
    airfoil[0, :, 1] = yc
    return airfoil

0.4 结构梁和气动网格¶

make_wing 同时生成结构梁和气动面：结构节点沿展向 $y$ 排列，气动网格通过 chord、elastic axis、twist 和 airfoil 附着到结构梁上。气动弹性案例需要这两部分同时存在，因为气动力和结构变形要互相传递。

In [6]:

Copied!





def make_wing(
    case_name,
    *,
    span=2.0,
    chord=0.4,
    num_node=9,
    num_chord_panels=3,
    mass_per_unit_length=0.75,
    EA=1.0e7,
    GAy=1.0e6,
    GAz=1.0e6,
    GJ=1.0e4,
    EIy=2.0e4,
    EIz=5.0e6,
    stiffness_modifier=None,
    twist_deg=0.0,
    airfoil_code="4415",
):
    node_r = np.zeros((num_node, 3))
    node_r[:, 1] = np.linspace(0.0, span, num_node)

    wing = gc.AeroelasticInformation()
    structure = wing.StructuralInformation
    structure.num_node = num_node
    structure.num_node_elem = 3
    structure.compute_basic_num_elem()
    structure.generate_uniform_beam(
        node_r,
        mass_per_unit_length,
        0.02,
        0.01,
        0.01,
        np.zeros(3),
        EA,
        GAy,
        GAz,
        GJ,
        EIy,
        EIz,
        num_node_elem=3,
        y_BFoR="x_AFoR",
    )
    structure.boundary_conditions[0] = 1
    structure.boundary_conditions[-1] = -1
    if stiffness_modifier is not None:
        stiffness_modifier(structure)

    wing.AerodynamicInformation.create_one_uniform_aerodynamics(
        structure,
        chord=chord,
        twist=np.deg2rad(twist_deg),
        sweep=0.0,
        num_chord_panels=num_chord_panels,
        m_distribution="uniform",
        elastic_axis=0.25,
        num_points_camber=80,
        airfoil=naca4_airfoil(airfoil_code),
    )
    return wing
def make_wing(
    case_name,
    *,
    span=2.0,
    chord=0.4,
    num_node=9,
    num_chord_panels=3,
    mass_per_unit_length=0.75,
    EA=1.0e7,
    GAy=1.0e6,
    GAz=1.0e6,
    GJ=1.0e4,
    EIy=2.0e4,
    EIz=5.0e6,
    stiffness_modifier=None,
    twist_deg=0.0,
    airfoil_code="4415",
):
    node_r = np.zeros((num_node, 3))
    node_r[:, 1] = np.linspace(0.0, span, num_node)

    wing = gc.AeroelasticInformation()
    structure = wing.StructuralInformation
    structure.num_node = num_node
    structure.num_node_elem = 3
    structure.compute_basic_num_elem()
    structure.generate_uniform_beam(
        node_r,
        mass_per_unit_length,
        0.02,
        0.01,
        0.01,
        np.zeros(3),
        EA,
        GAy,
        GAz,
        GJ,
        EIy,
        EIz,
        num_node_elem=3,
        y_BFoR="x_AFoR",
    )
    structure.boundary_conditions[0] = 1
    structure.boundary_conditions[-1] = -1
    if stiffness_modifier is not None:
        stiffness_modifier(structure)

    wing.AerodynamicInformation.create_one_uniform_aerodynamics(
        structure,
        chord=chord,
        twist=np.deg2rad(twist_deg),
        sweep=0.0,
        num_chord_panels=num_chord_panels,
        m_distribution="uniform",
        elastic_axis=0.25,
        num_points_camber=80,
        airfoil=naca4_airfoil(airfoil_code),
    )
    return wing

0.5 SHARPy 输入文件与结果提取¶

这一组函数负责三件事：写求解器配置、运行 SHARPy、从最后一个时间步提取结构节点位置和节点载荷。把这些包装起来后，后面的案例代码可以集中讨论“要跑什么物理问题”，而不是反复处理文件路径和日志。

In [7]:

Copied!





def solver_file(case_name, flow, *, route=None, u_inf=10.0, rho=1.225, dt=0.02, mstar=8):
    route = Path(route or RUN_ROOT / case_name)
    route.mkdir(parents=True, exist_ok=True)
    sim = gc.SimulationInformation()
    sim.set_default_values()
    sim.solvers["SHARPy"]["flow"] = flow
    sim.solvers["SHARPy"]["case"] = case_name
    sim.solvers["SHARPy"]["route"] = str(route) + "/"
    sim.solvers["SHARPy"]["write_screen"] = "off"
    sim.solvers["SHARPy"]["write_log"] = False
    for settings in sim.solvers.values():
        if isinstance(settings, dict) and "print_info" in settings:
            settings["print_info"] = False

    sim.solvers["BeamLoader"]["unsteady"] = "off"
    sim.solvers["AerogridLoader"]["unsteady"] = "off"
    sim.solvers["AerogridLoader"]["mstar"] = mstar
    sim.solvers["AerogridLoader"]["freestream_dir"] = np.array([1.0, 0.0, 0.0])
    sim.solvers["AerogridLoader"]["wake_shape_generator"] = "StraightWake"
    sim.solvers["AerogridLoader"]["wake_shape_generator_input"] = {
        "u_inf": u_inf,
        "u_inf_direction": np.array([1.0, 0.0, 0.0]),
        "dt": dt,
    }
    sim.solvers["StaticUvlm"]["rho"] = rho
    sim.solvers["StaticUvlm"]["velocity_field_generator"] = "SteadyVelocityField"
    sim.solvers["StaticUvlm"]["velocity_field_input"] = {
        "u_inf": u_inf,
        "u_inf_direction": np.array([1.0, 0.0, 0.0]),
    }
    sim.solvers["StepUvlm"]["rho"] = rho
    sim.solvers["StepUvlm"]["velocity_field_generator"] = "SteadyVelocityField"
    sim.solvers["StepUvlm"]["velocity_field_input"] = {
        "u_inf": u_inf,
        "u_inf_direction": np.array([1.0, 0.0, 0.0]),
    }
    sim.solvers["StepUvlm"]["dt"] = dt
    return sim, route
def solver_file(case_name, flow, *, route=None, u_inf=10.0, rho=1.225, dt=0.02, mstar=8):
    route = Path(route or RUN_ROOT / case_name)
    route.mkdir(parents=True, exist_ok=True)
    sim = gc.SimulationInformation()
    sim.set_default_values()
    sim.solvers["SHARPy"]["flow"] = flow
    sim.solvers["SHARPy"]["case"] = case_name
    sim.solvers["SHARPy"]["route"] = str(route) + "/"
    sim.solvers["SHARPy"]["write_screen"] = "off"
    sim.solvers["SHARPy"]["write_log"] = False
    for settings in sim.solvers.values():
        if isinstance(settings, dict) and "print_info" in settings:
            settings["print_info"] = False

    sim.solvers["BeamLoader"]["unsteady"] = "off"
    sim.solvers["AerogridLoader"]["unsteady"] = "off"
    sim.solvers["AerogridLoader"]["mstar"] = mstar
    sim.solvers["AerogridLoader"]["freestream_dir"] = np.array([1.0, 0.0, 0.0])
    sim.solvers["AerogridLoader"]["wake_shape_generator"] = "StraightWake"
    sim.solvers["AerogridLoader"]["wake_shape_generator_input"] = {
        "u_inf": u_inf,
        "u_inf_direction": np.array([1.0, 0.0, 0.0]),
        "dt": dt,
    }
    sim.solvers["StaticUvlm"]["rho"] = rho
    sim.solvers["StaticUvlm"]["velocity_field_generator"] = "SteadyVelocityField"
    sim.solvers["StaticUvlm"]["velocity_field_input"] = {
        "u_inf": u_inf,
        "u_inf_direction": np.array([1.0, 0.0, 0.0]),
    }
    sim.solvers["StepUvlm"]["rho"] = rho
    sim.solvers["StepUvlm"]["velocity_field_generator"] = "SteadyVelocityField"
    sim.solvers["StepUvlm"]["velocity_field_input"] = {
        "u_inf": u_inf,
        "u_inf_direction": np.array([1.0, 0.0, 0.0]),
    }
    sim.solvers["StepUvlm"]["dt"] = dt
    return sim, route

0.5.1 写文件并运行 SHARPy¶

run_case 是实际触发求解的地方。它清理旧目录、写入 h5 和 .sharpy 文件，然后调用 sharpy_main.main。

In [8]:

Copied!





def run_case(case_name, wing, sim, route):
    log = io.StringIO()
    with contextlib.redirect_stdout(log), contextlib.redirect_stderr(log), silence_process_output():
        if route.exists():
            shutil.rmtree(route)
        route.mkdir(parents=True, exist_ok=True)
        wing.generate_h5_files(str(route) + "/", case_name)
        sim.generate_solver_file()
        data = sharpy.sharpy_main.main(["", str(route / f"{case_name}.sharpy")])
    return data, log.getvalue()
def run_case(case_name, wing, sim, route):
    log = io.StringIO()
    with contextlib.redirect_stdout(log), contextlib.redirect_stderr(log), silence_process_output():
        if route.exists():
            shutil.rmtree(route)
        route.mkdir(parents=True, exist_ok=True)
        wing.generate_h5_files(str(route) + "/", case_name)
        sim.generate_solver_file()
        data = sharpy.sharpy_main.main(["", str(route / f"{case_name}.sharpy")])
    return data, log.getvalue()

0.5.2 结果提取¶

这些小函数只从 SHARPy 的 data 对象里提取常用量：翼尖位置、所有节点位置和节点载荷合力。

In [9]:

Copied!





def structural_tip(data):
    tstep = data.structure.timestep_info[-1]
    return tstep.pos[-1].copy()

def node_positions(data):
    tstep = data.structure.timestep_info[-1]
    return tstep.pos.copy()

def summed_structural_force(data):
    tstep = data.structure.timestep_info[-1]
    return np.sum(tstep.steady_applied_forces, axis=0)
def structural_tip(data):
    tstep = data.structure.timestep_info[-1]
    return tstep.pos[-1].copy()

def node_positions(data):
    tstep = data.structure.timestep_info[-1]
    return tstep.pos.copy()

def summed_structural_force(data):
    tstep = data.structure.timestep_info[-1]
    return np.sum(tstep.steady_applied_forces, axis=0)

0.6 三维结构/气动网格显示¶

show_model 用 Plotly 同时画结构梁、气动面板和 wake。它主要用于检查方向：来流、尾迹、展向、变形符号和气动面法向是否看起来合理。

0.6.1 结构梁轨迹¶

先画结构梁本身。每个结构单元用 lines+markers 表示，方便检查梁节点和单元连接顺序。

In [10]:

Copied!





def structural_beam_traces(data):
    if len(data.structure.timestep_info) == 0:
        struct_tstep = data.structure.ini_info
    else:
        struct_tstep = data.structure.timestep_info[-1]

    traces = []
    for ielem in range(data.structure.num_elem):
        nodes = data.structure.connectivities[ielem, :][[0, 2, 1]]
        xyz = struct_tstep.pos[nodes, :]
        traces.append({
            "type": "scatter3d",
            "mode": "lines+markers",
            "x": xyz[:, 0].tolist(),
            "y": xyz[:, 1].tolist(),
            "z": xyz[:, 2].tolist(),
            "line": {"color": "#1f77b4", "width": 6},
            "marker": {"size": 3, "color": "#1f77b4"},
            "name": "beam",
            "showlegend": ielem == 0,
        })
    return traces
def structural_beam_traces(data):
    if len(data.structure.timestep_info) == 0:
        struct_tstep = data.structure.ini_info
    else:
        struct_tstep = data.structure.timestep_info[-1]

    traces = []
    for ielem in range(data.structure.num_elem):
        nodes = data.structure.connectivities[ielem, :][[0, 2, 1]]
        xyz = struct_tstep.pos[nodes, :]
        traces.append({
            "type": "scatter3d",
            "mode": "lines+markers",
            "x": xyz[:, 0].tolist(),
            "y": xyz[:, 1].tolist(),
            "z": xyz[:, 2].tolist(),
            "line": {"color": "#1f77b4", "width": 6},
            "marker": {"size": 3, "color": "#1f77b4"},
            "name": "beam",
            "showlegend": ielem == 0,
        })
    return traces

0.6.2 气动面和 wake 轨迹¶

气动面板和 wake 是两个不同对象：气动面用于当前时刻的载荷计算，wake 表示尾迹历史。拆开后读者更容易看到 zeta 和 zeta_star 的区别。

In [11]:

Copied!





def wake_traces(aero_tstep, isurf, minus_mstar=6):
    traces = []
    zeta_star = aero_tstep.zeta_star[isurf]
    mstar, nstar = aero_tstep.dimensions_star[isurf]
    wake_stop = max(1, mstar + 1 - minus_mstar)

    for i_m in range(wake_stop):
        traces.append({
            "type": "scatter3d",
            "mode": "lines",
            "x": zeta_star[0, i_m, :].tolist(),
            "y": zeta_star[1, i_m, :].tolist(),
            "z": zeta_star[2, i_m, :].tolist(),
            "line": {"color": "#8ecae6", "width": 1},
            "name": "wake",
            "showlegend": isurf == 0 and i_m == 0,
        })
    for i_n in range(nstar + 1):
        traces.append({
            "type": "scatter3d",
            "mode": "lines",
            "x": zeta_star[0, :wake_stop, i_n].tolist(),
            "y": zeta_star[1, :wake_stop, i_n].tolist(),
            "z": zeta_star[2, :wake_stop, i_n].tolist(),
            "line": {"color": "#8ecae6", "width": 1},
            "name": "wake",
            "showlegend": False,
        })
    return traces
def wake_traces(aero_tstep, isurf, minus_mstar=6):
    traces = []
    zeta_star = aero_tstep.zeta_star[isurf]
    mstar, nstar = aero_tstep.dimensions_star[isurf]
    wake_stop = max(1, mstar + 1 - minus_mstar)

    for i_m in range(wake_stop):
        traces.append({
            "type": "scatter3d",
            "mode": "lines",
            "x": zeta_star[0, i_m, :].tolist(),
            "y": zeta_star[1, i_m, :].tolist(),
            "z": zeta_star[2, i_m, :].tolist(),
            "line": {"color": "#8ecae6", "width": 1},
            "name": "wake",
            "showlegend": isurf == 0 and i_m == 0,
        })
    for i_n in range(nstar + 1):
        traces.append({
            "type": "scatter3d",
            "mode": "lines",
            "x": zeta_star[0, :wake_stop, i_n].tolist(),
            "y": zeta_star[1, :wake_stop, i_n].tolist(),
            "z": zeta_star[2, :wake_stop, i_n].tolist(),
            "line": {"color": "#8ecae6", "width": 1},
            "name": "wake",
            "showlegend": False,
        })
    return traces

In [12]:

Copied!





def aero_grid_and_wake_traces(data, minus_mstar=6):
    if len(data.structure.timestep_info) == 0:
        aero_tstep = data.aero.ini_info
    else:
        aero_tstep = data.aero.timestep_info[-1]

    traces = []
    if aero_tstep is None:
        return traces

    for isurf in range(aero_tstep.n_surf):
        zeta = aero_tstep.zeta[isurf]
        m_dim, n_dim = aero_tstep.dimensions[isurf]
        for i_m in range(m_dim + 1):
            traces.append({
                "type": "scatter3d",
                "mode": "lines",
                "x": zeta[0, i_m, :].tolist(),
                "y": zeta[1, i_m, :].tolist(),
                "z": zeta[2, i_m, :].tolist(),
                "line": {"color": "#333333", "width": 2},
                "name": "aero grid",
                "showlegend": isurf == 0 and i_m == 0,
            })
        for i_n in range(n_dim + 1):
            traces.append({
                "type": "scatter3d",
                "mode": "lines",
                "x": zeta[0, :, i_n].tolist(),
                "y": zeta[1, :, i_n].tolist(),
                "z": zeta[2, :, i_n].tolist(),
                "line": {"color": "#333333", "width": 2},
                "name": "aero grid",
                "showlegend": False,
            })
        traces.extend(wake_traces(aero_tstep, isurf, minus_mstar))
    return traces
def aero_grid_and_wake_traces(data, minus_mstar=6):
    if len(data.structure.timestep_info) == 0:
        aero_tstep = data.aero.ini_info
    else:
        aero_tstep = data.aero.timestep_info[-1]

    traces = []
    if aero_tstep is None:
        return traces

    for isurf in range(aero_tstep.n_surf):
        zeta = aero_tstep.zeta[isurf]
        m_dim, n_dim = aero_tstep.dimensions[isurf]
        for i_m in range(m_dim + 1):
            traces.append({
                "type": "scatter3d",
                "mode": "lines",
                "x": zeta[0, i_m, :].tolist(),
                "y": zeta[1, i_m, :].tolist(),
                "z": zeta[2, i_m, :].tolist(),
                "line": {"color": "#333333", "width": 2},
                "name": "aero grid",
                "showlegend": isurf == 0 and i_m == 0,
            })
        for i_n in range(n_dim + 1):
            traces.append({
                "type": "scatter3d",
                "mode": "lines",
                "x": zeta[0, :, i_n].tolist(),
                "y": zeta[1, :, i_n].tolist(),
                "z": zeta[2, :, i_n].tolist(),
                "line": {"color": "#333333", "width": 2},
                "name": "aero grid",
                "showlegend": False,
            })
        traces.extend(wake_traces(aero_tstep, isurf, minus_mstar))
    return traces

0.6.3 组装 Plotly 图¶

最后这一格只负责把结构、气动面和 wake 的 traces 放进同一个 Plotly 3D 场景。

In [13]:

Copied!





def show_model(data, title, minus_mstar=6):
    traces = structural_beam_traces(data)
    traces.extend(aero_grid_and_wake_traces(data, minus_mstar=minus_mstar))

    div_id = "plotly-" + uuid.uuid4().hex
    layout = {
        "title": title,
        "height": 620,
        "margin": {"l": 0, "r": 0, "t": 42, "b": 0},
        "legend": {"x": 0.02, "y": 0.98},
        "scene": {
            "xaxis": {"title": "x [m]"},
            "yaxis": {"title": "y [m]"},
            "zaxis": {"title": "z [m]"},
            "aspectmode": "data",
        },
    }
    config = {"responsive": True, "displaylogo": False}
    html = f'''
    <div id="{div_id}" style="width: 100%; height: 620px;"></div>
    <script src="https://cdn.plot.ly/plotly-2.35.2.min.js"></script>
    <script>
    Plotly.newPlot(
      "{div_id}",
      {json.dumps(traces)},
      {json.dumps(layout)},
      {json.dumps(config)}
    );
    </script>
    '''
    display(HTML(html))
def show_model(data, title, minus_mstar=6):
    traces = structural_beam_traces(data)
    traces.extend(aero_grid_and_wake_traces(data, minus_mstar=minus_mstar))

    div_id = "plotly-" + uuid.uuid4().hex
    layout = {
        "title": title,
        "height": 620,
        "margin": {"l": 0, "r": 0, "t": 42, "b": 0},
        "legend": {"x": 0.02, "y": 0.98},
        "scene": {
            "xaxis": {"title": "x [m]"},
            "yaxis": {"title": "y [m]"},
            "zaxis": {"title": "z [m]"},
            "aspectmode": "data",
        },
    }
    config = {"responsive": True, "displaylogo": False}
    html = f'''
    <div id="{div_id}" style="width: 100%; height: 620px;"></div>
    <script src="https://cdn.plot.ly/plotly-2.35.2.min.js"></script>
    <script>
    Plotly.newPlot(
      "{div_id}",
      {json.dumps(traces)},
      {json.dumps(layout)},
      {json.dumps(config)}
    );
    </script>
    '''
    display(HTML(html))

1. 构造动态载荷¶

动态算例比静态算例多了时间轴。下面把时间离散、模型生成、动态力数组和输入载荷图拆开写，读者可以逐步确认：先有时间步，再有结构/气动模型，最后才有每个时间步上的外力。

In [14]:

Copied!





case_name = "case04_dynamic_tip_force"
frequency = 1.0
samples_per_period = 16
cycles = 2
dt = 1.0 / (frequency * samples_per_period)
n_time_steps = cycles * samples_per_period
time = np.arange(n_time_steps + 1) * dt

print(f"dt = {dt:.4f} s")
print(f"number of dynamic steps = {n_time_steps}")
case_name = "case04_dynamic_tip_force"
frequency = 1.0
samples_per_period = 16
cycles = 2
dt = 1.0 / (frequency * samples_per_period)
n_time_steps = cycles * samples_per_period
time = np.arange(n_time_steps + 1) * dt

print(f"dt = {dt:.4f} s")
print(f"number of dynamic steps = {n_time_steps}")

dt = 0.0625 s
number of dynamic steps = 32

1.1 生成结构/气动模型¶

这里仍然使用一个很小的悬臂翼模型。num_chord_panels=2 会让气动网格更粗，但能让教学页面运行更快。

In [15]:

Copied!





wing = make_wing(
    case_name,
    twist_deg=0.0,
    EIy=5.0e4,
    EIz=5.0e6,
    mass_per_unit_length=0.75,
    num_chord_panels=2,
)

print("number of structural nodes:", wing.StructuralInformation.num_node)
print("number of structural elements:", wing.StructuralInformation.num_elem)
wing = make_wing(
    case_name,
    twist_deg=0.0,
    EIy=5.0e4,
    EIz=5.0e6,
    mass_per_unit_length=0.75,
    num_chord_panels=2,
)

print("number of structural nodes:", wing.StructuralInformation.num_node)
print("number of structural elements:", wing.StructuralInformation.num_elem)

number of structural nodes: 9
number of structural elements: 4

1.2 写入周期端部力¶

dynamic_forces 的维度是 (time, node, 6)。本例只在最后一个结构节点的 $F_z$ 分量施加正弦力，其余节点和力矩分量保持为零。

In [16]:

Copied!

dynamic_forces = np.zeros((n_time_steps + 1, wing.StructuralInformation.num_node, 6))
dynamic_forces[:, -1, 2] = 8.0 * np.sin(2 * np.pi * frequency * time)

print("dynamic_forces shape:", dynamic_forces.shape)
print(f"tip-force amplitude: {np.max(np.abs(dynamic_forces[:, -1, 2])):.2f} N")
dynamic_forces = np.zeros((n_time_steps + 1, wing.StructuralInformation.num_node, 6))
dynamic_forces[:, -1, 2] = 8.0 * np.sin(2 * np.pi * frequency * time)

print("dynamic_forces shape:", dynamic_forces.shape)
print(f"tip-force amplitude: {np.max(np.abs(dynamic_forces[:, -1, 2])):.2f} N")

dynamic_forces shape: (33, 9, 6)
tip-force amplitude: 8.00 N

1.3 检查输入载荷¶

先画输入而不是直接跑求解器，是为了确认时间步、频率、幅值和符号都符合预期。

In [17]:

Copied!





fig, ax = plt.subplots(figsize=(5.8, 2.8))
ax.plot(time, dynamic_forces[:, -1, 2], "o-", color=NATURE_COLORS["blue"])
style_axes(
    ax,
    xlabel="time [s]",
    ylabel="tip $F_z$ [N]",
    title="Input Periodic Tip Force",
    zero_line=True,
)
finish_figure(fig)
fig, ax = plt.subplots(figsize=(5.8, 2.8))
ax.plot(time, dynamic_forces[:, -1, 2], "o-", color=NATURE_COLORS["blue"])
style_axes(
    ax,
    xlabel="time [s]",
    ylabel="tip $F_z$ [N]",
    title="Input Periodic Tip Force",
    zero_line=True,
)
finish_figure(fig)

No description has been provided for this image

2. 运行 DynamicCoupled¶

动态耦合前，先用 StaticCoupled 建一个初始平衡，再进入 DynamicCoupled 时间推进。这里把配置拆成四步：基础求解链、静态初始平衡、动态 FSI 设置、结构时间积分设置。

In [18]:

Copied!





sim, route = solver_file(
    case_name,
    ["BeamLoader", "AerogridLoader", "StaticCoupled", "DynamicCoupled"],
    u_inf=12.0,
    rho=1.225,
    dt=dt,
    mstar=6,
)
sim.solvers["BeamLoader"]["unsteady"] = "on"
sim.solvers["AerogridLoader"]["unsteady"] = "on"
sim.with_dynamic_forces = True
sim.dynamic_forces = dynamic_forces

print("solver flow:", sim.solvers["SHARPy"]["flow"])
sim, route = solver_file(
    case_name,
    ["BeamLoader", "AerogridLoader", "StaticCoupled", "DynamicCoupled"],
    u_inf=12.0,
    rho=1.225,
    dt=dt,
    mstar=6,
)
sim.solvers["BeamLoader"]["unsteady"] = "on"
sim.solvers["AerogridLoader"]["unsteady"] = "on"
sim.with_dynamic_forces = True
sim.dynamic_forces = dynamic_forces

print("solver flow:", sim.solvers["SHARPy"]["flow"])

solver flow: ['BeamLoader', 'AerogridLoader', 'StaticCoupled', 'DynamicCoupled']

2.1 静态初始平衡¶

这一步给动态计算一个稳态起点。结构静力求解器和 StaticUvlm 先做松弛迭代，避免动态第一步从一个不平衡状态突然开始。

In [19]:

Copied!





sim.solvers["StaticCoupled"]["structural_solver"] = "NonLinearStatic"
sim.solvers["StaticCoupled"]["structural_solver_settings"] = sim.solvers["NonLinearStatic"]
sim.solvers["StaticCoupled"]["aero_solver"] = "StaticUvlm"
sim.solvers["StaticCoupled"]["aero_solver_settings"] = sim.solvers["StaticUvlm"]
sim.solvers["StaticCoupled"]["max_iter"] = 30
sim.solvers["StaticCoupled"]["n_load_steps"] = 1
sim.solvers["StaticCoupled"]["tolerance"] = 1e-4
sim.solvers["StaticCoupled"]["relaxation_factor"] = 0.2
sim.solvers["StaticCoupled"]["structural_solver"] = "NonLinearStatic"
sim.solvers["StaticCoupled"]["structural_solver_settings"] = sim.solvers["NonLinearStatic"]
sim.solvers["StaticCoupled"]["aero_solver"] = "StaticUvlm"
sim.solvers["StaticCoupled"]["aero_solver_settings"] = sim.solvers["StaticUvlm"]
sim.solvers["StaticCoupled"]["max_iter"] = 30
sim.solvers["StaticCoupled"]["n_load_steps"] = 1
sim.solvers["StaticCoupled"]["tolerance"] = 1e-4
sim.solvers["StaticCoupled"]["relaxation_factor"] = 0.2

2.2 动态气动-结构耦合¶

DynamicCoupled 每个时间步都会调用结构动力学求解器和非定常 UVLM。fsi_substeps 与 fsi_tolerance 控制每个时间步内部的流固耦合收敛。

In [20]:

Copied!





sim.solvers["DynamicCoupled"]["structural_solver"] = "NonLinearDynamicPrescribedStep"
sim.solvers["DynamicCoupled"]["structural_solver_settings"] = sim.solvers["NonLinearDynamicPrescribedStep"]
sim.solvers["DynamicCoupled"]["aero_solver"] = "StepUvlm"
sim.solvers["DynamicCoupled"]["aero_solver_settings"] = sim.solvers["StepUvlm"]
sim.solvers["DynamicCoupled"]["n_time_steps"] = n_time_steps
sim.solvers["DynamicCoupled"]["dt"] = dt
sim.solvers["DynamicCoupled"]["fsi_substeps"] = 20
sim.solvers["DynamicCoupled"]["fsi_tolerance"] = 1e-4
sim.solvers["DynamicCoupled"]["relaxation_factor"] = 0.2
sim.solvers["StepUvlm"]["n_time_steps"] = n_time_steps + 1
sim.solvers["DynamicCoupled"]["structural_solver"] = "NonLinearDynamicPrescribedStep"
sim.solvers["DynamicCoupled"]["structural_solver_settings"] = sim.solvers["NonLinearDynamicPrescribedStep"]
sim.solvers["DynamicCoupled"]["aero_solver"] = "StepUvlm"
sim.solvers["DynamicCoupled"]["aero_solver_settings"] = sim.solvers["StepUvlm"]
sim.solvers["DynamicCoupled"]["n_time_steps"] = n_time_steps
sim.solvers["DynamicCoupled"]["dt"] = dt
sim.solvers["DynamicCoupled"]["fsi_substeps"] = 20
sim.solvers["DynamicCoupled"]["fsi_tolerance"] = 1e-4
sim.solvers["DynamicCoupled"]["relaxation_factor"] = 0.2
sim.solvers["StepUvlm"]["n_time_steps"] = n_time_steps + 1

2.3 结构动力学时间积分¶

结构求解器使用相同的 dt 和步数。这里关闭重力，只保留周期端部力和非定常气动力对翼尖响应的影响。

In [21]:

Copied!





sim.solvers["NonLinearDynamicPrescribedStep"]["dt"] = dt
sim.solvers["NonLinearDynamicPrescribedStep"]["num_steps"] = n_time_steps
sim.solvers["NonLinearDynamicPrescribedStep"]["gravity_on"] = False
sim.solvers["NonLinearDynamicPrescribedStep"]["newmark_damp"] = 1e-3
sim.solvers["NonLinearDynamicPrescribedStep"]["dt"] = dt
sim.solvers["NonLinearDynamicPrescribedStep"]["num_steps"] = n_time_steps
sim.solvers["NonLinearDynamicPrescribedStep"]["gravity_on"] = False
sim.solvers["NonLinearDynamicPrescribedStep"]["newmark_damp"] = 1e-3

2.4 写文件并运行¶

动态算例除了 .fem.h5、.aero.h5、.sharpy，还要写 .dyn.h5，因为时间相关外力保存在动态输入文件里。

In [22]:

Copied!





log = io.StringIO()
with contextlib.redirect_stdout(log), contextlib.redirect_stderr(log), silence_process_output():
    if route.exists():
        shutil.rmtree(route)
    route.mkdir(parents=True, exist_ok=True)
    wing.generate_h5_files(str(route) + "/", case_name)
    sim.generate_solver_file()
    sim.generate_dyn_file(n_time_steps)
    data = sharpy.sharpy_main.main(["", str(route / f"{case_name}.sharpy")])

print("saved time steps:", len(data.structure.timestep_info))
print("captured SHARPy log lines:", len(log.getvalue().splitlines()))
log = io.StringIO()
with contextlib.redirect_stdout(log), contextlib.redirect_stderr(log), silence_process_output():
    if route.exists():
        shutil.rmtree(route)
    route.mkdir(parents=True, exist_ok=True)
    wing.generate_h5_files(str(route) + "/", case_name)
    sim.generate_solver_file()
    sim.generate_dyn_file(n_time_steps)
    data = sharpy.sharpy_main.main(["", str(route / f"{case_name}.sharpy")])

print("saved time steps:", len(data.structure.timestep_info))
print("captured SHARPy log lines:", len(log.getvalue().splitlines()))

saved time steps: 33
captured SHARPy log lines: 0

3. 画出翼尖位移时间历程¶

动态响应先看时间序列：输入力长什么样，翼尖位移有没有同频变化，幅值有没有发散。

对简单线性单自由度受迫振动，响应常写成

$$ z(t)=Z\sin(2\pi f t-\phi), $$

$Z$ 是响应幅值，$\phi$ 是相位滞后。这个 Notebook 的模型当然不止一个自由度，但读图方法类似：先比周期，再看幅值和相位。如果位移一圈比一圈大，优先检查时间步、阻尼、FSI 收敛和载荷幅值。

In [23]:

Copied!





tip_z = np.array([step.pos[-1, 2] for step in data.structure.timestep_info])
t = np.arange(len(tip_z)) * dt

fig, axes = plt.subplots(2, 1, figsize=(6.4, 4.8), sharex=True)
axes[0].plot(time, dynamic_forces[:, -1, 2], "o-", color=NATURE_COLORS["blue"])
label_panel(axes[0], "a")
style_axes(axes[0], ylabel="input $F_z$ [N]", title="Tip Periodic Force", zero_line=True)

axes[1].plot(t, tip_z * 1000, "s-", color=NATURE_COLORS["red"])
label_panel(axes[1], "b")
style_axes(
    axes[1],
    xlabel="time [s]",
    ylabel="tip $z$ [mm]",
    title="Dynamic Response",
    zero_line=True,
)
finish_figure(fig)

print(f"tip displacement range: {tip_z.min()*1000:.3f} to {tip_z.max()*1000:.3f} mm")
tip_z = np.array([step.pos[-1, 2] for step in data.structure.timestep_info])
t = np.arange(len(tip_z)) * dt

fig, axes = plt.subplots(2, 1, figsize=(6.4, 4.8), sharex=True)
axes[0].plot(time, dynamic_forces[:, -1, 2], "o-", color=NATURE_COLORS["blue"])
label_panel(axes[0], "a")
style_axes(axes[0], ylabel="input $F_z$ [N]", title="Tip Periodic Force", zero_line=True)

axes[1].plot(t, tip_z * 1000, "s-", color=NATURE_COLORS["red"])
label_panel(axes[1], "b")
style_axes(
    axes[1],
    xlabel="time [s]",
    ylabel="tip $z$ [mm]",
    title="Dynamic Response",
    zero_line=True,
)
finish_figure(fig)

print(f"tip displacement range: {tip_z.min()*1000:.3f} to {tip_z.max()*1000:.3f} mm")

tip displacement range: -0.257 to -0.257 mm

4. 查看最后一个时间步的气动网格¶

动态耦合会把尾迹一步步往下游推进。最后一个时间步的交互图可以快速看两件事：尾迹是不是沿来流方向展开，结构变形有没有跑到不合理的范围。

非定常气动力不只看当前翼面姿态，也受尾迹历史影响。尾迹可以理解为前面时间步留下的涡量记录。若尾迹方向错了，或者很快折成奇怪的形状，后面解释气动力就没有太大意义。

In [24]:

Copied!

show_model(data, "DynamicCoupled: final time step", minus_mstar=3)
show_model(data, "DynamicCoupled: final time step", minus_mstar=3)

小结¶

动态耦合比静态耦合慢，因为每个时间步都要做气动、结构和 FSI 子迭代。教学时先用小网格、短时间序列确认流程；流程没问题后，再加周期数、采样率、迎角和载荷幅值。

这个例子只搭最小闭环：外载随时间变化，结构响应跟着推进，非定常气动力和尾迹历史进入下一步求解。真正做研究时，还要检查时间步收敛、初始瞬态、周期稳态是否建立，以及载荷频率和结构固有频率的关系。