案例二：三类悬臂梁刚度模型的静力响应对比¶

下载本案例 Jupyter Notebook

这一节继续使用悬臂梁。比较对象是三类结构刚度模型：Case 1 均匀刚度、Case 2 耦合刚度和 Case 3 沿展向变化的刚度。

0. 导入依赖¶

In [1]:

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]:

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,
})

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()

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, }) 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 运行目录与静默执行¶

In [3]:

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

@contextlib.contextmanager
def silence_process_output():
    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(): 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 纯结构 Beam 构造函数¶

In [4]:

def make_beam(
    case_name,
    *,
    span=2.0,
    num_node=17,
    mass_per_unit_length=0.75,
    EA=1.0e7,
    GAy=1.0e6,
    GAz=1.0e6,
    GJ=1.0e4,
    EIy=1.5e4,
    EIz=5.0e6,
    stiffness_modifier=None,
):
    node_r = np.zeros((num_node, 3))
    node_r[:, 1] = np.linspace(0.0, span, num_node)

    model = gc.AeroelasticInformation()
    structure = model.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)
    return model

def make_beam( case_name, *, span=2.0, num_node=17, mass_per_unit_length=0.75, EA=1.0e7, GAy=1.0e6, GAz=1.0e6, GJ=1.0e4, EIy=1.5e4, EIz=5.0e6, stiffness_modifier=None, ): node_r = np.zeros((num_node, 3)) node_r[:, 1] = np.linspace(0.0, span, num_node) model = gc.AeroelasticInformation() structure = model.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) return model

0.4 SHARPy 静力求解辅助函数¶

In [5]:

def solver_file(case_name, flow, *, route=None):
    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"
    return sim, route

def run_case(case_name, model, 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)
        model.StructuralInformation.generate_fem_file(str(route) + "/", case_name)
        sim.generate_solver_file()
        data = sharpy.sharpy_main.main(["", str(route / f"{case_name}.sharpy")])
    return data, log.getvalue()

def solver_file(case_name, flow, *, route=None): 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" return sim, route def run_case(case_name, model, 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) model.StructuralInformation.generate_fem_file(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 后处理辅助函数¶

In [6]:

def initial_positions(data):
    return data.structure.ini_info.pos.copy()

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

def nodal_displacements(data):
    return node_positions(data) - initial_positions(data)

def structural_tip(data):
    return node_positions(data)[-1].copy()

def tip_displacement(data):
    return nodal_displacements(data)[-1].copy()

def axial_strain_from_centerline(data):
    initial = initial_positions(data)
    deformed = node_positions(data)
    ds0 = np.linalg.norm(np.diff(initial, axis=0), axis=1)
    ds = np.linalg.norm(np.diff(deformed, axis=0), axis=1)
    y_mid = 0.5 * (initial[:-1, 1] + initial[1:, 1])
    return y_mid, ds / ds0 - 1.0

def inplane_curvature_z(data):
    y0 = initial_positions(data)[:, 1]
    deformed = node_positions(data)
    x = deformed[:, 0]
    y = deformed[:, 1]
    dx_dy = np.gradient(x, y0, edge_order=2)
    dy_dy = np.gradient(y, y0, edge_order=2)
    d2x_dy2 = np.gradient(dx_dy, y0, edge_order=2)
    d2y_dy2 = np.gradient(dy_dy, y0, edge_order=2)
    denominator = np.maximum((dx_dy**2 + dy_dy**2) ** 1.5, 1e-14)
    return y0, (dx_dy * d2y_dy2 - dy_dy * d2x_dy2) / denominator

def nominal_bending_stress(y, *, span, tip_force_z, weight_per_length=0.0, section_modulus=1.0e-5):
    lever = np.maximum(span - y, 0.0)
    moment = abs(tip_force_z) * lever + abs(weight_per_length) * lever**2 / 2.0
    return moment / section_modulus

def initial_positions(data): return data.structure.ini_info.pos.copy() def node_positions(data): return data.structure.timestep_info[-1].pos.copy() def nodal_displacements(data): return node_positions(data) - initial_positions(data) def structural_tip(data): return node_positions(data)[-1].copy() def tip_displacement(data): return nodal_displacements(data)[-1].copy() def axial_strain_from_centerline(data): initial = initial_positions(data) deformed = node_positions(data) ds0 = np.linalg.norm(np.diff(initial, axis=0), axis=1) ds = np.linalg.norm(np.diff(deformed, axis=0), axis=1) y_mid = 0.5 * (initial[:-1, 1] + initial[1:, 1]) return y_mid, ds / ds0 - 1.0 def inplane_curvature_z(data): y0 = initial_positions(data)[:, 1] deformed = node_positions(data) x = deformed[:, 0] y = deformed[:, 1] dx_dy = np.gradient(x, y0, edge_order=2) dy_dy = np.gradient(y, y0, edge_order=2) d2x_dy2 = np.gradient(dx_dy, y0, edge_order=2) d2y_dy2 = np.gradient(dy_dy, y0, edge_order=2) denominator = np.maximum((dx_dy**2 + dy_dy**2) ** 1.5, 1e-14) return y0, (dx_dy * d2y_dy2 - dy_dy * d2x_dy2) / denominator def nominal_bending_stress(y, *, span, tip_force_z, weight_per_length=0.0, section_modulus=1.0e-5): lever = np.maximum(span - y, 0.0) moment = abs(tip_force_z) * lever + abs(weight_per_length) * lever**2 / 2.0 return moment / section_modulus

0.6 只显示结构梁的 3D 视图¶

In [7]:

def show_beam_model(data, title):
    initial = initial_positions(data)
    deformed = node_positions(data)
    traces = []

    for positions, name, color, width in [
        (initial, "undeformed beam", "#9CA3AF", 4),
        (deformed, "deformed beam", "#3B6FB6", 7),
    ]:
        for ielem in range(data.structure.num_elem):
            nodes = data.structure.connectivities[ielem, :][[0, 2, 1]]
            xyz = positions[nodes, :]
            traces.append({
                "type": "scatter3d",
                "mode": "lines+markers",
                "x": xyz[:, 0].tolist(),
                "y": xyz[:, 1].tolist(),
                "z": xyz[:, 2].tolist(),
                "line": {"color": color, "width": width},
                "marker": {"size": 3, "color": color},
                "name": name,
                "showlegend": ielem == 0,
            })

    div_id = "plotly-" + uuid.uuid4().hex
    layout = {
        "title": {"text": title, "font": {"family": "Times New Roman, Times, serif", "size": 18}},
        "height": 560,
        "margin": {"l": 0, "r": 0, "t": 42, "b": 0},
        "legend": {"x": 0.02, "y": 0.98, "font": {"family": "Times New Roman, Times, serif", "size": 13}},
        "scene": {
            "xaxis": {"title": "x [m]"},
            "yaxis": {"title": "spanwise y [m]"},
            "zaxis": {"title": "z [m]"},
            "aspectmode": "data",
        },
        "font": {"family": "Times New Roman, Times, serif", "size": 13},
    }
    config = {"responsive": True, "displaylogo": False}
    html = f'''
    <div id="{div_id}" style="width: 100%; height: 560px;"></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_beam_model(data, title): initial = initial_positions(data) deformed = node_positions(data) traces = [] for positions, name, color, width in [ (initial, "undeformed beam", "#9CA3AF", 4), (deformed, "deformed beam", "#3B6FB6", 7), ]: for ielem in range(data.structure.num_elem): nodes = data.structure.connectivities[ielem, :][[0, 2, 1]] xyz = positions[nodes, :] traces.append({ "type": "scatter3d", "mode": "lines+markers", "x": xyz[:, 0].tolist(), "y": xyz[:, 1].tolist(), "z": xyz[:, 2].tolist(), "line": {"color": color, "width": width}, "marker": {"size": 3, "color": color}, "name": name, "showlegend": ielem == 0, }) div_id = "plotly-" + uuid.uuid4().hex layout = { "title": {"text": title, "font": {"family": "Times New Roman, Times, serif", "size": 18}}, "height": 560, "margin": {"l": 0, "r": 0, "t": 42, "b": 0}, "legend": {"x": 0.02, "y": 0.98, "font": {"family": "Times New Roman, Times, serif", "size": 13}}, "scene": { "xaxis": {"title": "x [m]"}, "yaxis": {"title": "spanwise y [m]"}, "zaxis": {"title": "z [m]"}, "aspectmode": "data", }, "font": {"family": "Times New Roman, Times, serif", "size": 13}, } config = {"responsive": True, "displaylogo": False} html = f''' <div id="{div_id}" style="width: 100%; height: 560px;"></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. 定义三种刚度模型¶

梁截面刚度矩阵的基本形式为

$$ \mathbf{K}=\begin{bmatrix} EA & 0 & 0 & 0 & 0 & 0\\ 0 & GA_y & 0 & 0 & 0 & 0\\ 0 & 0 & GA_z & 0 & 0 & 0\\ 0 & 0 & 0 & GJ & 0 & 0\\ 0 & 0 & 0 & 0 & EI_y & 0\\ 0 & 0 & 0 & 0 & 0 & EI_z \end{bmatrix}. $$

三种模型的含义如下：

• Case 1: uniform stiffness：全梁使用同一个对角刚度矩阵，作为基准；
• Case 2: coupled stiffness：加入对称耦合项，尤其是 $N_x$-$M_z$ 的轴向-面内弯曲耦合；
• Case 3: tapered stiffness：截面刚度沿展向逐渐降低，模拟从根部到端部变柔的梁。

对 Case 2，可把轴向拉伸和面内弯曲的主要耦合块写成

$$ \begin{bmatrix}N_x\\ M_z\end{bmatrix} = \begin{bmatrix}EA & C_{xz}\\ C_{xz} & EI_z\end{bmatrix} \begin{bmatrix}\epsilon_x\\ \kappa_z\end{bmatrix}. $$

因此在相同面内弯曲载荷下，Case 2 的 $\kappa_z$ 不再只由 $EI_z$ 决定，而会偏离 Case 1 的基准曲率。

In [8]:

def coupled_stiffness(structure):
    K = structure.stiffness_db[0].copy()
    pairs = {
        (0, 5): 0.28,
        (0, 3): 0.06,
        (3, 5): -0.08,
        (4, 5): 0.10,
    }
    for (i, j), rho in pairs.items():
        value = rho * np.sqrt(K[i, i] * K[j, j])
        K[i, j] = value
        K[j, i] = value
    structure.stiffness_db[0] = K

def tapered_stiffness(structure):
    base = structure.stiffness_db[0]
    span_pos = np.linspace(0.0625, 0.9375, structure.num_elem)
    scales = 1.25 + (0.55 - 1.25) * span_pos
    structure.stiffness_db = np.stack([scale * base for scale in scales])
    structure.elem_stiffness = np.arange(structure.num_elem)

def coupled_stiffness(structure): K = structure.stiffness_db[0].copy() pairs = { (0, 5): 0.28, (0, 3): 0.06, (3, 5): -0.08, (4, 5): 0.10, } for (i, j), rho in pairs.items(): value = rho * np.sqrt(K[i, i] * K[j, j]) K[i, j] = value K[j, i] = value structure.stiffness_db[0] = K def tapered_stiffness(structure): base = structure.stiffness_db[0] span_pos = np.linspace(0.0625, 0.9375, structure.num_elem) scales = 1.25 + (0.55 - 1.25) * span_pos structure.stiffness_db = np.stack([scale * base for scale in scales]) structure.elem_stiffness = np.arange(structure.num_elem)

1.1 组织算例标签¶

In [9]:

variants = {
    "Case 1: uniform stiffness": ("stiffness_case1_uniform", None),
    "Case 2: coupled stiffness": ("stiffness_case2_coupled", coupled_stiffness),
    "Case 3: tapered stiffness": ("stiffness_case3_tapered", tapered_stiffness),
}

for label, (case_id, modifier) in variants.items():
    beam = make_beam(case_id, stiffness_modifier=modifier)
    K = beam.StructuralInformation.stiffness_db[0]
    min_eig = np.linalg.eigvalsh(K).min()
    print(f"{label}: nodes={beam.StructuralInformation.num_node}, elements={beam.StructuralInformation.num_elem}, min eig(K)={min_eig:.3e}")

variants = { "Case 1: uniform stiffness": ("stiffness_case1_uniform", None), "Case 2: coupled stiffness": ("stiffness_case2_coupled", coupled_stiffness), "Case 3: tapered stiffness": ("stiffness_case3_tapered", tapered_stiffness), } for label, (case_id, modifier) in variants.items(): beam = make_beam(case_id, stiffness_modifier=modifier) K = beam.StructuralInformation.stiffness_db[0] min_eig = np.linalg.eigvalsh(K).min() print(f"{label}: nodes={beam.StructuralInformation.num_node}, elements={beam.StructuralInformation.num_elem}, min eig(K)={min_eig:.3e}")

Case 1: uniform stiffness: nodes=17, elements=8, min eig(K)=1.000e+04
Case 2: coupled stiffness: nodes=17, elements=8, min eig(K)=9.859e+03
Case 3: tapered stiffness: nodes=17, elements=8, min eig(K)=1.206e+04

1.2 比较刚度矩阵¶

色标越亮表示矩阵项量级越大。

In [10]:

fig, axes = plt.subplots(1, 3, figsize=(10.2, 3.4), constrained_layout=True)
tick_labels = ["EA", "GAy", "GAz", "GJ", "EIy", "EIz"]
for panel, ax, (label, (case_id, modifier)) in zip(["a", "b", "c"], axes, variants.items()):
    beam = make_beam(case_id, stiffness_modifier=modifier)
    matrix = beam.StructuralInformation.stiffness_db[0]
    im = ax.imshow(np.log10(np.abs(matrix) + 1.0), cmap="cividis")
    label_panel(ax, panel)
    ax.set_title(label, loc="left", pad=7)
    ax.set_xticks(range(6))
    ax.set_yticks(range(6))
    ax.set_xticklabels(tick_labels, rotation=45, ha="right")
    ax.set_yticklabels(tick_labels)
    ax.tick_params(length=0)
    for spine in ax.spines.values():
        spine.set_visible(False)
cbar = fig.colorbar(im, ax=axes, shrink=0.78, pad=0.02)
cbar.set_label("$\\log_{10}(|K_{ij}|+1)$")
plt.show()

fig, axes = plt.subplots(1, 3, figsize=(10.2, 3.4), constrained_layout=True) tick_labels = ["EA", "GAy", "GAz", "GJ", "EIy", "EIz"] for panel, ax, (label, (case_id, modifier)) in zip(["a", "b", "c"], axes, variants.items()): beam = make_beam(case_id, stiffness_modifier=modifier) matrix = beam.StructuralInformation.stiffness_db[0] im = ax.imshow(np.log10(np.abs(matrix) + 1.0), cmap="cividis") label_panel(ax, panel) ax.set_title(label, loc="left", pad=7) ax.set_xticks(range(6)) ax.set_yticks(range(6)) ax.set_xticklabels(tick_labels, rotation=45, ha="right") ax.set_yticklabels(tick_labels) ax.tick_params(length=0) for spine in ax.spines.values(): spine.set_visible(False) cbar = fig.colorbar(im, ax=axes, shrink=0.78, pad=0.02) cbar.set_label("$\\log_{10}(|K_{ij}|+1)$") plt.show()

2. 设置拉伸和面内弯曲两个静力工况¶

为了直接比较三类刚度模型的差别，这里不再只看一个端部竖向力，而是分别施加两个载荷。由于梁局部 $x$ 轴沿展向布置，app_forces 的第一个分量对应局部轴向拉伸，第二个分量对应 $x$-$y$ 平面内的横向弯曲载荷：

$$ \mathbf{F}_{\mathrm{tension}}=(F_x,0,0), \qquad F_x=2500\ \mathrm{N}, $$

$$ \mathbf{F}_{\mathrm{in\mbox{-}plane}}=(0,F_y,0), \qquad F_y=2500\ \mathrm{N}. $$

第一个工况主要考察拉伸载荷下沿梁长的轴向应变 $\epsilon_x$。第二个工况让梁在 $x$-$y$ 平面内弯曲，主要考察曲率 $\kappa_z$；这也是 Case 2 的 $N_x$-$M_z$ 耦合项最容易显现的地方。求解流程仍然只有结构加载和非线性静力求解。

In [11]:

load_cases = {
    "Tensile load": {
        "suffix": "tension",
        "force": np.array([2500.0, 0.0, 0.0]),
    },
    "In-plane bending load": {
        "suffix": "inplane_bending",
        "force": np.array([0.0, 2500.0, 0.0]),
    },
}

results = {load_label: {} for load_label in load_cases}

for load_label, load in load_cases.items():
    print(f"\n{load_label}: tip force = {load['force']} N")
    for label, (base_case_name, modifier) in variants.items():
        case_name = f"{base_case_name}_{load['suffix']}"
        beam = make_beam(case_name, stiffness_modifier=modifier)
        beam.StructuralInformation.app_forces[-1, :3] = load["force"]
        sim, route = solver_file(case_name, ["BeamLoader", "NonLinearStatic"])
        sim.solvers["NonLinearStatic"]["gravity_on"] = False
        sim.solvers["NonLinearStatic"]["num_load_steps"] = 5
        sim.solvers["NonLinearStatic"]["max_iterations"] = 100
        sim.solvers["NonLinearStatic"]["min_delta"] = 1e-7
        data, log = run_case(case_name, beam, sim, route)
        results[load_label][label] = data
        tip_u = tip_displacement(data)
        print(
            f"{label}: tip ux={tip_u[0]*1000:.3f} mm, "
            f"uy={tip_u[1]*1000:.3f} mm, uz={tip_u[2]*1000:.3f} mm, "
            f"log lines={len(log.splitlines())}"
        )

load_cases = { "Tensile load": { "suffix": "tension", "force": np.array([2500.0, 0.0, 0.0]), }, "In-plane bending load": { "suffix": "inplane_bending", "force": np.array([0.0, 2500.0, 0.0]), }, } results = {load_label: {} for load_label in load_cases} for load_label, load in load_cases.items(): print(f"\n{load_label}: tip force = {load['force']} N") for label, (base_case_name, modifier) in variants.items(): case_name = f"{base_case_name}_{load['suffix']}" beam = make_beam(case_name, stiffness_modifier=modifier) beam.StructuralInformation.app_forces[-1, :3] = load["force"] sim, route = solver_file(case_name, ["BeamLoader", "NonLinearStatic"]) sim.solvers["NonLinearStatic"]["gravity_on"] = False sim.solvers["NonLinearStatic"]["num_load_steps"] = 5 sim.solvers["NonLinearStatic"]["max_iterations"] = 100 sim.solvers["NonLinearStatic"]["min_delta"] = 1e-7 data, log = run_case(case_name, beam, sim, route) results[load_label][label] = data tip_u = tip_displacement(data) print( f"{label}: tip ux={tip_u[0]*1000:.3f} mm, " f"uy={tip_u[1]*1000:.3f} mm, uz={tip_u[2]*1000:.3f} mm, " f"log lines={len(log.splitlines())}" )

Tensile load: tip force = [2500.    0.    0.] N

Case 1: uniform stiffness: tip ux=-0.000 mm, uy=0.500 mm, uz=0.000 mm, log lines=0
Case 2: coupled stiffness: tip ux=-0.225 mm, uy=0.547 mm, uz=0.482 mm, log lines=0
Case 3: tapered stiffness: tip ux=-0.000 mm, uy=0.586 mm, uz=0.000 mm, log lines=0

In-plane bending load: tip force = [   0. 2500.    0.] N
Case 1: uniform stiffness: tip ux=6.333 mm, uy=-0.004 mm, uz=0.000 mm, log lines=0
Case 2: coupled stiffness: tip ux=6.476 mm, uy=-0.231 mm, uz=-3.280 mm, log lines=0
Case 3: tapered stiffness: tip ux=7.125 mm, uy=-0.005 mm, uz=0.000 mm, log lines=0

3. 比较 $\epsilon_x$ 与 $\kappa_z$¶

左图使用拉伸工况，从变形中心线相邻节点间距提取段平均轴向应变

$$ \epsilon_x \approx \frac{\Delta s-\Delta s_0}{\Delta s_0}. $$

右图使用面内弯曲工况，从变形中心线在 $x$-$y$ 平面内的几何曲率提取 $\kappa_z$

$$ \kappa_z = \frac{x' y'' - y' x''}{\left(x'^2+y'^2\right)^{3/2}}. $$

这样一来，Case 3 主要体现展向刚度变化带来的应变/曲率分布变化；Case 2 则会因为耦合项使 $\kappa_z$ 相对 Case 1 出现系统性偏离。图内文字保持英文。

In [12]:

fig, axes = plt.subplots(1, 2, figsize=(11.0, 4.0), constrained_layout=True)

ax = axes[0]
for color, (label, data) in zip(NATURE_SEQUENCE, results["Tensile load"].items()):
    y_mid, epsilon_x = axial_strain_from_centerline(data)
    ax.plot(y_mid, epsilon_x * 1e6, "o-", color=color, label=label)
style_axes(
    ax,
    xlabel="spanwise coordinate y [m]",
    ylabel=r"$\epsilon_x$ [$\mu\epsilon$]",
    title="Axial Strain under Tensile Load",
    zero_line=True,
)
label_panel(ax, "a")
ax.legend(loc="best")

ax = axes[1]
for color, (label, data) in zip(NATURE_SEQUENCE, results["In-plane bending load"].items()):
    y0, kappa_z = inplane_curvature_z(data)
    ax.plot(y0, kappa_z * 1e3, "o-", color=color, label=label)
style_axes(
    ax,
    xlabel="spanwise coordinate y [m]",
    ylabel=r"$\kappa_z$ [$10^{-3}$ 1/m]",
    title="In-plane Curvature under Bending Load",
    zero_line=True,
)
label_panel(ax, "b")
ax.legend(loc="best")

finish_figure(fig)

fig, axes = plt.subplots(1, 2, figsize=(11.0, 4.0), constrained_layout=True) ax = axes[0] for color, (label, data) in zip(NATURE_SEQUENCE, results["Tensile load"].items()): y_mid, epsilon_x = axial_strain_from_centerline(data) ax.plot(y_mid, epsilon_x * 1e6, "o-", color=color, label=label) style_axes( ax, xlabel="spanwise coordinate y [m]", ylabel=r"$\epsilon_x$ [$\mu\epsilon$]", title="Axial Strain under Tensile Load", zero_line=True, ) label_panel(ax, "a") ax.legend(loc="best") ax = axes[1] for color, (label, data) in zip(NATURE_SEQUENCE, results["In-plane bending load"].items()): y0, kappa_z = inplane_curvature_z(data) ax.plot(y0, kappa_z * 1e3, "o-", color=color, label=label) style_axes( ax, xlabel="spanwise coordinate y [m]", ylabel=r"$\kappa_z$ [$10^{-3}$ 1/m]", title="In-plane Curvature under Bending Load", zero_line=True, ) label_panel(ax, "b") ax.legend(loc="best") finish_figure(fig)

/var/folders/s3/vnf693hd3mg005d8ny9xwh900000gn/T/ipykernel_95696/2868343704.py:81: UserWarning: The figure layout has changed to tight
  fig.tight_layout()

小结¶

这个对比案例现在用两个更有辨识度的静力工况来读三类刚度模型。拉伸载荷下的 $\epsilon_x$ 图可以直接比较 Case 1、Case 2 和 Case 3 的轴向柔度差异；面内弯曲载荷下的 $\kappa_z$ 图则突出 Case 2 的轴向-弯曲耦合项如何让曲率偏离均匀刚度基准。Case 3 的主要特征不是耦合，而是沿展向逐渐变柔，因此它的应变和曲率分布会更明显地随 $y$ 改变。