otko/examples/moment_curvature.py
smillmorel 612936a00b
Some checks failed
CI / lint (push) Has been cancelled
CI / test (macos-latest, 3.10) (push) Has been cancelled
CI / test (macos-latest, 3.11) (push) Has been cancelled
CI / test (macos-latest, 3.12) (push) Has been cancelled
CI / test (ubuntu-latest, 3.10) (push) Has been cancelled
CI / test (ubuntu-latest, 3.11) (push) Has been cancelled
CI / test (ubuntu-latest, 3.12) (push) Has been cancelled
CI / test (windows-latest, 3.10) (push) Has been cancelled
CI / test (windows-latest, 3.11) (push) Has been cancelled
CI / test (windows-latest, 3.12) (push) Has been cancelled
feat: initial otko import
2026-09-08 02:12:15 -04:00

249 lines
8.8 KiB
Python
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

"""Moment-Curvature Example — OpenSees Examples Manual, Example 2.
Reinforced-concrete column cross-section, fibre discretisation, axial
preload + monotonic moment pushover. Mirrors the OpenSees Tcl script
at https://opensees.berkeley.edu/wiki/index.php?title=Moment_Curvature_Example
Model
-----
Two coincident nodes linked by a ``zeroLengthSection`` carrying the
RC fibre section. Node 1 is fully restrained; node 2 is free in Ux
(so axial can shorten) and Rz (the curvature DOF). A constant axial
load P = -180 kip is applied first via LoadControl(0); then a
linear reference moment pattern (Mz = 1 kip·in) is added and
DisplacementControl on DOF 3 ramps the curvature to μ·Ky where
μ = 15 and Ky is the elastic yield curvature estimate.
Units: kip, in, ksi (UnitSystem.US_IN_KIP).
GUI walkthrough: File → Open → moment_curvature.osmodel, Options →
Set Display Units → US (in, kip, kip·s²/in, s, ksi), Analyze → Run →
MK, Display → Show Pushover Curve.
"""
from __future__ import annotations
from pathlib import Path
from otko.core import (
Concrete01,
ConstantTimeSeries,
CoordinateGridSystem,
FiberSection,
GridSystem,
LinearTimeSeries,
NodalLoad,
Node,
PlainLoadPattern,
Project,
ProjectMeta,
PushoverCase,
RectangularPatch,
Steel01,
StraightLayer,
UnitSystem,
ZeroLengthSectionElement,
make_grid_lines,
)
from otko.services import load_project, save_project
# Cross-section parameters (kip-in-ksi, from the Tcl example).
COL_WIDTH = 15.0 # z-direction dimension
COL_DEPTH = 24.0 # y-direction dimension
COVER = 1.5
AS_BAR = 0.60 # area of one #7 rebar
FY = 60.0 # steel yield stress, ksi
E_STEEL = 30000.0 # steel Young's modulus, ksi
HARDENING = 0.01 # strain-hardening ratio
P_AXIAL = -180.0 # kip, compression
MU = 15 # target curvature ductility
NUM_INCR = 100 # DisplacementControl increments
def build_moment_curvature() -> Project:
"""Build the OpenSees Example 2 Moment-Curvature project."""
y1 = COL_DEPTH / 2.0 # 12
z1 = COL_WIDTH / 2.0 # 7.5
d = COL_DEPTH - COVER # 22.5
# Yield curvature estimate, assumed elastic + top/bottom steel only.
eps_y = FY / E_STEEL # 0.002
ky = eps_y / (0.7 * d) # ≈ 1.27e-4
max_k = ky * MU # ≈ 1.905e-3
d_k = max_k / NUM_INCR # per-step curvature increment
return Project(
meta=ProjectMeta(
name="Moment-Curvature (OpenSees Ex 2)",
author="OpenSees Examples Manual",
description=(
"RC column fibre section — constant axial P + DisplacementControl "
"curvature pushover. Kip-in-ksi units throughout."
),
units=UnitSystem.US_IN_KIP,
),
ndm=2, ndf=3,
coord_systems=[
# A tiny grid at the section origin so the two coincident
# nodes have a visual anchor in the canvas.
CoordinateGridSystem(
name="Global",
grid=GridSystem(
x_grid_lines=make_grid_lines("X", [0.0]),
y_grid_lines=make_grid_lines("Y", [0.0]),
z_grid_lines=make_grid_lines("Z", [0.0]),
),
),
],
nodes=[
# Node 1 — fully clamped.
Node(id=1, name="Support",
coords=(0.0, 0.0, 0.0),
restraint=(True, True, False, False, False, True)),
# Node 2 — free in Ux and Rz (axial + curvature).
Node(id=2, name="Crown",
coords=(0.0, 0.0, 0.0),
restraint=(False, True, False, False, False, False)),
],
materials=[
# Core concrete — confined (tag 1).
Concrete01(id=1, name="Core-Conc",
fpc=-6.0, epsc0=-0.004,
fpcu=-5.0, epsU=-0.014),
# Cover concrete — unconfined (tag 2).
Concrete01(id=2, name="Cover-Conc",
fpc=-5.0, epsc0=-0.002,
fpcu=0.0, epsU=-0.006),
# Reinforcing steel — bilinear hardening (tag 3).
Steel01(id=3, name="Steel-60",
Fy=FY, E0=E_STEEL, b=HARDENING),
],
sections=[
FiberSection(
id=1, name="RC-Column",
patches=[
# Core — confined concrete inside the rebar ring.
RectangularPatch(
material_id=1, n_fib_y=10, n_fib_z=1,
y_i=COVER - y1, z_i=COVER - z1,
y_j=y1 - COVER, z_j=z1 - COVER,
),
# Top cover (unconfined).
RectangularPatch(
material_id=2, n_fib_y=10, n_fib_z=1,
y_i=-y1, z_i=z1 - COVER,
y_j=y1, z_j=z1,
),
# Bottom cover.
RectangularPatch(
material_id=2, n_fib_y=10, n_fib_z=1,
y_i=-y1, z_i=-z1,
y_j=y1, z_j=COVER - z1,
),
# Left cover.
RectangularPatch(
material_id=2, n_fib_y=2, n_fib_z=1,
y_i=-y1, z_i=COVER - z1,
y_j=COVER - y1, z_j=z1 - COVER,
),
# Right cover.
RectangularPatch(
material_id=2, n_fib_y=2, n_fib_z=1,
y_i=y1 - COVER, z_i=COVER - z1,
y_j=y1, z_j=z1 - COVER,
),
],
layers=[
# Top rebar (3 × #7).
StraightLayer(
material_id=3, n_bars=3, bar_area=AS_BAR,
y_start=y1 - COVER, z_start=z1 - COVER,
y_end=y1 - COVER, z_end=COVER - z1,
),
# Middle rebar (2 × #7).
StraightLayer(
material_id=3, n_bars=2, bar_area=AS_BAR,
y_start=0.0, z_start=z1 - COVER,
y_end=0.0, z_end=COVER - z1,
),
# Bottom rebar (3 × #7).
StraightLayer(
material_id=3, n_bars=3, bar_area=AS_BAR,
y_start=COVER - y1, z_start=z1 - COVER,
y_end=COVER - y1, z_end=COVER - z1,
),
],
),
],
elements=[
ZeroLengthSectionElement(
id=1, name="MK-Link",
nodes=(1, 2), section_id=1,
),
],
time_series=[
ConstantTimeSeries(id=1, name="AxialP"),
LinearTimeSeries(id=2, name="RefMoment"),
],
load_patterns=[
# Constant axial preload at node 2 — Fx = P (compression).
PlainLoadPattern(
id=1, name="AxialP",
time_series_id=1,
nodal_loads=[
NodalLoad(node_id=2,
forces=(P_AXIAL, 0, 0, 0, 0, 0)),
],
),
# Linear reference moment — Mz = 1.0, DisplacementControl
# scales this as it ramps curvature.
PlainLoadPattern(
id=2, name="RefMoment",
time_series_id=2,
nodal_loads=[
NodalLoad(node_id=2,
forces=(0, 0, 0, 0, 0, 1.0)),
],
),
],
analyses=[
PushoverCase(
id=1, name="MK",
pattern_ids=[1, 2],
control_node=2, control_dof=3, # Rz = curvature
target_disp=max_k,
step_size=d_k,
base_nodes=[1],
test="NormUnbalance",
tolerance=1e-9, max_iter=25,
),
],
)
def main() -> None:
project = build_moment_curvature()
project.validate_references()
y1 = COL_DEPTH / 2.0
eps_y = FY / E_STEEL
ky = eps_y / (0.7 * (COL_DEPTH - COVER))
print(f"Built '{project.meta.name}'")
print(f" ndm={project.ndm}, ndf={project.ndf}, units={project.meta.units.value}")
print(f" Estimated yield curvature Ky = {ky:.4e} 1/in")
print(f" Target (mu * Ky) = {ky * MU:.4e} 1/in (mu = {MU})")
out_path = Path(__file__).with_suffix(".osmodel")
save_project(project, out_path)
print(f"Saved -> {out_path}")
restored = load_project(out_path)
restored.validate_references()
assert restored.model_dump(by_alias=True) == project.model_dump(by_alias=True)
print("Round-trip OK.")
if __name__ == "__main__":
main()