案例三：StaticCoupled 静气动弹性入门¶
这个 Notebook 用一个小型 NACA4415 悬臂翼跑 StaticCoupled。结构梁沿展向 $+y$，来流沿 $+x$。我们扫三个几何迎角，看总升力和变形后的气动网格。
静气动弹性麻烦在互相依赖：气动力要看结构变形，结构变形又由气动力决定。稳态下可以写成
$$ \mathbf{R}_s(\mathbf{q},\mathbf{f}_a)=\mathbf{0},\qquad \mathbf{f}_a=\mathcal{A}(\mathbf{q},\alpha,U_\infty,\rho), $$
$\mathbf{q}$ 是结构位形，$\mathbf{f}_a$ 是气动网格算出并映射到结构节点的气动力，$\alpha$ 是几何迎角，$U_\infty$ 是来流速度，$\rho$ 是空气密度。StaticCoupled 做的事很直接：结构求解器和 StaticUvlm 来回交换位形和载荷，直到残差够小。
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 四位数翼型的中弧线。
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. 静气动弹性求解器设置¶
StaticCoupled 在每个载荷步里交替调用结构静力求解器和 StaticUvlm。下面先把耦合参数写成一个小函数，再写单个迎角算例，最后才做迎角扫描。这样读者能分清楚三层逻辑：
- 耦合器怎么配置；
- 单个迎角如何生成模型并运行；
- 多个迎角如何批量比较。
In [14]:
Copied!
def configure_static_coupled(sim):
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"] = 60
sim.solvers["StaticCoupled"]["n_load_steps"] = 1
sim.solvers["StaticCoupled"]["tolerance"] = 1e-5
sim.solvers["StaticCoupled"]["relaxation_factor"] = 0.2
sim.solvers["NonLinearStatic"]["gravity_on"] = False
return sim
def configure_static_coupled(sim):
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"] = 60
sim.solvers["StaticCoupled"]["n_load_steps"] = 1
sim.solvers["StaticCoupled"]["tolerance"] = 1e-5
sim.solvers["StaticCoupled"]["relaxation_factor"] = 0.2
sim.solvers["NonLinearStatic"]["gravity_on"] = False
return sim
1.1 单个迎角算例¶
这个函数只负责一个迎角。twist_deg=alpha_deg 相当于把几何迎角写进气动面；flow 中只有 BeamLoader、AerogridLoader 和 StaticCoupled，说明这里求的是稳态气动弹性平衡。
In [15]:
Copied!
def run_static_aero(alpha_deg):
case_name = f"case03_static_alpha_{alpha_deg:+.0f}".replace("+", "p").replace("-", "m")
wing = make_wing(case_name, twist_deg=alpha_deg, EIy=4.0e4, EIz=5.0e6, num_chord_panels=3)
sim, route = solver_file(
case_name,
["BeamLoader", "AerogridLoader", "StaticCoupled"],
u_inf=15.0,
rho=1.225,
dt=0.03,
mstar=8,
)
configure_static_coupled(sim)
data, log = run_case(case_name, wing, sim, route)
force = summed_structural_force(data)
return data, force, log
def run_static_aero(alpha_deg):
case_name = f"case03_static_alpha_{alpha_deg:+.0f}".replace("+", "p").replace("-", "m")
wing = make_wing(case_name, twist_deg=alpha_deg, EIy=4.0e4, EIz=5.0e6, num_chord_panels=3)
sim, route = solver_file(
case_name,
["BeamLoader", "AerogridLoader", "StaticCoupled"],
u_inf=15.0,
rho=1.225,
dt=0.03,
mstar=8,
)
configure_static_coupled(sim)
data, log = run_case(case_name, wing, sim, route)
force = summed_structural_force(data)
return data, force, log
1.2 扫描三个迎角¶
这里故意只取 $-5^\circ$、$0^\circ$、$5^\circ$ 三个点。点数少，页面跑得快；趋势也足够说明迎角变化会同时改变气动力和结构变形。
In [16]:
Copied!
alphas = [-5, 0, 5]
aero_results = {}
for alpha in alphas:
data, force, log = run_static_aero(alpha)
aero_results[alpha] = {"data": data, "force": force, "log": log}
tip = structural_tip(data)
print(f"alpha={alpha:+} deg | sum force z={force[2]:.3f} N | tip z={tip[2]*1000:.3f} mm | log lines={len(log.splitlines())}")
alphas = [-5, 0, 5]
aero_results = {}
for alpha in alphas:
data, force, log = run_static_aero(alpha)
aero_results[alpha] = {"data": data, "force": force, "log": log}
tip = structural_tip(data)
print(f"alpha={alpha:+} deg | sum force z={force[2]:.3f} N | tip z={tip[2]*1000:.3f} mm | log lines={len(log.splitlines())}")
alpha=-5 deg | sum force z=65.497 N | tip z=-1.682 mm | log lines=0 alpha=+0 deg | sum force z=24.895 N | tip z=-0.641 mm | log lines=0 alpha=+5 deg | sum force z=-15.895 N | tip z=0.405 mm | log lines=0
2. 升力随迎角变化¶
这里把映射到结构节点上的稳态气动力求和，当作总升力的近似读数。小网格不能拿来评价气动精度，但可以看出迎角增大时升力上升的趋势。
从量纲上，升力常写成
$$ L=\frac{1}{2}\rho U_\infty^2 S C_L, $$
$S$ 是参考面积，$C_L$ 是升力系数。小迎角下常用的直觉是 $C_L\propto\alpha$。本页不严格换算 $C_L$，直接看结构节点上的 $F_z$ 总和，这样更容易和 SHARPy 的数据结构对上。
In [ ]:
Copied!
lift = np.array([aero_results[a]["force"][2] for a in alphas])
tip_z = np.array([structural_tip(aero_results[a]["data"])[2] * 1000 for a in alphas])
fig, axes = plt.subplots(1, 2, figsize=(7.4, 3.3))
axes[0].plot(alphas, lift, "o-", color=NATURE_COLORS["blue"])
label_panel(axes[0], "a")
style_axes(
axes[0],
xlabel="angle of attack alpha [deg]",
ylabel="total nodal $F_z$ [N]",
title="Total Lift Trend",
)
direct_label(axes[0], alphas[-1], lift[-1], f"{lift[-1]:.1f} N", NATURE_COLORS["blue"])
axes[1].plot(alphas, tip_z, "s-", color=NATURE_COLORS["orange"])
label_panel(axes[1], "b")
style_axes(
axes[1],
xlabel="angle of attack alpha [deg]",
ylabel="tip $z$ displacement [mm]",
title="Static Tip Deflection",
zero_line=True,
)
direct_label(axes[1], alphas[-1], tip_z[-1], f"{tip_z[-1]:.2f} mm", NATURE_COLORS["orange"])
finish_figure(fig)
#TODO 方向有问题，要 debug 下
lift = np.array([aero_results[a]["force"][2] for a in alphas])
tip_z = np.array([structural_tip(aero_results[a]["data"])[2] * 1000 for a in alphas])
fig, axes = plt.subplots(1, 2, figsize=(7.4, 3.3))
axes[0].plot(alphas, lift, "o-", color=NATURE_COLORS["blue"])
label_panel(axes[0], "a")
style_axes(
axes[0],
xlabel="angle of attack alpha [deg]",
ylabel="total nodal $F_z$ [N]",
title="Total Lift Trend",
)
direct_label(axes[0], alphas[-1], lift[-1], f"{lift[-1]:.1f} N", NATURE_COLORS["blue"])
axes[1].plot(alphas, tip_z, "s-", color=NATURE_COLORS["orange"])
label_panel(axes[1], "b")
style_axes(
axes[1],
xlabel="angle of attack alpha [deg]",
ylabel="tip $z$ displacement [mm]",
title="Static Tip Deflection",
zero_line=True,
)
direct_label(axes[1], alphas[-1], tip_z[-1], f"{tip_z[-1]:.2f} mm", NATURE_COLORS["orange"])
finish_figure(fig)
#TODO 方向有问题，要 debug 下
3. 查看 $\alpha=5^\circ$ 的三维网格¶
这张交互图同时显示结构梁、气动面和尾迹。读图时先看方向：来流和尾迹沿 $+x$，翼展沿 $+y$，升力对应 $z$ 方向。
静气动弹性结果不能只看升力曲线。气动力是在变形后的气动面上重新算出来的，所以变形网格本身也是结果。如果尾迹方向、翼面法向或结构变形符号不对，升力趋势再“合理”也不说明算例设对了。
In [18]:
Copied!
show_model(aero_results[5]["data"], "StaticCoupled: alpha = 5 deg", minus_mstar=4)
show_model(aero_results[5]["data"], "StaticCoupled: alpha = 5 deg", minus_mstar=4)
4. 沿展向的节点载荷分布¶
StaticCoupled 会把气动力插值到结构节点。下面画每个结构节点上的 $F_z$，用来解释根部弯矩为什么通常更大。
若把离散节点力近似看成连续分布 $q_z(y)$，截面处的弯矩贡献可写成
$$ M_y(y)\approx \int_y^L q_z(\eta)(\eta-y)\,\mathrm{d}\eta. $$
外侧载荷会带着力臂一起累积到根部，所以根部附近往往最吃力。这个图就是把气动力分布和结构弯曲响应接起来看。
In [19]:
Copied!
data = aero_results[5]["data"]
pos = node_positions(data)
loads = data.structure.timestep_info[-1].steady_applied_forces
fig, ax = plt.subplots(figsize=(5.6, 3.4))
color = NATURE_COLORS["teal"]
ax.plot(pos[:, 1], loads[:, 2], "o-", color=color)
style_axes(
ax,
xlabel="spanwise coordinate y [m]",
ylabel="nodal $F_z$ [N]",
title="Nodal Aerodynamic Force at alpha = 5 deg",
zero_line=True,
)
direct_label(ax, pos[np.argmax(loads[:, 2]), 1], np.max(loads[:, 2]), f"max {np.max(loads[:, 2]):.2f} N", color, dy=6)
finish_figure(fig)
data = aero_results[5]["data"]
pos = node_positions(data)
loads = data.structure.timestep_info[-1].steady_applied_forces
fig, ax = plt.subplots(figsize=(5.6, 3.4))
color = NATURE_COLORS["teal"]
ax.plot(pos[:, 1], loads[:, 2], "o-", color=color)
style_axes(
ax,
xlabel="spanwise coordinate y [m]",
ylabel="nodal $F_z$ [N]",
title="Nodal Aerodynamic Force at alpha = 5 deg",
zero_line=True,
)
direct_label(ax, pos[np.argmax(loads[:, 2]), 1], np.max(loads[:, 2]), f"max {np.max(loads[:, 2]):.2f} N", color, dy=6)
finish_figure(fig)
小结¶
静气动弹性算例先看三件事：坐标方向、耦合顺序、载荷符号。结构展向、来流方向、气动面法向对了，再去读升力、节点载荷和翼尖变形。这个 Notebook 的网格很小，适合教学和调试，不适合直接当高精度气动数据。
要把它推进到研究用例，下一步通常是网格收敛、参考面积和力系数定义、迎角范围检查，以及和解析结果或实验数据做对照。