From 71b9dd4bf2382ca48c0833e8e5a8d325e8baff98 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sat, 7 Nov 2026 22:38:54 +0900 Subject: [PATCH] fix(output): report the geometric curvature and the section moment of a stretched line The element energy is EI kappa^2/2 per unstretched length with kappa the curvature of the stretched line, so the recovered m, the conjugate of the tangent r', has the normal length EI kappa/(1 + e)^2. The node outputs read that length as the moment and divided it by EI, which put Curv 2e and BendMom e below the values they are defined by. The energy does not change under a rotation of the whole line, so the moment that balances the forces is r' x m, of length M = EI kappa/(1 + e). BendMom is now M, Curv is (1 + e) M/EI, the reciprocal of the bend radius, and StfMom the sleeve's share of M, with the stretch |r'| of the node's tangent. Tension, positions and forces do not change. Under OutputRecovery element the exact traces and EI x curvature stand on unsleeved elements. driver_curvature_stretch holds both against the elastica of this energy, strain swept to 1.2 %: the curvature error falls from 2.30 % to 0.006 % and the moment error from 1.15 % to under 0.001 %. driver_bend_stiffener takes the same elastica: the cable moment moves from 0.18 % to 0.01 % of the root value and the node positions from 0.6 to 0.03 mm. Signed-off-by: Jae Hoon Seo --- CMakeLists.txt | 9 + VALIDATION.md | 22 +-- doc/driver_format.md | 4 +- doc/options.rst | 3 +- doc/outputs.rst | 20 +- doc/theory.rst | 68 ++++--- python/tests/test_objects.py | 5 +- src/CableDyn_DeckDriver.f90 | 12 +- src/CableDyn_HermiteCableDynamic.f90 | 49 +++-- tests/check_bend_stiffener.py | 61 +++--- tests/check_curvature_stretch.py | 175 ++++++++++++++++++ tests/test_driver_channels.f90 | 12 +- tests/test_openfast_aggregate.f90 | 14 +- validation/scripts/orcaflex_bend_stiffener.py | 4 +- 14 files changed, 358 insertions(+), 100 deletions(-) create mode 100644 tests/check_curvature_stretch.py diff --git a/CMakeLists.txt b/CMakeLists.txt index be0cd833..fa0f0ea2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2763,6 +2763,15 @@ if(Python3_Interpreter_FOUND) LABELS "python;driver;integration;static;validation" SKIP_RETURN_CODE 77 ENVIRONMENT "${_cabledyn_python_test_environment}") + # The curvature of a stretched cable is the geometric one and the bend moment the section's: + # both follow the elastica of the element energy at strains up to 1.2 %. + add_test(NAME driver_curvature_stretch + COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/tests/check_curvature_stretch.py + $) + set_tests_properties(driver_curvature_stretch 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. diff --git a/VALIDATION.md b/VALIDATION.md index bf92b468..cd8f8fa6 100644 --- a/VALIDATION.md +++ b/VALIDATION.md @@ -78,7 +78,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.51% of root; cable moment 0.18%, section moment 0.18%, deflection angle 0.003°; local curvature 0.24% (0.1 m elements) | <1.5% / <2%, 0.05°, <0.5% | Pass | +| **L3-8** | Bend stiffener at a clamped hang-off (5 m tapered polyurethane) | OrcaFlex 11.6d (stiffener attachment) / planar elastica | curvature 0.45% of root; cable moment 0.01%, section moment 0.01%, deflection angle 0.0001°; local curvature 0.27% (0.1 m elements) | <1.5% / <2%, 0.05°, <0.5% | 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.11%; 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 | @@ -275,19 +275,19 @@ initialisation. 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.51% of the root curvature** (root 0.4523 against 0.4522 1/m; + within **0.45% of the root curvature** (root 0.4528 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.18%** 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°). Node by node, relative to the local value floored at 0.1% of the root - curvature, the curvature follows the elastica within **0.97 / 0.49 / 0.24 / 0.18%** on + Against an independent collocation solution of the planar extensible elastica of the element + energy with `EI(s)` and the sleeve's weight, the cable's own moment follows within **0.01%** of + its root value, cable and sleeve moments add up to the section moment within **0.01%**, the node + positions within **0.03 mm**, and the end force's deflection from the clamp axis within + **0.0001°** (32.37°). Node by node, relative to the local value floored at 0.1% of the root + curvature, the curvature follows the elastica within **0.99 / 0.52 / 0.27 / 0.06%** on 0.2 / 0.15 / 0.1 / 0.05 m elements (gates 1 / 1 / 0.5 / 0.5%), under either output recovery, and under the stiffener the element peaks of `.elements.out`, read through `bend_check`, - within **the same figures** on the same gates; 0.12% of it at the root is - the order-`T/EA` difference between the elastica's bending law and the element energy's, - not mesh error. + within **the same figures** on the same gates. On a cable strained up to 1.2% + (`driver_curvature_stretch`) the curvature and the section moment follow the same elastica + within **0.006%** and **0.001%** of the local value. - **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/driver_format.md b/doc/driver_format.md index 5a7f5447..5d584ca5 100644 --- a/doc/driver_format.md +++ b/doc/driver_format.md @@ -1553,8 +1553,8 @@ bending moment. None of these channels is mandatory. | `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) | 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. Under a bend stiffener: the cable's own share, the line type's EI × curvature | N·m | +| `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: the geometric curvature of the stretched line, (1 + ε) `BendMomN` / EI with ε the axial strain at the node, 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 / (1 + ε), the moment that balances the forces on the section, **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 section moment × the line type's EI / EI(s) | 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 | diff --git a/doc/options.rst b/doc/options.rst index abb903df..c2aa51d2 100644 --- a/doc/options.rst +++ b/doc/options.rst @@ -492,7 +492,8 @@ Finite-EI cable controls static-profile and range-graph columns, ``.t.out`` segment tensions). ``equilibrium`` recovers them from nodal equilibrium: the effective tension :math:`\mathbf F\cdot\mathbf t`, the moment the neighbouring elements carry, and - curvature = moment / ``EI`` (see :ref:`theory-output-recovery`). ``element`` reports the + the geometric curvature (1 + strain) × moment / ``EI`` (see + :ref:`theory-output-recovery`). ``element`` reports the length-weighted element-mean axial force :math:`EA(\lvert\mathbf r'\rvert - 1)`, the larger one-sided exact curvature (on an element with a bend stiffener sleeve, the recovered moment over ``EI(s)``), and the axial resultant at the element midpoint as the diff --git a/doc/outputs.rst b/doc/outputs.rst index 3ef03b1e..7770d0d8 100644 --- a/doc/outputs.rst +++ b/doc/outputs.rst @@ -695,9 +695,11 @@ Line-node channels * - ``CurvN`` - curvature at node ``J``. ``EI = 0`` lines: discrete 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 (under a bend stiffener the section moment over the stiffness ``EI(s)`` - at the node, :ref:`theory-bend-stiffeners`), and the exact geometric curvature where + Cubic-Hermite lines: the geometric curvature of the stretched line, the reciprocal of the + bend radius, ``(1 + ε) BendMomN`` / ``EI`` with ``ε`` the axial strain at the node, + the larger of the two sides at a section junction (under a bend stiffener ``(1 + ε)`` + times the section moment over the stiffness ``EI(s)`` at the node, + :ref:`theory-bend-stiffeners`), and the exact geometric curvature where ``EI`` is below ``max(|w| L0^2)/12`` times 1 m (:ref:`theory-output-recovery`) (``OutputRecovery element``: exact curvature of the continuous centreline, the larger of the two one-sided element values at an interior node; the @@ -717,15 +719,15 @@ Line-node channels returns them for every node - 1/m * - ``BendMomN`` - - bend moment ``EI × κ`` measured relative to the stress-free reference shape (the - cubic-Hermite reference is straight). Finite-EI lines: the component normal to the tangent - of the moment recovered from nodal equilibrium, the mean of the two neighbouring elements + - bend moment ``EI × κ / (1 + ε)``, the moment that balances the forces on the section (it + is ``EI × κ`` to within the axial strain ``ε``), measured relative to the stress-free + reference shape (the cubic-Hermite reference is straight). Finite-EI lines: 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. Under a bend stiffener (``STIFFENERS``) it is the cable's own share, - the line type's ``EI`` times the curvature, which is the section moment over the - stiffness ``EI(s)`` at the node; the stiffener's share is ``StfMomN`` and the two - add up to the section moment + the section moment times the line type's ``EI`` over the stiffness ``EI(s)`` at the node; + 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 diff --git a/doc/theory.rst b/doc/theory.rst index faf73946..9d5af7f5 100644 --- a/doc/theory.rst +++ b/doc/theory.rst @@ -183,10 +183,10 @@ steps at every node although the section does not, and the softer side would ove curvature by about :math:`h\lvert EI'\rvert/(2 EI)`, first order in the element length :math:`h`. A node takes :math:`EI(s)` at its own arc instead, from the stiffener rows, one limit from either side. The section moment :math:`M` (:ref:`theory-output-recovery`) is continuous and -second-order accurate, so the curvature :math:`\kappa = M/EI(s)` is too. Cable and sleeve take -that curvature and share the moment in proportion to the stiffnesses: the cable carries -:math:`EI_c\,\kappa`, reported as ``BendMom``, and the sleeve -:math:`(EI(s) - EI_c)\,\kappa`, reported as ``StfMom``. +second-order accurate, so the curvature :math:`\kappa = (1+\varepsilon) M/EI(s)` is too. Cable +and sleeve take that curvature and share the moment in proportion to the stiffnesses: the cable +carries :math:`EI_c\,M/EI(s)`, reported as ``BendMom``, and the sleeve the rest, +:math:`(EI(s) - EI_c)\,M/EI(s)`, reported as ``StfMom``. The two limits differ only at a real step of the section: at the bare tip, between stiffener rows whose end diameter, bore or modulus do not match, and between line types of different @@ -199,8 +199,8 @@ element, so at an element's end it is off by the same first-order step against : high where the curvature falls toward the tip and low where it rises, and :math:`EI_e` times it is low by :math:`h^2 M''/12` against :math:`M`. Two outputs that read the interpolation elsewhere read the section here. The peak of such an element in ``.elements.out`` is the largest -:math:`\lvert M(u)\rvert/EI(s(u))` along it, with :math:`M(u)` the quadratic through the -recovered moments of the end nodes that keeps the mean of the fit; a row end inside the element, +:math:`(1+\varepsilon)\lvert M(u)\rvert/EI(s(u))` along it, with :math:`M(u)` the quadratic +through the recovered moments of the end nodes that keeps the mean of the fit; a row end inside the element, or the tip, is a step of :math:`EI(s)` that the search keeps. Below the stiffness that resolves the recovered moment (:ref:`theory-output-recovery`) both the nodes and the peak fall back on the interpolation's curvature times :math:`EI_e/EI(s)`. ``OutputRecovery element`` reports @@ -217,19 +217,17 @@ 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:`\theta' = (1+\varepsilon)^2 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.2 % 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. Node by node, the curvature and the cable's moment follow the -elastica within 0.24 % of their local value on 0.1 m elements and 0.18 % on 0.05 m, 0.97 % on +3 m beyond follows the elastica within 0.01 % of its root value, cable and sleeve moments add up +to the section moment within 0.01 % at the root, the deflection angle agrees within 0.0001°, and +the node positions within 0.03 mm. Node by node, the curvature and the cable's moment follow the +elastica within 0.27 % of their local value on 0.1 m elements and 0.06 % on 0.05 m, 0.99 % on 0.2 m, under either output recovery, and the peaks of ``.elements.out`` under the stiffener -as closely. The local value is floored at 0.1 % of the root curvature, as the curvature changes sign -near the tip, and part of the difference is not mesh error: the elastica takes the moment -proportional to the rotation per unstretched length, the element energy the curvature of the -stretched line, a difference of order :math:`T/EA` that is 0.12 % at this root. On the same -hang-off, the curvature along the stiffener follows OrcaFlex 11.6d with its stiffener given as a -line attachment within 0.6 % of the root curvature. +as closely. The local value is floored at 0.1 % of the root curvature, as the curvature changes +sign near the tip. On the same hang-off, the curvature along the stiffener follows OrcaFlex +11.6d with its stiffener given as a line attachment within 0.5 % of the root curvature. .. _theory-output-recovery: @@ -247,11 +245,25 @@ and the conjugate of the material tangent are \mathbf{m} = \frac{\partial W}{\partial \mathbf{b}} = EI\,\frac{(\mathbf{a}\times\mathbf{b})\times\mathbf{a}}{\lvert\mathbf{a}\rvert^{6}} . +The energy is unchanged by a rotation of the whole line, so +:math:`\mathbf{r}\times\mathbf{F} + \mathbf{a}\times\mathbf{m}` changes along the line only by +the moment of the loads: the section moment is :math:`\mathbf{a}\times\mathbf{m}`, not +:math:`\mathbf{m}`, which is the conjugate of the tangent. With :math:`\kappa` the curvature of +the stretched line, the energy :math:`\tfrac12 EI\kappa^2` per unstretched metre gives + +.. math:: + + \lvert\mathbf{m}_\perp\rvert = \frac{EI\,\kappa}{(1+\varepsilon)^2}, + \qquad + M = \lvert\mathbf{a}\times\mathbf{m}\rvert = \frac{EI\,\kappa}{1+\varepsilon}, + \qquad + \kappa = (1+\varepsilon)\,\frac{M}{EI} . + On a stiff-``EA`` cable :math:`\lvert\mathbf{a}\rvert = 1 + \varepsilon` with :math:`\varepsilon \ll 1`, :math:`\mathbf{a}\cdot\mathbf{b} = O(\varepsilon)` and :math:`\mathbf{b} = \kappa\,\mathbf{n}`, so that :math:`\partial W/\partial\mathbf{a} = (EA\,\varepsilon - 2EI\kappa^2)\,\mathbf{t}` plus normal -terms, :math:`\mathbf{m} = EI\,\kappa\,\mathbf{n}`, and +terms, :math:`\mathbf{m} \approx EI\,\kappa\,\mathbf{n}`, and :math:`(\partial W/\partial\mathbf{b})'\cdot\mathbf{t} = EI\,\mathbf{b}'\cdot\mathbf{t} = -EI\kappa^2`. The axial internal force of the rod, the effective tension that balances the loads, and the bend moment are therefore @@ -260,7 +272,18 @@ bend moment are therefore N = \mathbf{F}\cdot\mathbf{t} = EA\,\varepsilon - EI\,\kappa^2, \qquad - M = \lvert \mathbf{m} - (\mathbf{m}\cdot\mathbf{t})\,\mathbf{t} \rvert = EI\,\kappa . + M = EI\,\kappa \quad\text{to } O(\varepsilon) . + +The curvature output is the geometric one, the reciprocal of the bend radius of the stretched +line: what the bend radius limit and the bending strain :math:`y\kappa` of a section point read. +``BendMom`` is :math:`M`, the moment that balances the forces, and at a line end the moment of +the end connection, whose potential depends on the direction of the tangent only. The check +``driver_curvature_stretch`` holds both against the elastica of this energy on a cable strained +up to 1.2 %: within 0.006 % and 0.001 % of the local value, against the 2.3 % and 1.2 % that +reading :math:`\lvert\mathbf{m}_\perp\rvert` as the moment gives. The stretch at a node is +:math:`\lvert\mathbf{r}'\rvert`, the length of its tangent, the one the discrete moment balance +holds with exactly; on a coarse stiff-``EA`` line it carries the pointwise oscillation above, +of order :math:`10^{-4}` of the curvature. The bending energy measures the curvature per deformed length, and its dependence on the stretch puts :math:`-EI\kappa^2` into the axial force: :math:`EA(\lvert\mathbf{r}'\rvert - 1)` @@ -306,12 +329,13 @@ it, and its tangent rows are :math:`\mathbf{m}` there. An interior node reports N_j = \mathbf{t}_j\cdot\tfrac{1}{2}\left(\mathbf{F}_j^- + \mathbf{F}_j^+\right), \qquad - M_j = \lvert \mathbf{m}_j - (\mathbf{m}_j\cdot\mathbf{t}_j)\,\mathbf{t}_j \rvert, + M_j = \lvert\mathbf{r}'_j\rvert\, + \lvert \mathbf{m}_j - (\mathbf{m}_j\cdot\mathbf{t}_j)\,\mathbf{t}_j \rvert, \qquad - \kappa_j = M_j / EI , + \kappa_j = \lvert\mathbf{r}'_j\rvert\, M_j / EI , with :math:`\mathbf{t}_j` the unit nodal tangent, :math:`\mathbf{m}_j` the mean of the two -one-sided values, and at a section junction the larger of :math:`M_j/EI` on the two sides (under +one-sided values, and at a section junction the larger of :math:`\kappa_j` on the two sides (under a bend stiffener :math:`EI` is the stiffness at the node, :ref:`theory-bend-stiffeners`). The recovered moment carries an error that does not depend on :math:`EI`, of the order of the largest consistent-load moment :math:`w L_0^2/12` of an element (a section junction leaves a jump of that diff --git a/python/tests/test_objects.py b/python/tests/test_objects.py index 1f728f14..07da0728 100644 --- a/python/tests/test_objects.py +++ b/python/tests/test_objects.py @@ -352,8 +352,9 @@ def test_finite_ei_views_report_bending(tmp_path): curv = line.curvature() moment = line.bend_moment() assert curv.max() > 0.01 - # EI is uniform (1.99e4 N.m^2) along this cable - np.testing.assert_allclose(moment, 1.99e4 * curv, rtol=1.0e-9, atol=1.0e-9) + # EI is uniform (1.99e4 N.m^2) along this cable; the moment is EI kappa over the stretch + stretch = 1.0 + (line.values("tension") + 1.99e4 * curv**2) / 4.69e8 + np.testing.assert_allclose(moment * stretch, 1.99e4 * curv, rtol=1.0e-6, atol=1.0e-9) assert line.bend_moment()[95] == model.channel("BendMom1N96") np.testing.assert_allclose(line.node_positions()[0], [5.0, 0.0, -14.0], atol=1e-9) np.testing.assert_allclose(line.node_positions()[-1], [125.0, 0.0, -80.0], atol=1e-9) diff --git a/src/CableDyn_DeckDriver.f90 b/src/CableDyn_DeckDriver.f90 index ffea06fb..c511ffb3 100644 --- a/src/CableDyn_DeckDriver.f90 +++ b/src/CableDyn_DeckDriver.f90 @@ -11151,7 +11151,8 @@ END FUNCTION sleeve_ei_at SUBROUTINE stiffened_element_peak(stf, qr, l0, ei_element, ei_line, ei_resolve, arc_a, arc_b, m0, m1, peak, u_peak, & ErrStat, ErrMsg) - !! Largest curvature of the section along one element under a bend stiffener, |M(u)|/EI(u). + !! Largest curvature of the section along one element under a bend stiffener, + !! (1 + e(u)) |M(u)|/EI(u), the stretch e running linearly between the end nodes. !! The element's interpolation bends under its one stiffness ei_element, so it gives the !! section neither its curvature nor, at its ends, its moment to second order: ei_element !! times the interpolation's curvature is the linear fit of the moment over the element. Its @@ -11177,9 +11178,9 @@ SUBROUTINE stiffened_element_peak(stf, qr, l0, ei_element, ei_line, ei_resolve, u_peak = CD_ZERO CALL CD_HermiteCable_Curvature(qr, l0, 0.5_wp, kappa, ErrStat, ErrMsg, kvec=kv) IF (ErrStat /= CD_HCABLE_OK) RETURN - ! the moment conjugate to the tangent is EI kappa n/|r'|^2, as the recovered ones are; the - ! curvature above has checked that the tangent r' at the midpoint does not vanish - bubble = ei_element*kv/SUM((1.5_wp*(qr(7:9) - qr(1:3))/l0 - 0.25_wp*(qr(4:6) + qr(10:12)))**2) - 0.5_wp*(m0 + m1) + ! the section moment is EI kappa n/|r'|, as the recovered ones are; the curvature above has + ! checked that the tangent r' at the midpoint does not vanish + bubble = ei_element*kv/NORM2(1.5_wp*(qr(7:9) - qr(1:3))/l0 - 0.25_wp*(qr(4:6) + qr(10:12))) - 0.5_wp*(m0 + m1) IF (.NOT. ALL(IEEE_IS_FINITE(bubble))) THEN ErrStat = CD_HCABLE_BADINPUT; ErrMsg = 'stiffened_element_peak: the element moments are not finite'; RETURN END IF @@ -11253,7 +11254,8 @@ REAL(wp) FUNCTION section_curvature(uu, from_b) ei = ei_line + sleeve_ei_at(stf, arc_a + (CD_ONE - uu)*l0, arc_b + uu*l0, from_b, STF_ARC_TOL*l0) IF (.NOT. ei > CD_ZERO) ei = ei_element IF (ei >= ei_resolve) THEN - section_curvature = NORM2((CD_ONE - uu)*m0 + uu*m1 + 6.0_wp*uu*(CD_ONE - uu)*bubble)/ei + section_curvature = ((CD_ONE - uu)*NORM2(qr(4:6)) + uu*NORM2(qr(10:12)))* & + NORM2((CD_ONE - uu)*m0 + uu*m1 + 6.0_wp*uu*(CD_ONE - uu)*bubble)/ei ELSE CALL CD_HermiteCable_Curvature(qr, l0, uu, section_curvature, ErrStat, ErrMsg) section_curvature = section_curvature*(ei_element/ei) diff --git a/src/CableDyn_HermiteCableDynamic.f90 b/src/CableDyn_HermiteCableDynamic.f90 index ca226fe0..5e0e0d7d 100644 --- a/src/CableDyn_HermiteCableDynamic.f90 +++ b/src/CableDyn_HermiteCableDynamic.f90 @@ -6075,8 +6075,8 @@ SUBROUTINE CD_HermiteCable_Dyn_Side_Curvature(model, e, node, curvature, ErrStat END IF CALL node_values(model, node, fm, mm, fp, mp, ten, kap, bend, ErrStat, ErrMsg, moment_vector=mv) IF (ErrStat /= CD_HCDYN_OK) RETURN - curvature = NORM2(mv)/ei - IF (PRESENT(kvec)) kvec = mv/ei + curvature = node_stretch(model, node)*NORM2(mv)/ei + IF (PRESENT(kvec)) kvec = node_stretch(model, node)*mv/ei END SUBROUTINE CD_HermiteCable_Dyn_Side_Curvature SUBROUTINE CD_HermiteCable_Dyn_Section_Curvature(model, ref_a, kx, ky, ErrStat, ErrMsg, frame) @@ -6679,8 +6679,11 @@ 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), 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 (under a bend + !! magnitude of CD_HermiteCable_Dyn_End_Force), the section moment + !! M_j = (1 + e_j) |m_j - (m_j . t_j) t_j|, 1 + e_j = |r'_j| the stretch at the node (m_j is the + !! conjugate of the tangent r', the section moment r' x m_j), + !! curvature(j) = (1 + e_j) M_j/EI, the geometric curvature of the stretched line, the larger + !! of the two sides at a section junction (under a bend !! stiffener EI is the stiffness at the node, not the element's; see node_values), 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 @@ -6690,10 +6693,10 @@ SUBROUTINE CD_HermiteCable_Dyn_Node_Recovery(model, ErrStat, ErrMsg, tension, cu !! mean of the axial components of element e's own two end forces. curvature_vector(:, j) !! is the curvature vector kappa n of node j in the solve frame, toward the centre of !! curvature, of length curvature(j): the section moment vector over EI on the side that - !! gives curvature(j), or the exact trace's vector. moment_vector(:, j) is the part of the - !! section moment normal to the tangent, of length M_j. Every output is optional (node arrays - !! nn, curvature_vector and moment_vector (3, nn), segment_tension ne); the routine allocates - !! nothing. + !! gives curvature(j), or the exact trace's vector. moment_vector(:, j) is the section moment + !! turned into the plane of the curvature, along kappa n, of length M_j. Every output is + !! optional (node arrays nn, curvature_vector and moment_vector (3, nn), segment_tension ne); + !! the routine allocates nothing. !! !! The condensed torque of a line restrained in torsion is not part of g_e: it adds a shear !! normal to the tangent and does no work on bending rotations, so tension and bend moment @@ -7113,14 +7116,19 @@ SUBROUTINE node_values(model, node, fm, mm, fp, mp, tension, curvature, bend_mom !! curvature vector kappa n (solve frame, toward the centre of curvature) of the side that !! gives curvature: the normal part of the section moment over EI, or the exact trace's. !! + !! The energy EI kappa^2/2 per unstretched length, with kappa the curvature of the stretched + !! line, gives m = EI kappa n/|r'|^2 as the conjugate of the tangent r'. Rotating the whole + !! line leaves the energy unchanged, so the moment that balances the forces is r' x m, of + !! length M = |r'| |m_perp| = EI kappa/|r'|, and kappa = |r'| M/EI. + !! !! The EI of a side is the element's, or on a stiffened line the limit of EI(s) at the node !! from that side (EI_node): the element's compliance average stands for its whole length !! and, under a taper, steps at every node although the section does not, so the softer side !! would overstate the curvature by about h |EI'|/(2 EI). The two limits differ only at a !! real step of the section (the bare tip, a junction of two line types, a step between !! stiffener rows), where the larger-side rule still applies. The exact trace of a side too - !! soft to recover takes the same stiffness (trace_scale). moment_vector is the normal part - !! of the section moment. + !! soft to recover takes the same stiffness (trace_scale). moment_vector is the section + !! moment turned into the plane of the curvature, along kappa n. TYPE(CD_HermiteCableDynType), INTENT(IN) :: model INTEGER, INTENT(IN) :: node REAL(wp), INTENT(IN) :: fm(3), mm(3), fp(3), mp(3) @@ -7129,6 +7137,7 @@ SUBROUTINE node_values(model, node, fm, mm, fp, mp, tension, curvature, bend_mom CHARACTER(*), INTENT(OUT) :: ErrMsg REAL(wp), INTENT(OUT), OPTIONAL :: sleeve_moment, curvature_vector(3), moment_vector(3) REAL(wp) :: t(3), f(3), m(3), fend(3), qe(12), trace, section_moment, cable_moment, m_perp(3), kv(3), ei_side + REAL(wp) :: stretch INTEGER :: e, side CHARACTER(200) :: em2 @@ -7148,7 +7157,8 @@ SUBROUTINE node_values(model, node, fm, mm, fp, mp, tension, curvature, bend_mom ELSE tension = DOT_PRODUCT(t, f) + 0.5_wp*(coupling_tension(model, node - 1) + coupling_tension(model, node)) END IF - m_perp = m - DOT_PRODUCT(m, t)*t + stretch = node_stretch(model, node) + m_perp = stretch*(m - DOT_PRODUCT(m, t)*t) section_moment = NORM2(m_perp) IF (PRESENT(moment_vector)) moment_vector = m_perp curvature = CD_ZERO @@ -7159,9 +7169,9 @@ SUBROUTINE node_values(model, node, fm, mm, fp, mp, tension, curvature, bend_mom IF (e < 1 .OR. e > model%ne) CYCLE ei_side = CD_HermiteCable_Dyn_Side_EI(model, e, node) IF (model%EI(e) > CD_ZERO .AND. ei_side >= model%ei_resolve) THEN - IF (PRESENT(curvature_vector) .AND. section_moment/ei_side > curvature) & - curvature_vector = m_perp/ei_side - curvature = MAX(curvature, section_moment/ei_side) + IF (PRESENT(curvature_vector) .AND. stretch*section_moment/ei_side > curvature) & + curvature_vector = stretch*m_perp/ei_side + curvature = MAX(curvature, stretch*section_moment/ei_side) IF (ABS(model%EI_line(e) - ei_side) <= CD_ZERO) THEN cable_moment = MAX(cable_moment, section_moment) ELSE @@ -7179,7 +7189,7 @@ SUBROUTINE node_values(model, node, fm, mm, fp, mp, tension, curvature, bend_mom IF (ABS(model%EI_line(e) - model%EI(e)) <= CD_ZERO) THEN cable_moment = MAX(cable_moment, section_moment) ELSE - cable_moment = MAX(cable_moment, model%EI_line(e)*trace) + cable_moment = MAX(cable_moment, model%EI_line(e)*trace/stretch) END IF END IF END DO @@ -7187,6 +7197,15 @@ SUBROUTINE node_values(model, node, fm, mm, fp, mp, tension, curvature, bend_mom IF (PRESENT(sleeve_moment)) sleeve_moment = MAX(CD_ZERO, section_moment - cable_moment) END SUBROUTINE node_values + PURE FUNCTION node_stretch(model, node) RESULT(stretch) + !! The stretch |r'| = 1 + e of a node, the length of its material tangent; 1 when degenerate. + TYPE(CD_HermiteCableDynType), INTENT(IN) :: model + INTEGER, INTENT(IN) :: node + REAL(wp) :: stretch + stretch = NORM2(model%q(6*(node - 1) + 4:6*(node - 1) + 6)) + IF (.NOT. stretch > CD_ZERO) stretch = CD_ONE + END FUNCTION node_stretch + PURE FUNCTION unit_tangent(model, node) RESULT(t) !! Unit material tangent of a node (internal arc direction); zero when degenerate. TYPE(CD_HermiteCableDynType), INTENT(IN) :: model diff --git a/tests/check_bend_stiffener.py b/tests/check_bend_stiffener.py index 60b494d2..13f123e6 100644 --- a/tests/check_bend_stiffener.py +++ b/tests/check_bend_stiffener.py @@ -4,11 +4,12 @@ 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, +table. The reference is an independent solution of the planar extensible elastica of the energy +EI kappa^2/2 + EA e^2/2 per unstretched length, kappa the curvature of the stretched line, - x' = (1 + e) cos(theta), z' = (1 + e) sin(theta), theta' = M / EI(s), + x' = (1 + e) cos(theta), z' = (1 + e) sin(theta), theta' = (1 + e)^2 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, + EA e - (1 + e) M^2 / EI(s) = n_x cos(theta) + n_z sin(theta), 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) @@ -142,6 +143,7 @@ "Curv1N41" "CurvX1N41" "CurvY1N41" +"Ten1N41" END """ @@ -170,6 +172,13 @@ def section(s: np.ndarray) -> tuple[np.ndarray, np.ndarray]: MOMENT_SCALE = 1.0e4 +def reference_stretch(y: np.ndarray, ei: np.ndarray) -> np.ndarray: + """1 + e of the elastica state y = (x, z, theta, M, n_x, n_z), scaled as in reference().""" + bending = (MOMENT_SCALE * y[3]) ** 2 / ei + axial = FORCE_SCALE * (y[4] * np.cos(y[2]) + y[5] * np.sin(y[2])) + return 1.0 + (axial + bending) / (EA - bending) + + 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: @@ -179,8 +188,8 @@ 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 = reference_stretch(y, ei) + return np.vstack([stretch * c, stretch * sn, stretch**2 * MOMENT_SCALE * m / ei, -stretch * FORCE_SCALE / MOMENT_SCALE * (c * nz - sn * nx), np.zeros_like(s), w / FORCE_SCALE]) @@ -279,7 +288,7 @@ def check_steps(exe: str, tmp: str) -> bool: the 5 m tip).""" cone = [(0.5, 0.60, 0.60), (4.5, 0.60, 0.22)] od_3m = 0.60 + 2.5 / 4.5 * (0.22 - 0.60) - # name, sections, stiffener rows, node, stiffness giving the curvature, cable stiffness of + # name, sections, stiffener rows, node, stiffness of the side of the curvature, cable stiffness of # the BendMom share (the larger share, on the side of the larger EI_cable) cases = [ ("bare tip", "1 cable 8.0 80\n", cone, 51, EI_CABLE, EI_CABLE), @@ -292,21 +301,22 @@ def check_steps(exe: str, tmp: str) -> bool: ] ok = True for name, sections, rows, node, ei, ei_cable in cases: - channels = [f"{c}1N{node}" for c in ("Curv", "BendMom", "StfMom")] + channels = [f"{c}1N{node}" for c in ("Curv", "BendMom", "StfMom", "Ten")] rc, log, _, out = solve(exe, tmp, "step", step_deck(sections, rows, channels)) if rc != 0: print(log[-3000:]) print(f"FAIL {name}: the driver did not solve the deck") return False - curvature, bend, stiffener = (out[c] for c in channels) + curvature, bend, stiffener, tension = (out[c] for c in channels) moment = bend + stiffener + stretch = 1.0 + tension / EA full = ei_cable == ei if name == "junction of two line types": # the stiffer type's own share: EI_cable M / (EI_cable + sleeve) share = ei_cable * moment / (ei_cable + sleeve(od_3m)) else: - share = moment if full else ei_cable * curvature - error = max(abs(curvature * ei / moment - 1.0), abs(bend / share - 1.0)) + share = moment if full else ei_cable * curvature / stretch + error = max(abs(curvature * ei / (stretch * moment) - 1.0), abs(bend / share - 1.0)) good = error <= IDENTITY_TOLERANCE print(f"{'ok ' if good else 'FAIL'} {name}: curvature {curvature:.5f} 1/m, " f"section moment {moment:.1f} N m, largest relative difference {error:.1e}") @@ -326,7 +336,7 @@ def check_refined(exe: str, tmp: str) -> bool: "(-) (-) (m) (m) (m) (m) (Pa) (kg/m^3)\n" f"4 A {length} {od0} {od1} {bore} {E_STF} {RHO_STF}\n" "--------------------- OPTIONS"), 1) - deck = deck.replace("FairTen4\n", "Curv4N3\nBendMom4N3\nStfMom4N3\n") + deck = deck.replace("FairTen4\n", "Curv4N3\nBendMom4N3\nStfMom4N3\nTen4N3\n") rc, log, rows, out = solve(exe, tmp, "refined", deck) if rc != 0 or "refined to" not in log: print(log[-3000:]) @@ -334,7 +344,8 @@ def check_refined(exe: str, tmp: str) -> bool: return False od = od0 + 40.0 / length * (od1 - od0) ei = 1.99e4 + E_STF * math.pi / 64.0 * (od**4 - bore**4) - error = abs(out["Curv4N3"] * ei / (out["BendMom4N3"] + out["StfMom4N3"]) - 1.0) + stretch = 1.0 + out["Ten4N3"] / 4.69e8 + error = abs(out["Curv4N3"] * ei / (stretch * (out["BendMom4N3"] + out["StfMom4N3"])) - 1.0) ok = error <= IDENTITY_TOLERANCE print(f"{'ok ' if ok else 'FAIL'} refined mesh of {len(rows) - 1} elements: curvature at the " f"deck node within {error:.1e} of the section moment over EI(s)") @@ -370,13 +381,15 @@ def check_soft_sleeve(exe: str, tmp: str) -> bool: return ok -def check_identities(out: dict[str, float], name: str = "node 41") -> bool: - """Curv, CurvX, CurvY, BendMom and StfMom of a node inside the cone (node 41, 2 m).""" +def check_identities(out: dict[str, float], name: str = "node 41", stretch: float | None = None) -> bool: + """Curv, CurvX, CurvY, BendMom and StfMom of a node inside the cone (node 41, 2 m): the moments + are the stiffnesses times the curvature over the stretch, 1 under OutputRecovery element.""" curvature = out["Curv1N41"] ei = float(section(np.array([2.0]))[0][0]) + stretch = 1.0 + out["Ten1N41"] / EA if stretch is None else stretch error = max(abs(math.hypot(out["CurvX1N41"], out["CurvY1N41"]) / curvature - 1.0), - abs(out["BendMom1N41"] / (EI_CABLE * curvature) - 1.0), - abs(out["StfMom1N41"] / ((ei - EI_CABLE) * curvature) - 1.0)) + abs(out["BendMom1N41"] * stretch / (EI_CABLE * curvature) - 1.0), + abs(out["StfMom1N41"] * stretch / ((ei - EI_CABLE) * curvature) - 1.0)) ok = error <= IDENTITY_TOLERANCE print(f"{'ok ' if ok else 'FAIL'} {name}: curvature components and moment shares agree " f"within {error:.1e}") @@ -392,8 +405,11 @@ def element_peaks(path: pathlib.Path, n: int) -> tuple[np.ndarray, np.ndarray, n return table[:, 2], table[:, 3], bend_check(read_output(path), 1.0, condition="storage", line_id=1).curvature[:n] -def reference_curvature(sol, s: np.ndarray) -> np.ndarray: - return np.abs(sol.sol(s)[3]) * MOMENT_SCALE / section(s)[0] +def reference_curvature(sol, s: np.ndarray, ei: np.ndarray | None = None) -> np.ndarray: + """Curvature of the stretched line, (1 + e) M / EI, at arc s.""" + ei = section(s)[0] if ei is None else ei + y = sol.sol(s) + return reference_stretch(y, ei) * np.abs(y[3]) * MOMENT_SCALE / ei def check_local(exe: str, tmp: str, sol) -> bool: @@ -409,7 +425,8 @@ def check_local(exe: str, tmp: str, sol) -> bool: s = np.arange(n + 1) * (WINDOW / n) ei = section(s)[0] ei[np.abs(s - tip) < 1.0e-9] = EI_CABLE - curvature = np.abs(sol.sol(s)[3]) * MOMENT_SCALE / ei + curvature = reference_curvature(sol, s, ei) + moment = np.abs(sol.sol(s)[3]) * MOMENT_SCALE / ei scale = np.maximum(curvature, LOCAL_FLOOR * root_curvature) for name, text in (("equilibrium nodes", deck), ("element nodes", deck.replace(*ELEMENT_RECOVERY))): rc, log, rows, _ = solve(exe, tmp, "local", text) @@ -425,7 +442,9 @@ def check_local(exe: str, tmp: str, sol) -> bool: # past the tip an element's own curvature is the section's, by definition covered = s < tip - 1.0e-9 if name == "element nodes" else s < np.inf error = np.where(covered, np.abs(data[:, 5] - curvature) / scale, 0.0) - moment_error = np.max(np.where(covered, np.abs(np.abs(data[:, 6]) / EI_CABLE - curvature) / scale, 0.0)) + # the element definition of the bend moment is EI times the curvature + shared = curvature if name == "element nodes" else moment + moment_error = np.max(np.where(covered, np.abs(np.abs(data[:, 6]) / EI_CABLE - shared) / scale, 0.0)) worst = int(np.argmax(error)) good = max(error[worst], moment_error) <= tolerance print(f"{'ok ' if good else 'FAIL'} {WINDOW / n:.3f} m elements, {name}: curvature within " @@ -466,7 +485,7 @@ def main() -> int: print(log[-3000:]) print("FAIL: the driver did not solve the deck under OutputRecovery element") return 1 - notes_ok = check_identities(element_out, "node 41, OutputRecovery element") and notes_ok + notes_ok = check_identities(element_out, "node 41, OutputRecovery element", 1.0) and notes_ok 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)") diff --git a/tests/check_curvature_stretch.py b/tests/check_curvature_stretch.py new file mode 100644 index 00000000..683a1714 --- /dev/null +++ b/tests/check_curvature_stretch.py @@ -0,0 +1,175 @@ +# SPDX-License-Identifier: Apache-2.0 +"""Curvature and section moment of a stretched finite-EI line (driver_curvature_stretch in +CMakeLists.txt). + +The element energy is EI kappa^2/2 + EA e^2/2 per unstretched length, with kappa the curvature of +the stretched line. A soft cable (axial strain up to 1.2 %) is clamped at End A and pinned at End +B, and the driver's static profile is compared with an independent solution of the planar +elastica of that energy, + + x' = (1 + e) cos(theta), z' = (1 + e) sin(theta), theta' = (1 + e)^2 M / EI, + M' = -(1 + e) (cos(theta) n_z - sin(theta) n_x), n' = (0, w), + EA e - (1 + e) M^2 / EI = n_x cos(theta) + n_z sin(theta), + +where M is the moment that balances the forces and kappa = theta'/(1 + e) = (1 + e) M / EI. The +strain is swept with EA, and at every node, relative to the local value floored at 1 % of the +peak, Curvature must follow kappa and BendMoment must follow M within 0.1 % (the old definitions +were low by 2e and e: 2.3 % and 1.2 % at the root). The exact trace of OutputRecovery element must +agree with the recovered curvature to the mesh's second order, within 0.5 %. The channels must +keep CurvX^2 + CurvY^2 = Curv^2 and BendMom (1 + e) = EI Curv. + +Usage: check_curvature_stretch.py +Exit 77 (skipped) when scipy is not installed. +""" + +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 + +G = 9.80665 +RHO_W = 1025.0 +D_CABLE = 0.20 +M_CABLE = 60.0 +EI = 2.0e4 +LENGTH = 60.0 +END_B = (45.0, -30.0) +W = (M_CABLE - RHO_W * math.pi / 4.0 * D_CABLE**2) * G +SWEEP = (2.0e6, 1.0e7, 5.0e8) +ELEMENTS = 480 +TOLERANCE = 1.0e-3 +TRACE_TOLERANCE = 5.0e-3 +IDENTITY_TOLERANCE = 2.0e-5 +FLOOR = 1.0e-2 +NODE = 10 +MOMENT_SCALE = 1.0e3 +FORCE_SCALE = 1.0e4 + + +def deck(ea: float, recovery: str = "") -> str: + return f"""--------------------- CableDyn Input File ------------------------------------ +Curvature of a stretched cable: clamped at End A, 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} 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 {LENGTH} {ELEMENTS} +--------------------- END CONNECTIONS ---------------------------------------- +LineID End Stiffness EzX EzY EzZ +(-) (-) (N-m/rad) (-) (-) (-) +1 A Rigid 1.0 0.0 0.0 +--------------------- OPTIONS ------------------------------------------------ +{recovery}9.80665 g +1025.0 rhoW +0.05 dtM +0.0 TMax +--------------------- OUTPUTS ------------------------------------------------ +"Curv1N{NODE}" +"CurvX1N{NODE}" +"CurvY1N{NODE}" +"BendMom1N{NODE}" +"Ten1N{NODE}" +END +""" + + +def solve(exe: str, tmp: str, name: str, text: str) -> tuple[np.ndarray, dict[str, float]]: + """The static profile (columns from ArcLength) and the output channels at t = 0.""" + path = pathlib.Path(tmp) / f"{name}.dat" + path.write_text(text) + root = pathlib.Path(tmp) / name + run = subprocess.run([exe, str(path), str(root)], capture_output=True, text=True) + if run.returncode != 0: + raise SystemExit(f"FAIL: the driver did not solve {name}:\n{(run.stdout + run.stderr)[-2000:]}") + rows = pathlib.Path(f"{root}.static.out").read_text().splitlines()[3:] + lines = pathlib.Path(f"{root}.out").read_text().splitlines() + channels = dict(zip(lines[1].split(), (float(v) for v in lines[2].split()))) + return np.array([[float(v) for v in r.split()[2:9]] for r in rows if r.strip()]), channels + + +def stretch_of(y: np.ndarray, ea: float) -> np.ndarray: + bending = (MOMENT_SCALE * y[3]) ** 2 / EI + axial = FORCE_SCALE * (y[4] * np.cos(y[2]) + y[5] * np.sin(y[2])) + return 1.0 + (axial + bending) / (ea - bending) + + +def reference(arc: np.ndarray, x_seed: np.ndarray, z_seed: np.ndarray, ea: float): + """The elastica, y = (x, z, theta, M, n_x, n_z), the moment and force scaled, seeded by the + driver's shape.""" + + def rhs(s, y): + stretch = stretch_of(y, ea) + c, sn = np.cos(y[2]), np.sin(y[2]) + return np.vstack([stretch * c, stretch * sn, stretch**2 * MOMENT_SCALE * y[3] / EI, + -stretch * FORCE_SCALE / MOMENT_SCALE * (c * y[5] - sn * y[4]), + np.zeros_like(s), np.full_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]]) + + mesh = np.linspace(0.0, LENGTH, 3001) + xs, zs = np.interp(mesh, arc, x_seed), np.interp(mesh, arc, 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.full_like(mesh, np.cos(th[-1])), np.sin(th)]) + sol = solve_bvp(rhs, bc, mesh, y0, tol=1.0e-7, max_nodes=400000) + if not sol.success: + raise SystemExit(f"FAIL: the reference elastica did not converge: {sol.message}") + return sol + + +def check(exe: str, tmp: str, ea: float) -> bool: + data, out = solve(exe, tmp, "stretch", deck(ea)) + element, _ = solve(exe, tmp, "trace", deck(ea, "element OutputRecovery\n")) + arc = np.arange(ELEMENTS + 1) * (LENGTH / ELEMENTS) + y = reference(arc, data[:, 1], data[:, 3] + 20.0, ea).sol(arc) + stretch = stretch_of(y, ea) + moment = np.abs(y[3]) * MOMENT_SCALE + kappa = stretch * moment / EI + scale = np.maximum(kappa, FLOOR * kappa.max()) + curvature_error = np.max(np.abs(data[:, 5] - kappa) / scale) + moment_error = np.max(np.abs(data[:, 6] - moment) / (EI * scale)) + trace_error = np.max(np.abs(element[:, 5] - data[:, 5]) / scale) + curvature = out[f"Curv1N{NODE}"] + identity_error = max( + abs(math.hypot(out[f"CurvX1N{NODE}"], out[f"CurvY1N{NODE}"]) / curvature - 1.0), + abs(out[f"BendMom1N{NODE}"] * stretch[NODE - 1] / (EI * curvature) - 1.0)) + ok = (curvature_error <= TOLERANCE and moment_error <= TOLERANCE and trace_error <= TRACE_TOLERANCE + and identity_error <= IDENTITY_TOLERANCE) + print(f"{'ok ' if ok else 'FAIL'} EA {ea:.1e} N, strain up to {100 * (stretch.max() - 1.0):.3f} %: " + f"curvature within {100 * curvature_error:.3f} %, moment within {100 * moment_error:.3f} %, " + f"exact trace within {100 * trace_error:.3f} %, identities within {identity_error:.1e}") + return ok + + +def main() -> int: + if solve_bvp is None: + print("SKIP: scipy is not installed") + return 77 + with tempfile.TemporaryDirectory() as tmp: + ok = all([check(sys.argv[1], tmp, ea) for ea in SWEEP]) + print("PASS" if ok else "FAIL") + return 0 if ok else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/test_driver_channels.f90 b/tests/test_driver_channels.f90 index 41bba17e..a5758b3e 100644 --- a/tests/test_driver_channels.f90 +++ b/tests/test_driver_channels.f90 @@ -279,7 +279,7 @@ SUBROUTINE test_static_config_file() WRITE (*, '(A,I0,A,F10.5)') 'static.out rows=', nrow, ' node3 curvature=', cv3 END SUBROUTINE test_static_config_file - ! ---- Part 4: finite-EI initialised deck; non-zero bend moment = EI*curvature ---- + ! ---- Part 4: finite-EI initialised deck; non-zero bend moment = EI*curvature/stretch ---- SUBROUTINE test_finite_dynamic_channels() INTEGER :: es, ok @@ -310,10 +310,12 @@ SUBROUTINE test_finite_dynamic_channels() ! element-interior maximum and is retained only for EI=0 nodal lines. CALL require(ABS(v(4) - ref_curv(n2, n3, n4)) > 5.0e-6_wp, & 'finite-EI Curv1N3 uses exact Hermite traces rather than Menger curvature') - ! BendMom = EI(node) * |kvec_current - kvec_reference|; straight reference => EI * geometric, - ! EI = 2.0e5 for this homogeneous line (the required straight-reference agreement check). - CALL require(ABS(v(5) - EI*v(4)) <= 1.0e-6_wp*ABS(v(5)) + 1.0e-9_wp, & - 'BendMom1N3 = EI * geometric curvature for the straight reference') + ! BendMom = EI(node) * |kvec_current - kvec_reference|/|r'|; straight reference => EI * geometric + ! over the stretch, EI = 2.0e5 for this homogeneous line (the required straight-reference + ! agreement check). The channels do not give |r'|, within 1e-4 of 1 on this stiff line: the + ! exact relation is gated in test_openfast_aggregate and check_curvature_stretch.py. + CALL require(ABS(v(5) - EI*v(4)) <= 1.0e-4_wp*ABS(v(5)) + 1.0e-9_wp, & + 'BendMom1N3 = EI * geometric curvature, over the stretch, for the straight reference') CALL require(v(5) > 1.0_wp, 'finite-EI BendMom1N3 is non-zero (bent vs straight reference)') CALL check_continuous_element_file('chan_fin.elements.out', 20, EI) ! tensile_safety True judges the element-mean axial force (the pointwise resultant of a diff --git a/tests/test_openfast_aggregate.f90 b/tests/test_openfast_aggregate.f90 index 04546052..8109673a 100644 --- a/tests/test_openfast_aggregate.f90 +++ b/tests/test_openfast_aggregate.f90 @@ -831,7 +831,7 @@ SUBROUTINE case_output_channels() CALL require(vFairTen2 > 1.0e3_wp, 'chan:FairTen2-physical') ! (b) cable Curv1N5 and the public HFMF nodal curvature share one definition at deck node 5: - ! the bend moment recovered from nodal equilibrium over EI by default, and under + ! the recovered bend moment times the stretch over EI by default, and under ! OutputRecovery element the exact one-sided Hermite curvature, the larger adjacent-element ! trace (CD_HFMF_Curvature_Exact). The channel follows the cable's setting. nnode = agg%cables(1)%line%nn @@ -864,10 +864,14 @@ SUBROUTINE case_output_channels() CALL require(ABS(vFairIncl1 - (decexp - 90.0_wp)) <= 1.0e-9_wp + 1.0e-6_wp*ABS(decexp), & 'chan:FairIncl1-signed-below-horizontal') - ! (d) cable BendMom1N5 == EI(node) * geometric curvature; the cable EI is uniform (bare/buoy - ! both 1.99e4), so nodal EI = EI_BARE and BendMom reduces to EI_BARE * Curv1N5. - CALL require(ABS(vBendMom1N5 - EI_BARE*vCurv1N5) <= 1.0e-6_wp*MAX(1.0_wp, ABS(EI_BARE*vCurv1N5)), & - 'chan:BendMom1N5-equals-EI-times-curv') + ! (d) cable BendMom1N5 == EI(node) * geometric curvature/|r'|, the stretch of the node's tangent; + ! the cable EI is uniform (bare/buoy both 1.99e4), so nodal EI = EI_BARE. + BLOCK + INTEGER :: i0 + i0 = 6*(nnode - 5) + CALL require(ABS(vBendMom1N5*NORM2(agg%cables(1)%line%q(i0 + 4:i0 + 6)) - EI_BARE*vCurv1N5) <= & + 1.0e-9_wp*MAX(1.0_wp, ABS(EI_BARE*vCurv1N5)), 'chan:BendMom1N5-equals-EI-times-curv-over-stretch') + END BLOCK WRITE (*, '(A,ES12.4,A,ES12.4,A,ES12.4)') 'channels: FairTen2 = ', vFairTen2, ' N; Curv1N5 = ', & vCurv1N5, ' 1/m; FairAngle1 = ', vFairAngle1 diff --git a/validation/scripts/orcaflex_bend_stiffener.py b/validation/scripts/orcaflex_bend_stiffener.py index 627a0bcc..bb78762d 100644 --- a/validation/scripts/orcaflex_bend_stiffener.py +++ b/validation/scripts/orcaflex_bend_stiffener.py @@ -18,8 +18,8 @@ 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.6 % of the root curvature. +Result (OrcFxAPI 11.6d): root curvature 0.4522 /m (CableDyn 0.4528), cable bend moment at the +root 9.04 kN.m; CableDyn's curvature follows the profile within 0.5 % of the root curvature. OrcaFlex SI (m, te, kN). Usage: python validation/scripts/orcaflex_bend_stiffener.py [output.csv]