From b2f07d46c8ba556314e7d62fa3da2ed232763470 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 26 Oct 2026 15:12:42 +0900 Subject: [PATCH 1/3] feat(stiffeners): model tapered bend stiffeners on finite-EI lines A STIFFENERS table describes a polyurethane bend stiffener over End A or End B of a finite-EI line as conical segments from the root outward (length, outer diameter at each side, bore, modulus, density). The sleeve bends with the cable: every element it covers takes the compliance average of EI_cable + E pi/64 (OD^4 - ID^4) over its length, the sleeve's mean mass, a displaced area of cable plus solid sleeve, and the mean outer diameter for drag and added mass. Elements it does not reach keep their line type's values exactly, so decks without stiffeners are unchanged. The overlay sits in the two routines every finite-EI route builds its element properties from, so it applies standalone, on the OpenFAST aggregate, through the C API and in Python. Under a stiffener BendMomN (and the BendMoment columns) report the cable's own moment, the line type's EI times the curvature; the new StfMomN reports the sleeve's share. FairDefl/AnchDefl give the angle between the end force and the direction a held END CONNECTIONS row declares, the stiffener's load case with the end tension. Initialisation notes name a stiffener spanning fewer than 10 elements and one at a pinned end. The Python reader and builder mirror the native table and channel rules. On a hang-off clamped through a 5 m stiffener (16 times the cable's EI at the root), the curvature follows OrcaFlex 11.6d with its stiffener as a line attachment within 0.74 % of the root value, and an independent collocation solution of the planar extensible elastica within 0.74 % (cable moment), 0.18 % (section moment) and 0.003 deg (deflection). Example lazy_wave_bend_stiffener.dat. Signed-off-by: Jae Hoon Seo --- CHANGELOG.md | 20 + CMakeLists.txt | 19 +- VALIDATION.md | 13 + doc/capabilities.rst | 6 + doc/concepts.rst | 4 +- doc/driver_format.md | 80 +++- doc/examples.rst | 6 + doc/migrating.rst | 4 +- doc/modeling.rst | 38 +- doc/outputs.rst | 23 +- doc/references.rst | 6 + doc/theory.rst | 53 +++ examples/README.md | 2 +- examples/lazy_wave_bend_stiffener.dat | 72 ++++ python/cabledyn/builder.py | 93 +++- python/cabledyn/deck_file.py | 133 +++++- python/tests/test_builder.py | 21 + python/tests/test_deck_native_parity.py | 57 +++ src/CableDyn_DeckDriver.f90 | 401 +++++++++++++++++- src/CableDyn_HermiteCableDynamic.f90 | 48 ++- tests/check_bend_stiffener.py | 264 ++++++++++++ .../orcaflex_bend_stiffener_curvature.csv | 163 +++++++ tests/test_bend_stiffener.f90 | 332 +++++++++++++++ validation/scripts/orcaflex_bend_stiffener.py | 101 +++++ 24 files changed, 1902 insertions(+), 57 deletions(-) create mode 100644 examples/lazy_wave_bend_stiffener.dat create mode 100644 tests/check_bend_stiffener.py create mode 100644 tests/data/orcaflex_bend_stiffener_curvature.csv create mode 100644 tests/test_bend_stiffener.f90 create mode 100644 validation/scripts/orcaflex_bend_stiffener.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 37ef50c..41a39e8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,26 @@ All notable changes to CableDyn are recorded in this file. The format follows ## [Unreleased] +### Added + +- Bend stiffeners (`STIFFENERS` table): a tapered polyurethane sleeve over an end of a finite-EI + line, given as conical segments from the root outward (outer diameters, bore, modulus, + density). Each element it covers takes the compliance average of the cable-plus-sleeve + bending stiffness and the sleeve's mass, buoyancy and drag; standalone and coupled, with the + Python reader and builder. A hang-off clamped through a 5 m stiffener follows an independent + solution of the planar elastica within 0.2 % of the root moment, and its curvature follows + OrcaFlex 11.6d within 0.8 %. Initialisation notes name a stiffener spanning fewer than 10 + elements and one at a pinned end. Example `lazy_wave_bend_stiffener.dat`. +- Output channels `StfMomN`, the share of the section moment a bend stiffener carries, + and `FairDefl`/`AnchDefl`, the deflection of the end force from the direction a held + END CONNECTIONS row declares (with the end tension, a stiffener's load case). + +### Changed + +- Under a bend stiffener, `BendMomN` (and the `BendMoment` columns of the static profile, + range graphs and element file) report the cable's own moment, the line type's EI times the + curvature; lines without stiffeners are unchanged. + ## [0.3.0] - 2026-10-23 ### Added diff --git a/CMakeLists.txt b/CMakeLists.txt index 27126a4..6c66403 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1087,6 +1087,14 @@ file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/line_attachments) add_test(NAME line_attachments COMMAND test_line_attachments WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/line_attachments) set_tests_properties(line_attachments PROPERTIES LABELS "fortran;static;dynamic;validation") +# Bend stiffeners (STIFFENERS): element properties against closed forms and a fine integration, +# the deck grammar, and the stiffener's end on a stock-order LINES row. +add_executable(test_bend_stiffener tests/test_bend_stiffener.f90) +target_link_libraries(test_bend_stiffener PRIVATE cabledyn_core) +file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bend_stiffener) +add_test(NAME bend_stiffener COMMAND test_bend_stiffener WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/bend_stiffener) +set_tests_properties(bend_stiffener PROPERTIES LABELS "fortran;driver;static") + # OpenFAST CompMooring=5 aggregate: one deck carrying BOTH EI=0 mooring lines AND finite-EI # Hermite cables composed behind the same FMF moving-point facade. Proves pure-mooring # equivalence to CD_FMF (bit-for-bit), mixed-deck composition + column routing, pure-cable @@ -2293,6 +2301,15 @@ if(Python3_Interpreter_FOUND) set_tests_properties(driver_touchdown_recovery PROPERTIES LABELS "python;driver;integration;static" ENVIRONMENT "${_cabledyn_python_test_environment}") + # A hang-off through a tapered bend stiffener follows the planar elastica (scipy solve_bvp) + # within 2 % of the root moment and 1 cm; skipped where scipy is not installed. + add_test(NAME driver_bend_stiffener + COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/tests/check_bend_stiffener.py + $) + set_tests_properties(driver_bend_stiffener PROPERTIES + LABELS "python;driver;integration;static;validation" + SKIP_RETURN_CODE 77 + ENVIRONMENT "${_cabledyn_python_test_environment}") endif() # A clamped finite-EI cable with a buoyant section on a Rigid6 body: the body static solve # holds the objects at their deck poses, warns once, and reports no balanced clamped end loads. @@ -2772,7 +2789,7 @@ endforeach() # spectra, stream waves, anisotropic friction, modal analysis, clamped ends on bodies, MoorDyn-C # kinematics files, bodies and rods, line failure): every shipped deck must run to completion. foreach(_feature_example IN ITEMS - lazy_wave_vessel_motion lazy_wave_vessel_rao lazy_wave_buoyancy_modules lazy_wave_buoyancy_smeared + lazy_wave_vessel_motion lazy_wave_vessel_rao lazy_wave_bend_stiffener lazy_wave_buoyancy_modules lazy_wave_buoyancy_smeared chain_range_tdp chain_torsethaugen_spread chain_two_train_sea stream_wave_shallow_chain airy_wave_shallow_chain laid_cable_cross_current chain_modes lazy_wave_modes buoy_clamped_cable moordynC_wavekin/chain_wavekin7_currents1 spar_pinned_rods mixed_body_rods_points diff --git a/VALIDATION.md b/VALIDATION.md index 99243cd..db4eb82 100644 --- a/VALIDATION.md +++ b/VALIDATION.md @@ -75,6 +75,7 @@ article's methods, are compared in | **L3-6d** | Moving-touchdown dynamic contact | OrcaFlex 11.6d installed | curv 0.48/0.58%, tension 0.74/0.90/0.80%, TDP 0.08 m, migration 2–3 m/cycle | curv <3–5%, tension <2–3%, TDP excursion bounds <1.5 m | Pass | | **L3-6e** | Out-of-plane / 3D whirl response | OrcaFlex 11.6d | 3D curv 0.52%, interior Y env ≤1.0 mm, sag ≤22 mm, tension 0.14/0.17/0.12% | curv <3%, env <5 cm, tension <2–3% | Pass | | **L3-7** | 1-hour VolturnUS-S mooring dynamic | OrcaFlex 11.6d | `FairTen` mean 0.27%, std 0.26%, min/max 0.50/0.43%, Fourier 0.27/1.14% | ≤5% | Pass | +| **L3-8** | Bend stiffener at a clamped hang-off (5 m tapered polyurethane) | OrcaFlex 11.6d (stiffener attachment) / planar elastica | curvature 0.74% of root; cable moment 0.74%, section moment 0.18%, deflection angle 0.003° | <1.5% / <2%, 0.05° | Pass | | **L3-H** | Holcombe static lazy-wave experiment | measured marker geometry | position RMSE 1.07–1.89% of suspended length; peak marker-curvature error 0.34–5.54% | end-restraint bracket | Pass | | **L3-B** | Bergdahl chain experiment, 30 orbit conditions | measured upper-end force | RMSE 1.58 N; MAPE 3.13%; Pearson 0.993; Lin concordance 0.989 | full published matrix | Pass | | **L-fig** | Along-arc curvature overlay (CableDyn / OrcaFlex / MoorDyn-F ±4×) | 3-code committed record | CableDyn peak 0.11% from OrcaFlex | <0.5% | Pass | @@ -261,6 +262,18 @@ initialisation. offset of L2-V in its mean and extremes. The force-blended `EI = 0` step moved the minimum by +7 kN (0.15% before) and the other statistics by at most 0.01%; no step of this drive needs the rotation- and snap-resolved sub-steps, so they leave it unchanged bit for bit. +- **L3-8 bend stiffener at a clamped hang-off** (`driver_bend_stiffener`): a 60 m power + cable (EI 20 kN·m²) clamped horizontal and pulled 30° off by its weight through a 5 m + polyurethane stiffener (0.5 m cylinder of 0.60 m, cone to 0.22 m, bore 0.20 m, E 50 MPa; + 16 times the cable's EI at the root), on 0.05 m elements. The static curvature over the + stiffener and 3 m beyond follows OrcaFlex 11.6d with the stiffener given as a line attachment + within **0.74% of the root curvature** (root 0.4523 against 0.4522 1/m; + [`orcaflex_bend_stiffener.py`](https://github.com/SMI-Lab-Inha/CableDyn/blob/main/validation/scripts/orcaflex_bend_stiffener.py)). + Against an independent collocation solution of the planar extensible elastica with + `EI(s)` and the sleeve's weight, the cable's own moment follows within **0.74%** of its root + value, cable and sleeve moments add up to the section moment within **0.18%**, the node + positions within **0.6 mm**, and the end force's deflection from the clamp axis within + **0.003°** (32.37°). - **L-fig curvature overlay**: `l3_lazywave_dynamic` writes each site's along-arc curvature profile. [`validation/scripts/curvature_profile_overlay.py`](https://github.com/SMI-Lab-Inha/CableDyn/blob/main/validation/scripts/curvature_profile_overlay.py) diff --git a/doc/capabilities.rst b/doc/capabilities.rst index e6fcac7..2225ca7 100644 --- a/doc/capabilities.rst +++ b/doc/capabilities.rst @@ -102,6 +102,12 @@ Line and constitutive capability - Supported - lumped at nodes of a finite-EI cable, standalone and coupled; rejected on the two-moving-end compatibility solver + * - Bend stiffeners (``STIFFENERS``) + - Supported with limits + - tapered linear-elastic sleeve over a finite-EI line end, adding bending stiffness, mass, + buoyancy and drag per element, standalone and coupled; nonlinear (moment–curvature) + materials and bend restrictors are not provided; rejected on the two-moving-end + compatibility solver * - Modal analysis (``nModes``) - Standalone only - natural frequencies and mode shapes of each line about its static equilibrium (all diff --git a/doc/concepts.rst b/doc/concepts.rst index a44f0b4..e627790 100644 --- a/doc/concepts.rst +++ b/doc/concepts.rst @@ -61,8 +61,8 @@ line type and its own mesh density: 1 buoy 50.0 16 1 bare 55.0 18 <- End B (anchor) side -A bare cable + a buoyancy stretch + a bend stiffener is **three sections of one line**, not -three lines. A single-material chain is **one section**. The unstretched length is the sum of +A bare cable + a buoyancy stretch + more bare cable is **three sections of one line**, not +three lines; a bend stiffener over its end is a ``STIFFENERS`` row on that line. A single-material chain is **one section**. The unstretched length is the sum of the section lengths; the mesh is refined per section. .. admonition:: Porting from MoorDyn diff --git a/doc/driver_format.md b/doc/driver_format.md index 7e33617..d039530 100644 --- a/doc/driver_format.md +++ b/doc/driver_format.md @@ -133,16 +133,17 @@ except that `SYROPE IC` must come **after** `LINES`, because its rows name lines exist. A conventional layout is: LINE TYPES → BODIES → EXTERNAL LOADS → ROD TYPES → RODS → POINTS → TURBINES → LINES → SYROPE IC → -SECTIONS → END CONNECTIONS / EQUIVALENT BUOYANCY / ATTACHMENTS / FAILURE / CONTROL → OPTIONS → -OUTPUTS +SECTIONS → END CONNECTIONS / EQUIVALENT BUOYANCY / ATTACHMENTS / STIFFENERS / FAILURE / CONTROL → +OPTIONS → OUTPUTS Within LINE TYPES, the header row must precede the data rows because it selects their column order. A **line is one object spanning two end points** (End A = `NodeA`, End B = `NodeB`), built from an ordered list of **sections**, each with its own line type and mesh density (the same -line-and-section arrangement OrcaFlex uses): a bare cable, a bend stiffener, and a buoyancy -stretch are sections of one line. A single-material line is one section. MoorDyn composites written +line-and-section arrangement OrcaFlex uses): a bare cable and a buoyancy stretch are sections of +one line, and a bend stiffener over its end is a [STIFFENERS](#stiffeners-optional-finite-ei-lines) +row. A single-material line is one section. MoorDyn composites written as separate lines joined at a `Connect` point are also valid. @@ -208,6 +209,7 @@ are collapsed. It must match one of the spellings below; any other name fails wi | END CONNECTIONS | `END CONNECTIONS`, `END CONNECTION` | | EQUIVALENT BUOYANCY | `EQUIVALENT BUOYANCY`, `EQUIVALENT SECTIONS`, `BUOYANCY SECTIONS` | | ATTACHMENTS | `ATTACHMENTS`, `LINE ATTACHMENTS`, `CLUMPS` | +| STIFFENERS | `STIFFENERS`, `STIFFENER`, `BEND STIFFENERS`, `BEND STIFFENER` | | FAILURE | `FAILURE`, `FAILURES` | | CONTROL | `CONTROL`, `CONTROLS` | | TURBINES | `TURBINES`, `TURBINE` | @@ -687,8 +689,9 @@ as written. | `Length` | section unstretched length [m]; > 0 and at most 1e6, with `Length/NumSegs` at least 1e-6 | | `NumSegs` | section element count, 1 to 1 000 000. A finite-EI line needs far fewer: its static solve is verified up to 20 480 elements on the 80 m reference cable (elements of about 8 mm) and does not converge at 40 960 (about 4 mm). Keep finite-EI elements longer than about `sqrt(EI/EA)`, and refine only where the curvature needs it | -A single-material line is one `SECTIONS` row; a bare cable with a buoyancy stretch and a bend -stiffener is three sections of one line. +A single-material line is one `SECTIONS` row; a bare cable with a buoyancy stretch is three sections +of one line (bare, buoyant, bare). A bend stiffener over a line end is a +[STIFFENERS](#stiffeners-optional-finite-ei-lines) row, not a section. - A standalone deck with a finite-EI section needs `dtM`/`TMax`, unless it is a mixed `EI = 0` + finite-EI deck, which may be static-only. @@ -760,7 +763,8 @@ Behaviour: with `vesselMotion` or `vesselRAO` it is stored in vessel axes (equal to global axes at the reference pose) and turns with the vessel. - The connection moment equals the bending moment of the line at its end node, so - `BendMomN1` (End A) and `BendMomN` (End B) report it. + `BendMomN1` (End A) and `BendMomN` (End B) report it; under a bend stiffener + it is `BendMom` plus `StfMom` at that node. - On a `Body` point the direction is given in the body frame, like the point offset, and turns with the body. On a rod end it must be parallel to the rod axis (either sense, in the frame of the rod's deck coordinates); another direction is rejected by name. The connection @@ -931,6 +935,64 @@ Behaviour: attachments at theirs; a checkpoint restart is bit-identical. They are rejected on the two-moving-end compatibility route. +### STIFFENERS (optional, finite-EI lines) + +A bend stiffener: a tapered polyurethane sleeve over a line end, usually at a hang-off, that +spreads the bending of the cable over its length. The sleeve is not a separate object; it adds +its bending stiffness, mass, buoyancy and drag to the elements it covers. Each row is one +conical segment, and the rows of one line end follow each other outward from the stiffener +root, each starting where the previous one ends. A cylindrical root, a cone and a tip are three +rows. + +```text +--------------------- STIFFENERS ------------------------------------ +LineID End Length RootOD TipOD InnerD E Density +(-) (-) (m) (m) (m) (m) (Pa) (kg/m^3) +1 A 0.4 0.45 0.45 0.17 5.0e7 1150 +1 A 3.6 0.45 0.18 0.17 5.0e7 1150 +``` + +| Column | Meaning | +|--------|---------| +| `LineID` | existing `LINES` id of a finite-EI line | +| `End` | the end the stiffener sleeves: `A` or `B` (also `EndA`, `EndB`) | +| `Length` | length of the segment along the line [m] | +| `RootOD` | outer diameter at the segment's root side [m] | +| `TipOD` | outer diameter at the segment's tip side [m]; the diameter varies linearly in between | +| `InnerD` | inner (bore) diameter [m], smaller than both outer diameters | +| `E` | Young's modulus of the sleeve material [Pa] | +| `Density` | density of the sleeve material [kg/m³] | + +Behaviour: + +- Over the stiffener the section bends as cable and sleeve together: + `EI(s) = EI_cable + E π/64 (OD(s)⁴ − InnerD⁴)`. The sleeve carries no axial load and no + torque, so `EA` and `GJ` are the cable's. Each element takes the compliance average of + `EI(s)` over its length, `L / ∫ ds/EI(s)`, the cable mass plus the mean sleeve mass, a + displaced area of cable plus solid sleeve, and the mean outer diameter for drag and added mass + (see [theory](theory.rst), *Bend stiffeners*). Elements the stiffener does not reach keep + their line type's values exactly. +- The element properties are piecewise constant, one value per deck element, so the deck mesh + must resolve the taper: an initialisation note names a stiffener that spans fewer than 10 + elements and the element length that resolves it. Split the section at the stiffener tip and + mesh the part under it finely (about 0.1–0.2 m elements for a typical 3–8 m stiffener). +- The stiffener root takes the moment only where the end is held in rotation: declare the end + `Rigid`, or a rotational spring, in [END CONNECTIONS](#end-connections-optional-finite-ei-lines). + A stiffener at a pinned end prints a note. +- Outputs: `CurvN` is the curvature cable and sleeve share; `BendMomN` the cable's + own moment (the line type's `EI` × curvature, the value to check against the cable's limits); + `StfMomN` the sleeve's share, so that the two add up to the section moment the root + connection carries. `FairDefl`/`AnchDefl` give the angle between the end force and the + held direction, which with `FairTen`/`AnchTen` is the stiffener's load case. Keep the fine + mesh a couple of metres past the stiffener tip, where the cable's curvature can peak. +- The material is linear elastic. A row fails if the line is not finite-EI, a value is + non-finite, above 1e15 or out of range (`Length`, `RootOD`, `TipOD` and `E` must be positive, + `InnerD` and `Density` non-negative, `InnerD` below both outer diameters), or the stiffeners + of End A and End B together are longer than the line. On a stock-order line (anchor as + `NodeA`) the `End` column follows the written row, so the stiffener stays on the same physical + end. Stiffeners run wherever the cubic-Hermite cable does, standalone and coupled; they are + rejected on the two-moving-end compatibility route. + ### FAILURE (optional, EI = 0 decks) A FAILURE row detaches line ends from a point during a dynamic run, with MoorDyn's trigger rules. @@ -1295,13 +1357,15 @@ bending moment. None of these channels is mandatory. | `FairDecl` | line L fairlead (End A) declination from +GZ (0 = up, 90 = horizontal, 180 = down) | deg | | `AnchDecl` | line L anchor (End B) declination from +GZ (0 = up, 90 = horizontal, 180 = down) | deg | | `FairAngle` / `AnchAngle` | aliases of `FairDecl` / `AnchDecl` | deg | +| `FairDefl` / `AnchDefl` | line L deflection of the End-A / End-B end force from the direction a `Rigid` or rotational-spring END CONNECTIONS row holds there (0 when the line leaves along it); only at such an end | deg | | `Point

p{x,y,z}` | point P position component | m | | `Con

p{x,y,z}` | alias of `Point

p{x,y,z}` (MoorDyn v1 spelling) | m | | `Body…`, `Rod…` | Rigid6 body and rod channels in the MoorDyn-F names: position, attitude, velocity and acceleration (`P`, `R`, `V`, `RV`, `A`, `RA`), net force and moment (`F`, `M`), rod end tensions (`TenA`, `TenB`), submerged fraction (`Sub`) and rod node positions (`RodNP`). Available on every route that carries bodies or rods, OpenFAST included; the full list and units are in [Outputs](outputs.rst) | see Outputs | | `Point

F{x,y,z}`, `Point

FH` | resultant of the line-end and cable end forces on point P (each the `FairTen`/`AnchTen` end force): the anchor load of a Fixed point shared by several lines, the line load a free point balances; `FH` is its horizontal magnitude. Available on every route (static, `EI = 0`, finite-EI, multibody and OpenFAST) except the two-moving-end finite-EI route, which rejects it by name | N | | `TenN` | line L effective tension at node J. EI = 0 lines: the mean of the two adjacent segment tensions. Finite-EI lines: the axial component of the section force recovered from nodal equilibrium (`OutputRecovery element`: the length-weighted mean of the adjacent element-mean axial forces). End nodes report the line-end force of `FairTen`/`AnchTen` | N | | `CurvN` | line L curvature at node J. `EI = 0` lines: geometric curvature of the circle through the node and its two neighbours (end nodes take the adjacent interior value). Cubic-Hermite lines: `BendMomN` / EI, the larger of the two sides at a section junction (`OutputRecovery element`: exact curvature of the continuous centreline, the larger one-sided element value at an interior node) | 1/m | -| `BendMomN` | line L bend moment = EI × curvature **relative to the stress-free reference shape** (≈ 0 at the reference; the finite-EI reference is straight). Finite-EI lines: recovered from nodal equilibrium, the moment the neighbouring elements carry (`OutputRecovery element`: EI × the larger one-sided exact curvature); 0 on EI = 0 lines | N·m | +| `BendMomN` | line L bend moment = EI × curvature **relative to the stress-free reference shape** (≈ 0 at the reference; the finite-EI reference is straight). Finite-EI lines: recovered from nodal equilibrium, the moment the neighbouring elements carry (`OutputRecovery element`: EI × the larger one-sided exact curvature); 0 on EI = 0 lines. Under a bend stiffener: the cable's own share, the line type's EI × curvature | N·m | +| `StfMomN` | line L share of the section moment a bend stiffener carries at node J (the section moment less `BendMomN`; 0 outside the stiffener); only on a line with STIFFENERS | N·m | | `LNp{x,y,z}` | line L node J position component | m | | `LNv{x,y,z}` | line L node J velocity component (0 in static-only runs) | m/s | | `LNa{x,y,z}` | line L node J acceleration component (0 in static-only runs) | m/s² | diff --git a/doc/examples.rst b/doc/examples.rst index b3a6f5d..342e0d7 100644 --- a/doc/examples.rst +++ b/doc/examples.rst @@ -335,6 +335,12 @@ touchdown and lazy-wave shape from the geometry alone; there is no initial-shape cable — :doc:`tutorial_motion` - Hermite, dynamic - ~9 s + * - ``lazy_wave_bend_stiffener.dat`` + - the RAO-driven cable with a 4 m polyurethane bend stiffener (a cylindrical root and a cone, + ``STIFFENERS``) at its clamped hang-off; ``BendMom1N1`` is the stiffener root moment and the + curvature channels sample the stiffener and the cable past its tip — :doc:`modeling` + - Hermite, dynamic + - ~9 s * - ``lazy_wave_buoyancy_modules.dat`` - the 80 m cable in a 0.5 m/s current with its buoyant stretch built from ten discrete buoyancy modules (``ATTACHMENTS`` series ``first:pitch:last``) — :doc:`tutorial_lazywave` diff --git a/doc/migrating.rst b/doc/migrating.rst index 0474872..de3e5fc 100644 --- a/doc/migrating.rst +++ b/doc/migrating.rst @@ -62,8 +62,8 @@ When to use SECTIONS instead ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The 7-column row can describe only a single-type line. A composite (chain–polyester–chain, or a -power cable with bare, buoyancy, and bend-stiffener zones) is, in MoorDyn, several lines joined at -``Connect`` points; in CableDyn it is **one line with several ordered sections**: +power cable with bare and buoyancy zones) is, in MoorDyn, several lines joined at ``Connect`` +points; in CableDyn it is **one line with several ordered sections**: .. code-block:: text diff --git a/doc/modeling.rst b/doc/modeling.rst index 25eb2eb..46d7e51 100644 --- a/doc/modeling.rst +++ b/doc/modeling.rst @@ -78,6 +78,37 @@ viscoelastic model is a four-parameter series-Kelvin solid that reduces to a sta solid only when ``BA_D = 0``). Preserve the units shown in the ``LINE TYPES`` table; the parser does not infer or convert vendor units. +.. _modeling-bend-stiffeners: + +**Bend stiffeners.** Give a bend stiffener as a ``STIFFENERS`` table on the cable, not as a +section of its own: one row per conical segment from the root outward (outer diameter at each +side, bore, modulus and density of the polyurethane), at End A or End B. The driver adds the +sleeve's bending stiffness, mass, buoyancy and drag to the elements it covers +(:ref:`theory-bend-stiffeners`). Three choices decide the result: + +- **the root connection.** The stiffener root is clamped to the hang-off or I-tube: declare that + end ``Rigid`` (or a rotational spring for a compliant connection) in ``END CONNECTIONS``, along + the hang-off direction. A pinned end carries no moment, so the stiffener then takes none; +- **the mesh under the stiffener.** One stiffness per element: split the section at the + stiffener tip and mesh the part under it with elements of about 0.1–0.2 m, at least ten over + the stiffener (the driver notes fewer); +- **the material.** The sleeve is linear elastic. Polyurethane stiffness depends on temperature + and strain rate; use the modulus for the design temperature and check the result at the + ends of its range. + +The design checks are the curvature of the cable against its minimum bend radius along the +stiffener and just past its tip (``CurvN`` channels and the range graphs; keep the fine +mesh a couple of metres beyond the tip), the cable's own moment (``BendMomN``), and the +stiffener's load case, the end force with its deflection from the hang-off direction +(``FairTen`` with ``FairDefl``). ``StfMomN`` is the sleeve's share of the moment; +at the root, ``BendMomN1`` plus ``StfMomN1`` is the moment the hang-off takes. + +A sleeve that can slide or rattle on the cable, a steel root section and friction between +sleeve and cable are not modelled: the sleeve is bonded to the cable in bending. + +``lazy_wave_bend_stiffener.dat`` shows a 4 m stiffener at a clamped hang-off on a vessel in +waves. + 5. Choose and converge the mesh ------------------------------- @@ -125,9 +156,10 @@ two- or four-fold local refinement. When end curvature matters: - report the curvature 1–2 m from a pinned end, never the value at the pinned node: a pinned end carries no moment, so the exact curvature there is zero, and the computed nodal value is a discretisation residue (at ``dtM = 0.05 s`` up to 37 times the converged value 1 m away); -- for bend-stiffener or hang-off design, model the stiffener as a section and the end as a - clamped or rotational-spring ``END CONNECTIONS`` row (see :doc:`driver_format`); a pinned end - does not represent the moment a stiffener carries. +- for bend-stiffener or hang-off design, give the stiffener as a ``STIFFENERS`` table (see + :ref:`bend stiffeners `) and the end as a clamped or rotational-spring + ``END CONNECTIONS`` row (see :doc:`driver_format`); a pinned end does not represent the moment + a stiffener carries. **Snap loads on slack lines.** When a slack line snaps taut it excites axial waves up to the segment scale, near :math:`2c/l` with :math:`c = \sqrt{EA/m}` and :math:`l` the segment length. diff --git a/doc/outputs.rst b/doc/outputs.rst index fc24a61..cfe4e6f 100644 --- a/doc/outputs.rst +++ b/doc/outputs.rst @@ -223,8 +223,8 @@ finite-EI line of the run: - nodal curvature (same definition as ``CurvN``) * - ``BendMoment`` - ``(N.m)`` - - bend moment (same definition as ``BendMomN``); written as 0 by the ``EI = 0`` - static route + - bend moment (same definition as ``BendMomN``: the cable's own share under a bend + stiffener); written as 0 by the ``EI = 0`` static route * - ``Declination`` - ``(deg)`` - axial-tangent angle from +Z (0 = up, 90 = horizontal, 180 = down) @@ -301,7 +301,8 @@ B)``, header row, units row, data (12 columns): - maximum curvature over the element * - ``BendMomentAtPeak`` - ``(N.m)`` - - ``EI × PeakCurvature`` + - ``EI × PeakCurvature``, with the line type's ``EI`` (the cable's share under a bend + stiffener) * - ``MinimumAxialResultant`` - ``(N)`` - minimum axial force resultant over the element (negative = compression) @@ -558,6 +559,13 @@ Line-end channels * - ``FairAngle`` / ``AnchAngle`` - aliases of ``FairDecl`` / ``AnchDecl`` - deg + * - ``FairDefl`` / ``AnchDefl`` + - deflection of the End-A / End-B end force (``FairTen`` / ``AnchTen``) from the direction + a ``Rigid`` or rotational-spring ``END CONNECTIONS`` row holds at that end, turning with + the body or vessel it is on: 0 when the line leaves the connection along that direction. + With the end force, the load case of a bend stiffener or hang-off design. Only at an end + with such a row + - deg A grounded anchor segment reads ``AnchIncl ≈ 0``; the small non-zero value reports the actual orientation of the last element, which CableDyn does not force to horizontal. @@ -648,7 +656,14 @@ Line-node channels of the moment recovered from nodal equilibrium, the mean of the two neighbouring elements (one-sided at a line end, where it is the end connection moment); ``OutputRecovery element``: ``EI`` times the larger one-sided exact curvature. 0 on - ``EI = 0`` lines + ``EI = 0`` lines. Under a bend stiffener (``STIFFENERS``) it is the cable's own share, + the line type's ``EI`` times the curvature; the stiffener's share is ``StfMomN`` + and the two add up to the section moment + - N·m + * - ``StfMomN`` + - the share of the section moment a bend stiffener carries at node ``J``, the section + moment less ``BendMomN``; 0 where no stiffener covers the node. Only on a line with + ``STIFFENERS`` - N·m * - ``LNpx`` / ``py`` / ``pz`` - node position component diff --git a/doc/references.rst b/doc/references.rst index 0fc43e6..371c264 100644 --- a/doc/references.rst +++ b/doc/references.rst @@ -28,6 +28,12 @@ Element formulation element formulation for geometrically exact curved Kirchhoff rods. *Computer Methods in Applied Mechanics and Engineering* **278**, 445–478. https://doi.org/10.1016/j.cma.2014.05.017 +.. _ref-boef1990: + +- **Boef, W. J. C. & Out, J. M. M. (1990).** Analysis of a flexible-riser top connection with + bend restrictor. *Offshore Technology Conference*, Houston, OTC-6436-MS. + https://doi.org/10.4043/6436-MS + .. _ref-bergou2008: - **Bergou, M., Wardetzky, M., Robinson, S., Audoly, B. & Grinspun, E. (2008).** Discrete diff --git a/doc/theory.rst b/doc/theory.rst index 312d462..5ae0426 100644 --- a/doc/theory.rst +++ b/doc/theory.rst @@ -140,6 +140,59 @@ The element follows the geometrically exact Kirchhoff-rod formulation of interpolation of :ref:`Meier, Popp & Wall (2015) `, reduced to a torsion-free cable. Torsion enters only through the condensed line term below. +.. _theory-bend-stiffeners: + +Bend stiffeners +~~~~~~~~~~~~~~~ + +A bend stiffener (``STIFFENERS``) is a sleeve of outer diameter :math:`D(s)`, bore :math:`d_i`, +modulus :math:`E_s` and density :math:`\rho_s` over a line end, made of conical segments in +which :math:`D` varies linearly with the arc length. Cable and sleeve bend together, as in the +stiffener design model of :ref:`Boef & Out (1990) `, so over the stiffener + +.. math:: + + EI(s) = EI_c + E_s \frac{\pi}{64}\left(D(s)^4 - d_i^4\right), + +with :math:`EI_c` the cable's bending stiffness. The sleeve carries no axial load and no torque +(:math:`EA` and :math:`GJ` are the cable's), and its material is linear elastic. The element +energy above takes one :math:`EI` per element, so an element over the arc :math:`[a, b]` of +length :math:`L` takes the compliance average + +.. math:: + + EI_e = L \left( \int_a^b \frac{\mathrm{d}s}{EI(s)} \right)^{-1}, + +which gives the exact relative end rotation of the element under a uniform moment: the stiff +root of a tapered sleeve is not allowed to stiffen a whole element that only partly lies on it. +The element mass per length is the cable's plus :math:`\rho_s \tfrac{\pi}{4}\, +\overline{D^2 - d_i^2}`; buoyancy acts on the displaced area of cable and solid sleeve, +:math:`\tfrac{\pi}{4}\,(d^2 + \overline{D^2 - d_i^2})`, and drag and added mass on the mean outer +diameter :math:`\overline{\max(D, d)}`, where the bar is the mean over the element and +:math:`d` the cable diameter. The integrals use six-point Gauss–Legendre quadrature on each +part of the element a segment covers. With these one-per-element values the taper is resolved +by the mesh, as in a segmented line model; the driver notes a stiffener that spans fewer than +ten elements. + +The recovered section moment :math:`M` (:ref:`theory-output-recovery`) is shared in proportion +to the stiffnesses, since cable and sleeve take the same curvature :math:`\kappa = M/EI_e`: the +cable carries :math:`EI_c\,\kappa`, reported as ``BendMom``, and the sleeve +:math:`(EI_e - EI_c)\,\kappa`, reported as ``StfMom``. At a node between a covered and a bare +element the side of the larger curvature, the bare one, gives the values. The deflection angle +``FairDefl``/``AnchDefl`` is :math:`\arccos(\hat{\mathbf F}\cdot\hat{\mathbf d})` between the +end force :math:`\mathbf F` and the direction :math:`\mathbf d` the end connection holds (into +the line), the angle through which the stiffener turns the cable. + +The planar check ``driver_bend_stiffener`` compares a hang-off clamped through a 5 m stiffener +(a cylindrical root and a cone, 16 times the cable's stiffness at the root) with an independent +collocation solution of the extensible elastica +:math:`\theta' = M/EI(s)`, :math:`M' = -(1 + \varepsilon)(\cos\theta\,n_z - \sin\theta\,n_x)`, +:math:`\mathbf n' = (0, w(s))`: on 0.05 m elements the cable's moment over the stiffener and +3 m beyond follows the elastica within 0.8 % of its root value, cable and sleeve moments add up +to the section moment within 0.2 % at the root, the deflection angle agrees within 0.003°, and +the node positions within 1 mm. On the same hang-off, the curvature along the stiffener follows +OrcaFlex 11.6d with its stiffener given as a line attachment within 0.8 % of the root curvature. + .. _theory-output-recovery: Node tension, curvature and bend moment diff --git a/examples/README.md b/examples/README.md index cacac9a..416de1b 100644 --- a/examples/README.md +++ b/examples/README.md @@ -45,7 +45,7 @@ Copy a deck under a new name and keep every file it references at the same relat | `openfast/` | IEA-15MW VolturnUS-S `.fst` pair (`CompMooring = 5` and `3`) with their mooring decks; see [openfast/README.md](openfast/README.md). | | `data/syrope/` | Syrope settings and original working-curve table used by `syrope_polyester_mooring.dat`. | | `data/lozon/` | Prescribed hang-off heave history for `lozon_gomex80_power_cable_motion.dat` and the script that generates it. | -| `data/vessel/` | 6-DOF vessel-motion record for `lazy_wave_vessel_motion.dat`, the script that generates it, and the illustrative RAO table for `lazy_wave_vessel_rao.dat`. | +| `data/vessel/` | 6-DOF vessel-motion record for `lazy_wave_vessel_motion.dat`, the script that generates it, and the illustrative RAO table for `lazy_wave_vessel_rao.dat` and `lazy_wave_bend_stiffener.dat`. | | `data/torsion/` | Hang-off roll history (motionFile roll column) for `torsion_lazy_wave_hangoff_twist.dat` and the script that generates it. | | `data/range_tdp/` | Fairlead surge history for `chain_range_tdp.dat` and the script that generates it. | | `moordynC_wavekin/` | A deck with the MoorDyn-C `wave_frequencies.txt` and `current_profile.txt` files it reads from its own folder. | diff --git a/examples/lazy_wave_bend_stiffener.dat b/examples/lazy_wave_bend_stiffener.dat new file mode 100644 index 0000000..86e9c99 --- /dev/null +++ b/examples/lazy_wave_bend_stiffener.dat @@ -0,0 +1,72 @@ +--------------------- CableDyn Input File ------------------------------------ +Lazy-wave power cable through a bend stiffener at a clamped hang-off on a vessel in a JONSWAP sea +-- The cable, vessel motion and sea of lazy_wave_vessel_rao.dat, with a 4 m polyurethane bend +-- stiffener at the hang-off: a 0.4 m cylindrical root of 0.45 m outer diameter and a 3.6 m cone +-- to 0.18 m, bored to 0.17 m (E = 50 MPa, 1150 kg/m^3). At the root the sleeve is about five +-- times as stiff in bending as the cable (EI 9.9e4 against 1.99e4 N m^2). The root is clamped +-- (Rigid) along the static hang-off tangent; BendMom1N1 is the root moment, the curvature +-- channels sample the stiffener (nodes 14 and 27, 2 m and 4 m from the root) and the bare cable +-- just past its tip (node 34). The 0.152 m elements give 27 elements over the stiffener. +--------------------- LINE TYPES --------------------------------------------- +TypeName Diam MassDenInAir EA BA EI Cd_n Cd_t Ca_n Ca_t +(-) (m) (kg/m) (N) (-) (N-m2) (-) (-) (-) (-) +bare 0.160 36.70 4.69e8 0.0 1.99e4 1.2 0.1 1.0 0.0 +buoy 0.290 59.53 4.69e8 0.0 1.99e4 1.2 0.1 1.0 0.0 +--------------------- POINTS ------------------------------------------------ +ID Type X Y Z +(-) (-) (m) (m) (m) +1 Coupled 5.0 0.0 -14.0 +2 Fixed 125.0 0.0 -80.0 +--------------------- LINES ------------------------------------------------- +ID NodeA NodeB Outputs +(-) (-) (-) (-) +1 1 2 r +--------------------- SECTIONS ---------------------------------------------- +LineID LineType Length NumSegs +(-) (-) (m) (-) +1 bare 68.114 448 +1 buoy 50.000 320 +1 bare 52.101 256 +--------------------- STIFFENERS -------------------------------------------- +-- Rows run outward from the root at End A: the cylindrical root, then the cone. +LineID End Length RootOD TipOD InnerD E Density +(-) (-) (m) (m) (m) (m) (Pa) (kg/m^3) +1 A 0.4 0.45 0.45 0.17 5.0e7 1150 +1 A 3.6 0.45 0.18 0.17 5.0e7 1150 +--------------------- END CONNECTIONS --------------------------------------- +-- End A is clamped along the static hang-off tangent (81.6 deg below horizontal, toward End B). +LineID End Stiffness EzX EzY EzZ +(-) (-) (N-m/rad) (-) (-) (-) +1 A Rigid 0.1461 0.0 -0.9893 +--------------------- OPTIONS ----------------------------------------------- +9.80665 g - Gravitational acceleration (m/s^2) [default: 9.80665] +1025.0 rhoW - Water density (kg/m^3) [default: 1025] +80.0 WtrDpth - Water depth for the flat seabed (m) [default: absent standalone; host-owned OpenFAST] +1.0e5 kBot - Seabed penalty stiffness base (Pa/m) [default: 1.0e5] +1.0e4 cBot - Seabed normal damping base (Pa-s/m) [default: 1.0e4] +0.4 rhoInf - Generalised-alpha spectral radius {0<=rhoInf<=1} [default: 0.4] +False modified_newton - Tangent-reuse mode {False=full Newton; True=guarded reuse} [default: False] +0.05 dtM - CableDyn internal time step (s) [default: omitted standalone; the glue DT in OpenFAST] +40.0 TMax - Standalone simulation duration (s), including the 10 s wave ramp +data/vessel/sample_rao.txt vesselRAO - Displacement RAO table: vessel motion from the deck waves +0.0|0.0|0.0 vesselRef - Vessel reference point, rotation centre and RAO origin (m) [default: 0|0|0] +10.0 RangeStart - First output time in the range graphs (s) [default: 0] +False adaptive_mesh - Automatic finite-EI refinement {False; True} [default: False] +True tensile_safety - Refine and reject local axial compression in this tensile cable +dynamic_solver 1.0e-4 1.0e-14 100 12 - Relative tolerance, absolute tolerance, max iterations, backtracks +none frictionMu - Seabed friction coefficient {none; value>=0} [default: none] +none current - Current model and velocity (m/s) {none; uniform; profile} [default: none] +jonswap 4.0 10.0 3.3 0.0 waves - Wave model, JONSWAP Hs, Tp, gamma, and direction (m, s, deg) {none; airy; jonswap} +1 WaveSeed - Seed of the random sea [default: 1] +10.0 rampTime - Half-cosine start-up ramp of the wave amplitudes (s) [default: 0] +0 WaterKin - Water-kinematics source {0/none; file; SEASTATE coupled only} [default: no policy] +-- Complete option/default reference: cabledyn_options_reference.dat and doc/options.rst. +--------------------- OUTPUTS ----------------------------------------------- +"FairTen1" +"AnchTen1" +"BendMom1N1" +"Curv1N14" +"Curv1N27" +"Curv1N34" +"TDP1s" +--------------------- need this line ---------------------------------------- diff --git a/python/cabledyn/builder.py b/python/cabledyn/builder.py index d7a7463..8bf7c1a 100644 --- a/python/cabledyn/builder.py +++ b/python/cabledyn/builder.py @@ -84,6 +84,7 @@ "RodEnd", "RodType", "Section", + "Stiffener", "SyropeEA", "SyropeIC", "Turbine", @@ -104,13 +105,14 @@ ( "line", re.compile( - r"(fairten|anchten|anchfz|anchfh|fairangle|anchangle|fairdecl|anchdecl|fairincl|anchincl)" + r"(fairten|anchten|anchfz|anchfh|fairangle|anchangle|fairdecl|anchdecl|fairincl|anchincl" + r"|fairdefl|anchdefl)" r"([0-9]+)()", re.IGNORECASE, ), ), ("line", re.compile(r"(tdp)([0-9]+)(s|x|y|z|lay|exc|grd)", re.IGNORECASE)), - ("line", re.compile(r"(ten|curv|bendmom|torq)([0-9]+)(n[0-9]+)", re.IGNORECASE)), + ("line", re.compile(r"(ten|curv|bendmom|stfmom|torq)([0-9]+)(n[0-9]+)", re.IGNORECASE)), ("line", re.compile(r"(twist)([0-9]+)(n[0-9]+|)", re.IGNORECASE)), ("line", re.compile(r"(l)([0-9]+)(n[0-9]+(?:[pva][xyz]|dec|azi))", re.IGNORECASE)), ("point", re.compile(r"(point|con)([0-9]+)(p[xyz])", re.IGNORECASE)), @@ -615,6 +617,43 @@ class Attachment: cdax: float | None = None +@dataclass(eq=False) +class Stiffener: + """One ``STIFFENERS`` row: a conical segment of a bend stiffener on a finite-EI line. + + The rows of one line end follow each other outward from the stiffener root, each starting + where the previous one ends; a cylindrical root, a cone and a tip are three rows. + + Attributes + ---------- + line : Line + The finite-EI line. + end : str + The end the stiffener sleeves, ``A`` or ``B``. + length : float + Length of the segment, in m. + root_od : float + Outer diameter at the segment's root side, in m. + tip_od : float + Outer diameter at the segment's tip side, in m. + inner_d : float + Inner diameter, in m. + modulus : float + Young's modulus of the stiffener material, in Pa. + density : float + Density of the stiffener material, in kg/m^3. + """ + + line: Line + end: str + length: float + root_od: float + tip_od: float + inner_d: float + modulus: float + density: float + + @dataclass(eq=False) class SyropeIC: """One ``SYROPE IC`` row: the prior load history of Syrope lines. @@ -1257,6 +1296,7 @@ def __init__( self._end_connections: list[EndConnection] = [] self._equivalent_buoyancy: list[EquivalentBuoyancy] = [] self._attachments: list[Attachment] = [] + self._stiffeners: list[Stiffener] = [] self._syrope_ic: list[SyropeIC] = [] self._failures: list[Failure] = [] self._controls: list[Control] = [] @@ -1426,6 +1466,11 @@ def attachments(self) -> tuple[Attachment, ...]: """``ATTACHMENTS`` rows.""" return tuple(self._attachments) + @property + def stiffeners(self) -> tuple[Stiffener, ...]: + """``STIFFENERS`` rows.""" + return tuple(self._stiffeners) + @property def syrope_ic(self) -> tuple[SyropeIC, ...]: """``SYROPE IC`` rows.""" @@ -2005,6 +2050,28 @@ def add_attachment( self._attachments.append(item) return item + def add_stiffener( + self, + line: Line | int, + end: str, + length: float, + root_od: float, + tip_od: float, + inner_d: float, + modulus: float, + density: float, + ) -> Stiffener: + """Add a ``STIFFENERS`` row (one conical segment); see :class:`Stiffener`. + + Returns + ------- + Stiffener + The new row. + """ + item = Stiffener(self._line(line), end, length, root_od, tip_od, inner_d, modulus, density) + self._stiffeners.append(item) + return item + def add_syrope_ic( self, lines: Line | int | Iterable[Line | int], tmax0: float, tmean0: float ) -> SyropeIC: @@ -2151,6 +2218,7 @@ def _registry(self, obj: object) -> list[Any] | None: (EndConnection, self._end_connections), (EquivalentBuoyancy, self._equivalent_buoyancy), (Attachment, self._attachments), + (Stiffener, self._stiffeners), (SyropeIC, self._syrope_ic), (Failure, self._failures), (Control, self._controls), @@ -2234,6 +2302,7 @@ def references(self, obj: object) -> tuple[object, ...]: elif isinstance(obj, Line): found.extend(row for row in self._end_connections if row.line is obj) found.extend(row for row in self._attachments if row.line is obj) + found.extend(row for row in self._stiffeners if row.line is obj) for rows in (self._syrope_ic, self._failures, self._controls): found.extend(row for row in rows if any(line is obj for line in row.lines)) found.extend(self._line_options(obj.id)) @@ -2397,7 +2466,11 @@ def _check_references(self) -> None: self._member(point.body, self._bodies, "body") if point.rod_end is not None: self._member(point.rod_end.rod, self._rods, "rod") - owners: list[EndConnection | Attachment] = [*self._end_connections, *self._attachments] + owners: list[EndConnection | Attachment | Stiffener] = [ + *self._end_connections, + *self._attachments, + *self._stiffeners, + ] for owned in owners: self._member(owned.line, self._lines, "line") for eq in self._equivalent_buoyancy: @@ -2453,6 +2526,7 @@ def to_text(self, *, validate: bool = True) -> str: _emit(out, "END CONNECTIONS", *self._end_connection_table()) _emit(out, "EQUIVALENT BUOYANCY", *self._equivalent_table()) _emit(out, "ATTACHMENTS", *self._attachment_table()) + _emit(out, "STIFFENERS", *self._stiffener_table()) _emit(out, "FAILURE", *self._failure_table()) _emit(out, "CONTROL", *self._control_table()) _emit(out, "EXTERNAL LOADS", *self._external_load_table()) @@ -2785,6 +2859,16 @@ def _attachment_table(self) -> tuple[str, str, list[list[str]]]: header = "LineID ArcLength Mass Volume CdA Ca CdAx" return header, "(-) (m) (kg) (m^3) (m^2) (-) (m^2)", rows + def _stiffener_table(self) -> tuple[str, str, list[list[str]]]: + rows: list[list[str]] = [] + for row in self._stiffeners: + what = f"line {row.line.id} stiffener" + values = (row.length, row.root_od, row.tip_od, row.inner_d, row.modulus, row.density) + cells = [_int(row.line.id, what), _end_letter(row.end)] + rows.append([*cells, *(_num(value, what) for value in values)]) + header = "LineID End Length RootOD TipOD InnerD E Density" + return header, "(-) (-) (m) (m) (m) (m) (Pa) (kg/m^3)", rows + def _syrope_table(self) -> tuple[str, str, list[list[str]]]: rows = [ [_line_list(row.lines), _num(row.tmax0, "SYROPE IC"), _num(row.tmean0, "SYROPE IC")] @@ -2980,6 +3064,9 @@ def load(self) -> None: for row in self.rows("EQUIVALENT BUOYANCY"): t = row.tokens model.add_equivalent_buoyancy(t[0], _f(t[1]), _f(t[2])) + for row in self.rows("STIFFENERS"): + t = row.tokens + model.add_stiffener(_native_int(t[0]), t[1][-1], *(_f(token) for token in t[2:8])) for row in self.rows("ATTACHMENTS"): t = row.tokens arc: float | tuple[float, float, float] diff --git a/python/cabledyn/deck_file.py b/python/cabledyn/deck_file.py index 3ff8dd6..9f7f1ea 100644 --- a/python/cabledyn/deck_file.py +++ b/python/cabledyn/deck_file.py @@ -70,6 +70,10 @@ "ATTACHMENTS": "ATTACHMENTS", "LINE ATTACHMENTS": "ATTACHMENTS", "CLUMPS": "ATTACHMENTS", + "STIFFENERS": "STIFFENERS", + "STIFFENER": "STIFFENERS", + "BEND STIFFENERS": "STIFFENERS", + "BEND STIFFENER": "STIFFENERS", } # Table sections whose rows the native reader takes through a Fortran # list-directed READ. Their descriptive and units rows are recognized by the @@ -91,7 +95,9 @@ ) # Table sections the native reader splits on whitespace and parses column by # column (their comma-separated line lists are valid there). -_WHITESPACE_TABLE_SECTIONS = frozenset({"FAILURE", "CONTROL", "SYROPE IC", "ATTACHMENTS"}) +_WHITESPACE_TABLE_SECTIONS = frozenset( + {"FAILURE", "CONTROL", "SYROPE IC", "ATTACHMENTS", "STIFFENERS"} +) _TABLE_SECTIONS = _LIST_DIRECTED_SECTIONS | _WHITESPACE_TABLE_SECTIONS _NATIVE_RECORD_LENGTH = 512 _INT32_MAX = 2**31 - 1 @@ -691,7 +697,14 @@ def _parse_line_node_channel(lowered: str) -> tuple[int, int, int, int]: # Twist: the line's total twist (node 1) line_id = _native_positive_int(lowered, 6, last) return (line_id, 1, 11, 0) if line_id > 0 else (line_id, 0, 0, 0) - for keyword, code in (("torq", 9), ("twist", 10), ("ten", 3), ("curv", 5), ("bendmom", 6)): + for keyword, code in ( + ("torq", 9), + ("twist", 10), + ("ten", 3), + ("curv", 5), + ("bendmom", 6), + ("stfmom", 12), + ): if lowered.startswith(keyword): prefix = len(keyword) separator = lowered.find("n", prefix) @@ -745,6 +758,9 @@ def _channel_identity_key(channel: str) -> str: if lowered.startswith(("fairincl", "anchincl")): return f"{lowered[:8]}:{line_id}" return f"{lowered[:4]}decl:{line_id}" + if lowered.startswith(("fairdefl", "anchdefl")): + line_id = _native_id_suffix(lowered[8:]) + return f"{lowered[:8]}:{line_id}" if line_id > 0 else lowered if lowered.startswith("point") or ( lowered.startswith("con") and lowered[3:4].isascii() and lowered[3:4].isdigit() ): @@ -761,7 +777,7 @@ def _channel_identity_key(channel: str) -> str: if touchdown is not None and int(touchdown.group(1)) > 0: component = ("s", "x", "y", "z", "lay", "exc", "grd").index(touchdown.group(2)) + 1 return f"tdp:{int(touchdown.group(1))}:{component}" - if lowered.startswith(("ten", "curv", "bendmom", "l", "torq", "twist")): + if lowered.startswith(("ten", "curv", "bendmom", "stfmom", "l", "torq", "twist")): line_id, node_id, kind, component = _parse_line_node_channel(lowered) if line_id > 0 and node_id > 0 and kind > 0: return f"node:{line_id}:{node_id}:{kind}:{component}" @@ -2034,6 +2050,83 @@ def _check_attachments( f"(unstretched length {line_lengths[line_id]:.10g} m)" ) + def _check_stiffeners( + self, + line_ids: set[str], + section_types_by_line: dict[str, list[str]], + line_type_meta: dict[str, tuple[float, float, float, float]], + line_lengths: dict[str, float], + ) -> None: + """Native ``append_line_stiffener`` and its finalize checks for STIFFENERS rows. + + A row ``LineID End Length RootOD TipOD InnerD E Density`` is one conical segment of the + bend stiffener sleeving End A or End B; the rows of one line end follow each other + outward from the root. The stiffeners of both ends together must fit on the line. + """ + extent: dict[tuple[str, str], float] = {} + for row in self._rows("STIFFENERS"): + where = self._at(row) + if len(row.tokens) != 8: + raise DeckFormatError( + f"{where}: STIFFENERS rows must be: " + "LineID End Length RootOD TipOD InnerD E Density" + ) + try: + line_id = str(_native_int(row.tokens[0])) + except ValueError as exc: + raise DeckFormatError(f"{where}: STIFFENERS LineID must be an integer") from exc + if line_id not in line_ids: + raise DeckFormatError( + f"{where}: STIFFENERS row references undefined LINE id {line_id}" + ) + end = {"a": "a", "enda": "a", "b": "b", "endb": "b"}.get(row.tokens[1].lower()) + if end is None: + raise DeckFormatError(f"{where}: STIFFENERS End must be A or B") + try: + length, root_od, tip_od, inner_d, modulus, density = ( + _native_float(token) for token in row.tokens[2:] + ) + except ValueError as exc: + raise DeckFormatError( + f"{where}: STIFFENERS Length, RootOD, TipOD, InnerD, E and Density must be " + "plain numbers" + ) from exc + values = (length, root_od, tip_od, inner_d, modulus, density) + if not all( + math.isfinite(value) and value <= _NATIVE_MAX_COEFFICIENT for value in values + ): + raise DeckFormatError( + f"{where}: STIFFENERS Length, RootOD, TipOD, InnerD, E and Density must be " + "finite and at most 1e15" + ) + if not (length > 0.0 and root_od > 0.0 and tip_od > 0.0 and modulus > 0.0): + raise DeckFormatError( + f"{where}: STIFFENERS Length, RootOD, TipOD and E must be > 0" + ) + if inner_d < 0.0 or density < 0.0: + raise DeckFormatError(f"{where}: STIFFENERS InnerD and Density must be >= 0") + if not inner_d < min(root_od, tip_od): + raise DeckFormatError( + f"{where}: STIFFENERS InnerD must be smaller than RootOD and TipOD" + ) + extent[(line_id, end)] = extent.get((line_id, end), 0.0) + length + for line_id in sorted({key[0] for key in extent}, key=int): + if not any( + line_type_meta[name][3] > 0.0 + for name in section_types_by_line.get(line_id, []) + if name in line_type_meta + ): + raise DeckFormatError( + f"{self._label}: STIFFENERS name line {line_id}, which is not a finite-EI " + "line: a stiffener adds bending stiffness, which an EI = 0 line does not have" + ) + total = extent.get((line_id, "a"), 0.0) + extent.get((line_id, "b"), 0.0) + if total > line_lengths[line_id] * (1.0 + 1.0e-12): + raise DeckFormatError( + f"{self._label}: the STIFFENERS of line {line_id} are longer than the line " + "(End A and End B stiffeners together)" + ) + def _legacy_rod_type_columns(self, record: DeckRecord) -> bool: """Return whether ``record`` sits under a legacy CableDyn ROD TYPES header. @@ -2763,6 +2856,8 @@ def end_key(token: str) -> str: line_nelems: dict[str, int] = dict.fromkeys(line_ids, 0) line_lengths: dict[str, float] = dict.fromkeys(line_ids, 0.0) normalized_line_endpoints: dict[str, tuple[str, str]] = {} + # lines whose stock anchor-first row the native reader turns around (its ends swap) + swapped_lines: set[str] = set() range_lines = 0 series_lines = 0 # lines with per-line time-series files (flags p, t) for row in self.lines: @@ -2845,6 +2940,7 @@ def end_key(token: str) -> str: if type_a == "fixed" and upper_b: node_a, node_b = node_b, node_a type_a, type_b = type_b, type_a + swapped_lines.add(str(_native_int(row.tokens[0]))) normalized_line_endpoints[str(_native_int(row.tokens[0]))] = (node_a, node_b) if not has_dynamic_points: type_a_is_coupled = type_a in {"coupled", "vessel"} or _is_turbine_point_type( @@ -2978,9 +3074,12 @@ def end_key(token: str) -> str: "to run the dynamic working-curve path" ) + self._check_stiffeners(line_ids, section_types_by_line, line_type_meta, line_lengths) self._check_attachments(line_ids, section_types_by_line, line_type_meta, line_lengths) assigned_end_connections: set[tuple[str, str]] = set() + # (line, end) held by a Rigid or rotational-spring row, in native End A/B order + held_ends: set[tuple[str, str]] = set() # line id -> {end: torsionally restrained} from the optional torsion columns torsion_ends: dict[str, dict[str, bool]] = {} torsion_pretwist: dict[str, dict[str, float]] = {} @@ -3076,6 +3175,8 @@ def end_key(token: str) -> str: self._check_torsion_normal(where, row, normal, direction, norm) torsion_ends.setdefault(line_id, {})[end] = restrained if non_pinned: + swapped = line_id in swapped_lines + held_ends.add((line_id, {"a": "b", "b": "a"}[end] if swapped else end)) has_finite_ei = any( line_type_meta[line_type][3] > 0.0 for line_type in section_types_by_line[line_id] @@ -3762,6 +3863,7 @@ def is_body_kind(kind: str, pinned: bool) -> bool: # Native check_channel: ids are digit strings compared by value, so leading zeros # are accepted (FairTen01 is FairTen1). + stiffened_lines = {str(_native_int(row.tokens[0])) for row in self._rows("STIFFENERS")} for channel in channels: lowered = channel.lower() endpoint_channel = re.fullmatch( @@ -3818,6 +3920,22 @@ def is_body_kind(kind: str, pinned: bool) -> bool: f"{self._label}: OUTPUT {channel!r} references an unknown line" ) continue + # FairDefl/AnchDefl: deflection of the end force from a held end's direction + deflection_channel = re.fullmatch(r"(fair|anch)defl([0-9]+)", lowered) + if deflection_channel: + line_id = _channel_id(deflection_channel.group(2)) + if line_id not in line_ids: + raise DeckFormatError( + f"{self._label}: OUTPUT {channel!r} references an unknown line" + ) + end = "a" if deflection_channel.group(1) == "fair" else "b" + if (line_id, end) not in held_ends: + raise DeckFormatError( + f"{self._label}: OUTPUT {channel!r}: the deflection angle needs a Rigid or " + "rotational-spring END CONNECTIONS row at End " + f"{end.upper()} of line {line_id}" + ) + continue if point_channel: # (a POINT row merged into a zero-length rod's connector no longer exists) point_id = _channel_id(point_channel.group(1)) @@ -3864,8 +3982,8 @@ def is_body_kind(kind: str, pinned: bool) -> bool: "no torque" ) continue - if lowered.startswith(("ten", "curv", "bendmom", "l")): - # native parse_line_node_channel: Ten/Curv/BendMomN, LN... + if lowered.startswith(("ten", "curv", "bendmom", "stfmom", "l")): + # native parse_line_node_channel: Ten/Curv/BendMom/StfMomN, LN... line_number, node_id, kind, _ = _parse_line_node_channel(lowered) if ( node_id == 0 @@ -3886,6 +4004,11 @@ def is_body_kind(kind: str, pinned: bool) -> bool: raise DeckFormatError( f"{self._label}: OUTPUT {channel!r} node exceeds the line node count" ) + if kind == 12 and line_id not in stiffened_lines: + raise DeckFormatError( + f"{self._label}: OUTPUT {channel!r}: line {line_id} has no STIFFENERS, " + "so no stiffener moment" + ) continue raise DeckFormatError(f"{self._label}: unsupported OUTPUT channel {channel!r}") diff --git a/python/tests/test_builder.py b/python/tests/test_builder.py index c24acc9..c45b3f4 100644 --- a/python/tests/test_builder.py +++ b/python/tests/test_builder.py @@ -107,6 +107,7 @@ def _semantic(deck: DeckFile) -> dict[str, object]: "TURBINES", "EQUIVALENT BUOYANCY", "ATTACHMENTS", + "STIFFENERS", "CONTROL", "EXTERNAL LOADS", ) @@ -425,6 +426,26 @@ def test_end_connections_attachments_and_equivalent_buoyancy() -> None: round_trip = DeckModel.from_deck_file(deck) assert round_trip.attachments[0].arc_length == (42.5, 5.0, 87.5) assert round_trip.end_connections[1].stiffness == "Rigid" + # a bend stiffener: a cylindrical root and a cone at End A + model.add_stiffener(1, "a", 0.4, 0.45, 0.45, 0.17, 5.0e7, 1150.0) + model.add_stiffener(model.lines[1], "EndA", 3.6, 0.45, 0.18, 0.17, 5.0e7, 1150.0) + deck = model.to_deck_file() + assert [row.tokens[:3] for row in deck._rows("STIFFENERS")] == [ + ("1", "A", "0.4"), + ("1", "A", "3.6"), + ] + round_trip = DeckModel.from_deck_file(deck) + assert [(row.end, row.length, row.tip_od) for row in round_trip.stiffeners] == [ + ("A", 0.4, 0.45), + ("A", 3.6, 0.18), + ] + refs = round_trip.references(round_trip.lines[1]) + assert all(any(ref is row for ref in refs) for row in round_trip.stiffeners) + with pytest.raises(ValueError, match="end must be A or B"): + model.add_stiffener(1, "C", 1.0, 0.4, 0.4, 0.2, 5.0e7, 1150.0) + model.to_text(validate=False) + model.remove(model.stiffeners[-1]) + assert len(model.stiffeners) == 2 # an end connection on an EI = 0 line is a native error model.line_types["cable"].ei = 0.0 with pytest.raises(DeckFormatError): diff --git a/python/tests/test_deck_native_parity.py b/python/tests/test_deck_native_parity.py index bc1f2d4..21d15d4 100644 --- a/python/tests/test_deck_native_parity.py +++ b/python/tests/test_deck_native_parity.py @@ -202,6 +202,63 @@ def test_attachment_arc_length_lies_on_the_line(arc, accepted): _rejects(text, "ATTACHMENTS ArcLength lies beyond the end of line 1") +_STIFFENER_ROOT = "1 A 0.4 0.45 0.45 0.17 5.0e7 1150" + + +@pytest.mark.parametrize( + ("row", "message"), + [ + ("1 A 0.4 0.45 0.45 0.17 5.0e7", "STIFFENERS rows must be: LineID End Length"), + ("x A 0.4 0.45 0.45 0.17 5.0e7 1150", "STIFFENERS LineID must be an integer"), + ("2 A 0.4 0.45 0.45 0.17 5.0e7 1150", "STIFFENERS row references undefined LINE id 2"), + ("1 C 0.4 0.45 0.45 0.17 5.0e7 1150", "STIFFENERS End must be A or B"), + ("1 A 0.4 0.45 0.45 0.17 5.0e7 dense", "must be plain numbers"), + ("1 A 0.4 0.45 0.45 0.17 2.0e15 1150", "must be finite and at most 1e15"), + ("1 A 0.0 0.45 0.45 0.17 5.0e7 1150", "Length, RootOD, TipOD and E must be > 0"), + ("1 A 0.4 0.45 0.45 -0.1 5.0e7 1150", "InnerD and Density must be >= 0"), + ("1 A 0.4 0.45 0.15 0.17 5.0e7 1150", "InnerD must be smaller than RootOD and TipOD"), + ("1 B 167.0 0.45 0.45 0.17 5.0e7 1150", "are longer than the line"), + ], +) +def test_stiffener_rows_are_native(row, message): + _rejects(_example("lazy_wave_bend_stiffener.dat", (_STIFFENER_ROOT, row)), message) + + +def test_stiffener_needs_a_finite_ei_line(): + text = _section(_HELD, "STIFFENERS", "1 A 1.0 0.5 0.5 0.3 5.0e7 1150") + _rejects(text, "STIFFENERS name line 1, which is not a finite-EI line") + + +@pytest.mark.parametrize( + ("channel", "message"), + [ + ("StfMom1N3", None), + ("FairDefl1", None), + ("AnchDefl1", "needs a Rigid or rotational-spring END CONNECTIONS row at End B of line 1"), + ], +) +def test_stiffener_and_deflection_channels_are_native(channel, message): + text = _example("lazy_wave_bend_stiffener.dat", ('"TDP1s"', f'"TDP1s"\n"{channel}"')) + if message is None: + DeckFile.from_text(text) + else: + _rejects(text, message) + + +def test_stiffener_moment_needs_a_stiffener(): + text = _example("lazy_wave_vessel_rao.dat", ('"TDP1s"', '"TDP1s"\n"StfMom1N3"')) + _rejects(text, "line 1 has no STIFFENERS, so no stiffener moment") + + +def test_stiffener_aliases_and_end_spellings_are_native(): + text = _example( + "lazy_wave_bend_stiffener.dat", + ("STIFFENERS ---", "BEND STIFFENERS ---"), + (_STIFFENER_ROOT, "1 enda 0.4 0.45 0.45 0.17 5.0e7 1150"), + ) + DeckFile.from_text(text) + + # --------------------------------------------------------------------------- output channels diff --git a/src/CableDyn_DeckDriver.f90 b/src/CableDyn_DeckDriver.f90 index c43aef6..ac0f6ac 100644 --- a/src/CableDyn_DeckDriver.f90 +++ b/src/CableDyn_DeckDriver.f90 @@ -77,7 +77,8 @@ MODULE CableDyn_DeckDriver USE CableDyn_HermiteArch, ONLY: CD_FiniteEI_Hermite_Arch_Seed_LengthMatched, CD_HARCH_OK USE CableDyn_HermiteCable, ONLY: CD_HermiteCable_Curvature, & CD_HermiteCable_Peak_Curvature, CD_HermiteCable_Axial_Resultant, & - CD_HermiteCable_Axial_Resultant_Range, CD_HermiteCable_Element, CD_HCABLE_OK + CD_HermiteCable_Axial_Resultant_Range, CD_HermiteCable_Element, CD_HCABLE_OK, & + CD_HermiteCable_Gauss_Rule USE CableDyn_HermiteCableStatic, ONLY: CD_HermiteCable_Static_Solve_AutoMesh, CD_HermiteCable_Static_Solve, & CD_HermiteCable_Prolong_Arc, & CD_HC_OSCILLATION_REL_TOL, & @@ -184,6 +185,7 @@ MODULE CableDyn_DeckDriver PUBLIC :: CD_Init_Deck_Models PUBLIC :: CD_Init_Deck_System PUBLIC :: CD_Init_Deck_HermiteCable + PUBLIC :: CD_Stiffener_Overlay PUBLIC :: CD_Init_Deck_Aggregate PUBLIC :: CD_Deck_Query_dtM PUBLIC :: CD_End_Deck_Aggregate @@ -475,6 +477,11 @@ MODULE CableDyn_DeckDriver ! End A [m], mass [kg], displaced volume [m^3], drag area CdA [m^2], added-mass Ca [-], ! axial drag area CdAx [m^2] (rows 1-6). REAL(wp), ALLOCATABLE :: att(:, :) + ! STIFFENERS (finite-EI lines): one column per conical segment of a bend stiffener -- the end + ! it sleeves (1 = End A, 2 = End B), the segment's arc from that end to its root side and to + ! its tip side [m], outer diameter at the root side and at the tip side [m], inner diameter + ! [m], Young's modulus [Pa] and density [kg/m^3] (rows 1-8). + REAL(wp), ALLOCATABLE :: stf(:, :) END TYPE TYPE :: DeckEndConnection @@ -1504,6 +1511,10 @@ SUBROUTINE init_one_finite_model(ln, sections, types, points, opts, deck_model, CALL fail(ErrStat, ErrMsg, 'finite-EI two-moving-end compatibility solver does not support '// & 'ATTACHMENTS; use a Fixed End B'); RETURN END IF + IF (ALLOCATED(ln%stf)) THEN + CALL fail(ErrStat, ErrMsg, 'finite-EI two-moving-end compatibility solver does not support '// & + 'STIFFENERS; use a Fixed End B'); RETURN + END IF IF (opts%friction_mu_axial > CD_ZERO) THEN CALL fail(ErrStat, ErrMsg, 'the finite-EI two-moving-end compatibility solver has isotropic seabed '// & 'friction only; remove frictionMuAxial/frictionMuLateral') @@ -1725,6 +1736,9 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, LOGICAL :: has_att REAL(wp), ALLOCATABLE :: w_bare_deck(:), l0_deck_att(:), w_bare(:) INTEGER, ALLOCATABLE :: att_node(:) + ! STIFFENERS: the cable's own bending stiffness per element (internal order), without the + ! sleeve, for the cable's share of the bend moment (allocated only on a stiffened line) + REAL(wp), ALLOCATABLE :: ei_cable(:) REAL(wp) :: att_offset ! Condensed torsion (both END CONNECTIONS rows torsionally restrained): solved as the last ! static load stage on the final mesh, then installed in the cable. @@ -1891,6 +1905,10 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, CALL build_finite_section_mesh(ln, sections, types, conn, l0, ea, mpl, diam, ei, & gas_in, gj_in, i_rt_in, i_rn_in, ErrStat, ErrMsg) IF (ErrStat /= CD_DECKDRV_OK) RETURN + IF (ALLOCATED(ln%stf)) THEN + ALLOCATE (ei_cable(SIZE(l0))) + CALL finite_line_type_ei(ln, sections, types, ei_cable) + END IF ! A seabed sloping across that plane under the as-laid run pushes the grounded nodes out ! of it against their friction: solved in 3D as well, so the march starts at rest. IF (.NOT. require_3d_static) require_3d_static = as_laid_cross_slope() @@ -2837,6 +2855,7 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, 'dynamic tolerance ', static_tol_target, '; equilibrium accepted at the relative residual ', res CALL report_note(note_text) END IF + IF (ALLOCATED(ln%stf)) CALL stiffener_notes() IF (saddle_kept) THEN WRITE (note_text, '(A,I0,A)') 'line ', ln%id, ': the equilibrium in the current is held in the '// & 'vertical plane of the line; it is unstable out of that plane and no out-of-plane equilibrium was '// & @@ -2900,6 +2919,7 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, CALL remap_elem_by_reference_arc(can, l0, l0_am) CALL remap_elem_by_reference_arc(cat, l0, l0_am) CALL remap_elem_by_reference_arc(ba_raw, l0, l0_am) + IF (ALLOCATED(ei_cable)) CALL remap_elem_by_reference_arc(ei_cable, l0, l0_am) END IF ! The node output channels read deck nodes: map each onto the final mesh (set on the cable ! once it is built below). @@ -3003,6 +3023,8 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, CALL fail(ErrStat, ErrMsg, 'finite-EI cable coupled-shell init failed: '//TRIM(em)); RETURN END IF IF (ALLOCATED(deck_map)) CALL MOVE_ALLOC(deck_map, cable%deck_node) + ! a stiffened line: the outputs report the cable's own share of the bend moment + IF (ALLOCATED(ei_cable)) cable%line%EI_line = ei_cable IF (ANY(ba > CD_ZERO)) THEN CALL CD_HFMF_Set_Axial_Damping(cable, ba, es, em) IF (es /= CD_HFMF_OK) THEN @@ -4375,6 +4397,50 @@ SUBROUTINE report_note(text) IF (PRESENT(init_note)) CALL append_note_line(init_note, text) END SUBROUTINE report_note + SUBROUTINE stiffener_notes() + !! Notes on each bend stiffener of the line, nothing changed: a stiffener spanning fewer than + !! STIFFENER_MIN_ELEMENTS deck elements, whose taper the element properties (one value per + !! deck element) cannot follow, with the element length that resolves it; and a stiffener at + !! a pinned end, where the root moment it is built to carry cannot develop. + INTEGER, PARAMETER :: STIFFENER_MIN_ELEMENTS = 10 + CHARACTER(1), PARAMETER :: END_NAME(2) = ['A', 'B'] + REAL(wp) :: s_end, len_stf, lo, hi + INTEGER :: iend, n_cov, e + CHARACTER(32) :: b1, b2 + CHARACTER(512) :: text + DO iend = 1, 2 + IF (.NOT. ANY(NINT(ln%stf(1, :)) == iend)) CYCLE + len_stf = MAXVAL(ln%stf(3, :), MASK=NINT(ln%stf(1, :)) == iend) + ! deck elements in internal order run from End B; count those the stiffener reaches + n_cov = 0 + s_end = CD_ZERO + DO e = 1, SIZE(l0) + IF (iend == 2) THEN + lo = s_end; hi = s_end + l0(e) + ELSE + lo = SUM(l0) - s_end - l0(e); hi = SUM(l0) - s_end + END IF + IF (lo < len_stf*(CD_ONE - 1.0e-12_wp)) n_cov = n_cov + 1 + s_end = s_end + l0(e) + END DO + WRITE (b1, '(F12.2)') len_stf + IF (n_cov < STIFFENER_MIN_ELEMENTS) THEN + WRITE (b2, '(F12.3)') len_stf/REAL(STIFFENER_MIN_ELEMENTS, wp) + text = 'line '//TRIM(int_to_str(ln%id))//': the '//TRIM(ADJUSTL(b1))//' m stiffener at End '// & + END_NAME(iend)//' spans '//TRIM(int_to_str(n_cov))//' elements, too few to follow its taper '// & + '(one stiffness per element): give the section there elements of at most '// & + TRIM(ADJUSTL(b2))//' m over the stiffener' + CALL report_note(text) + END IF + IF (ln%endconn_mode(iend) == CD_ENDCONN_PINNED) THEN + text = 'line '//TRIM(int_to_str(ln%id))//': the stiffener at End '//END_NAME(iend)//' sleeves a '// & + 'pinned end, which carries no moment, so the stiffener root takes none; declare that end '// & + 'Rigid or a rotational spring in END CONNECTIONS to model a stiffener clamped at its root' + CALL report_note(text) + END IF + END DO + END SUBROUTINE stiffener_notes + SUBROUTINE friction_resolution_note() !! Seabed friction acts at the nodes: as the line lands and lifts at touchdown, a grounded !! element locks or releases up to mu w L of friction (mu the axial friction coefficient, @@ -9082,6 +9148,7 @@ SUBROUTINE build_finite_section_mesh(ln, sections, types, conn, l0, ea, mpl, dia i_rn(e) = types(tidx)%i_rn END DO END DO + IF (ALLOCATED(ln%stf)) CALL CD_Stiffener_Overlay(ln%stf, l0, diam, ei=ei, mpl=mpl, volume_diameter=.TRUE.) l0 = l0(n_elem:1:-1) ea = ea(n_elem:1:-1) mpl = mpl(n_elem:1:-1) @@ -9115,12 +9182,33 @@ SUBROUTINE finite_raw_ba_array(ln, sections, types, ba_raw) ba_raw = ba_raw(SIZE(ba_raw):1:-1) END SUBROUTINE finite_raw_ba_array + SUBROUTINE finite_line_type_ei(ln, sections, types, ei) + !! The line types' own bending stiffness on the deck elements (internal order), without + !! any STIFFENERS sleeve. + TYPE(DeckLine), INTENT(IN) :: ln + TYPE(DeckSection), INTENT(IN) :: sections(:) + TYPE(DeckLineType), INTENT(IN) :: types(:) + REAL(wp), INTENT(OUT) :: ei(:) + INTEGER :: i, j, e + + e = 0 + DO i = 1, SIZE(sections) + IF (sections(i)%line_id /= ln%id) CYCLE + DO j = 1, sections(i)%nsegs + e = e + 1 + ei(e) = types(sections(i)%ltype)%ei + END DO + END DO + ei = ei(SIZE(ei):1:-1) + END SUBROUTINE finite_line_type_ei + SUBROUTINE finite_hydro_arrays(ln, sections, types, diam, cdn, cdt, can, cat) !! Expand section hydro coefficients to finite-EI element arrays. TYPE(DeckLine), INTENT(IN) :: ln TYPE(DeckSection), INTENT(IN) :: sections(:) TYPE(DeckLineType), INTENT(IN) :: types(:) REAL(wp), INTENT(OUT) :: diam(:), cdn(:), cdt(:), can(:), cat(:) + REAL(wp) :: l0(SIZE(diam)) INTEGER :: i, j, e, tidx e = 0 @@ -9129,6 +9217,7 @@ SUBROUTINE finite_hydro_arrays(ln, sections, types, diam, cdn, cdt, can, cat) tidx = sections(i)%ltype DO j = 1, sections(i)%nsegs e = e + 1 + l0(e) = sections(i)%length/REAL(sections(i)%nsegs, wp) diam(e) = types(tidx)%diam cdn(e) = types(tidx)%cdn cdt(e) = types(tidx)%cdt @@ -9136,6 +9225,8 @@ SUBROUTINE finite_hydro_arrays(ln, sections, types, diam, cdn, cdt, can, cat) cat(e) = types(tidx)%cat END DO END DO + ! a stiffener widens the line to its outer diameter for drag and added mass + IF (ALLOCATED(ln%stf)) CALL CD_Stiffener_Overlay(ln%stf, l0, diam, volume_diameter=.FALSE.) diam = diam(SIZE(diam):1:-1) cdn = cdn(SIZE(cdn):1:-1) cdt = cdt(SIZE(cdt):1:-1) @@ -9143,6 +9234,83 @@ SUBROUTINE finite_hydro_arrays(ln, sections, types, diam, cdn, cdt, can, cat) cat = cat(SIZE(cat):1:-1) END SUBROUTINE finite_hydro_arrays + SUBROUTINE CD_Stiffener_Overlay(stf, l0, diam, ei, mpl, volume_diameter) + !! Add the bend stiffeners of a finite-EI line (ln%stf columns) to its element properties, + !! in public element order (End A first). A stiffener is a polyurethane sleeve bonded to the + !! cable in bending: over its length the section has the bending stiffness of the cable plus + !! that of the sleeve, EI_c + E pi/64 (OD^4 - ID^4), with the outer diameter OD varying + !! linearly along each conical segment; it carries no axial load and no torque. An element it + !! covers, wholly or in part, takes + !! ei the compliance average l / int ds/EI(s), exact for the element's rotation under a + !! uniform moment, so a stiff root next to a soft tip is not over-stiffened; + !! mpl the cable mass plus the mean sleeve mass, rho pi/4 (OD^2 - ID^2); + !! diam the diameter of the mean displaced area, cable plus solid sleeve, sqrt(d^2 + mean + !! (OD^2 - ID^2)) (volume_diameter: weight and buoyancy), or the mean outer diameter, + !! mean max(OD, d) (otherwise: drag and added mass). + !! The integrals use 6-point Gauss-Legendre on each part of the element a segment covers. + !! An element no stiffener reaches is left untouched. + REAL(wp), INTENT(IN) :: stf(:, :), l0(:) + REAL(wp), INTENT(INOUT) :: diam(:) + REAL(wp), INTENT(INOUT), OPTIONAL :: ei(:), mpl(:) + LOGICAL, INTENT(IN) :: volume_diameter + REAL(wp), PARAMETER :: PI_64 = 0.0490873852123405193_wp ! pi/64 + REAL(wp) :: gp(6), gw(6), s_end(0:SIZE(l0)), line_len, a, b, lo, hi, x, t, od, d_line + REAL(wp) :: compliance, sleeve_mass, sleeve_area, outer, covered, ei_line + INTEGER :: e, k, g + + CALL CD_HermiteCable_Gauss_Rule(6, gp, gw) + s_end(0) = CD_ZERO + DO e = 1, SIZE(l0) + s_end(e) = s_end(e - 1) + l0(e) + END DO + line_len = s_end(SIZE(l0)) + DO e = 1, SIZE(l0) + a = s_end(e - 1) + b = s_end(e) + d_line = diam(e) + ei_line = CD_ZERO + IF (PRESENT(ei)) ei_line = ei(e) + compliance = CD_ZERO; sleeve_mass = CD_ZERO; sleeve_area = CD_ZERO; outer = CD_ZERO; covered = CD_ZERO + DO k = 1, SIZE(stf, 2) + ! the segment's span on the public arc + IF (NINT(stf(1, k)) == 1) THEN + lo = stf(2, k); hi = stf(3, k) + ELSE + lo = line_len - stf(3, k); hi = line_len - stf(2, k) + END IF + lo = MAX(lo, a); hi = MIN(hi, b) + IF (.NOT. (hi > lo)) CYCLE + covered = covered + (hi - lo) + DO g = 1, 6 + x = lo + gp(g)*(hi - lo) + ! fraction of the way from the segment's root side to its tip side + IF (NINT(stf(1, k)) == 1) THEN + t = (x - stf(2, k))/(stf(3, k) - stf(2, k)) + ELSE + t = (line_len - x - stf(2, k))/(stf(3, k) - stf(2, k)) + END IF + t = MIN(CD_ONE, MAX(CD_ZERO, t)) + od = stf(4, k) + t*(stf(5, k) - stf(4, k)) + compliance = compliance + gw(g)*(hi - lo)/(ei_line + stf(7, k)*PI_64*(od**4 - stf(6, k)**4)) + sleeve_mass = sleeve_mass + gw(g)*(hi - lo)*stf(8, k)*QUARTER_PI*(od**2 - stf(6, k)**2) + sleeve_area = sleeve_area + gw(g)*(hi - lo)*(od**2 - stf(6, k)**2) + outer = outer + gw(g)*(hi - lo)*MAX(od, d_line) + END DO + END DO + IF (.NOT. (covered > CD_ZERO)) CYCLE + IF (PRESENT(ei)) THEN + IF (ei_line > CD_ZERO) compliance = compliance + (l0(e) - covered)/ei_line + ei(e) = l0(e)/compliance + END IF + IF (PRESENT(mpl)) mpl(e) = mpl(e) + sleeve_mass/l0(e) + IF (volume_diameter) THEN + diam(e) = SQRT(d_line**2 + sleeve_area/l0(e)) + ELSE + diam(e) = (outer + (l0(e) - covered)*d_line)/l0(e) + END IF + END DO + END SUBROUTINE CD_Stiffener_Overlay + SUBROUTINE remap_elem_by_reference_arc(a, l0_parent, l0_child) !! Map a piecewise-constant per-element field from a parent reference mesh to !! any nested child mesh. Child midpoints select the parent, so interfaces are @@ -10871,6 +11039,9 @@ SUBROUTINE parse_deck(path, types, bodies, rod_types, rods, points, lines, secti CHARACTER(DECK_RECLEN), ALLOCATABLE :: att_rows(:) INTEGER, ALLOCATABLE :: att_lines(:) INTEGER :: n_att + CHARACTER(DECK_RECLEN), ALLOCATABLE :: stf_rows(:) + INTEGER, ALLOCATABLE :: stf_lines(:) + INTEGER :: n_stf CHARACTER(2*NAMELEN), ALLOCATABLE :: channel_keys(:) LOGICAL :: ext_march, types_stock_cols, rodtypes_legacy_cols, turb_ok, static_finite_ok INTEGER :: unit, ios, section, nt, nb, nrt, nr, np, nlin, nsec, neq, nch, nfl, nct, nec, lineno @@ -10887,7 +11058,8 @@ SUBROUTINE parse_deck(path, types, bodies, rod_types, rods, points, lines, secti INTEGER, PARAMETER :: SEC_NONE = 0, SEC_TYPES = 1, SEC_BODIES = 2, SEC_POINTS = 3, SEC_LINES = 4, & SEC_SECTIONS = 5, SEC_OPTIONS = 6, SEC_OUTPUTS = 7, SEC_RODTYPES = 8, SEC_RODS = 9, & SEC_EQUIV = 10, SEC_FAILURE = 11, SEC_CONTROL = 12, SEC_SYROPE_IC = 13, & - SEC_ENDCONN = 14, SEC_TURBINES = 15, SEC_EXTLOADS = 16, SEC_ATTACH = 17 + SEC_ENDCONN = 14, SEC_TURBINES = 15, SEC_EXTLOADS = 16, SEC_ATTACH = 17, & + SEC_STIFF = 18 ErrStat = 0; ErrMsg = '' ext_march = .FALSE. @@ -10907,6 +11079,8 @@ SUBROUTINE parse_deck(path, types, bodies, rod_types, rods, points, lines, secti next_ext = 0 ALLOCATE (att_rows(0), att_lines(0)) n_att = 0 + ALLOCATE (stf_rows(0), stf_lines(0)) + n_stf = 0 nt = 0; nb = 0; nrt = 0; nr = 0; np = 0; nlin = 0; nsec = 0; neq = 0; nch = 0; nfl = 0; nct = 0; nec = 0 section = SEC_NONE ALLOCATE (channel_lines(0)) @@ -11023,6 +11197,15 @@ SUBROUTINE parse_deck(path, types, bodies, rod_types, rods, points, lines, secti END IF att_rows(n_att) = line(1:MIN(LEN(line), DECK_RECLEN)) att_lines(n_att) = lineno + CASE (SEC_STIFF) + IF (.NOT. has_numeric_token(line)) CYCLE + n_stf = n_stf + 1 + IF (n_stf > SIZE(stf_rows)) THEN + stf_rows = [stf_rows, [(REPEAT(' ', DECK_RECLEN), i_ch=1, 8)]] + stf_lines = [stf_lines, [(0, i_ch=1, 8)]] + END IF + stf_rows(n_stf) = line(1:MIN(LEN(line), DECK_RECLEN)) + stf_lines(n_stf) = lineno CASE (SEC_OPTIONS) CALL apply_option(line, opts, ErrStat, ErrMsg) CASE (SEC_OUTPUTS) @@ -11113,6 +11296,13 @@ SUBROUTINE parse_deck(path, types, bodies, rod_types, rods, points, lines, secti RETURN END IF END DO + DO i_ch = 1, n_stf + CALL append_line_stiffener(stf_rows(i_ch), lines, ErrStat, ErrMsg) + IF (ErrStat /= 0) THEN + CALL add_row_context('deck', stf_lines(i_ch), stf_rows(i_ch), ErrMsg) + RETURN + END IF + END DO ! Caller-supplied environment overrides land HERE -- before equivalent-buoyancy ! conversion and validation -- so both run under the environment the model is @@ -12302,6 +12492,8 @@ SUBROUTINE set_section(secname, section, ErrStat, ErrMsg) section = 16 CASE ('ATTACHMENTS', 'LINE ATTACHMENTS', 'CLUMPS') section = 17 + CASE ('STIFFENERS', 'STIFFENER', 'BEND STIFFENERS', 'BEND STIFFENER') + section = 18 CASE DEFAULT IF (LEN_TRIM(secname) > NAMELEN) THEN CALL fail(ErrStat, ErrMsg, 'unknown deck section "'//secname(1:NAMELEN)//'..."') @@ -13637,6 +13829,81 @@ SUBROUTINE append_line_attachment(row, lines, ErrStat, ErrMsg) CALL MOVE_ALLOC(grown, lines(li)%att) END SUBROUTINE append_line_attachment + SUBROUTINE append_line_stiffener(row, lines, ErrStat, ErrMsg) + !! One STIFFENERS row: LineID End Length RootOD TipOD InnerD E Density -- a conical segment of + !! the bend stiffener sleeving End A or End B of a finite-EI line. The rows of one line end + !! follow each other outward from the root, each starting where the previous one ends, so a + !! cylindrical root, a cone and a tip are three rows. RootOD and TipOD are the outer diameters + !! of the segment at its root side and at its tip side. Stored columns: end (1 = A, 2 = B), + !! arc of the root side and of the tip side from that end, RootOD, TipOD, InnerD, E, Density. + CHARACTER(*), INTENT(IN) :: row + TYPE(DeckLine), INTENT(INOUT) :: lines(:) + INTEGER, INTENT(OUT) :: ErrStat + CHARACTER(*), INTENT(OUT) :: ErrMsg + CHARACTER(DECK_RECLEN) :: tok + INTEGER :: li, lid, k, iend, nold + REAL(wp) :: vals(6), s0 + REAL(wp), ALLOCATABLE :: grown(:, :) + LOGICAL :: found, ok + + ErrStat = 0; ErrMsg = '' + IF (count_tokens(row) /= 8) THEN + CALL fail(ErrStat, ErrMsg, 'STIFFENERS rows must be: LineID End Length RootOD TipOD InnerD E Density'); RETURN + END IF + CALL get_token(row, 1, tok, found) + CALL parse_int_token(tok, lid, ok) + IF (.NOT. ok) THEN + CALL fail(ErrStat, ErrMsg, 'STIFFENERS LineID must be an integer'); RETURN + END IF + li = find_line(lines, lid) + IF (li == 0) THEN + CALL fail(ErrStat, ErrMsg, 'STIFFENERS row references undefined LINE id '//TRIM(int_to_str(lid))); RETURN + END IF + CALL get_token(row, 2, tok, found) + SELECT CASE (to_lower(TRIM(tok))) + CASE ('a', 'enda') + iend = 1 + CASE ('b', 'endb') + iend = 2 + CASE DEFAULT + CALL fail(ErrStat, ErrMsg, 'STIFFENERS End must be A or B'); RETURN + END SELECT + DO k = 1, 6 + CALL get_token(row, k + 2, tok, found) + CALL parse_real_token(tok, vals(k), ok) + IF (.NOT. ok) THEN + CALL fail(ErrStat, ErrMsg, 'STIFFENERS Length, RootOD, TipOD, InnerD, E and Density must be plain numbers') + RETURN + END IF + END DO + IF (.NOT. ALL(IEEE_IS_FINITE(vals)) .OR. ANY(vals > MAX_COEFFICIENT)) THEN + CALL fail(ErrStat, ErrMsg, 'STIFFENERS Length, RootOD, TipOD, InnerD, E and Density must be finite '// & + 'and at most 1e15'); RETURN + END IF + IF (.NOT. (vals(1) > CD_ZERO .AND. vals(2) > CD_ZERO .AND. vals(3) > CD_ZERO .AND. vals(5) > CD_ZERO)) THEN + CALL fail(ErrStat, ErrMsg, 'STIFFENERS Length, RootOD, TipOD and E must be > 0'); RETURN + END IF + IF (vals(4) < CD_ZERO .OR. vals(6) < CD_ZERO) THEN + CALL fail(ErrStat, ErrMsg, 'STIFFENERS InnerD and Density must be >= 0'); RETURN + END IF + IF (.NOT. (vals(4) < MIN(vals(2), vals(3)))) THEN + CALL fail(ErrStat, ErrMsg, 'STIFFENERS InnerD must be smaller than RootOD and TipOD'); RETURN + END IF + ! a segment starts where the previous segment of the same line end stops + s0 = CD_ZERO + nold = 0 + IF (ALLOCATED(lines(li)%stf)) THEN + nold = SIZE(lines(li)%stf, 2) + DO k = 1, nold + IF (NINT(lines(li)%stf(1, k)) == iend) s0 = MAX(s0, lines(li)%stf(3, k)) + END DO + END IF + ALLOCATE (grown(8, nold + 1)) + IF (nold > 0) grown(:, 1:nold) = lines(li)%stf + grown(:, nold + 1) = [REAL(iend, wp), s0, s0 + vals(1), vals(2:6)] + CALL MOVE_ALLOC(grown, lines(li)%stf) + END SUBROUTINE append_line_stiffener + SUBROUTINE trim_end_connections(endconns, nec) TYPE(DeckEndConnection), ALLOCATABLE, INTENT(INOUT) :: endconns(:) INTEGER, INTENT(IN) :: nec @@ -14804,6 +15071,9 @@ FUNCTION channel_identity_key(ch) RESULT(key) ELSE key = lo(1:4)//'decl:'//TRIM(int_to_str(a)) END IF + ELSE IF (is_end_deflection_channel(lo)) THEN + a = int_suffix(ch, 9) + IF (a > 0) key = lo(1:8)//':'//TRIM(int_to_str(a)) ELSE IF (is_point_force_channel(lo)) THEN CALL parse_point_force_channel(ch, a, b) IF (a > 0 .AND. b > 0) key = 'pointf:'//TRIM(int_to_str(a))//':'//TRIM(int_to_str(b)) @@ -14814,7 +15084,8 @@ FUNCTION channel_identity_key(ch) RESULT(key) CALL parse_point_channel(ch, a, b) IF (a > 0 .AND. b > 0) key = 'point:'//TRIM(int_to_str(a))//':'//TRIM(int_to_str(b)) ELSE IF (starts_with(lo, 'ten') .OR. starts_with(lo, 'curv') .OR. starts_with(lo, 'bendmom') .OR. & - starts_with(lo, 'l') .OR. starts_with(lo, 'torq') .OR. starts_with(lo, 'twist')) THEN + starts_with(lo, 'stfmom') .OR. starts_with(lo, 'l') .OR. starts_with(lo, 'torq') .OR. & + starts_with(lo, 'twist')) THEN CALL parse_line_node_channel(ch, a, b, c, d) IF (a > 0 .AND. b > 0 .AND. c > 0) key = 'node:'//TRIM(int_to_str(a))//':'//TRIM(int_to_str(b))// & ':'//TRIM(int_to_str(c))//':'//TRIM(int_to_str(d)) @@ -14935,6 +15206,8 @@ SUBROUTINE normalize_stock_endpoint_order(points, lines, sections) ! attachment arcs are measured from End A, which is now the other end IF (ALLOCATED(lines(i)%att)) & lines(i)%att(1, :) = SUM(sections%length, MASK=sections%line_id == lines(i)%id) - lines(i)%att(1, :) + ! a stiffener sleeves the same physical end, which is now the other end index + IF (ALLOCATED(lines(i)%stf)) lines(i)%stf(1, :) = 3.0_wp - lines(i)%stf(1, :) CALL reverse_line_sections(sections, lines(i)%id) END IF END DO @@ -15919,6 +16192,23 @@ SUBROUTINE finalize_checks(types, bodies, rod_types, rods, points, lines, sectio END IF END DO + ! STIFFENERS: bend stiffeners sleeving the ends of a finite-EI line + DO j = 1, SIZE(lines) + IF (.NOT. ALLOCATED(lines(j)%stf)) CYCLE + at_line = lines(j)%deck_line + IF (.NOT. line_is_finite_ei(lines(j), types, sections)) THEN + CALL fail(ErrStat, ErrMsg, 'STIFFENERS name line '//TRIM(int_to_str(lines(j)%id))// & + ', which is not a finite-EI line: a stiffener adds bending stiffness, which an EI = 0 line '// & + 'does not have'); RETURN + END IF + IF (MAX(CD_ZERO, MAXVAL(lines(j)%stf(3, :), MASK=NINT(lines(j)%stf(1, :)) == 1)) + & + MAX(CD_ZERO, MAXVAL(lines(j)%stf(3, :), MASK=NINT(lines(j)%stf(1, :)) == 2)) > & + SUM(sections%length, MASK=sections%line_id == lines(j)%id)*(CD_ONE + 1.0e-12_wp)) THEN + CALL fail(ErrStat, ErrMsg, 'the STIFFENERS of line '//TRIM(int_to_str(lines(j)%id))// & + ' are longer than the line (End A and End B stiffeners together)'); RETURN + END IF + END DO + ! ATTACHMENTS: discrete clumps/buoyancy modules on a finite-EI line DO j = 1, SIZE(lines) IF (.NOT. ALLOCATED(lines(j)%att)) CYCLE @@ -27602,7 +27892,7 @@ SUBROUTINE write_hermite_static_element_profile(out_root, lines, cables, ErrStat minimum_u_public = CD_ONE - minimum_u_internal maximum_u_public = CD_ONE - maximum_u_internal arc_peak = arc_start + u_curvature_public*l0 - moment_peak = cables(li)%line%EI(internal_element)*curvature_peak + moment_peak = cables(li)%line%EI_line(internal_element)*curvature_peak WRITE (unit, '(I0,A,I0,10(A,ES15.7E3))') lines(li)%id, tab, public_element, & tab, arc_start, tab, arc_end, tab, u_curvature_public, tab, arc_peak, tab, curvature_peak, & tab, moment_peak, tab, minimum_resultant, tab, minimum_u_public, tab, maximum_resultant, & @@ -27787,6 +28077,8 @@ SUBROUTINE write_hermite_model_row(unit, time, channels, lines, points, cables, line_id = anchor_force_line_id(channels(i)) ELSE IF (is_endpoint_angle_channel(lo)) THEN line_id = endpoint_angle_line_id(channels(i)) + ELSE IF (is_end_deflection_channel(lo)) THEN + line_id = int_suffix(channels(i), 9) ELSE CALL parse_line_node_channel(channels(i), line_id, pidx, pp, comp) END IF @@ -31647,6 +31939,13 @@ LOGICAL FUNCTION is_endpoint_angle_channel(lo) RESULT(is_angle) starts_with(lo, 'fairincl') .OR. starts_with(lo, 'anchincl') END FUNCTION is_endpoint_angle_channel + LOGICAL FUNCTION is_end_deflection_channel(lo) RESULT(is_defl) + !! FairDefl/AnchDefl: the angle between the end force and the direction a Rigid or + !! rotational-spring END CONNECTIONS row holds at End A / End B [deg]. + CHARACTER(*), INTENT(IN) :: lo + is_defl = starts_with(lo, 'fairdefl') .OR. starts_with(lo, 'anchdefl') + END FUNCTION is_end_deflection_channel + LOGICAL FUNCTION is_fair_endpoint_angle(lo) RESULT(is_fair) CHARACTER(*), INTENT(IN) :: lo is_fair = starts_with(lo, 'fairangle') .OR. starts_with(lo, 'fairdecl') .OR. & @@ -31814,6 +32113,17 @@ SUBROUTINE check_channel(ch, lines, points, sections, ErrStat, ErrMsg, bodies, r IF (find_line(lines, endpoint_angle_line_id(ch)) == 0) THEN CALL fail(ErrStat, ErrMsg, 'OUTPUT "'//TRIM(ch)//'": unknown line id'); RETURN END IF + ELSE IF (is_end_deflection_channel(lo)) THEN + li = find_line(lines, int_suffix(ch, 9)) + IF (li == 0) THEN + CALL fail(ErrStat, ErrMsg, 'OUTPUT "'//TRIM(ch)//'": unknown line id'); RETURN + END IF + ! the angle is measured from the direction a held end connection declares + IF (lines(li)%endconn_mode(MERGE(1, 2, starts_with(lo, 'fairdefl'))) == CD_ENDCONN_PINNED) THEN + CALL fail(ErrStat, ErrMsg, 'OUTPUT "'//TRIM(ch)//'": the deflection angle needs a Rigid or '// & + 'rotational-spring END CONNECTIONS row at End '//MERGE('A', 'B', starts_with(lo, 'fairdefl'))// & + ' of line '//TRIM(int_to_str(lines(li)%id))); RETURN + END IF ELSE IF (is_point_force_channel(lo)) THEN CALL parse_point_force_channel(ch, pidx, comp) IF (comp < 1 .OR. point_index(points, pidx) == 0) THEN @@ -31854,7 +32164,7 @@ SUBROUTINE check_channel(ch, lines, points, sections, ErrStat, ErrMsg, bodies, r 'torque'); RETURN END IF ELSE IF (starts_with(lo, 'ten') .OR. starts_with(lo, 'curv') .OR. & - starts_with(lo, 'bendmom') .OR. starts_with(lo, 'l')) THEN + starts_with(lo, 'bendmom') .OR. starts_with(lo, 'stfmom') .OR. starts_with(lo, 'l')) THEN CALL parse_line_node_channel(ch, line_id, node, kind, comp) li = find_line(lines, line_id) IF (li /= 0 .AND. node == 0 .AND. names_node_zero(lo)) THEN @@ -31868,12 +32178,16 @@ SUBROUTINE check_channel(ch, lines, points, sections, ErrStat, ErrMsg, bodies, r IF (node > nelem + 1) THEN CALL fail(ErrStat, ErrMsg, 'OUTPUT "'//TRIM(ch)//'": node index exceeds line node count'); RETURN END IF + IF (kind == 12 .AND. .NOT. ALLOCATED(lines(li)%stf)) THEN + CALL fail(ErrStat, ErrMsg, 'OUTPUT "'//TRIM(ch)//'": line '//TRIM(int_to_str(line_id))//' has no '// & + 'STIFFENERS, so no stiffener moment'); RETURN + END IF ELSE CALL fail(ErrStat, ErrMsg, 'OUTPUT channel "'//TRIM(ch)// & '" not supported (use FairTen, AnchTen, AnchFz, AnchFH, FairIncl, AnchIncl, '// & - 'FairDecl, AnchDecl (FairAngle/AnchAngle aliases), '// & + 'FairDecl, AnchDecl (FairAngle/AnchAngle aliases), FairDefl, AnchDefl, '// & 'Point

p{x,y,z} (alias Con

p{x,y,z}), Point

F{x,y,z,H}, '// & - 'TenN, CurvN, BendMomN, TorqN, TwistN, Twist, '// & + 'TenN, CurvN, BendMomN, StfMomN, TorqN, TwistN, Twist, '// & 'LNp{x,y,z}, LNv{x,y,z}, LNa{x,y,z}, LNDec, LNAzi, '// & 'TDP{s,x,y,z,Lay,Exc,Grd}, or the Body/Rod channels)'); RETURN END IF @@ -32107,7 +32421,8 @@ SUBROUTINE parse_line_node_channel(ch, line_id, node_id, kind, comp) !! Parse supported line-node output channels (kind codes shared with the writers): !! TenN -> kind=3, comp=0 nodal tension [N] !! CurvN -> kind=5, comp=0 curvature [1/m] - !! BendMomN -> kind=6, comp=0 bend moment [N.m] + !! BendMomN -> kind=6, comp=0 cable bend moment [N.m] + !! StfMomN -> kind=12, comp=0 stiffener moment [N.m] !! LNp{x,y,z} -> kind=1, comp=1..3 position [m] !! LNv{x,y,z} -> kind=2, comp=1..3 velocity [m/s] !! LNa{x,y,z} -> kind=4, comp=1..3 acceleration [m/s^2] @@ -32160,6 +32475,10 @@ SUBROUTINE parse_line_node_channel(ch, line_id, node_id, kind, comp) CALL parse_keyword_line_node(lo, 7, last, line_id, node_id) IF (line_id > 0 .AND. node_id > 0) kind = 6 RETURN + ELSE IF (starts_with(lo, 'stfmom')) THEN + CALL parse_keyword_line_node(lo, 6, last, line_id, node_id) + IF (line_id > 0 .AND. node_id > 0) kind = 12 + RETURN END IF ! LN... channels: p/v/a{x,y,z} (component) or Dec/Azi (scalar). IF (.NOT. starts_with(lo, 'l')) RETURN @@ -33099,11 +33418,11 @@ SUBROUTINE range_sample_cable_element(ErrStat, ErrMsg) bend = CD_ZERO IF (inode > 1) THEN curv = MAX(curv, rl%as(inode - 1)) - bend = MAX(bend, cable%line%EI(inode - 1)*rl%as(inode - 1)) + bend = MAX(bend, cable%line%EI_line(inode - 1)*rl%as(inode - 1)) END IF IF (inode <= ne) THEN curv = MAX(curv, rl%vs(inode)) - bend = MAX(bend, cable%line%EI(inode)*rl%vs(inode)) + bend = MAX(bend, cable%line%EI_line(inode)*rl%vs(inode)) END IF rl%val(CD_RQ_CURVATURE, j) = curv rl%val(CD_RQ_BEND, j) = bend @@ -33465,6 +33784,8 @@ SUBROUTINE CD_Eval_Aggregate_Channel(ch, has_sys, system, cables, line_is_cable, line_id = anchor_force_line_id(ch) ELSE IF (is_endpoint_angle_channel(lo)) THEN line_id = endpoint_angle_line_id(ch) + ELSE IF (is_end_deflection_channel(lo)) THEN + line_id = int_suffix(ch, 9) ELSE CALL parse_line_node_channel(ch, line_id, node, kind, comp) END IF @@ -33632,6 +33953,8 @@ LOGICAL FUNCTION CD_Channel_Token_Parses(ch) RESULT(ok) ok = anchor_force_line_id(ch) > 0 ELSE IF (is_endpoint_angle_channel(lo)) THEN ok = endpoint_angle_line_id(ch) > 0 + ELSE IF (is_end_deflection_channel(lo)) THEN + ok = int_suffix(ch, 9) > 0 ELSE CALL parse_line_node_channel(ch, line_id, node, kind, comp) ok = line_id > 0 .AND. node > 0 .AND. kind > 0 @@ -33844,6 +34167,27 @@ SUBROUTINE eval_cable_channel(cable, ch, val, ErrStat, ErrMsg, mesh_node) END IF val = endpoint_angle_from_declination(lo, val) RETURN + ELSE IF (is_end_deflection_channel(lo)) THEN + ! the angle between the force on the end's support and the held direction: internally End + ! A is node nn and End B node 1, and the held direction follows the internal tangent (out + ! of the line at End A, into it at End B), while the end force points into the line + BLOCK + REAL(wp) :: force(3), axis(3), cosine + INTEGER :: inode_end + inode_end = MERGE(nnode, 1, starts_with(lo, 'fairdefl')) + CALL CD_HermiteCable_Dyn_End_Force(cable%line, inode_end, force, es, em) + IF (es /= CD_HCDYN_OK) THEN + CALL fail_solve(ErrStat, ErrMsg, 'OUTPUT "'//TRIM(ch)//'": finite-EI end-force recovery failed: '// & + TRIM(em)); RETURN + END IF + axis = MERGE(-cable%line%endconn_d0(:, 2), cable%line%endconn_d0(:, 1), inode_end == nnode) + val = CD_ZERO + IF (NORM2(force) > CD_ZERO .AND. NORM2(axis) > CD_ZERO) THEN + cosine = DOT_PRODUCT(force, axis)/(NORM2(force)*NORM2(axis)) + val = ACOS(MAX(-CD_ONE, MIN(CD_ONE, cosine)))*CD_RAD2DEG + END IF + END BLOCK + RETURN END IF CALL parse_line_node_channel(ch, line_id, node, kind, comp) deck_numbering = ALLOCATED(cable%deck_node) @@ -33905,6 +34249,15 @@ SUBROUTINE eval_cable_channel(cable, ch, val, ErrStat, ErrMsg, mesh_node) IF (es /= CD_HCABLE_OK) THEN CALL fail_solve(ErrStat, ErrMsg, 'OUTPUT "'//TRIM(ch)//'": '//TRIM(em)); RETURN END IF + CASE (12) ! the bend stiffener's share of the section moment (zero outside a stiffener) + IF (.NOT. cable%line%output_element_recovery) THEN + CALL cable_node_values(cable, inode, sec_ten, sec_curv, sec_bend, ErrStat, ErrMsg, sleeve_moment=val) + RETURN + END IF + CALL cable_nodal_bend_moment(cable, inode, sec_bend, es, em, sleeve_moment=val) + IF (es /= CD_HCABLE_OK) THEN + CALL fail_solve(ErrStat, ErrMsg, 'OUTPUT "'//TRIM(ch)//'": '//TRIM(em)); RETURN + END IF CASE (9, 10, 11) ! condensed torsion: torque [N m], twist from End A to node J, total twist [deg] BLOCK REAL(wp) :: th, mt @@ -34080,18 +34433,20 @@ SUBROUTINE cable_end_tension(cable, inode, tension, ErrStat, ErrMsg) tension = NORM2(force) END SUBROUTINE cable_end_tension - SUBROUTINE cable_node_values(cable, inode, tension, curvature, bend_moment, ErrStat, ErrMsg) + SUBROUTINE cable_node_values(cable, inode, tension, curvature, bend_moment, ErrStat, ErrMsg, sleeve_moment) !! Equilibrium node values of a finite-EI cable at internal node inode - !! (CD_HermiteCable_Dyn_Node_Values): effective tension, curvature and bend moment. + !! (CD_HermiteCable_Dyn_Node_Values): effective tension, curvature, the cable's bend moment and + !! (optional) the share a bend stiffener carries. TYPE(CD_HFMF_ModuleType), INTENT(INOUT) :: cable INTEGER, INTENT(IN) :: inode REAL(wp), INTENT(OUT) :: tension, curvature, bend_moment INTEGER, INTENT(OUT) :: ErrStat CHARACTER(*), INTENT(OUT) :: ErrMsg + REAL(wp), INTENT(OUT), OPTIONAL :: sleeve_moment INTEGER :: es CHARACTER(1024) :: em ErrStat = CD_DECKDRV_OK; ErrMsg = '' - CALL CD_HermiteCable_Dyn_Node_Values(cable%line, inode, tension, curvature, bend_moment, es, em) + CALL CD_HermiteCable_Dyn_Node_Values(cable%line, inode, tension, curvature, bend_moment, es, em, sleeve_moment) IF (es /= CD_HCDYN_OK) CALL fail_solve(ErrStat, ErrMsg, 'finite-EI node-value recovery failed: '//TRIM(em)) END SUBROUTINE cable_node_values @@ -34147,31 +34502,39 @@ SUBROUTINE cable_nodal_curvature(cable, inode, curvature, ErrStat, ErrMsg) END IF END SUBROUTINE cable_nodal_curvature - SUBROUTINE cable_nodal_bend_moment(cable, inode, bend_moment, ErrStat, ErrMsg) + SUBROUTINE cable_nodal_bend_moment(cable, inode, bend_moment, ErrStat, ErrMsg, sleeve_moment) !! Magnitudes of the two one-sided EI*kappa traces at a node. Taking their - !! maximum avoids inventing an averaged EI at a stepped section interface. + !! maximum avoids inventing an averaged EI at a stepped section interface. bend_moment is the + !! cable's own share EI_line kappa; sleeve_moment (optional) the share of a bend stiffener, + !! (EI - EI_line) kappa. TYPE(CD_HFMF_ModuleType), INTENT(IN) :: cable INTEGER, INTENT(IN) :: inode REAL(wp), INTENT(OUT) :: bend_moment INTEGER, INTENT(OUT) :: ErrStat CHARACTER(*), INTENT(OUT) :: ErrMsg + REAL(wp), INTENT(OUT), OPTIONAL :: sleeve_moment REAL(wp) :: curvature ErrStat = CD_HCABLE_OK ErrMsg = '' bend_moment = CD_ZERO + IF (PRESENT(sleeve_moment)) sleeve_moment = CD_ZERO IF (inode < 1 .OR. inode > cable%line%nn) THEN ErrStat = 1; ErrMsg = 'cable_nodal_bend_moment: node index out of range'; RETURN END IF IF (inode > 1) THEN CALL cable_element_curvature(cable, inode - 1, CD_ONE, curvature, ErrStat, ErrMsg) IF (ErrStat /= CD_HCABLE_OK) RETURN - bend_moment = MAX(bend_moment, cable%line%EI(inode - 1)*curvature) + bend_moment = MAX(bend_moment, cable%line%EI_line(inode - 1)*curvature) + IF (PRESENT(sleeve_moment)) sleeve_moment = MAX(sleeve_moment, & + (cable%line%EI(inode - 1) - cable%line%EI_line(inode - 1))*curvature) END IF IF (inode <= cable%line%ne) THEN CALL cable_element_curvature(cable, inode, CD_ZERO, curvature, ErrStat, ErrMsg) IF (ErrStat /= CD_HCABLE_OK) RETURN - bend_moment = MAX(bend_moment, cable%line%EI(inode)*curvature) + bend_moment = MAX(bend_moment, cable%line%EI_line(inode)*curvature) + IF (PRESENT(sleeve_moment)) sleeve_moment = MAX(sleeve_moment, & + (cable%line%EI(inode) - cable%line%EI_line(inode))*curvature) END IF END SUBROUTINE cable_nodal_bend_moment @@ -34219,7 +34582,7 @@ FUNCTION CD_Channel_Unit(ch) RESULT(unit) END SELECT ELSE IF (starts_with(lo, 'fairten') .OR. starts_with(lo, 'anchten') .OR. is_anchor_force_channel(lo)) THEN unit = '(N)' - ELSE IF (is_endpoint_angle_channel(lo)) THEN + ELSE IF (is_endpoint_angle_channel(lo) .OR. is_end_deflection_channel(lo)) THEN unit = '(deg)' ELSE IF (is_point_force_channel(lo)) THEN unit = '(N)' @@ -34230,7 +34593,7 @@ FUNCTION CD_Channel_Unit(ch) RESULT(unit) SELECT CASE (kind) CASE (3); unit = '(N)' CASE (5); unit = '(1/m)' - CASE (6); unit = '(N.m)' + CASE (6, 12); unit = '(N.m)' CASE (7, 8); unit = '(deg)' CASE (9); unit = '(N.m)' CASE (10, 11); unit = '(deg)' diff --git a/src/CableDyn_HermiteCableDynamic.f90 b/src/CableDyn_HermiteCableDynamic.f90 index edf6b65..05df1cc 100644 --- a/src/CableDyn_HermiteCableDynamic.f90 +++ b/src/CableDyn_HermiteCableDynamic.f90 @@ -216,6 +216,10 @@ MODULE CableDyn_HermiteCableDynamic INTEGER :: ne = 0, nn = 0, ndof = 0 INTEGER :: axial_quadrature_order = 4, bending_quadrature_order = 4 REAL(wp), ALLOCATABLE :: l0(:), EA(:), EI(:), rho_a(:), w(:) ! per element + ! The cable's own bending stiffness per element, for the outputs: EI less the sleeve of a + ! bend stiffener (STIFFENERS), and equal to EI on every element without one. The bend + ! moment outputs report the cable's share EI_line kappa and the sleeve's the rest. + REAL(wp), ALLOCATABLE :: EI_line(:) ! Optional axial Kelvin-Voigt damping. BA is the axial constitutive dashpot ! coefficient [N s] after resolving the deck's direct-BA / negative-zeta convention. ! It is integrated over reference arc length with the same Hermite kinematics and @@ -671,8 +675,9 @@ SUBROUTINE CD_HermiteCable_Dyn_Init(model, l0, EA, EI, rho_a, w, seed, fixed_dof model%ne = ne; model%nn = nn; model%ndof = ndof model%seabed_z = seabed_z; model%kn = kn - ALLOCATE (model%l0(ne), model%EA(ne), model%EI(ne), model%rho_a(ne), model%w(ne)) + ALLOCATE (model%l0(ne), model%EA(ne), model%EI(ne), model%EI_line(ne), model%rho_a(ne), model%w(ne)) model%l0 = l0; model%EA = EA; model%EI = EI; model%rho_a = rho_a; model%w = w + model%EI_line = EI model%ei_resolve = MAXVAL(ABS(w)*l0*l0)/12.0_wp ! gen-alpha constants @@ -5464,7 +5469,8 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_Output_Recovery(model, element_recovery) model%output_element_recovery = element_recovery END SUBROUTINE CD_HermiteCable_Dyn_Set_Output_Recovery - SUBROUTINE CD_HermiteCable_Dyn_Node_Values(model, node, tension, curvature, bend_moment, ErrStat, ErrMsg) + SUBROUTINE CD_HermiteCable_Dyn_Node_Values(model, node, tension, curvature, bend_moment, ErrStat, ErrMsg, & + sleeve_moment) !! Equilibrium node values at one node (internal order), from the two adjacent elements; !! the same values as CD_HermiteCable_Dyn_Node_Recovery returns for that node. The element end !! residuals are kept in the model for the element state they were evaluated at (see the rc_ @@ -5477,10 +5483,12 @@ SUBROUTINE CD_HermiteCable_Dyn_Node_Values(model, node, tension, curvature, bend REAL(wp), INTENT(OUT) :: tension, curvature, bend_moment INTEGER, INTENT(OUT) :: ErrStat CHARACTER(*), INTENT(OUT) :: ErrMsg + REAL(wp), INTENT(OUT), OPTIONAL :: sleeve_moment REAL(wp) :: g(12), fm(3), mm(3), fp(3), mp(3) ErrStat = CD_HCDYN_OK; ErrMsg = '' tension = CD_ZERO; curvature = CD_ZERO; bend_moment = CD_ZERO + IF (PRESENT(sleeve_moment)) sleeve_moment = CD_ZERO IF (.NOT. model%initialized) THEN ErrStat = CD_HCDYN_BADINPUT; ErrMsg = 'CD_HermiteCable_Dyn_Node_Values: model not initialised'; RETURN END IF @@ -5505,7 +5513,8 @@ SUBROUTINE CD_HermiteCable_Dyn_Node_Values(model, node, tension, curvature, bend IF (ErrStat /= CD_HCDYN_OK) RETURN fp = -g(1:3); mp = -g(4:6) END IF - CALL node_values(model, node, fm, mm, fp, mp, tension, curvature, bend_moment, ErrStat, ErrMsg) + CALL node_values(model, node, fm, mm, fp, mp, tension, curvature, bend_moment, ErrStat, ErrMsg, & + sleeve_moment) CONTAINS @@ -5561,8 +5570,10 @@ SUBROUTINE CD_HermiteCable_Dyn_Node_Recovery(model, ErrStat, ErrMsg, tension, cu !! oscillation of the stiff-EA axial strain. !! !! tension(j) = t_j . F_j (t_j the unit nodal tangent; a line end reports the end force - !! magnitude of CD_HermiteCable_Dyn_End_Force), bend_moment(j) = |m_j - (m_j . t_j) t_j|, - !! curvature(j) = bend_moment(j)/EI, the larger of the two sides at a section junction (the + !! magnitude of CD_HermiteCable_Dyn_End_Force), the section moment M_j = |m_j - (m_j . t_j) t_j|, + !! curvature(j) = M_j/EI, the larger of the two sides at a section junction, and + !! bend_moment(j) the cable's own share EI_line curvature (M_j itself where no bend stiffener + !! sleeves the element; see node_values) (the !! exact geometric trace on an element without bending stiffness, or with EI below ei_resolve, !! where the recovered moment error w L0^2/12 exceeds the elastic moment of a 1 m bend radius). !! segment_tension(e) is the @@ -5629,16 +5640,21 @@ PURE LOGICAL FUNCTION size_ok(x, n) RESULT(ok) END FUNCTION size_ok END SUBROUTINE CD_HermiteCable_Dyn_Node_Recovery - SUBROUTINE node_values(model, node, fm, mm, fp, mp, tension, curvature, bend_moment, ErrStat, ErrMsg) + SUBROUTINE node_values(model, node, fm, mm, fp, mp, tension, curvature, bend_moment, ErrStat, ErrMsg, & + sleeve_moment) !! Node values from the one-sided section force and moment before (fm, mm) and after !! (fp, mp) the node; a line end uses its only side (see CD_HermiteCable_Dyn_Node_Recovery). + !! bend_moment is the cable's own share of the section moment, EI_line kappa on the side of + !! the larger value (the whole section moment on a side without a bend stiffener), and + !! sleeve_moment the rest, the share a stiffener sleeve carries. TYPE(CD_HermiteCableDynType), INTENT(IN) :: model INTEGER, INTENT(IN) :: node REAL(wp), INTENT(IN) :: fm(3), mm(3), fp(3), mp(3) REAL(wp), INTENT(OUT) :: tension, curvature, bend_moment INTEGER, INTENT(OUT) :: ErrStat CHARACTER(*), INTENT(OUT) :: ErrMsg - REAL(wp) :: t(3), f(3), m(3), fend(3), qe(12), trace + REAL(wp), INTENT(OUT), OPTIONAL :: sleeve_moment + REAL(wp) :: t(3), f(3), m(3), fend(3), qe(12), trace, section_moment, cable_moment INTEGER :: e, side CHARACTER(200) :: em2 @@ -5658,13 +5674,19 @@ SUBROUTINE node_values(model, node, fm, mm, fp, mp, tension, curvature, bend_mom ELSE tension = DOT_PRODUCT(t, f) END IF - bend_moment = NORM2(m - DOT_PRODUCT(m, t)*t) + section_moment = NORM2(m - DOT_PRODUCT(m, t)*t) curvature = CD_ZERO + cable_moment = CD_ZERO DO side = 0, 1 e = node - 1 + side IF (e < 1 .OR. e > model%ne) CYCLE IF (model%EI(e) > CD_ZERO .AND. model%EI(e) >= model%ei_resolve) THEN - curvature = MAX(curvature, bend_moment/model%EI(e)) + curvature = MAX(curvature, section_moment/model%EI(e)) + IF (model%EI_line(e) == model%EI(e)) THEN + cable_moment = MAX(cable_moment, section_moment) + ELSE + cable_moment = MAX(cable_moment, model%EI_line(e)*(section_moment/model%EI(e))) + END IF ELSE CALL gather_elem(model, e, qe) CALL CD_HermiteCable_Curvature(qe, model%l0(e), REAL(1 - side, wp), trace, ErrStat, em2) @@ -5672,8 +5694,15 @@ SUBROUTINE node_values(model, node, fm, mm, fp, mp, tension, curvature, bend_mom ErrStat = CD_HCDYN_NOCONVERGE; ErrMsg = 'CD_HermiteCable_Dyn_Node_Recovery: '//TRIM(em2); RETURN END IF curvature = MAX(curvature, trace) + IF (model%EI_line(e) == model%EI(e)) THEN + cable_moment = MAX(cable_moment, section_moment) + ELSE + cable_moment = MAX(cable_moment, model%EI_line(e)*trace) + END IF END IF END DO + bend_moment = cable_moment + IF (PRESENT(sleeve_moment)) sleeve_moment = MAX(CD_ZERO, section_moment - cable_moment) END SUBROUTINE node_values PURE FUNCTION unit_tangent(model, node) RESULT(t) @@ -6569,6 +6598,7 @@ SUBROUTINE CD_HermiteCable_Dyn_End(model) IF (ALLOCATED(model%l0)) DEALLOCATE (model%l0) IF (ALLOCATED(model%EA)) DEALLOCATE (model%EA) IF (ALLOCATED(model%EI)) DEALLOCATE (model%EI) + IF (ALLOCATED(model%EI_line)) DEALLOCATE (model%EI_line) IF (ALLOCATED(model%rho_a)) DEALLOCATE (model%rho_a) IF (ALLOCATED(model%w)) DEALLOCATE (model%w) IF (ALLOCATED(model%BA)) DEALLOCATE (model%BA) diff --git a/tests/check_bend_stiffener.py b/tests/check_bend_stiffener.py new file mode 100644 index 0000000..c1e077a --- /dev/null +++ b/tests/check_bend_stiffener.py @@ -0,0 +1,264 @@ +# SPDX-License-Identifier: Apache-2.0 +"""Bend stiffener against the planar elastica (driver_bend_stiffener in CMakeLists.txt). + +A power cable leaves a hang-off clamped horizontal at End A through a 5 m polyurethane bend +stiffener (a 0.5 m cylindrical root and a 4.5 m cone) and hangs to a fixed End B 52 m away and +30 m lower. The driver solves the static equilibrium with the stiffener given as a STIFFENERS +table. The reference is an independent solution of the planar extensible elastica, + + x' = (1 + e) cos(theta), z' = (1 + e) sin(theta), theta' = M / EI(s), + M' = -(1 + e) (cos(theta) n_z - sin(theta) n_x), n' = (0, w(s)), + e = (n_x cos(theta) + n_z sin(theta)) / EA, + +with EI(s) the cable's EI plus E pi/64 (OD(s)^4 - ID^4) of the sleeve and w(s) the submerged +weight of cable and sleeve per unstretched length, solved by collocation (scipy solve_bvp) +with the clamped direction and position at End A, and the position and zero moment of the +pinned End B. Over the stiffener and the 3 m beyond it the cable's own bend moment (the +static profile's BendMoment, EI_cable kappa) must follow the reference within 2 % of its root +value, and the node positions within 1 cm; at the root and 2 m along, the cable's and the +stiffener's moments (BendMom1N + StfMom1N) must add up to the reference section moment +within 2 % of the root moment; and FairDefl1, the angle between the end force and the clamp +axis, must follow the reference within 0.05 deg. + +The curvature over the same 8 m must also follow OrcaFlex 11.6d, with the stiffener given as a +line attachment (tests/data/orcaflex_bend_stiffener_curvature.csv), within 1.5 % of the root +curvature; this part needs no scipy. + +The initialisation notes are checked too: none on this deck; the resolution note when the 5 m +stiffener spans 3 elements (2 m elements); the pinned-end note when End A has no END +CONNECTIONS row. + +Usage: check_bend_stiffener.py +Exit 77 (skipped) when scipy is not installed and the notes pass. +""" + +import math +import pathlib +import subprocess +import sys +import tempfile + +import numpy as np + +try: + from scipy.integrate import solve_bvp +except ImportError: # pragma: no cover - environment dependent + solve_bvp = None + +MOMENT_TOLERANCE = 0.02 +POSITION_TOLERANCE = 0.01 +ANGLE_TOLERANCE = 0.05 +ORCAFLEX_TOLERANCE = 0.015 +ORCAFLEX_CSV = pathlib.Path(__file__).resolve().parent / "data" / "orcaflex_bend_stiffener_curvature.csv" +WINDOW = 8.0 + +G = 9.80665 +RHO_W = 1025.0 +D_CABLE = 0.20 +M_CABLE = 60.0 +EA = 5.0e8 +EI_CABLE = 2.0e4 +LENGTH = 60.0 +END_B = (51.90, -30.0) +# root side arc, tip side arc, root OD, tip OD (outward from End A) +SEGMENTS = [(0.0, 0.5, 0.60, 0.60), (0.5, 5.0, 0.60, 0.22)] +ID_STF = 0.20 +E_STF = 5.0e7 +RHO_STF = 1150.0 + +DECK = f"""--------------------- CableDyn Input File ------------------------------------ +Bend stiffener check: hang-off clamped horizontal, static +--------------------- LINE TYPES --------------------------------------------- +TypeName Diam MassDenInAir EA BA/-zeta EI Cd_n Cd_t Ca_n Ca_t +(-) (m) (kg/m) (N) (N-s/-) (N-m^2) (-) (-) (-) (-) +cable {D_CABLE} {M_CABLE} {EA} -1.0 {EI_CABLE} 1.2 0.0 1.0 0.0 +--------------------- POINTS ------------------------------------------------- +ID Type X Y Z Mass Vol CdA Ca +(-) (-) (m) (m) (m) (kg) (m^3) (m^2) (-) +1 Coupled 0.0 0.0 -20.0 0 0 0 0 +2 Fixed {END_B[0]} 0.0 {-20.0 + END_B[1]} 0 0 0 0 +--------------------- LINES -------------------------------------------------- +ID NodeA NodeB Outputs +(-) (-) (-) (-) +1 1 2 - +--------------------- SECTIONS ----------------------------------------------- +LineID LineType Length NumSegs +(-) (-) (m) (-) +1 cable 8.0 160 +1 cable 52.0 52 +--------------------- STIFFENERS --------------------------------------------- +LineID End Length RootOD TipOD InnerD E Density +(-) (-) (m) (m) (m) (m) (Pa) (kg/m^3) +""" + "".join( + f"1 A {s1 - s0} {od0} {od1} {ID_STF} {E_STF} {RHO_STF}\n" for s0, s1, od0, od1 in SEGMENTS +) + """--------------------- END CONNECTIONS ---------------------------------------- +LineID End Stiffness EzX EzY EzZ +(-) (-) (N-m/rad) (-) (-) (-) +1 A Rigid 1.0 0.0 0.0 +--------------------- OPTIONS ------------------------------------------------ +9.80665 g +1025.0 rhoW +0.05 dtM +0.05 TMax +--------------------- OUTPUTS ------------------------------------------------ +"FairTen1" +"FairDefl1" +"BendMom1N1" +"StfMom1N1" +"BendMom1N41" +"StfMom1N41" +END +""" + + +def outer_diameter(s: np.ndarray) -> np.ndarray: + od = np.full_like(s, np.nan) + for s0, s1, od0, od1 in SEGMENTS: + inside = (s >= s0) & (s <= s1) + od[inside] = od0 + (s[inside] - s0) / (s1 - s0) * (od1 - od0) + return od + + +def section(s: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Bending stiffness and submerged weight per unstretched length at arc s from End A.""" + od = outer_diameter(s) + sleeve = ~np.isnan(od) + ei = np.full_like(s, EI_CABLE) + w = np.full_like(s, (M_CABLE - RHO_W * math.pi / 4.0 * D_CABLE**2) * G) + area = math.pi / 4.0 * (od[sleeve] ** 2 - ID_STF**2) + ei[sleeve] += E_STF * math.pi / 64.0 * (od[sleeve] ** 4 - ID_STF**4) + w[sleeve] += (RHO_STF - RHO_W) * area * G + return ei, w + + +FORCE_SCALE = 1.0e5 +MOMENT_SCALE = 1.0e4 + + +def reference(s_nodes: np.ndarray, x_seed: np.ndarray, z_seed: np.ndarray): + """Planar extensible elastica along the unstretched arc, y = (x, z, theta, M, n_x, n_z), with + the moment and the force scaled by MOMENT_SCALE and FORCE_SCALE. The seed is the shape only: + the driver's node positions; moment zero and a tension of FORCE_SCALE along the line.""" + + def rhs(s, y): + x, z, th, m, nx, nz = y + ei, w = section(s) + c, sn = np.cos(th), np.sin(th) + stretch = 1.0 + FORCE_SCALE * (nx * c + nz * sn) / EA + return np.vstack([stretch * c, stretch * sn, MOMENT_SCALE * m / ei, + -stretch * FORCE_SCALE / MOMENT_SCALE * (c * nz - sn * nx), + np.zeros_like(s), w / FORCE_SCALE]) + + def bc(ya, yb): + return np.array([ya[0], ya[1], ya[2], yb[0] - END_B[0], yb[1] - END_B[1], yb[3]]) + + # a dense start mesh: the collocation residual is met at 1e-4 without refinement, which + # would split the mesh at the steps of the sleeve's weight and stiffness at its tip + mesh = np.unique(np.round(np.concatenate([s_nodes, np.linspace(0.0, WINDOW, 1601), + np.linspace(0.0, LENGTH, 1201)]), 9)) + xs = np.interp(mesh, s_nodes, x_seed) + zs = np.interp(mesh, s_nodes, z_seed) + th = np.unwrap(np.arctan2(np.gradient(zs, mesh), np.gradient(xs, mesh))) + y0 = np.vstack([xs, zs, th, np.zeros_like(mesh), np.cos(th[-1]) * np.ones_like(mesh), np.sin(th)]) + sol = solve_bvp(rhs, bc, mesh, y0, tol=1.0e-4, max_nodes=50000) + if not sol.success: + raise RuntimeError(f"reference elastica did not converge: {sol.message}") + return sol + + +RESOLUTION_NOTE = "stiffener at End A spans" +PINNED_NOTE = "the stiffener at End A sleeves a pinned end" + + +def solve(exe: str, tmp: str, name: str, deck_text: str): + """Run the driver on a deck; return (exit code, stdout + stderr, static profile rows, the + output channels at t = 0 by name).""" + deck = pathlib.Path(tmp) / f"{name}.dat" + deck.write_text(deck_text) + root = str(pathlib.Path(tmp) / name) + run = subprocess.run([exe, str(deck), root], capture_output=True, text=True) + static = pathlib.Path(root + ".static.out") + rows = static.read_text().splitlines()[3:] if static.exists() else [] + channels = {} + series = pathlib.Path(root + ".out") + if series.exists(): + lines = series.read_text().splitlines() + channels = dict(zip(lines[1].split(), (float(v) for v in lines[2].split()))) + return run.returncode, run.stdout + run.stderr, rows, channels + + +def check_notes(exe: str, tmp: str, log: str) -> bool: + ok = RESOLUTION_NOTE not in log and PINNED_NOTE not in log + print(f"{'ok ' if ok else 'FAIL'} no stiffener note on the resolved, clamped deck") + coarse = DECK.replace("1 cable 8.0 160", "1 cable 8.0 4") + rc, out, _, _ = solve(exe, tmp, "coarse", coarse) + good = rc == 0 and RESOLUTION_NOTE in out and "spans 3 elements" in out + print(f"{'ok ' if good else 'FAIL'} resolution note with 2 m elements") + ok = ok and good + start = DECK.index("--------------------- END CONNECTIONS") + pinned = DECK[:start] + DECK[DECK.index("--------------------- OPTIONS"):] + # (the deflection angle needs a held end) + rc, out, _, _ = solve(exe, tmp, "pinned", pinned.replace('"FairDefl1"\n', "")) + good = rc == 0 and PINNED_NOTE in out + print(f"{'ok ' if good else 'FAIL'} pinned-end note without END CONNECTIONS") + return ok and good + + +def check_orcaflex(arc: np.ndarray, curvature: np.ndarray) -> bool: + """The curvature over the first WINDOW m against OrcaFlex 11.6d with its stiffener given as a + line attachment (tests/data/orcaflex_bend_stiffener_curvature.csv, from + validation/scripts/orcaflex_bend_stiffener.py), within ORCAFLEX_TOLERANCE of the root value.""" + table = np.loadtxt(ORCAFLEX_CSV, delimiter=",", comments="#", skiprows=2) + window = arc <= WINDOW + reference = np.interp(arc[window], table[:, 0], table[:, 1]) + error = np.max(np.abs(curvature[window] - reference)) / table[0, 1] + ok = error <= ORCAFLEX_TOLERANCE + print(f"{'ok ' if ok else 'FAIL'} curvature against OrcaFlex 11.6d: largest difference " + f"{100 * error:.2f} % of the root curvature {table[0, 1]:.4f} 1/m") + return ok + + +def main() -> int: + exe = sys.argv[1] + with tempfile.TemporaryDirectory() as tmp: + rc, log, rows, out = solve(exe, tmp, "stiffener", DECK) + if rc != 0: + print(log[-3000:]) + print("FAIL: the driver did not solve the stiffener deck") + return 1 + notes_ok = check_notes(exe, tmp, log) + data = np.array([[float(v) for v in r.split()[2:9]] for r in rows if r.strip()]) + orcaflex_ok = check_orcaflex(data[:, 0], data[:, 5]) + if solve_bvp is None: + print("SKIP: scipy is not installed (the elastica comparison needs solve_bvp)") + return 77 if notes_ok and orcaflex_ok else 1 + arc, x, z, moment = data[:, 0], data[:, 1], data[:, 3] + 20.0, data[:, 6] + s_nodes = np.cumsum(np.r_[0.0, np.full(160, 8.0 / 160), np.full(52, 1.0)]) + sol = reference(s_nodes, x, z) + ref = sol.sol(s_nodes) + ref[3] *= MOMENT_SCALE + section_moment = np.abs(ref[3]) + cable_moment = EI_CABLE * section_moment / section(s_nodes)[0] + window = s_nodes <= WINDOW + cable_error = np.max(np.abs(np.abs(moment[window]) - cable_moment[window])) / cable_moment[0] + position_error = np.max(np.hypot(x - ref[0], z - ref[1])) + print(f"root cable moment: driver {moment[0]:.1f} N m, elastica {cable_moment[0]:.1f} N m") + print(f"largest cable-moment difference over the first {WINDOW} m: {100 * cable_error:.2f} % of the root value") + print(f"largest node position difference: {1000 * position_error:.2f} mm") + total_error = 0.0 + for node in (1, 41): + total = out[f"BendMom1N{node}"] + out[f"StfMom1N{node}"] + total_error = max(total_error, abs(total - section_moment[node - 1]) / section_moment[0]) + print(f"node {node}: cable + stiffener {total:.1f} N m, elastica section moment {section_moment[node - 1]:.1f} N m") + end_force = sol.sol(np.array([0.0]))[4:6, 0] + deflection = math.degrees(math.atan2(abs(end_force[1]), end_force[0])) + angle_error = abs(out["FairDefl1"] - deflection) + print(f"FairDefl1: driver {out['FairDefl1']:.4f} deg, elastica {deflection:.4f} deg") + ok = (cable_error <= MOMENT_TOLERANCE and position_error <= POSITION_TOLERANCE and notes_ok + and total_error <= MOMENT_TOLERANCE and angle_error <= ANGLE_TOLERANCE and orcaflex_ok) + print("PASS" if ok else "FAIL") + return 0 if ok else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/data/orcaflex_bend_stiffener_curvature.csv b/tests/data/orcaflex_bend_stiffener_curvature.csv new file mode 100644 index 0000000..8204179 --- /dev/null +++ b/tests/data/orcaflex_bend_stiffener_curvature.csv @@ -0,0 +1,163 @@ +# OrcFxAPI 11.6d static curvature along the hang-off (validation/scripts/orcaflex_bend_stiffener.py) +arc_m,curvature_per_m +0.000000,4.52202935e-01 +0.025000,4.42374765e-01 +0.075000,4.23063500e-01 +0.125000,4.04431762e-01 +0.175000,3.86458130e-01 +0.225000,3.69120886e-01 +0.275000,3.52398148e-01 +0.325000,3.36267975e-01 +0.375000,3.20708463e-01 +0.425000,3.05697832e-01 +0.475000,2.91214188e-01 +0.525000,2.80981913e-01 +0.575000,2.74616508e-01 +0.625000,2.68251090e-01 +0.675000,2.61890066e-01 +0.725000,2.55537806e-01 +0.775000,2.49198641e-01 +0.825000,2.42876863e-01 +0.875000,2.36576718e-01 +0.925000,2.30302408e-01 +0.975000,2.24058084e-01 +1.025000,2.17847846e-01 +1.075000,2.11675739e-01 +1.125000,2.05545752e-01 +1.175000,1.99461812e-01 +1.225000,1.93427784e-01 +1.275000,1.87447466e-01 +1.325000,1.81524586e-01 +1.375000,1.75662800e-01 +1.425000,1.69865690e-01 +1.475000,1.64136755e-01 +1.525000,1.58479413e-01 +1.575000,1.52896996e-01 +1.625000,1.47392744e-01 +1.675000,1.41969803e-01 +1.725000,1.36631219e-01 +1.775000,1.31379938e-01 +1.825000,1.26218795e-01 +1.875000,1.21150517e-01 +1.925000,1.16177711e-01 +1.975000,1.11302864e-01 +2.025000,1.06528339e-01 +2.075000,1.01856364e-01 +2.125000,9.72890354e-02 +2.175000,9.28283055e-02 +2.225000,8.84759825e-02 +2.275000,8.42337235e-02 +2.325000,8.01030302e-02 +2.375000,7.60852439e-02 +2.425000,7.21815410e-02 +2.475000,6.83929286e-02 +2.525000,6.47202403e-02 +2.575000,6.11641320e-02 +2.625000,5.77250786e-02 +2.675000,5.44033709e-02 +2.725000,5.11991123e-02 +2.775000,4.81122169e-02 +2.825000,4.51424075e-02 +2.875000,4.22892147e-02 +2.925000,3.95519762e-02 +2.975000,3.69298370e-02 +3.025000,3.44217505e-02 +3.075000,3.20264801e-02 +3.125000,2.97426016e-02 +3.175000,2.75685072e-02 +3.225000,2.55024093e-02 +3.275000,2.35423460e-02 +3.325000,2.16861872e-02 +3.375000,1.99316415e-02 +3.425000,1.82762645e-02 +3.475000,1.67174676e-02 +3.525000,1.52525274e-02 +3.575000,1.38785968e-02 +3.625000,1.25927159e-02 +3.675000,1.13918242e-02 +3.725000,1.02727733e-02 +3.775000,9.23234047e-03 +3.825000,8.26724185e-03 +3.875000,7.37414703e-03 +3.925000,6.54969332e-03 +3.975000,5.79050032e-03 +4.025000,5.09318450e-03 +4.075000,4.45437373e-03 +4.125000,3.87072148e-03 +4.175000,3.33892068e-03 +4.225000,2.85571710e-03 +4.275000,2.41792197e-03 +4.325000,2.02242392e-03 +4.375000,1.66619989e-03 +4.425000,1.34632515e-03 +4.475000,1.05998206e-03 +4.525000,8.04467736e-04 +4.575000,5.77200370e-04 +4.625000,3.75724273e-04 +4.675000,1.97713570e-04 +4.725000,4.09745402e-05 +4.775000,9.65533700e-05 +4.825000,2.16797819e-04 +4.875000,3.21555530e-04 +4.925000,4.12480920e-04 +4.975000,4.91672140e-04 +5.025000,6.05407857e-04 +5.075000,6.57380692e-04 +5.125000,7.00952817e-04 +5.175000,7.37483590e-04 +5.225000,7.68112385e-04 +5.275000,7.93794202e-04 +5.325000,8.15329500e-04 +5.375000,8.33389211e-04 +5.425000,8.48535698e-04 +5.475000,8.61240324e-04 +5.525000,8.71898179e-04 +5.575000,8.80840429e-04 +5.625000,8.88344649e-04 +5.675000,8.94643502e-04 +5.725000,8.99932005e-04 +5.775000,9.04373624e-04 +5.825000,9.08105373e-04 +5.875000,9.11242098e-04 +5.925000,9.13880065e-04 +5.975000,9.16099965e-04 +6.025000,9.17969431e-04 +6.075000,9.19545159e-04 +6.125000,9.20874668e-04 +6.175000,9.21997791e-04 +6.225000,9.22947919e-04 +6.275000,9.23753040e-04 +6.325000,9.24436612e-04 +6.375000,9.25018304e-04 +6.425000,9.25514600e-04 +6.475000,9.25939319e-04 +6.525000,9.26304046e-04 +6.575000,9.26618494e-04 +6.625000,9.26890805e-04 +6.675000,9.27127807e-04 +6.725000,9.27335224e-04 +6.775000,9.27517858e-04 +6.825000,9.27679737e-04 +6.875000,9.27824243e-04 +6.925000,9.27954213e-04 +6.975000,9.28072031e-04 +7.025000,9.28179701e-04 +7.075000,9.28278914e-04 +7.125000,9.28371095e-04 +7.175000,9.28457446e-04 +7.225000,9.28538991e-04 +7.275000,9.28616604e-04 +7.325000,9.28691036e-04 +7.375000,9.28762937e-04 +7.425000,9.28832879e-04 +7.475000,9.28901374e-04 +7.525000,9.28968886e-04 +7.575000,9.29035852e-04 +7.625000,9.29102688e-04 +7.675000,9.29169809e-04 +7.725000,9.29237639e-04 +7.775000,9.29306619e-04 +7.825000,9.29377234e-04 +7.875000,9.29450012e-04 +7.925000,9.29525554e-04 +7.975000,9.29604543e-04 diff --git a/tests/test_bend_stiffener.f90 b/tests/test_bend_stiffener.f90 new file mode 100644 index 0000000..ae03150 --- /dev/null +++ b/tests/test_bend_stiffener.f90 @@ -0,0 +1,332 @@ +! File: tests/test_bend_stiffener.f90 +! SPDX-License-Identifier: Apache-2.0 +! Copyright (c) 2026 Jae Hoon Seo, SMI Lab, Inha University +PROGRAM test_bend_stiffener + !! Gates for bend stiffeners (STIFFENERS) on a finite-EI line: + !! 1. Element properties (CD_Stiffener_Overlay) against closed forms: a cylindrical sleeve + !! over half an element gives the compliance average of the bending stiffness, the mean + !! sleeve mass, the displaced-area diameter and the mean outer diameter; the same sleeve + !! at End B lands on the mirrored element; elements it does not reach are unchanged bit + !! for bit. + !! 2. A cone: the element values follow a 20 000-interval midpoint integration of the same + !! averages to 1e-10. + !! 3. Deck grammar: malformed rows, a stiffener on an EI = 0 line and stiffeners longer than + !! the line fail closed by name; on a stock-order LINES row the stiffener sleeves the same + !! physical end, and the static profile is the native deck's; StfMom on a line without a + !! stiffener and AnchDefl at a pinned end fail closed, FairDefl at the clamped end runs. + !! 4. Coupled route (the OpenFAST aggregate): the cable carries the stiffened root element, + !! and its initial fairlead load is the standalone static solve's. + !! The static equilibrium against the planar elastica is tests/check_bend_stiffener.py. + USE CableDyn_Precision, ONLY: wp + USE CableDyn_DeckDriver, ONLY: CD_Run_Deck_Driver, CD_DECKDRV_OK, CD_Stiffener_Overlay + USE CableDyn_OpenFAST_Aggregate, ONLY: CD_AGG_ModuleType, CD_AGG_Init_From_Deck, CD_AGG_GetMovingPointMesh, & + CD_AGG_End, CD_AGG_OK + IMPLICIT NONE + + REAL(wp), PARAMETER :: PI = 3.14159265358979323846_wp + INTEGER :: nfail + nfail = 0 + + CALL case_cylinder() + CALL case_cone() + CALL case_grammar() + CALL case_coupled() + + IF (nfail > 0) THEN + WRITE (*, '(A,I0,A)') 'test_bend_stiffener: ', nfail, ' check(s) failed' + ERROR STOP 1 + END IF + WRITE (*, '(A)') 'test_bend_stiffener: all checks passed' + +CONTAINS + + SUBROUTINE require(cond, msg) + LOGICAL, INTENT(IN) :: cond + CHARACTER(*), INTENT(IN) :: msg + IF (.NOT. cond) THEN + nfail = nfail + 1 + WRITE (*, '(A)') 'FAIL: '//TRIM(msg) + END IF + END SUBROUTINE require + + LOGICAL FUNCTION close_rel(a, b, tol) + REAL(wp), INTENT(IN) :: a, b, tol + close_rel = ABS(a - b) <= tol*MAX(ABS(a), ABS(b)) + END FUNCTION close_rel + + SUBROUTINE case_cylinder() + !! Four 1 m elements; a 1.5 m cylindrical sleeve (OD 0.5 m, ID 0.2 m) covers element 1 and + !! half of element 2 from End A (or from End B: elements 4 and 3). + REAL(wp), PARAMETER :: EI_C = 2.0e4_wp, M_C = 60.0_wp, D_C = 0.2_wp + REAL(wp), PARAMETER :: OD = 0.5_wp, ID = 0.2_wp, EMOD = 5.0e7_wp, RHO = 1150.0_wp + REAL(wp) :: stf(8, 1), l0(4), ei(4), mpl(4), dv(4), dh(4), ei_s, ei_half, m_s + INTEGER :: iend, e_full, e_half, e_far(2) + l0 = 1.0_wp + ei_s = EMOD*PI/64.0_wp*(OD**4 - ID**4) + m_s = RHO*PI/4.0_wp*(OD**2 - ID**2) + ei_half = 1.0_wp/(0.5_wp/EI_C + 0.5_wp/(EI_C + ei_s)) + DO iend = 1, 2 + stf(:, 1) = [REAL(iend, wp), 0.0_wp, 1.5_wp, OD, OD, ID, EMOD, RHO] + ei = EI_C; mpl = M_C; dv = D_C; dh = D_C + CALL CD_Stiffener_Overlay(stf, l0, dv, ei=ei, mpl=mpl, volume_diameter=.TRUE.) + CALL CD_Stiffener_Overlay(stf, l0, dh, volume_diameter=.FALSE.) + IF (iend == 1) THEN + e_full = 1; e_half = 2; e_far = [3, 4] + ELSE + e_full = 4; e_half = 3; e_far = [1, 2] + END IF + CALL require(close_rel(ei(e_full), EI_C + ei_s, 1.0e-13_wp), 'cylinder: covered element EI') + CALL require(close_rel(ei(e_half), ei_half, 1.0e-13_wp), 'cylinder: half-covered element takes the '// & + 'compliance average') + CALL require(close_rel(mpl(e_full), M_C + m_s, 1.0e-13_wp), 'cylinder: covered element mass') + CALL require(close_rel(mpl(e_half), M_C + 0.5_wp*m_s, 1.0e-13_wp), 'cylinder: half-covered mass') + CALL require(close_rel(dv(e_full), SQRT(D_C**2 + OD**2 - ID**2), 1.0e-13_wp), & + 'cylinder: displaced-area diameter of cable plus solid sleeve') + CALL require(close_rel(dv(e_half), SQRT(D_C**2 + 0.5_wp*(OD**2 - ID**2)), 1.0e-13_wp), & + 'cylinder: half-covered displaced-area diameter') + CALL require(close_rel(dh(e_full), OD, 1.0e-13_wp), 'cylinder: outer diameter for drag') + CALL require(close_rel(dh(e_half), 0.5_wp*(OD + D_C), 1.0e-13_wp), 'cylinder: half-covered outer diameter') + CALL require(ALL(ei(e_far) == EI_C) .AND. ALL(mpl(e_far) == M_C) .AND. ALL(dv(e_far) == D_C) .AND. & + ALL(dh(e_far) == D_C), 'cylinder: elements the sleeve does not reach are unchanged') + END DO + END SUBROUTINE case_cylinder + + SUBROUTINE case_cone() + !! A 2.5 m cone (OD 0.6 -> 0.22 m) from End A over five 0.7 m elements: the fourth element + !! (2.1 - 2.8 m) holds the cone's tip; the fifth is bare. + REAL(wp), PARAMETER :: EI_C = 2.0e4_wp, M_C = 60.0_wp, D_C = 0.2_wp + REAL(wp), PARAMETER :: OD0 = 0.6_wp, OD1 = 0.22_wp, ID = 0.2_wp, EMOD = 5.0e7_wp, RHO = 1150.0_wp + INTEGER, PARAMETER :: NQ = 20000 + REAL(wp) :: stf(8, 1), l0(5), ei(5), mpl(5), dv(5), dh(5), a, s, od, comp, mass, area, outer, h, tol + INTEGER :: e, k + l0 = 0.7_wp + stf(:, 1) = [1.0_wp, 0.0_wp, 2.5_wp, OD0, OD1, ID, EMOD, RHO] + ei = EI_C; mpl = M_C; dv = D_C; dh = D_C + CALL CD_Stiffener_Overlay(stf, l0, dv, ei=ei, mpl=mpl, volume_diameter=.TRUE.) + CALL CD_Stiffener_Overlay(stf, l0, dh, volume_diameter=.FALSE.) + DO e = 1, 4 + a = 0.7_wp*REAL(e - 1, wp) + h = 0.7_wp/REAL(NQ, wp) + comp = 0.0_wp; mass = 0.0_wp; area = 0.0_wp; outer = 0.0_wp + DO k = 1, NQ + s = a + (REAL(k, wp) - 0.5_wp)*h + IF (s < 2.5_wp) THEN + od = OD0 + s/2.5_wp*(OD1 - OD0) + comp = comp + h/(EI_C + EMOD*PI/64.0_wp*(od**4 - ID**4)) + mass = mass + h*RHO*PI/4.0_wp*(od**2 - ID**2) + area = area + h*(od**2 - ID**2) + outer = outer + h*MAX(od, D_C) + ELSE + comp = comp + h/EI_C + outer = outer + h*D_C + END IF + END DO + ! the cone tip inside element 4 limits the midpoint rule there to O(h/l) = 5e-5 + tol = MERGE(2.0e-4_wp, 1.0e-6_wp, e == 4) + CALL require(close_rel(ei(e), 0.7_wp/comp, tol), 'cone: compliance-average EI of element') + CALL require(close_rel(mpl(e), M_C + mass/0.7_wp, tol), 'cone: mean mass of element') + CALL require(close_rel(dv(e), SQRT(D_C**2 + area/0.7_wp), tol), 'cone: displaced-area diameter') + CALL require(close_rel(dh(e), outer/0.7_wp, tol), 'cone: mean outer diameter') + IF (e < 4) CALL require(close_rel(ei(e), 0.7_wp/comp, 1.0e-10_wp), & + 'cone: Gauss average matches the fine integration inside the cone') + END DO + CALL require(ei(5) == EI_C .AND. mpl(5) == M_C .AND. dv(5) == D_C .AND. dh(5) == D_C, & + 'cone: bare elements unchanged') + END SUBROUTINE case_cone + + SUBROUTINE write_deck(path, stiffener_rows, ei, stock, channel) + !! A 60 m cable clamped horizontal at a 20 m deep hang-off and fixed 30 m lower; stock: the + !! 7-column LINES row written anchor first, so the parser swaps the ends. + CHARACTER(*), INTENT(IN) :: path + CHARACTER(*), INTENT(IN) :: stiffener_rows(:) + REAL(wp), INTENT(IN) :: ei + LOGICAL, INTENT(IN) :: stock + CHARACTER(*), INTENT(IN), OPTIONAL :: channel + INTEGER :: u, k + OPEN (NEWUNIT=u, FILE=path, STATUS='REPLACE', ACTION='WRITE') + WRITE (u, '(A)') '--------------------- CableDyn Input File ------------------------------------' + WRITE (u, '(A)') 'Bend stiffener test deck' + WRITE (u, '(A)') '--------------------- LINE TYPES ---------------------------------------------' + WRITE (u, '(A)') 'TypeName Diam MassDenInAir EA BA/-zeta EI Cd_n Cd_t Ca_n Ca_t' + WRITE (u, '(A)') '(-) (m) (kg/m) (N) (N-s/-) (N-m^2) (-) (-) (-) (-)' + WRITE (u, '(A,ES12.4,A)') 'cable 0.20 60.0 5.0e8 -1.0 ', ei, ' 1.2 0.0 1.0 0.0' + WRITE (u, '(A)') '--------------------- POINTS -------------------------------------------------' + WRITE (u, '(A)') 'ID Type X Y Z Mass Vol CdA Ca' + WRITE (u, '(A)') '(-) (-) (m) (m) (m) (kg) (m^3) (m^2) (-)' + WRITE (u, '(A)') '1 Coupled 0.0 0.0 -20.0 0 0 0 0' + WRITE (u, '(A)') '2 Fixed 51.90 0.0 -50.0 0 0 0 0' + WRITE (u, '(A)') '--------------------- LINES --------------------------------------------------' + IF (stock) THEN + WRITE (u, '(A)') 'ID LineType AttachA AttachB UnstrLen NumSegs Outputs' + WRITE (u, '(A)') '(-) (-) (-) (-) (m) (-) (-)' + WRITE (u, '(A)') '1 cable 2 1 60.0 120 -' + ELSE + WRITE (u, '(A)') 'ID NodeA NodeB Outputs' + WRITE (u, '(A)') '(-) (-) (-) (-)' + WRITE (u, '(A)') '1 1 2 -' + WRITE (u, '(A)') '--------------------- SECTIONS -----------------------------------------------' + WRITE (u, '(A)') 'LineID LineType Length NumSegs' + WRITE (u, '(A)') '(-) (-) (m) (-)' + WRITE (u, '(A)') '1 cable 60.0 120' + END IF + IF (SIZE(stiffener_rows) > 0) THEN + WRITE (u, '(A)') '--------------------- STIFFENERS ---------------------------------------------' + WRITE (u, '(A)') 'LineID End Length RootOD TipOD InnerD E Density' + WRITE (u, '(A)') '(-) (-) (m) (m) (m) (m) (Pa) (kg/m^3)' + DO k = 1, SIZE(stiffener_rows) + WRITE (u, '(A)') TRIM(stiffener_rows(k)) + END DO + END IF + IF (ei > 0.0_wp) THEN + WRITE (u, '(A)') '--------------------- END CONNECTIONS ----------------------------------------' + WRITE (u, '(A)') 'LineID End Stiffness EzX EzY EzZ' + WRITE (u, '(A)') '(-) (-) (N-m/rad) (-) (-) (-)' + IF (stock) THEN + WRITE (u, '(A)') '1 B Rigid -1.0 0.0 0.0' + ELSE + WRITE (u, '(A)') '1 A Rigid 1.0 0.0 0.0' + END IF + END IF + WRITE (u, '(A)') '--------------------- OPTIONS ------------------------------------------------' + WRITE (u, '(A)') '9.80665 g' + WRITE (u, '(A)') '1025.0 rhoW' + WRITE (u, '(A)') '0.05 dtM' + WRITE (u, '(A)') '0.0 TMax' + WRITE (u, '(A)') '--------------------- OUTPUTS ------------------------------------------------' + WRITE (u, '(A)') '"FairTen1"' + IF (PRESENT(channel)) WRITE (u, '(A)') '"'//channel//'"' + WRITE (u, '(A)') 'END' + CLOSE (u) + END SUBROUTINE write_deck + + SUBROUTINE run(deck, root, ok, msg) + CHARACTER(*), INTENT(IN) :: deck, root + LOGICAL, INTENT(OUT) :: ok + CHARACTER(*), INTENT(OUT) :: msg + LOGICAL :: conv + INTEGER :: es + CALL CD_Run_Deck_Driver(deck, root, conv, es, msg) + ok = es == CD_DECKDRV_OK .AND. conv + END SUBROUTINE run + + SUBROUTINE read_static(path, prof, n) + !! Columns 3-9 of a .static.out (arc, X, Y, Z, tension, curvature, moment), node by node. + CHARACTER(*), INTENT(IN) :: path + REAL(wp), ALLOCATABLE, INTENT(OUT) :: prof(:, :) + INTEGER, INTENT(OUT) :: n + CHARACTER(1024) :: rec + INTEGER :: u, ios, idum(2) + ALLOCATE (prof(7, 2000)) + n = 0 + OPEN (NEWUNIT=u, FILE=path, STATUS='OLD', ACTION='READ', IOSTAT=ios) + IF (ios /= 0) RETURN + READ (u, '(A)') rec + READ (u, '(A)') rec + READ (u, '(A)') rec + DO + READ (u, '(A)', IOSTAT=ios) rec + IF (ios /= 0) EXIT + IF (LEN_TRIM(rec) == 0) CYCLE + n = n + 1 + READ (rec, *) idum, prof(:, n) + END DO + CLOSE (u) + END SUBROUTINE read_static + + SUBROUTINE case_grammar() + CHARACTER(80) :: rows(2), none(0) + CHARACTER(512) :: msg + LOGICAL :: ok + REAL(wp), ALLOCATABLE :: pa(:, :), pb(:, :) + INTEGER :: na, nb + rows(1) = '1 A 0.5 0.60 0.60 0.20 5.0e7 1150' + rows(2) = '1 A 4.5 0.60 0.22 0.20 5.0e7 1150' + CALL write_deck('stf_native.dat', rows, 2.0e4_wp, .FALSE.) + CALL run('stf_native.dat', 'stf_native', ok, msg) + CALL require(ok, 'grammar: native deck solves: '//TRIM(msg)) + ! stock order: the deck's End B is the hang-off; the stiffener sleeves the same physical end + rows(1) = '1 B 0.5 0.60 0.60 0.20 5.0e7 1150' + rows(2) = '1 B 4.5 0.60 0.22 0.20 5.0e7 1150' + CALL write_deck('stf_stock.dat', rows, 2.0e4_wp, .TRUE.) + CALL run('stf_stock.dat', 'stf_stock', ok, msg) + CALL require(ok, 'grammar: stock-order deck solves: '//TRIM(msg)) + CALL read_static('stf_native.static.out', pa, na) + CALL read_static('stf_stock.static.out', pb, nb) + CALL require(na == nb .AND. na > 0, 'grammar: both decks write a static profile') + IF (na == nb .AND. na > 0) THEN + CALL require(MAXVAL(ABS(pa(2:4, :na) - pb(2:4, :nb))) < 1.0e-6_wp .AND. & + MAXVAL(ABS(pa(7, :na) - pb(7, :nb))) < 1.0e-6_wp*MAXVAL(ABS(pa(7, :na))), & + 'grammar: on a stock-order row the stiffener sleeves the same physical end') + END IF + rows(1) = '1 A 0.5 0.60 0.60 0.20 5.0e7' + CALL write_deck('stf_bad0.dat', rows(1:1), 2.0e4_wp, .FALSE.) + CALL run('stf_bad0.dat', 'stf_bad0', ok, msg) + CALL require(.NOT. ok .AND. INDEX(msg, 'LineID End Length') > 0, 'grammar: a seven-column row: '//TRIM(msg)) + rows(1) = '1 C 0.5 0.60 0.60 0.20 5.0e7 1150' + CALL write_deck('stf_bad1.dat', rows(1:1), 2.0e4_wp, .FALSE.) + CALL run('stf_bad1.dat', 'stf_bad1', ok, msg) + CALL require(.NOT. ok .AND. INDEX(msg, 'End must be A or B') > 0, 'grammar: an unknown end: '//TRIM(msg)) + rows(1) = '1 A 0.5 0.60 0.60 0.70 5.0e7 1150' + CALL write_deck('stf_bad2.dat', rows(1:1), 2.0e4_wp, .FALSE.) + CALL run('stf_bad2.dat', 'stf_bad2', ok, msg) + CALL require(.NOT. ok .AND. INDEX(msg, 'InnerD must be smaller') > 0, 'grammar: InnerD above OD: '//TRIM(msg)) + rows(1) = '1 A 0.5 0.60 0.60 0.20 0.0 1150' + CALL write_deck('stf_bad3.dat', rows(1:1), 2.0e4_wp, .FALSE.) + CALL run('stf_bad3.dat', 'stf_bad3', ok, msg) + CALL require(.NOT. ok .AND. INDEX(msg, 'must be > 0') > 0, 'grammar: a zero modulus: '//TRIM(msg)) + rows(1) = '2 A 0.5 0.60 0.60 0.20 5.0e7 1150' + CALL write_deck('stf_bad4.dat', rows(1:1), 2.0e4_wp, .FALSE.) + CALL run('stf_bad4.dat', 'stf_bad4', ok, msg) + CALL require(.NOT. ok .AND. INDEX(msg, 'undefined LINE') > 0, 'grammar: an unknown line: '//TRIM(msg)) + rows(1) = '1 A 40.0 0.60 0.60 0.20 5.0e7 1150' + rows(2) = '1 B 25.0 0.60 0.60 0.20 5.0e7 1150' + CALL write_deck('stf_bad5.dat', rows, 2.0e4_wp, .FALSE.) + CALL run('stf_bad5.dat', 'stf_bad5', ok, msg) + CALL require(.NOT. ok .AND. INDEX(msg, 'longer than the line') > 0, 'grammar: stiffeners overlap: '//TRIM(msg)) + rows(1) = '1 A 0.5 0.60 0.60 0.20 5.0e7 1150' + CALL write_deck('stf_bad6.dat', rows(1:1), 0.0_wp, .FALSE.) + CALL run('stf_bad6.dat', 'stf_bad6', ok, msg) + CALL require(.NOT. ok .AND. INDEX(msg, 'not a finite-EI line') > 0, 'grammar: an EI = 0 line: '//TRIM(msg)) + ! no STIFFENERS table: the same deck still solves (the table is optional) + CALL write_deck('stf_none.dat', none, 2.0e4_wp, .FALSE.) + CALL run('stf_none.dat', 'stf_none', ok, msg) + CALL require(ok, 'grammar: the deck without stiffeners solves: '//TRIM(msg)) + ! the stiffener moment needs a stiffener; the deflection angle a held end + CALL write_deck('stf_ch1.dat', none, 2.0e4_wp, .FALSE., channel='StfMom1N3') + CALL run('stf_ch1.dat', 'stf_ch1', ok, msg) + CALL require(.NOT. ok .AND. INDEX(msg, 'has no STIFFENERS') > 0, 'grammar: StfMom without a stiffener: '//TRIM(msg)) + rows(1) = '1 A 0.5 0.60 0.60 0.20 5.0e7 1150' + CALL write_deck('stf_ch2.dat', rows(1:1), 2.0e4_wp, .FALSE., channel='AnchDefl1') + CALL run('stf_ch2.dat', 'stf_ch2', ok, msg) + CALL require(.NOT. ok .AND. INDEX(msg, 'END CONNECTIONS row at End B') > 0, & + 'grammar: AnchDefl at a pinned End B: '//TRIM(msg)) + CALL write_deck('stf_ch3.dat', rows(1:1), 2.0e4_wp, .FALSE., channel='FairDefl1') + CALL run('stf_ch3.dat', 'stf_ch3', ok, msg) + CALL require(ok, 'grammar: FairDefl at the clamped End A: '//TRIM(msg)) + END SUBROUTINE case_grammar + + SUBROUTINE case_coupled() + !! The native deck of case_grammar on the coupled route: its root element, wholly inside the + !! 0.5 m cylindrical root (OD 0.6 m, ID 0.2 m), has EI_c + E pi/64 (OD^4 - ID^4), and the + !! initial fairlead load matches the standalone static fairlead tension. + TYPE(CD_AGG_ModuleType), SAVE :: agg + REAL(wp), PARAMETER :: EI_ROOT = 2.0e4_wp + 5.0e7_wp*PI/64.0_wp*(0.6_wp**4 - 0.2_wp**4) + REAL(wp) :: r0(3, 1), vel(3, 1), acc(3, 1), load(3, 1) + REAL(wp), ALLOCATABLE :: prof(:, :) + CHARACTER(512) :: em + INTEGER :: es, ne, n + CALL CD_AGG_Init_From_Deck(agg, 'stf_native.dat', 0.05_wp, es, em) + CALL require(es == CD_AGG_OK, 'coupled: the stiffener deck initializes: '//TRIM(em)) + IF (es /= CD_AGG_OK) RETURN + ! internal element order runs from End B: the root element at End A is the last + ne = SIZE(agg%cables(1)%line%EI) + CALL require(close_rel(agg%cables(1)%line%EI(ne), EI_ROOT, 1.0e-12_wp), & + 'coupled: the root element carries the stiffener') + CALL CD_AGG_GetMovingPointMesh(agg, r0, vel, acc, load, es, em) + CALL read_static('stf_native.static.out', prof, n) + IF (n > 0) CALL require(close_rel(NORM2(load(:, 1)), prof(5, 1), 1.0e-6_wp), & + 'coupled: the initial fairlead load is the standalone static tension') + CALL CD_AGG_End(agg, es, em) + END SUBROUTINE case_coupled + +END PROGRAM test_bend_stiffener diff --git a/validation/scripts/orcaflex_bend_stiffener.py b/validation/scripts/orcaflex_bend_stiffener.py new file mode 100644 index 0000000..99d9fa3 --- /dev/null +++ b/validation/scripts/orcaflex_bend_stiffener.py @@ -0,0 +1,101 @@ +# File: validation/scripts/orcaflex_bend_stiffener.py +# SPDX-License-Identifier: Apache-2.0 +# Copyright (c) 2026 Jae Hoon Seo, SMI Lab, Inha University +"""OrcaFlex static reference for a hang-off clamped through a tapered bend stiffener. + +Provenance tool (NOT built by CMake): regenerates +`tests/data/orcaflex_bend_stiffener_curvature.csv`, the OrcaFlex column of the +`driver_bend_stiffener` check (tests/check_bend_stiffener.py) and of the L3-7 row of +VALIDATION.md. The problem is the check's deck: a 60 m power cable (OD 0.20 m, 60 kg/m, +EA 500 MN, EI 20 kN.m^2) clamped horizontal at (0, 0, -20) m and fixed (pinned) at +(51.9, 0, -50) m, in still water 100 m deep, through a 5 m polyurethane stiffener at End A: a +0.5 m cylinder of 0.60 m outer diameter and a 4.5 m cone to 0.22 m, bored to 0.20 m +(E 50 MPa, 1150 kg/m^3). + +In OrcaFlex the stiffener is a line attachment (Stiffener type) whose line type is a +homogeneous pipe of that material, with its outer diameter given as a line-type-diameter +profile; the attachment and the cable bend together. The cable is meshed at 0.05 m over its +first 8 m and 1 m beyond, as in the check. The CSV holds the static curvature of the cable +along its first 8 m (arc from End A, curvature in 1/m). + +Result (OrcFxAPI 11.6d): root curvature 0.4522 /m (CableDyn 0.4523), cable bend moment at the +root 9.04 kN.m; CableDyn's curvature follows the profile within 0.8 % of the root curvature. +OrcaFlex SI (m, te, kN). + +Usage: python validation/scripts/orcaflex_bend_stiffener.py [output.csv] +""" + +import pathlib +import sys + +import OrcFxAPI as ofx + +ROOT = pathlib.Path(__file__).resolve().parents[2] +OUTPUT = ROOT / "tests" / "data" / "orcaflex_bend_stiffener_curvature.csv" +WINDOW = 8.0 + + +def build() -> tuple[ofx.Model, ofx.OrcaFlexObject]: + m = ofx.Model() + m.environment.WaterDepth = 100.0 + m.environment.Density = 1.025 + profile = m.CreateObject(ofx.otLineTypeDiameter, "StiffenerOD") + profile.SetDataRowCount("IndependentValue", 3) + for i, (arc, od) in enumerate([(0.0, 0.60), (0.5, 0.60), (5.0, 0.22)]): + profile.SetData("IndependentValue", i, arc) + profile.SetData("DependentValue", i, od) + cable = m.CreateObject(ofx.otLineType, "cable") + cable.OD, cable.ID = 0.20, 0.0 + cable.MassPerUnitLength = 0.060 + cable.EIx = 20.0 + cable.EA = 5.0e5 + sleeve = m.CreateObject(ofx.otLineType, "stiffener material") + sleeve.Category = "Homogeneous pipe" + sleeve.OD = "StiffenerOD" + sleeve.ID = 0.20 + sleeve.E = 5.0e4 # kPa + sleeve.MaterialDensity = 1.15 + stiffener = m.CreateObject(ofx.otStiffenerType, "stiffener") + stiffener.LineType = "stiffener material" # its length (5 m) follows the diameter profile + line = m.CreateObject(ofx.otLine, "hangoff") + line.EndAConnection, line.EndBConnection = "Fixed", "Fixed" + line.EndAX, line.EndAY, line.EndAZ = 0.0, 0.0, -20.0 + line.EndBX, line.EndBY, line.EndBZ = 51.90, 0.0, -50.0 + line.EndAAzimuth, line.EndADeclination = 0.0, 90.0 + line.EndAxBendingStiffness = ofx.OrcinaInfinity() + line.EndBxBendingStiffness = 0.0 + line.NumberOfSections = 2 + line.LineType = ("cable", "cable") + line.Length = (8.0, 52.0) + line.TargetSegmentLength = (0.05, 1.0) + line.NumberOfAttachments = 1 + line.AttachmentType = ("stiffener",) + line.Attachmentz = (0.0,) + line.AttachmentzRelativeTo = ("End A",) + return m, line + + +def main() -> int: + output = pathlib.Path(sys.argv[1]) if len(sys.argv) > 1 else OUTPUT + m, line = build() + m.CalculateStatics() + curvature = line.RangeGraph("Curvature", ofx.pnStaticState) + moment = line.RangeGraph("Bend moment", ofx.pnStaticState) + rows = [ + (arc, k) + for arc, k in zip(curvature.X, curvature.Mean) + if arc <= WINDOW + 1.0e-9 + ] + with output.open("w", encoding="utf-8", newline="\n") as handle: + handle.write("# OrcFxAPI 11.6d static curvature along the hang-off (validation/scripts/" + "orcaflex_bend_stiffener.py)\n") + handle.write("arc_m,curvature_per_m\n") + for arc, k in rows: + handle.write(f"{arc:.6f},{k:.8e}\n") + print(f"root curvature {curvature.Mean[0]:.5f} /m, cable bend moment {moment.Mean[0]:.3f} kN.m") + print(f"wrote {len(rows)} rows to {output}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) From ddbfbbc6ac690f8a31110cab357a7e7d63531f2c Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 26 Oct 2026 15:19:59 +0900 Subject: [PATCH 2/3] fix(stiffeners): keep EI = 0 on a partly sleeved EI = 0 element A stiffener ending partway through an element of an EI = 0 section gave the element the sleeve's bending stiffness over its whole length, since the bare part's compliance was skipped. That part has an infinite compliance, so the compliance average stays zero: only an element the sleeve covers completely now gains bending stiffness. Signed-off-by: Jae Hoon Seo --- doc/theory.rst | 2 ++ src/CableDyn_DeckDriver.f90 | 10 ++++++++-- tests/test_bend_stiffener.f90 | 6 ++++++ 3 files changed, 16 insertions(+), 2 deletions(-) diff --git a/doc/theory.rst b/doc/theory.rst index 5ae0426..98866a7 100644 --- a/doc/theory.rst +++ b/doc/theory.rst @@ -165,6 +165,8 @@ length :math:`L` takes the compliance average which gives the exact relative end rotation of the element under a uniform moment: the stiff root of a tapered sleeve is not allowed to stiffen a whole element that only partly lies on it. +On an ``EI = 0`` section of the line the bare part of an element has an infinite compliance, +so only an element the sleeve covers completely gains bending stiffness. The element mass per length is the cable's plus :math:`\rho_s \tfrac{\pi}{4}\, \overline{D^2 - d_i^2}`; buoyancy acts on the displaced area of cable and solid sleeve, :math:`\tfrac{\pi}{4}\,(d^2 + \overline{D^2 - d_i^2})`, and drag and added mass on the mean outer diff --git a/src/CableDyn_DeckDriver.f90 b/src/CableDyn_DeckDriver.f90 index ac0f6ac..e43c5c5 100644 --- a/src/CableDyn_DeckDriver.f90 +++ b/src/CableDyn_DeckDriver.f90 @@ -9299,8 +9299,14 @@ SUBROUTINE CD_Stiffener_Overlay(stf, l0, diam, ei, mpl, volume_diameter) END DO IF (.NOT. (covered > CD_ZERO)) CYCLE IF (PRESENT(ei)) THEN - IF (ei_line > CD_ZERO) compliance = compliance + (l0(e) - covered)/ei_line - ei(e) = l0(e)/compliance + IF (ei_line > CD_ZERO) THEN + compliance = compliance + (l0(e) - covered)/ei_line + ei(e) = l0(e)/compliance + ELSE IF (covered >= l0(e)*(CD_ONE - 1.0e-12_wp)) THEN + ! an EI = 0 element bends only where the sleeve covers all of it; a bare part has an + ! infinite compliance and leaves the element without bending stiffness + ei(e) = l0(e)/compliance + END IF END IF IF (PRESENT(mpl)) mpl(e) = mpl(e) + sleeve_mass/l0(e) IF (volume_diameter) THEN diff --git a/tests/test_bend_stiffener.f90 b/tests/test_bend_stiffener.f90 index ae03150..15fabbe 100644 --- a/tests/test_bend_stiffener.f90 +++ b/tests/test_bend_stiffener.f90 @@ -88,6 +88,12 @@ SUBROUTINE case_cylinder() CALL require(close_rel(dh(e_half), 0.5_wp*(OD + D_C), 1.0e-13_wp), 'cylinder: half-covered outer diameter') CALL require(ALL(ei(e_far) == EI_C) .AND. ALL(mpl(e_far) == M_C) .AND. ALL(dv(e_far) == D_C) .AND. & ALL(dh(e_far) == D_C), 'cylinder: elements the sleeve does not reach are unchanged') + ! an EI = 0 section: the covered element bends with the sleeve alone; the half-covered one + ! keeps EI = 0, its bare half having an infinite compliance + ei = 0.0_wp; mpl = M_C; dv = D_C + CALL CD_Stiffener_Overlay(stf, l0, dv, ei=ei, mpl=mpl, volume_diameter=.TRUE.) + CALL require(close_rel(ei(e_full), ei_s, 1.0e-13_wp), 'cylinder: covered EI = 0 element takes the sleeve') + CALL require(ei(e_half) == 0.0_wp, 'cylinder: half-covered EI = 0 element keeps EI = 0') END DO END SUBROUTINE case_cylinder From bd571673a84c5cfdb54e13ce457d404cd9fe5ee2 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 26 Oct 2026 15:31:02 +0900 Subject: [PATCH 3/3] fix(stiffeners): give a sleeved element one diameter for all its loads The drag diameter of a sleeved element was the mean outer diameter while its weight used the displaced-area diameter, but the dynamic model also recovers the dry-part buoyancy above the surface from the drag diameter. A sleeved element at or above the waterline then lost a different buoyancy than the static solve assumed. Like a line type, an element now has one diameter, that of the displaced area of cable and solid sleeve, for its buoyancy, added mass, Froude-Krylov force and drag; with the bore equal to the cable diameter it is the sleeve's outer diameter. Signed-off-by: Jae Hoon Seo --- doc/driver_format.md | 5 +++-- doc/theory.rst | 9 +++++---- src/CableDyn_DeckDriver.f90 | 25 ++++++++++--------------- tests/test_bend_stiffener.f90 | 34 +++++++++++++++------------------- 4 files changed, 33 insertions(+), 40 deletions(-) diff --git a/doc/driver_format.md b/doc/driver_format.md index d039530..14bf2e2 100644 --- a/doc/driver_format.md +++ b/doc/driver_format.md @@ -969,8 +969,9 @@ Behaviour: `EI(s) = EI_cable + E π/64 (OD(s)⁴ − InnerD⁴)`. The sleeve carries no axial load and no torque, so `EA` and `GJ` are the cable's. Each element takes the compliance average of `EI(s)` over its length, `L / ∫ ds/EI(s)`, the cable mass plus the mean sleeve mass, a - displaced area of cable plus solid sleeve, and the mean outer diameter for drag and added mass - (see [theory](theory.rst), *Bend stiffeners*). Elements the stiffener does not reach keep + diameter of the displaced area of cable plus solid sleeve, `√(d² + OD² − InnerD²)`, for its + buoyancy, added mass and drag alike, as a line type has one diameter (the sleeve's outer + diameter when `InnerD` is the cable diameter; see [theory](theory.rst), *Bend stiffeners*). Elements the stiffener does not reach keep their line type's values exactly. - The element properties are piecewise constant, one value per deck element, so the deck mesh must resolve the taper: an initialisation note names a stiffener that spans fewer than 10 diff --git a/doc/theory.rst b/doc/theory.rst index 98866a7..046aea2 100644 --- a/doc/theory.rst +++ b/doc/theory.rst @@ -168,10 +168,11 @@ root of a tapered sleeve is not allowed to stiffen a whole element that only par On an ``EI = 0`` section of the line the bare part of an element has an infinite compliance, so only an element the sleeve covers completely gains bending stiffness. The element mass per length is the cable's plus :math:`\rho_s \tfrac{\pi}{4}\, -\overline{D^2 - d_i^2}`; buoyancy acts on the displaced area of cable and solid sleeve, -:math:`\tfrac{\pi}{4}\,(d^2 + \overline{D^2 - d_i^2})`, and drag and added mass on the mean outer -diameter :math:`\overline{\max(D, d)}`, where the bar is the mean over the element and -:math:`d` the cable diameter. The integrals use six-point Gauss–Legendre quadrature on each +\overline{D^2 - d_i^2}`, where the bar is the mean over the element. Like a line type, the +element has one diameter, that of the displaced area of cable and solid sleeve, +:math:`\sqrt{d^2 + \overline{D^2 - d_i^2}}` with :math:`d` the cable diameter, for its +buoyancy (the dry part above the surface included), added mass, Froude–Krylov force and drag; +with the bore equal to the cable diameter it is the sleeve's outer diameter. The integrals use six-point Gauss–Legendre quadrature on each part of the element a segment covers. With these one-per-element values the taper is resolved by the mesh, as in a segmented line model; the driver notes a stiffener that spans fewer than ten elements. diff --git a/src/CableDyn_DeckDriver.f90 b/src/CableDyn_DeckDriver.f90 index e43c5c5..5e2b639 100644 --- a/src/CableDyn_DeckDriver.f90 +++ b/src/CableDyn_DeckDriver.f90 @@ -9148,7 +9148,7 @@ SUBROUTINE build_finite_section_mesh(ln, sections, types, conn, l0, ea, mpl, dia i_rn(e) = types(tidx)%i_rn END DO END DO - IF (ALLOCATED(ln%stf)) CALL CD_Stiffener_Overlay(ln%stf, l0, diam, ei=ei, mpl=mpl, volume_diameter=.TRUE.) + IF (ALLOCATED(ln%stf)) CALL CD_Stiffener_Overlay(ln%stf, l0, diam, ei=ei, mpl=mpl) l0 = l0(n_elem:1:-1) ea = ea(n_elem:1:-1) mpl = mpl(n_elem:1:-1) @@ -9225,8 +9225,8 @@ SUBROUTINE finite_hydro_arrays(ln, sections, types, diam, cdn, cdt, can, cat) cat(e) = types(tidx)%cat END DO END DO - ! a stiffener widens the line to its outer diameter for drag and added mass - IF (ALLOCATED(ln%stf)) CALL CD_Stiffener_Overlay(ln%stf, l0, diam, volume_diameter=.FALSE.) + ! a stiffener widens the line: the same displaced-area diameter as its weight and buoyancy + IF (ALLOCATED(ln%stf)) CALL CD_Stiffener_Overlay(ln%stf, l0, diam) diam = diam(SIZE(diam):1:-1) cdn = cdn(SIZE(cdn):1:-1) cdt = cdt(SIZE(cdt):1:-1) @@ -9234,7 +9234,7 @@ SUBROUTINE finite_hydro_arrays(ln, sections, types, diam, cdn, cdt, can, cat) cat = cat(SIZE(cat):1:-1) END SUBROUTINE finite_hydro_arrays - SUBROUTINE CD_Stiffener_Overlay(stf, l0, diam, ei, mpl, volume_diameter) + SUBROUTINE CD_Stiffener_Overlay(stf, l0, diam, ei, mpl) !! Add the bend stiffeners of a finite-EI line (ln%stf columns) to its element properties, !! in public element order (End A first). A stiffener is a polyurethane sleeve bonded to the !! cable in bending: over its length the section has the bending stiffness of the cable plus @@ -9245,17 +9245,17 @@ SUBROUTINE CD_Stiffener_Overlay(stf, l0, diam, ei, mpl, volume_diameter) !! uniform moment, so a stiff root next to a soft tip is not over-stiffened; !! mpl the cable mass plus the mean sleeve mass, rho pi/4 (OD^2 - ID^2); !! diam the diameter of the mean displaced area, cable plus solid sleeve, sqrt(d^2 + mean - !! (OD^2 - ID^2)) (volume_diameter: weight and buoyancy), or the mean outer diameter, - !! mean max(OD, d) (otherwise: drag and added mass). + !! (OD^2 - ID^2)): one diameter, like a line type's, for the weight and buoyancy (dry + !! part above the surface included), the added mass, Froude-Krylov force and drag. With + !! the bore equal to the cable diameter it is the sleeve's outer diameter. !! The integrals use 6-point Gauss-Legendre on each part of the element a segment covers. !! An element no stiffener reaches is left untouched. REAL(wp), INTENT(IN) :: stf(:, :), l0(:) REAL(wp), INTENT(INOUT) :: diam(:) REAL(wp), INTENT(INOUT), OPTIONAL :: ei(:), mpl(:) - LOGICAL, INTENT(IN) :: volume_diameter REAL(wp), PARAMETER :: PI_64 = 0.0490873852123405193_wp ! pi/64 REAL(wp) :: gp(6), gw(6), s_end(0:SIZE(l0)), line_len, a, b, lo, hi, x, t, od, d_line - REAL(wp) :: compliance, sleeve_mass, sleeve_area, outer, covered, ei_line + REAL(wp) :: compliance, sleeve_mass, sleeve_area, covered, ei_line INTEGER :: e, k, g CALL CD_HermiteCable_Gauss_Rule(6, gp, gw) @@ -9270,7 +9270,7 @@ SUBROUTINE CD_Stiffener_Overlay(stf, l0, diam, ei, mpl, volume_diameter) d_line = diam(e) ei_line = CD_ZERO IF (PRESENT(ei)) ei_line = ei(e) - compliance = CD_ZERO; sleeve_mass = CD_ZERO; sleeve_area = CD_ZERO; outer = CD_ZERO; covered = CD_ZERO + compliance = CD_ZERO; sleeve_mass = CD_ZERO; sleeve_area = CD_ZERO; covered = CD_ZERO DO k = 1, SIZE(stf, 2) ! the segment's span on the public arc IF (NINT(stf(1, k)) == 1) THEN @@ -9294,7 +9294,6 @@ SUBROUTINE CD_Stiffener_Overlay(stf, l0, diam, ei, mpl, volume_diameter) compliance = compliance + gw(g)*(hi - lo)/(ei_line + stf(7, k)*PI_64*(od**4 - stf(6, k)**4)) sleeve_mass = sleeve_mass + gw(g)*(hi - lo)*stf(8, k)*QUARTER_PI*(od**2 - stf(6, k)**2) sleeve_area = sleeve_area + gw(g)*(hi - lo)*(od**2 - stf(6, k)**2) - outer = outer + gw(g)*(hi - lo)*MAX(od, d_line) END DO END DO IF (.NOT. (covered > CD_ZERO)) CYCLE @@ -9309,11 +9308,7 @@ SUBROUTINE CD_Stiffener_Overlay(stf, l0, diam, ei, mpl, volume_diameter) END IF END IF IF (PRESENT(mpl)) mpl(e) = mpl(e) + sleeve_mass/l0(e) - IF (volume_diameter) THEN - diam(e) = SQRT(d_line**2 + sleeve_area/l0(e)) - ELSE - diam(e) = (outer + (l0(e) - covered)*d_line)/l0(e) - END IF + diam(e) = SQRT(d_line**2 + sleeve_area/l0(e)) END DO END SUBROUTINE CD_Stiffener_Overlay diff --git a/tests/test_bend_stiffener.f90 b/tests/test_bend_stiffener.f90 index 15fabbe..f4fcc37 100644 --- a/tests/test_bend_stiffener.f90 +++ b/tests/test_bend_stiffener.f90 @@ -5,7 +5,8 @@ PROGRAM test_bend_stiffener !! Gates for bend stiffeners (STIFFENERS) on a finite-EI line: !! 1. Element properties (CD_Stiffener_Overlay) against closed forms: a cylindrical sleeve !! over half an element gives the compliance average of the bending stiffness, the mean - !! sleeve mass, the displaced-area diameter and the mean outer diameter; the same sleeve + !! sleeve mass and the displaced-area diameter (the outer one when the bore is the + !! cable's); the same sleeve !! at End B lands on the mirrored element; elements it does not reach are unchanged bit !! for bit. !! 2. A cone: the element values follow a 20 000-interval midpoint integration of the same @@ -59,7 +60,7 @@ SUBROUTINE case_cylinder() !! half of element 2 from End A (or from End B: elements 4 and 3). REAL(wp), PARAMETER :: EI_C = 2.0e4_wp, M_C = 60.0_wp, D_C = 0.2_wp REAL(wp), PARAMETER :: OD = 0.5_wp, ID = 0.2_wp, EMOD = 5.0e7_wp, RHO = 1150.0_wp - REAL(wp) :: stf(8, 1), l0(4), ei(4), mpl(4), dv(4), dh(4), ei_s, ei_half, m_s + REAL(wp) :: stf(8, 1), l0(4), ei(4), mpl(4), dv(4), ei_s, ei_half, m_s INTEGER :: iend, e_full, e_half, e_far(2) l0 = 1.0_wp ei_s = EMOD*PI/64.0_wp*(OD**4 - ID**4) @@ -67,9 +68,8 @@ SUBROUTINE case_cylinder() ei_half = 1.0_wp/(0.5_wp/EI_C + 0.5_wp/(EI_C + ei_s)) DO iend = 1, 2 stf(:, 1) = [REAL(iend, wp), 0.0_wp, 1.5_wp, OD, OD, ID, EMOD, RHO] - ei = EI_C; mpl = M_C; dv = D_C; dh = D_C - CALL CD_Stiffener_Overlay(stf, l0, dv, ei=ei, mpl=mpl, volume_diameter=.TRUE.) - CALL CD_Stiffener_Overlay(stf, l0, dh, volume_diameter=.FALSE.) + ei = EI_C; mpl = M_C; dv = D_C + CALL CD_Stiffener_Overlay(stf, l0, dv, ei=ei, mpl=mpl) IF (iend == 1) THEN e_full = 1; e_half = 2; e_far = [3, 4] ELSE @@ -84,14 +84,14 @@ SUBROUTINE case_cylinder() 'cylinder: displaced-area diameter of cable plus solid sleeve') CALL require(close_rel(dv(e_half), SQRT(D_C**2 + 0.5_wp*(OD**2 - ID**2)), 1.0e-13_wp), & 'cylinder: half-covered displaced-area diameter') - CALL require(close_rel(dh(e_full), OD, 1.0e-13_wp), 'cylinder: outer diameter for drag') - CALL require(close_rel(dh(e_half), 0.5_wp*(OD + D_C), 1.0e-13_wp), 'cylinder: half-covered outer diameter') - CALL require(ALL(ei(e_far) == EI_C) .AND. ALL(mpl(e_far) == M_C) .AND. ALL(dv(e_far) == D_C) .AND. & - ALL(dh(e_far) == D_C), 'cylinder: elements the sleeve does not reach are unchanged') + CALL require(close_rel(dv(e_full), OD, 1.0e-13_wp), & + 'cylinder: with the bore equal to the cable diameter the diameter is the outer one') + CALL require(ALL(ei(e_far) == EI_C) .AND. ALL(mpl(e_far) == M_C) .AND. ALL(dv(e_far) == D_C), & + 'cylinder: elements the sleeve does not reach are unchanged') ! an EI = 0 section: the covered element bends with the sleeve alone; the half-covered one ! keeps EI = 0, its bare half having an infinite compliance ei = 0.0_wp; mpl = M_C; dv = D_C - CALL CD_Stiffener_Overlay(stf, l0, dv, ei=ei, mpl=mpl, volume_diameter=.TRUE.) + CALL CD_Stiffener_Overlay(stf, l0, dv, ei=ei, mpl=mpl) CALL require(close_rel(ei(e_full), ei_s, 1.0e-13_wp), 'cylinder: covered EI = 0 element takes the sleeve') CALL require(ei(e_half) == 0.0_wp, 'cylinder: half-covered EI = 0 element keeps EI = 0') END DO @@ -103,17 +103,16 @@ SUBROUTINE case_cone() REAL(wp), PARAMETER :: EI_C = 2.0e4_wp, M_C = 60.0_wp, D_C = 0.2_wp REAL(wp), PARAMETER :: OD0 = 0.6_wp, OD1 = 0.22_wp, ID = 0.2_wp, EMOD = 5.0e7_wp, RHO = 1150.0_wp INTEGER, PARAMETER :: NQ = 20000 - REAL(wp) :: stf(8, 1), l0(5), ei(5), mpl(5), dv(5), dh(5), a, s, od, comp, mass, area, outer, h, tol + REAL(wp) :: stf(8, 1), l0(5), ei(5), mpl(5), dv(5), a, s, od, comp, mass, area, h, tol INTEGER :: e, k l0 = 0.7_wp stf(:, 1) = [1.0_wp, 0.0_wp, 2.5_wp, OD0, OD1, ID, EMOD, RHO] - ei = EI_C; mpl = M_C; dv = D_C; dh = D_C - CALL CD_Stiffener_Overlay(stf, l0, dv, ei=ei, mpl=mpl, volume_diameter=.TRUE.) - CALL CD_Stiffener_Overlay(stf, l0, dh, volume_diameter=.FALSE.) + ei = EI_C; mpl = M_C; dv = D_C + CALL CD_Stiffener_Overlay(stf, l0, dv, ei=ei, mpl=mpl) DO e = 1, 4 a = 0.7_wp*REAL(e - 1, wp) h = 0.7_wp/REAL(NQ, wp) - comp = 0.0_wp; mass = 0.0_wp; area = 0.0_wp; outer = 0.0_wp + comp = 0.0_wp; mass = 0.0_wp; area = 0.0_wp DO k = 1, NQ s = a + (REAL(k, wp) - 0.5_wp)*h IF (s < 2.5_wp) THEN @@ -121,10 +120,8 @@ SUBROUTINE case_cone() comp = comp + h/(EI_C + EMOD*PI/64.0_wp*(od**4 - ID**4)) mass = mass + h*RHO*PI/4.0_wp*(od**2 - ID**2) area = area + h*(od**2 - ID**2) - outer = outer + h*MAX(od, D_C) ELSE comp = comp + h/EI_C - outer = outer + h*D_C END IF END DO ! the cone tip inside element 4 limits the midpoint rule there to O(h/l) = 5e-5 @@ -132,11 +129,10 @@ SUBROUTINE case_cone() CALL require(close_rel(ei(e), 0.7_wp/comp, tol), 'cone: compliance-average EI of element') CALL require(close_rel(mpl(e), M_C + mass/0.7_wp, tol), 'cone: mean mass of element') CALL require(close_rel(dv(e), SQRT(D_C**2 + area/0.7_wp), tol), 'cone: displaced-area diameter') - CALL require(close_rel(dh(e), outer/0.7_wp, tol), 'cone: mean outer diameter') IF (e < 4) CALL require(close_rel(ei(e), 0.7_wp/comp, 1.0e-10_wp), & 'cone: Gauss average matches the fine integration inside the cone') END DO - CALL require(ei(5) == EI_C .AND. mpl(5) == M_C .AND. dv(5) == D_C .AND. dh(5) == D_C, & + CALL require(ei(5) == EI_C .AND. mpl(5) == M_C .AND. dv(5) == D_C, & 'cone: bare elements unchanged') END SUBROUTINE case_cone