案例一：悬臂梁基础建模与重力/端部载荷静力分析¶

下载本案例 Jupyter Notebook

这一节只研究一根单悬臂梁。梁根部固定，梁轴沿展向坐标 $y$ 从 $0$ 延伸到 $L=2\ \mathrm{m}$。

本案例将实现三件事：

1. 用 SHARPy 的结构梁数据结构生成一个悬臂梁模型；
2. 同时施加重力和端部向下载荷，使端部竖向位移达到清楚可见的量级；
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

Matplotlib is building the font cache; this may take a moment.

0.1 图表样式¶

图表内文字统一使用英文，并优先采用 Times New Roman。

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

每次运行会把 SHARPy 输入文件写到 _sharpy_runs/ 下，避免污染源码目录。

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 梁结构构造函数¶

这里的 make_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 静力求解辅助函数¶

求解流程只包含 BeamLoader 和 NonLinearStatic。

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 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 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. 建立最小悬臂梁¶

梁的截面刚度矩阵采用 SHARPy 常用的 6 自由度顺序

$$ \mathbf{K}=\mathrm{diag}(EA,GA_y,GA_z,GJ,EI_y,EI_z). $$

本例将竖向弯曲刚度设置得相对柔一些，使 2 m 梁在较大端部力下产生肉眼可见的静力位移。根部边界条件为 1，表示固支；端部边界条件为 -1，表示无约束。

In [8]:

case_name = "case01_gravity_tip_load"
span = 2.0
section_modulus = 1.0e-5
model = make_beam(case_name, span=span, EIy=1.5e4, EIz=5.0e6, mass_per_unit_length=0.75)
structure = model.StructuralInformation
K = structure.stiffness_db[0]

print("number of nodes:", structure.num_node)
print("number of elements:", structure.num_elem)
print("root boundary condition:", structure.boundary_conditions[0])
print("tip boundary condition:", structure.boundary_conditions[-1])
print("stiffness diagonal [EA, GAy, GAz, GJ, EIy, EIz]:")
print(np.diag(K))

case_name = "case01_gravity_tip_load" span = 2.0 section_modulus = 1.0e-5 model = make_beam(case_name, span=span, EIy=1.5e4, EIz=5.0e6, mass_per_unit_length=0.75) structure = model.StructuralInformation K = structure.stiffness_db[0] print("number of nodes:", structure.num_node) print("number of elements:", structure.num_elem) print("root boundary condition:", structure.boundary_conditions[0]) print("tip boundary condition:", structure.boundary_conditions[-1]) print("stiffness diagonal [EA, GAy, GAz, GJ, EIy, EIz]:") print(np.diag(K))

number of nodes: 17
number of elements: 8
root boundary condition: 1
tip boundary condition: -1
stiffness diagonal [EA, GAy, GAz, GJ, EIy, EIz]:
[10000000.  1000000.  1000000.    10000.    15000.  5000000.]

2. 施加重力和端部向下载荷¶

端部集中力直接写入最后一个结构节点的 app_forces。这里取

$$ \mathbf{F}_{\mathrm{tip}}=(0,0,-1600)\ \mathrm{N}, $$

在本例的 SHARPy 梁输入约定下，app_forces[-1, 2] 取正值会得到负的 $u_z$ 位移，也就是图中的向下弯曲。对长度 $L=2\ \mathrm{m}$、弯曲刚度约 $EI=1.5\times10^4\ \mathrm{N\,m^2}$ 的线性梁，端部力本身给出的量级估计为

$$ |u_z(L)|\approx \frac{|F_z|L^3}{3EI}\approx 0.28\ \mathrm{m}, $$

In [ ]:

tip_force_z = 1600.0
model.StructuralInformation.app_forces[-1, 2] = tip_force_z

# 以下就是在定义求解流
sim, route = solver_file(case_name, ["BeamLoader", "NonLinearStatic"])
sim.solvers["NonLinearStatic"]["gravity_on"] = True
sim.solvers["NonLinearStatic"]["gravity"] = 9.81
sim.solvers["NonLinearStatic"]["gravity_dir"] = np.array([0.0, 0.0, -1.0])
sim.solvers["NonLinearStatic"]["num_load_steps"] = 5
sim.solvers["NonLinearStatic"]["max_iterations"] = 100
sim.solvers["NonLinearStatic"]["min_delta"] = 1e-7

print(f"tip force z: {tip_force_z:.1f} N")
print("solver flow:", sim.solvers["SHARPy"]["flow"])

tip_force_z = 1600.0 model.StructuralInformation.app_forces[-1, 2] = tip_force_z # 以下就是在定义求解流 sim, route = solver_file(case_name, ["BeamLoader", "NonLinearStatic"]) sim.solvers["NonLinearStatic"]["gravity_on"] = True sim.solvers["NonLinearStatic"]["gravity"] = 9.81 sim.solvers["NonLinearStatic"]["gravity_dir"] = np.array([0.0, 0.0, -1.0]) sim.solvers["NonLinearStatic"]["num_load_steps"] = 5 sim.solvers["NonLinearStatic"]["max_iterations"] = 100 sim.solvers["NonLinearStatic"]["min_delta"] = 1e-7 print(f"tip force z: {tip_force_z:.1f} N") print("solver flow:", sim.solvers["SHARPy"]["flow"])

tip force z: 1600.0 N
solver flow: ['BeamLoader', 'NonLinearStatic']

3. 运行静力求解¶

In [10]:

data, log = run_case(case_name, model, sim, route)
tip_pos = structural_tip(data)
tip_u = tip_displacement(data)

print(f"tip position: x={tip_pos[0]:.4f} m, y={tip_pos[1]:.4f} m, z={tip_pos[2]:.4f} m")
print(f"tip displacement: ux={tip_u[0]*1000:.2f} mm, uy={tip_u[1]*1000:.2f} mm, uz={tip_u[2]*1000:.2f} mm")
print(f"tip displacement magnitude: {np.linalg.norm(tip_u)*1000:.2f} mm")
print(f"SHARPy log lines captured: {len(log.splitlines())}")

data, log = run_case(case_name, model, sim, route) tip_pos = structural_tip(data) tip_u = tip_displacement(data) print(f"tip position: x={tip_pos[0]:.4f} m, y={tip_pos[1]:.4f} m, z={tip_pos[2]:.4f} m") print(f"tip displacement: ux={tip_u[0]*1000:.2f} mm, uy={tip_u[1]*1000:.2f} mm, uz={tip_u[2]*1000:.2f} mm") print(f"tip displacement magnitude: {np.linalg.norm(tip_u)*1000:.2f} mm") print(f"SHARPy log lines captured: {len(log.splitlines())}")

tip position: x=0.0000 m, y=1.9756 m, z=-0.2846 m
tip displacement: ux=0.00 mm, uy=-24.43 mm, uz=-284.62 mm
tip displacement magnitude: 285.66 mm
SHARPy log lines captured: 0

4. 查看纯结构 Beam 网格和变形¶

In [11]:

show_beam_model(data, "Cantilever Beam: Undeformed and Deformed Shape")

show_beam_model(data, "Cantilever Beam: Undeformed and Deformed Shape")

5. 位移曲线¶

In [12]:

initial = initial_positions(data)
deformed = node_positions(data)

y_initial = initial[:, 1]
z_initial = initial[:, 2]
y_deformed = deformed[:, 1]
z_deformed = deformed[:, 2]

fig, ax = plt.subplots(figsize=(6.4, 4.0))
ax.plot(
    y_initial,
    z_initial,
    "o-",
    color=NATURE_COLORS["grey"],
    label="undeformed nodes $(Y_0,Z_0)$",
)
ax.plot(
    y_deformed,
    z_deformed,
    "o-",
    color=NATURE_COLORS["blue"],
    label="deformed nodes $(Y,Z)$",
)
style_axes(
    ax,
    xlabel="spanwise coordinate y [m]",
    ylabel="z coordinate z [m]",
    title="Cantilever Beam Node Positions in the y-z Plane",
    zero_line=True,
)
ax.legend()
finish_figure(fig)

print(f"tip initial position: y={y_initial[-1]:.4f} m, z={z_initial[-1]:.4f} m")
print(f"tip deformed position: y={y_deformed[-1]:.4f} m, z={z_deformed[-1]:.4f} m")
print(f"tip y displacement: {(y_deformed[-1] - y_initial[-1])*1000:.2f} mm")
print(f"tip z displacement: {(z_deformed[-1] - z_initial[-1])*1000:.2f} mm")
assert abs(z_deformed[-1] - z_initial[-1]) >= 0.20, "Tip z displacement did not reach the 20 cm target."

initial = initial_positions(data) deformed = node_positions(data) y_initial = initial[:, 1] z_initial = initial[:, 2] y_deformed = deformed[:, 1] z_deformed = deformed[:, 2] fig, ax = plt.subplots(figsize=(6.4, 4.0)) ax.plot( y_initial, z_initial, "o-", color=NATURE_COLORS["grey"], label="undeformed nodes $(Y_0,Z_0)$", ) ax.plot( y_deformed, z_deformed, "o-", color=NATURE_COLORS["blue"], label="deformed nodes $(Y,Z)$", ) style_axes( ax, xlabel="spanwise coordinate y [m]", ylabel="z coordinate z [m]", title="Cantilever Beam Node Positions in the y-z Plane", zero_line=True, ) ax.legend() finish_figure(fig) print(f"tip initial position: y={y_initial[-1]:.4f} m, z={z_initial[-1]:.4f} m") print(f"tip deformed position: y={y_deformed[-1]:.4f} m, z={z_deformed[-1]:.4f} m") print(f"tip y displacement: {(y_deformed[-1] - y_initial[-1])*1000:.2f} mm") print(f"tip z displacement: {(z_deformed[-1] - z_initial[-1])*1000:.2f} mm") assert abs(z_deformed[-1] - z_initial[-1]) >= 0.20, "Tip z displacement did not reach the 20 cm target."

tip initial position: y=2.0000 m, z=0.0000 m
tip deformed position: y=1.9756 m, z=-0.2846 m
tip y displacement: -24.43 mm
tip z displacement: -284.62 mm