diff --git a/CHANGELOG.md b/CHANGELOG.md index 8ea3bb5..ee1ed8d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,16 @@ All notable changes to CableDyn are recorded in this file. The format follows ### Added +- Documentation of the negative `BA` (damping ratio) form: it damps each element's own axial + mode at the given ratio, so the line's axial modes are damped at `zeta n pi / (2N)` and less on + a finer mesh (`BA = -1`: 1.6 % for the first mode on 100 elements, 0.8 % on 200); the manual + gives the `BA` in N·s for a mesh-independent material damping ratio. +- Seabed-friction resolution note: friction acts at the nodes, so on a coarse grounded mesh the + tension of the laid part and near touchdown steps by up to mu w L per node as the line lands + and lifts (a 10 m lazy-wave mesh swung a grounded segment tension by 1600 N where 1 m elements + give 120 N). At initialisation a finite-EI line whose grounded elements lock more than 15 % of + the touchdown tension is noted with the grounded element length that resolves it (flat + seabed; standalone and OpenFAST). Results do not change. - The run-time profile (`CABLEDYN_PROFILE=1`, standalone driver and OpenFAST) splits the steps of `EI = 0` lines into residual-only evaluations, residual and tangent evaluations, LU factorisations and back-substitutions, with their counts and times, and separates the @@ -49,6 +59,11 @@ All notable changes to CableDyn are recorded in this file. The format follows ### Changed +- A finite-EI static solve in a current that fails because seabed friction holds the grounded + run while the current pushes the suspended span into the touchdown now stops after the + catenary-continuation route has diagnosed it, as the frictionless grounded-run case already + did, instead of also running the mesh-sequenced route to its iteration budget (a 1000 m lazy + wave in 440 m of water: 22 s to 17.7 s, same message). Converging decks do not change. - The Windows release executables (`CableDyn_driver.exe`, `openfast.exe`) are compiled for processors with AVX2 (x86-64-v3: Intel Haswell and AMD Excavator or newer; CMake option `CABLEDYN_TARGET_AVX2`), which makes the driver 13-33% faster on the benchmark decks @@ -75,6 +90,13 @@ All notable changes to CableDyn are recorded in this file. The format follows from the step and a fresh evaluation differ slightly; the restore discarded the cache and took the fresh path. The snapshot now keeps the cached force with the state, and the restore brings it back unless a load setter invalidated it since. Runs without restores do not change. +- Finite-EI node tension, curvature and bend moment at touchdown (node channels, range graphs): + the recovery spread a grounded node's floor reaction over the touchdown element, which is on + the floor at one end only, and understated the tension of the last suspended node by the + floor reaction of half an element (5.3 % on a 10 m lazy-wave mesh). A node's floor reaction is + now shared by the grounded elements beside it; the touchdown element takes none. The same + coarse mesh now follows a four times finer one within 3 % through touchdown. Line end forces + (`FairTen`, `AnchTen`) and the solution do not change. - OpenFAST with `dtM` larger than `DT`: an `openfast -restart` re-sampled the SeaState fields of the committed step at its end positions instead of restoring the fields the step held (sampled at its start positions), so the loads served until the next solve, and the run after them, diff --git a/CMakeLists.txt b/CMakeLists.txt index 28b84a8..f1b2ba5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2284,6 +2284,15 @@ if(Python3_Interpreter_FOUND) set_tests_properties(driver_ei0_fold_layouts PROPERTIES LABELS "python;driver;integration;static" ENVIRONMENT "${_cabledyn_python_test_environment}") + # Node tension through touchdown on a coarse mesh (10 m grounded elements) follows a four + # times finer reference within 3.5 %: the touchdown element, grounded at one end only, takes + # no share of its grounded node's floor reaction. + add_test(NAME driver_touchdown_recovery + COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/tests/check_touchdown_recovery.py + $) + set_tests_properties(driver_touchdown_recovery PROPERTIES + LABELS "python;driver;integration;static" + 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/doc/conventions.rst b/doc/conventions.rst index 5b5b8b8..811baa1 100644 --- a/doc/conventions.rst +++ b/doc/conventions.rst @@ -48,7 +48,9 @@ are in degrees. - ``0`` selects the ``EI = 0`` cable path * - Axial damping ``BA`` - N·s - - a negative value is a damping ratio :math:`-\zeta` (dimensionless); see :doc:`theory` + - a negative value is a damping ratio :math:`-\zeta` (dimensionless) of each element's axial + mode, so the damping of the line's own axial modes depends on the mesh; give ``BA`` in + N·s for mesh-independent damping (see :doc:`theory`) * - End-connection stiffness - N·m/rad - ``END CONNECTIONS`` table diff --git a/doc/options.rst b/doc/options.rst index fd1f4a3..3d6090d 100644 --- a/doc/options.rst +++ b/doc/options.rst @@ -109,7 +109,10 @@ Environment and seabed ``dtM``/``TMax`` in the standalone driver. Line nodes carry stick-slip friction springs; in a current the static solve holds the line with friction springs from its still-water laid shape, and the march continues them (``staticsFriction asLaid`` holds it toward its - as-laid position instead, in still water too). The isotropic shorthand of + as-laid position instead, in still water too). Friction acts at the nodes, so a coarse + grounded mesh steps the tension of the laid part; the initialisation report notes a + finite-EI line whose grounded elements are too long (*Mesh resolution* in :doc:`theory`). + The isotropic shorthand of ``frictionMuAxial``/``frictionMuLateral``. Not supported (the deck is rejected) with a ``FAILURE`` section, on a deck with ``Connect``/``Free`` points or ``Point3`` buoys, or on a deck whose lines attach to ``Rigid6`` bodies or rods unless it runs on the default @@ -346,7 +349,12 @@ rejected as unknown. The dynamic Newton controls are set with ``dynamic_solver`` - logical - ``False`` - ``True`` enables guarded within-step tangent reuse on every dynamic solver the deck builds. - It does not relax convergence tolerances and should be qualified against full Newton + It does not relax convergence tolerances and should be qualified against full Newton. + Measured against full Newton: 8-12 % faster on free-hanging finite-EI cables (lazy wave, + buoyed cable), about 1.5 times slower on finite-EI lines on the seabed (contact fails the + reuse contraction test and rebuilds the tangent), 3 % slower to 5 % faster on ``EI = 0`` + systems; + results move at the solver tolerance (Bergdahl mean cycle maxima within 0.16 %) - both - ``False modified_newton`` * - ``recovery_max_substeps`` / ``recoverymaxsubsteps`` diff --git a/doc/theory.rst b/doc/theory.rst index b696d11..312d462 100644 --- a/doc/theory.rst +++ b/doc/theory.rst @@ -456,6 +456,24 @@ element as using the final element length after any mesh refinement. The damping share is included in the reported tension. +**The damping ratio form depends on the mesh.** ``BA`` resolved this way is proportional to the +element length, so it is a property of the discretisation, not of the material. For a uniform +line of :math:`N` elements the axial mode :math:`n` (wavenumber :math:`n\pi/L`) has the damping +ratio + +.. math:: + + \zeta_n = \frac{BA\,\omega_n}{2\,EA} = \zeta\,\frac{n\pi}{2N} , + +so the element-scale modes (:math:`n \approx N`) are damped at up to :math:`(\pi/2)\,\zeta`, as intended +(the convention suppresses segment-scale axial noise, the same one as MoorDyn's negative +``BA``), while the line's global axial modes are damped much less, and less again on a finer +mesh: ``BA = -1`` on 100 elements gives the first axial mode 1.6 %, and 0.8 % on 200 elements. +Results that depend on the axial damping (tension peaks of snatch loads, the decay of axial +oscillations) therefore change with the mesh. For mesh-independent damping give ``BA`` in N·s: +a material damping ratio :math:`\zeta_1` of the first axial mode of a line of length :math:`L` is +:math:`BA = 2\,\zeta_1\,EA/\omega_1` with :math:`\omega_1 = (\pi/L)\sqrt{EA/\rho_A}`. + Viscoelastic rope (MoorDyn ``ElasticMod`` 2 and 3) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -884,6 +902,16 @@ failed or subdivided step, join the aggregate snapshot, and are part of the Open checkpoint. The force and its position and velocity Jacobians (through :math:`N`) are exact. Unlike a velocity-regularised Coulomb law, the spring holds its force at rest. +**Mesh resolution.** Because friction acts at the nodes, the tension along the laid part changes +in steps of up to :math:`\mu w L` (:math:`w` the submerged weight per length, :math:`L` the +element length) instead of continuously, and as the line lands and lifts at touchdown each node +locks or releases that much. On a coarse grounded mesh the tension of the laid part and near +touchdown then swings in time by about that amount. On a 600 m lazy-wave cable in 200 m of +water (:math:`\mu = 0.35`) the standard deviation of a grounded segment tension was 491, 136, 42 +and 33 N with 10, 5, 2.5 and 1 m elements, where :math:`\mu w L` is 46, 22, 11 and 4.5 % of the +touchdown tension. When a finite-EI line's grounded elements exceed 15 %, the initialisation +report notes it with the grounded element length that resolves it. + .. _static-friction: **Friction in statics.** With seabed friction a static equilibrium is not unique: a grounded diff --git a/src/CableDyn_DeckDriver.f90 b/src/CableDyn_DeckDriver.f90 index ce6604e..c43aef6 100644 --- a/src/CableDyn_DeckDriver.f90 +++ b/src/CableDyn_DeckDriver.f90 @@ -2235,12 +2235,15 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, IF (pass == 3) THEN IF (.NOT. symmetric_current_retry()) EXIT route_pass END IF - ! A catenary-continuation route that stopped on the frictionless grounded-run - ! diagnosis (the current drives the run on the seabed into compression) is not - ! rescued by the mesh-sequenced route, which reaches the same planar limit after - ! up to a minute; only the symmetry-breaking passes can. + ! A catenary-continuation route that stopped on a grounded-run diagnosis (the current + ! drives the run on the seabed into compression; with friction, the friction holds the + ! grounded run while the current pushes the suspended span into the touchdown) is not + ! rescued by the mesh-sequenced route, which reaches the same limit only after its + ! iteration budget (17 s on a 1000 m lazy wave in 440 m of water); only the + ! symmetry-breaking passes can. IF (pass == 2 .AND. .NOT. opts%cable_statics_sequenced .AND. & - INDEX(primary_msg, 'the run on the seabed') > 0) THEN + (INDEX(primary_msg, 'the run on the seabed') > 0 .OR. & + INDEX(primary_msg, 'the friction holds the grounded run') > 0)) THEN sequence_msg = 'not run: the catenary-continuation route diagnosed the grounded run' CYCLE route_pass END IF @@ -3056,6 +3059,8 @@ RECURSIVE SUBROUTINE build_one_hermite_cable(ln, sections, types, points, opts, IF (es /= CD_HCDYN_OK) THEN CALL fail(ErrStat, ErrMsg, 'finite-EI Hermite seabed recompute wiring failed: '//TRIM(em)); RETURN END IF + IF (opts%friction_mu > CD_ZERO .AND. .NOT. PRESENT(bathymetry) .AND. & + (opts%tmax > CD_ZERO .OR. .NOT. standalone_mode)) CALL friction_resolution_note() ! The static friction force of every node becomes the initial force of its ! elastic-plastic spring, so the march starts in the static equilibrium. IF (static_friction) THEN @@ -4370,6 +4375,72 @@ SUBROUTINE report_note(text) IF (PRESENT(init_note)) CALL append_note_line(init_note, text) END SUBROUTINE report_note + 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, + !! w its submerged weight per length, L its length), and the tension of the laid part and + !! near touchdown steps by that much. A note, on a flat seabed, when a grounded element's mu w L + !! exceeds FRICTION_TENSION_FRACTION of the static tension at touchdown, with the element + !! length that keeps it below; nothing is changed. On a 600 m lazy wave in 200 m of water + !! (mu 0.35) the standard deviation of a grounded segment tension was 491, 136, 42 and 33 N + !! with 10, 5, 2.5 and 1 m elements, where mu w L is 46, 22, 11 and 4.5 % of the touchdown + !! tension: converged below about 15 %. + REAL(wp), PARAMETER :: FRICTION_TENSION_FRACTION = 0.15_wp + REAL(wp) :: mu, quantum, q_max, l_max, w_max, t_td, gap_tol, l_res + REAL(wp), ALLOCATABLE :: t_nodes(:) + INTEGER :: k, nnd, ned, i_td, es_rec + LOGICAL, ALLOCATABLE :: grounded(:) + CHARACTER(32) :: b1, b2, b3, b4, b5 + CHARACTER(512) :: text + CHARACTER(1024) :: em_rec + mu = opts%friction_mu + IF (opts%friction_mu_axial > CD_ZERO) mu = opts%friction_mu_axial + nnd = cable%line%nn + ned = cable%line%ne + IF (nnd < 2) RETURN + ALLOCATE (grounded(nnd)) + gap_tol = 0.01_wp*MAXVAL(diam) + DO k = 1, nnd + grounded(k) = cable%line%q(6*(k - 1) + 3) <= contact_z + gap_tol + END DO + q_max = CD_ZERO; l_max = CD_ZERO; w_max = CD_ZERO + DO k = 1, ned + IF (.NOT. (grounded(k) .AND. grounded(k + 1))) CYCLE + quantum = mu*MAX(cable%line%w(k), CD_ZERO)*cable%line%l0(k) + IF (quantum > q_max) THEN + q_max = quantum; l_max = cable%line%l0(k); w_max = cable%line%w(k) + END IF + END DO + IF (.NOT. (q_max > CD_ZERO)) RETURN + ! the touchdown node: a grounded node beside a suspended one (the first from End B, + ! internal node 1), and its tension from nodal equilibrium, as the tension channels report + ! it (the elastic EA (|r'| - 1) of a node is dominated by the discretisation of a stiff EA) + i_td = 0 + DO k = 1, nnd - 1 + IF (grounded(k) .NEQV. grounded(k + 1)) THEN + i_td = MERGE(k, k + 1, grounded(k)); EXIT + END IF + END DO + IF (i_td < 1) RETURN + ALLOCATE (t_nodes(nnd)) + CALL CD_HermiteCable_Dyn_Node_Recovery(cable%line, es_rec, em_rec, tension=t_nodes) + IF (es_rec /= CD_HCDYN_OK) RETURN + t_td = t_nodes(i_td) + IF (.NOT. (t_td > CD_ZERO) .OR. q_max <= FRICTION_TENSION_FRACTION*t_td) RETURN + l_res = FRICTION_TENSION_FRACTION*t_td/(mu*w_max) + WRITE (b1, '(F12.1)') l_max + WRITE (b2, '(I0)') NINT(q_max) + WRITE (b3, '(I0)') NINT(100.0_wp*q_max/t_td) + WRITE (b4, '(I0)') NINT(t_td) + WRITE (b5, '(F12.2)') l_res + text = 'line '//TRIM(int_to_str(ln%id))//': seabed friction acts at the nodes, so a grounded element of '// & + TRIM(ADJUSTL(b1))//' m locks or releases up to '//TRIM(ADJUSTL(b2))//' N (friction coefficient x '// & + 'submerged weight x length) as the line lands and lifts, '//TRIM(ADJUSTL(b3))//' % of the '// & + TRIM(ADJUSTL(b4))//' N tension at touchdown: on this mesh the tensions of the laid part and near '// & + 'touchdown step by that much. Grounded elements of at most '//TRIM(ADJUSTL(b5))//' m resolve it' + CALL report_note(text) + END SUBROUTINE friction_resolution_note + END SUBROUTINE build_one_hermite_cable SUBROUTINE attachment_nodes(l0, arc_public, node, max_offset) diff --git a/src/CableDyn_HermiteCableDynamic.f90 b/src/CableDyn_HermiteCableDynamic.f90 index 3e5b24b..edf6b65 100644 --- a/src/CableDyn_HermiteCableDynamic.f90 +++ b/src/CableDyn_HermiteCableDynamic.f90 @@ -494,14 +494,16 @@ MODULE CableDyn_HermiteCableDynamic ! Cache of the element end residuals behind the node outputs (CD_HermiteCable_Dyn_Node_Values). ! Element e's entry rc_g(:, e) is valid for the time and configuration epoch the cache was ! keyed at (rc_t, rc_key_epoch) and for the element's own 12 q, v and a DOFs it was evaluated - ! at (rc_q/rc_v/rc_a(:, e)), the only state its residual depends on, so a node query costs two - ! element comparisons and the channels of one output row share one evaluation per element. + ! at (rc_q/rc_v/rc_a(:, e)), and for the positions of the neighbouring nodes e - 1 and e + 2 + ! (rc_qnb(:, e); their floor contact decides the share of the floor reaction the element + ! carries), the only state its residual depends on, so a node query costs two element + ! comparisons and the channels of one output row share one evaluation per element. ! Allocated once at initialisation; rc_epoch counts the loads configured after it (drag, ! axial damping, waves, added mass, held fluid, contact). INTEGER :: rc_epoch = 0, rc_key_epoch = -1 LOGICAL :: rc_key_ok = .FALSE. REAL(wp) :: rc_t = CD_ZERO - REAL(wp), ALLOCATABLE :: rc_q(:, :), rc_v(:, :), rc_a(:, :), rc_g(:, :) + REAL(wp), ALLOCATABLE :: rc_q(:, :), rc_v(:, :), rc_a(:, :), rc_g(:, :), rc_qnb(:, :) LOGICAL, ALLOCATABLE :: rc_have(:) LOGICAL :: initialized = .FALSE. END TYPE CD_HermiteCableDynType @@ -737,7 +739,8 @@ SUBROUTINE CD_HermiteCable_Dyn_Init(model, l0, EA, EI, rho_a, w, seed, fixed_dof ! state: at-rest seed, v = 0, consistent a0 ALLOCATE (model%q(ndof), model%v(ndof), model%a(ndof)) - ALLOCATE (model%rc_q(12, ne), model%rc_v(12, ne), model%rc_a(12, ne), model%rc_g(12, ne), model%rc_have(ne)) + ALLOCATE (model%rc_q(12, ne), model%rc_v(12, ne), model%rc_a(12, ne), model%rc_g(12, ne), model%rc_have(ne), & + model%rc_qnb(6, ne)) model%rc_have = .FALSE. model%rc_key_ok = .FALSE. model%q = seed @@ -5509,14 +5512,19 @@ SUBROUTINE CD_HermiteCable_Dyn_Node_Values(model, node, tension, curvature, bend SUBROUTINE cached_end_residual(e, ge) INTEGER, INTENT(IN) :: e REAL(wp), INTENT(OUT) :: ge(12) - REAL(wp) :: gl(12) + REAL(wp) :: gl(12), qnb(6) INTEGER :: i0 - ! the element's residual depends on its own two nodes' q, v and a only (element_end_residual) + ! the element's residual depends on its own two nodes' q, v and a and on the floor contact + ! of the neighbouring nodes, a function of their positions (element_end_residual) i0 = 6*(e - 1) + qnb = CD_ZERO + IF (e > 1) qnb(1:3) = model%q(i0 - 5:i0 - 3) + IF (e + 2 <= model%nn) qnb(4:6) = model%q(i0 + 13:i0 + 15) IF (model%rc_have(e)) THEN IF (ALL(reals_equal(model%rc_q(:, e), model%q(i0 + 1:i0 + 12))) .AND. & ALL(reals_equal(model%rc_v(:, e), model%v(i0 + 1:i0 + 12))) .AND. & - ALL(reals_equal(model%rc_a(:, e), model%a(i0 + 1:i0 + 12)))) THEN + ALL(reals_equal(model%rc_a(:, e), model%a(i0 + 1:i0 + 12))) .AND. & + ALL(reals_equal(model%rc_qnb(:, e), qnb))) THEN ge = model%rc_g(:, e); RETURN END IF END IF @@ -5528,6 +5536,7 @@ SUBROUTINE cached_end_residual(e, ge) model%rc_q(:, e) = model%q(i0 + 1:i0 + 12) model%rc_v(:, e) = model%v(i0 + 1:i0 + 12) model%rc_a(:, e) = model%a(i0 + 1:i0 + 12) + model%rc_qnb(:, e) = qnb model%rc_g(:, e) = gl model%rc_have(e) = .TRUE. ge = gl @@ -5881,16 +5890,24 @@ SUBROUTINE element_end_residual(model, e, g, ErrStat, ErrMsg) IF (es /= CD_HCDYN_OK) GOTO 900 g = g - fx END IF - ! Seabed contact acts at the nodes; spread over each node's tributary length it is the + ! Seabed contact acts at the nodes; spread over the grounded line next to each node it is the ! distributed floor reaction it stands for. Its consistent share is taken out of the element, - ! so a line resting on the floor carries no fixed-end moment of its own weight. + ! so a line resting on the floor carries no fixed-end moment of its own weight. A node's + ! reaction is shared by the half-elements beside it whose far node is on the floor as well + ! (grounded elements), by length. The touchdown element, on the floor at one end only, is + ! mostly suspended and takes none of it; spreading the reaction onto it understated the + ! tension of the last suspended node by the floor reaction of half an element (5 % on a + ! 10 m lazy-wave mesh). A node with no grounded element beside it (a node just touching + ! down, or an isolated contact) shares its reaction by tributary length. IF (model%has_contact .OR. model%kn > CD_ZERO) THEN CALL node_contact_normal(model, e, fca, es, em2) IF (es /= CD_HCDYN_OK) GOTO 900 CALL node_contact_normal(model, e + 1, fcb, es, em2) IF (es /= CD_HCDYN_OK) GOTO 900 - IF (model%trib(e) > CD_ZERO) fca = fca*l0/model%trib(e) - IF (model%trib(e + 1) > CD_ZERO) fcb = fcb*l0/model%trib(e + 1) + CALL grounded_share(e, e, fca, es) + IF (es /= CD_HCDYN_OK) GOTO 900 + CALL grounded_share(e + 1, e, fcb, es) + IF (es /= CD_HCDYN_OK) GOTO 900 ! A held end whose end element lies on the floor (assemble_residual): the seabed carries the ! element's weight with its consistent end moment, which the held end node's own contact ! (a point reaction at an anchor sitting on the floor) cannot. The floor reaction of the @@ -5937,6 +5954,75 @@ SUBROUTINE element_end_residual(model, e, g, ErrStat, ErrMsg) RETURN 900 ErrStat = CD_HCDYN_NOCONVERGE ErrMsg = 'CD_HermiteCable_Dyn_Node_Recovery: element '//TRIM(em2) + + CONTAINS + + SUBROUTINE grounded_share(node, el, fc, esx) + !! Scale the floor reaction fc of node (an end of element el) to el's share of it, as the + !! factor l0(el)/(grounded length beside the node) that the half-element shape integrals + !! below expect: zero when el is not grounded but another element beside the node is, the + !! tributary-length share when no element beside the node is grounded. + INTEGER, INTENT(IN) :: node, el + REAL(wp), INTENT(INOUT) :: fc(3) + INTEGER, INTENT(OUT) :: esx + REAL(wp) :: grounded_len + LOGICAL :: on_el, grounded + INTEGER :: k + esx = CD_HCDYN_OK + IF (.NOT. ANY(ABS(fc) > CD_ZERO)) RETURN + grounded_len = CD_ZERO + on_el = .FALSE. + DO k = MAX(1, node - 1), MIN(model%ne, node) + CALL element_on_floor(k, grounded, esx) + IF (esx /= CD_HCDYN_OK) RETURN + IF (grounded) THEN + grounded_len = grounded_len + 0.5_wp*model%l0(k) + IF (k == el) on_el = .TRUE. + END IF + END DO + IF (grounded_len > CD_ZERO) THEN + IF (on_el) THEN + fc = fc*model%l0(el)/grounded_len + ELSE + fc = CD_ZERO + END IF + ELSE IF (model%trib(node) > CD_ZERO) THEN + fc = fc*model%l0(el)/model%trib(node) + END IF + END SUBROUTINE grounded_share + + SUBROUTINE element_on_floor(k, grounded, esx) + !! Element k is grounded when both its nodes are on the floor (node_on_floor). + INTEGER, INTENT(IN) :: k + LOGICAL, INTENT(OUT) :: grounded + INTEGER, INTENT(OUT) :: esx + LOGICAL :: on1 + grounded = .FALSE. + CALL node_on_floor(k, on1, esx) + IF (esx /= CD_HCDYN_OK .OR. .NOT. on1) RETURN + CALL node_on_floor(k + 1, grounded, esx) + END SUBROUTINE element_on_floor + + SUBROUTINE node_on_floor(node, on, esx) + !! A node is on the floor when it carries a floor reaction, or when it is a held line end + !! whose end element the floor supports (held_end_floor_support): an anchor held exactly + !! on the floor has no contact force of its own, yet its laid end element is grounded. + INTEGER, INTENT(IN) :: node + LOGICAL, INTENT(OUT) :: on + INTEGER, INTENT(OUT) :: esx + REAL(wp) :: fn(3), sup(3), dsup(3), nn_in(3) + INTEGER :: in_node + on = .FALSE. + CALL node_contact_normal(model, node, fn, esx, em2) + IF (esx /= CD_HCDYN_OK) RETURN + on = ANY(ABS(fn) > CD_ZERO) + IF (on) RETURN + IF (node == 1) THEN + on = held_end_floor_support(model, model%q, 1, sup, dsup, in_node, nn_in) + ELSE IF (node == model%nn) THEN + on = held_end_floor_support(model, model%q, 2, sup, dsup, in_node, nn_in) + END IF + END SUBROUTINE node_on_floor END SUBROUTINE element_end_residual PURE REAL(wp) FUNCTION CD_HermiteCable_Dyn_Max_Step_Rotation(model) RESULT(angle) @@ -6509,7 +6595,7 @@ SUBROUTINE CD_HermiteCable_Dyn_End(model) model%ws_elem_Kwave, model%ws_elem_Kheld, & model%ws_elem_es, model%ws_elem_em) END IF - IF (ALLOCATED(model%rc_q)) DEALLOCATE (model%rc_q, model%rc_v, model%rc_a, model%rc_g, model%rc_have) + IF (ALLOCATED(model%rc_q)) DEALLOCATE (model%rc_q, model%rc_v, model%rc_a, model%rc_g, model%rc_have, model%rc_qnb) model%rc_key_ok = .FALSE. IF (ALLOCATED(model%freemask)) DEALLOCATE (model%freemask) IF (ALLOCATED(model%solve_mask)) DEALLOCATE (model%solve_mask) diff --git a/tests/check_driver_advisories.cmake b/tests/check_driver_advisories.cmake index d34e8a4..4b653d4 100644 --- a/tests/check_driver_advisories.cmake +++ b/tests/check_driver_advisories.cmake @@ -5,7 +5,9 @@ # * a slack anchored line on a deck with no seabed (an OpenFAST MooringFile run # standalone) prints a no-seabed warning; a deck with a water depth does not; # * a dtM that resolves a line's axial period (2 x 410 m / sqrt(1.674e9 / 390) = 0.396 s) in -# fewer than 6 steps prints the dtM advisory: 4.0 steps at 0.1 s, none at 0.05 s (7.9). +# fewer than 6 steps prints the dtM advisory: 4.0 steps at 0.1 s, none at 0.05 s (7.9); +# * a finite-EI cable whose grounded elements lock more than 15 % of the touchdown tension in +# seabed friction prints the friction resolution note; a finer grounded mesh does not. foreach(_required DRIVER WORK NO_SEABED_DECK SEABED_DECK) if(NOT DEFINED ${_required} OR "${${_required}}" STREQUAL "") message(FATAL_ERROR "check_driver_advisories: missing ${_required}") @@ -111,3 +113,57 @@ require_text("no seabed" "${_no_seabed_log}" "below its Fixed anchor and the dec require_text("no seabed" "${_no_seabed_log}" "for a standalone run add WtrDpth") run_driver("seabed" "${SEABED_DECK}" "${WORK}/seabed" _seabed_log) forbid_text("seabed" "${_seabed_log}" "WARNING") + +# The seabed-friction resolution note: a lazy-wave cable whose grounded run has 10 m elements +# locks up to mu w L = 0.35 x 157.8 N/m x 10 m = 552 N per node, 46 % of the 1204 N tension at +# touchdown, and is told the grounded element length that resolves it; with 2.5 m elements +# (11 %, below the 15 % limit) it is not. +set(_friction_deck [=[--------------------- CableDyn Input File ------------------------------------ +Driver advisory check: lazy-wave cable with seabed friction +--------------------- LINE TYPES --------------------------------------------- +Name Diam Mass EA BA EI GAs GJ Irt Irn Cdn Cdt Can Cat +(-) (m) (kg/m) (N) (-) (Nm2) (N) (Nm2) (kgm) (kgm) (-) (-) (-) (-) +bare 0.16 36.70 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +buoy 0.29 59.53 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +--------------------- POINTS ------------------------------------------------- +ID Type X Y Z +(-) (-) (m) (m) (m) +7 Coupled 30.0 0.0 -20.0 +8 Fixed 500.0 0.0 -200.0 +--------------------- LINES -------------------------------------------------- +ID NodeA NodeB Outputs +(-) (-) (-) (-) +4 7 8 - +--------------------- SECTIONS ----------------------------------------------- +LineID LineType Length NumSegs +(-) (-) (m) (-) +4 bare 80.0 20 +4 buoy 100.0 25 +4 bare 420.0 42 +--------------------- END CONNECTIONS ---------------------------------------- +LineID End Stiffness EzX EzY EzZ TorsStiffness NxX NxY NxZ Pretwist +(-) (-) (N-m/rad) (-) (-) (-) (N-m/rad) (-) (-) (-) (deg) +4 A Rigid 0.1058 0 -0.9944 Free 0 1 0 0 +4 B Rigid 1 0 0 Free 0 1 0 0 +--------------------- OPTIONS ------------------------------------------------ +9.80665 g +1025.0 rhoW +3.0e6 kBot +3.0e5 cBot +0.35 frictionMu +200.0 WtrDpth +0.05 dtM +0.1 TMax +--------------------- OUTPUTS ------------------------------------------------ +FairTen4 +--------------------- END ---------------------------------------------------- +]=]) +file(WRITE "${WORK}/friction_coarse.dat" "${_friction_deck}") +string(REPLACE "4 bare 420.0 42" "4 bare 420.0 168" _friction_fine "${_friction_deck}") +file(WRITE "${WORK}/friction_fine.dat" "${_friction_fine}") +run_driver("friction coarse" "${WORK}/friction_coarse.dat" "${WORK}/friction_coarse" _fc_log) +require_text("friction coarse" "${_fc_log}" + "Note: line 4: seabed friction acts at the nodes, so a grounded element of 10.0 m locks or releases up to 552 N") +require_text("friction coarse" "${_fc_log}" "Grounded elements of at most 3.27 m resolve it.") +run_driver("friction fine" "${WORK}/friction_fine.dat" "${WORK}/friction_fine" _ff_log) +forbid_text("friction fine" "${_ff_log}" "seabed friction acts at the nodes") diff --git a/tests/check_touchdown_recovery.py b/tests/check_touchdown_recovery.py new file mode 100644 index 0000000..7d59371 --- /dev/null +++ b/tests/check_touchdown_recovery.py @@ -0,0 +1,121 @@ +# SPDX-License-Identifier: Apache-2.0 +"""Node tension at touchdown on a coarse mesh (driver_touchdown_recovery in CMakeLists.txt). + +A lazy-wave cable in 200 m of water is solved statically on a coarse mesh (10 m elements along +the grounded run) and on a reference mesh four times finer. The node tensions of the coarse +mesh, recovered from nodal equilibrium, must follow the reference within 3.5 % at every node +within 30 m of touchdown. The touchdown element is on the floor at one end only; spreading the +floor reaction of its grounded node over it understated the tension of the last suspended node +by 5.3 % (the rest of the profile is within 1.5 %). + +Usage: check_touchdown_recovery.py +""" + +import pathlib +import subprocess +import sys +import tempfile + +TOLERANCE = 0.035 +WINDOW = 30.0 +FLOOR_Z = -200.0 + +DECK = """--------------------- CableDyn Input File ------------------------------------ +Touchdown recovery check: lazy-wave cable, static +--------------------- LINE TYPES --------------------------------------------- +Name Diam Mass EA BA EI GAs GJ Irt Irn Cdn Cdt Can Cat +(-) (m) (kg/m) (N) (-) (Nm2) (N) (Nm2) (kgm) (kgm) (-) (-) (-) (-) +bare 0.16 36.70 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +buoy 0.29 59.53 4.69e8 0.0 1.99e4 1.0e8 1.5e4 1.0 1.0 1.2 0.1 1.0 0.0 +--------------------- POINTS ------------------------------------------------- +ID Type X Y Z +(-) (-) (m) (m) (m) +7 Coupled 30.0 0.0 -20.0 +8 Fixed 500.0 0.0 -200.0 +--------------------- LINES -------------------------------------------------- +ID NodeA NodeB Outputs +(-) (-) (-) (-) +4 7 8 - +--------------------- SECTIONS ----------------------------------------------- +LineID LineType Length NumSegs +(-) (-) (m) (-) +4 bare 80.0 {n1} +4 buoy 100.0 {n2} +4 bare 420.0 {n3} +--------------------- END CONNECTIONS ---------------------------------------- +LineID End Stiffness EzX EzY EzZ TorsStiffness NxX NxY NxZ Pretwist +(-) (-) (N-m/rad) (-) (-) (-) (N-m/rad) (-) (-) (-) (deg) +4 A Rigid 0.1058 0 -0.9944 Free 0 1 0 0 +4 B Rigid 1 0 0 Free 0 1 0 0 +--------------------- OPTIONS ------------------------------------------------ +9.80665 g +1025.0 rhoW +3.0e6 kBot +3.0e5 cBot +200.0 WtrDpth +0.025 dtM +0.0 TMax +--------------------- OUTPUTS ------------------------------------------------ +FairTen4 +--------------------- END ---------------------------------------------------- +""" + + +def static_profile(exe, work, name, segments): + """Arc length, z and node tension of every node of the static configuration.""" + deck = work / f"{name}.dat" + n1, n2, n3 = segments + deck.write_text(DECK.format(n1=n1, n2=n2, n3=n3), encoding="utf-8") + run = subprocess.run([exe, str(deck), str(work / name)], capture_output=True, text=True, check=False) + if run.returncode != 0: + raise SystemExit(f"FAIL: the {name} deck did not solve (exit {run.returncode}):\n{run.stdout}{run.stderr}") + rows = [] + for line in (work / f"{name}.static.out").read_text(encoding="utf-8").splitlines(): + fields = line.split() + try: + values = [float(x) for x in fields] + except ValueError: + continue + rows.append(values) + # columns: LineID Node ArcLength X Y Z Tension ... + return [(r[2], r[5], r[6]) for r in rows] + + +def interpolate(arcs, values, s): + for k in range(1, len(arcs)): + if arcs[k] >= s: + f = (s - arcs[k - 1]) / (arcs[k] - arcs[k - 1]) + return values[k - 1] + f * (values[k] - values[k - 1]) + return values[-1] + + +def main(): + if len(sys.argv) != 2: + raise SystemExit(__doc__) + exe = sys.argv[1] + with tempfile.TemporaryDirectory() as tmp: + work = pathlib.Path(tmp) + coarse = static_profile(exe, work, "coarse", (20, 25, 42)) + fine = static_profile(exe, work, "fine", (80, 100, 168)) + touchdown = next(s for s, z, _ in fine if z <= FLOOR_Z + 0.01) + arcs = [s for s, _, _ in fine] + tensions = [t for _, _, t in fine] + worst = 0.0 + failures = [] + for s, z, t in coarse: + if abs(s - touchdown) > WINDOW: + continue + ref = interpolate(arcs, tensions, s) + err = (t - ref) / ref + worst = max(worst, abs(err)) + print(f" s = {s:7.1f} m z = {z:8.2f} m T = {t:8.1f} N reference {ref:8.1f} N {100 * err:+6.2f} %") + if abs(err) > TOLERANCE: + failures.append(f"s = {s:.1f} m: {100 * err:+.2f} %") + print(f"touchdown at s = {touchdown:.1f} m; largest deviation {100 * worst:.2f} % (limit {100 * TOLERANCE:.1f} %)") + if failures: + raise SystemExit("FAIL: coarse-mesh node tension near touchdown off the reference: " + "; ".join(failures)) + print("PASS: coarse-mesh node tension follows the reference through touchdown") + + +if __name__ == "__main__": + main() diff --git a/tests/test_hermite_cable_static.f90 b/tests/test_hermite_cable_static.f90 index c1b6d44..0dd1b2a 100644 --- a/tests/test_hermite_cable_static.f90 +++ b/tests/test_hermite_cable_static.f90 @@ -18,7 +18,8 @@ PROGRAM test_hermite_cable_static USE CableDyn_HermiteCableDynamic, ONLY: CD_HermiteCableDynType, CD_HermiteCable_Dyn_Init, & CD_HermiteCable_Dyn_Set_Drag, CD_HermiteCable_Dyn_Reaction, & CD_HermiteCable_Dyn_End, CD_HCDYN_OK, CD_HermiteCable_Dyn_End_Force, & - CD_HermiteCable_Drag_Element, CD_HermiteCable_Dyn_Step + CD_HermiteCable_Drag_Element, CD_HermiteCable_Dyn_Step, & + CD_HermiteCable_Dyn_Node_Recovery USE CableDyn_SeabedContact, ONLY: CD_HermiteCable_Friction_Spring => CD_Seabed_Friction_Spring, & CD_Seabed_End_Element_Share, CD_Seabed_Normal_Law USE CableDyn_Line, ONLY: CD_Nodal_Seabed_Stiffness, CD_LINE_OK @@ -41,6 +42,7 @@ PROGRAM test_hermite_cable_static CALL check_seabed_mesh_jump() CALL check_scalar_contact_on_slope() CALL check_taut_line_anchored_on_floor() + CALL check_laid_line_node_recovery() CALL check_end_support_energy_gradient() CALL check_frictionless_slope_run() CALL check_long_sequence_hierarchy() @@ -385,6 +387,55 @@ SUBROUTINE check_taut_line_anchored_on_floor() END DO END SUBROUTINE check_taut_line_anchored_on_floor + SUBROUTINE check_laid_line_node_recovery() + !! A straight line laid on the scalar penalty floor of CD_HermiteCable_Dyn_Init (no contact + !! model), both ends held exactly on the floor: the floor carries the weight, so the + !! recovered bend moment is only that of the line settling into the floor, at the held + !! ends too. An end node on the + !! floor has no contact force of its own; its laid end element is still grounded (the + !! held-end floor support), so the reaction of the next node stays shared with it. Counted + !! as suspended, the end element lost that share and the anchor and its neighbour showed a + !! weight moment of order w l0^2/12 (521 N.m on this line, 31 %; the bending of the line + !! settling 4 mm into the penalty floor next to its held ends is 4.6 N.m, 0.3 %). + INTEGER, PARAMETER :: NE = 20, NN = NE + 1 + REAL(wp) :: l0(NE), EA(NE), EI(NE), w(NE), rho_a(NE), seed(6*NN), q(6*NN), curv(NN), res, kn + REAL(wp) :: trib(NN), tension(NN), bend(NN), scale + INTEGER :: fixed(6 + 2*NN), nfix, i, es, it + TYPE(CD_HermiteCableDynType) :: dyn + CHARACTER(512) :: em + l0 = 5.0_wp; EA = 5.0e8_wp; EI = 2.0e4_wp; w = 800.0_wp; rho_a = 120.0_wp; kn = 2.0e5_wp + trib = 0.0_wp + DO i = 1, NE + trib(i) = trib(i) + 0.5_wp*l0(i); trib(i + 1) = trib(i + 1) + 0.5_wp*l0(i) + END DO + ! both ends on the 100 m floor, 100 m apart: a 100 m line laid straight + CALL CD_HermiteCable_Trivial_Seed([0.0_wp, 0.0_wp, -100.0_wp], [100.0_wp, 0.0_wp, -100.0_wp], l0, -100.0_wp, & + seed, es, em) + CALL require(es == CD_HCSTAT_OK, 'laid line seed: '//TRIM(em)) + fixed(1:6) = [1, 2, 3, 6*NE + 1, 6*NE + 2, 6*NE + 3] + nfix = 6 + DO i = 1, NN + fixed(nfix + 1:nfix + 2) = [6*i - 4, 6*i - 1] + nfix = nfix + 2 + END DO + CALL CD_HermiteCable_Static_Solve(l0, EA, EI, w, seed, fixed(1:nfix), -100.0_wp, 0.0_wp, 4, 200, 1.0e-10_wp, & + 1.0_wp, q, curv, res, it, es, em, contact_kn=kn*trib) + CALL require(es == CD_HCSTAT_OK, 'laid line on the floor: '//TRIM(em)) + IF (es /= CD_HCSTAT_OK) RETURN + CALL CD_HermiteCable_Dyn_Init(dyn, l0, EA, EI, rho_a, w, q, fixed(1:nfix), -100.0_wp, kn, 0.8_wp, es, em) + CALL require(es == CD_HCDYN_OK, 'laid line dynamics: '//TRIM(em)) + IF (es /= CD_HCDYN_OK) RETURN + CALL CD_HermiteCable_Dyn_Node_Recovery(dyn, es, em, tension=tension, bend_moment=bend) + CALL require(es == CD_HCDYN_OK, 'laid line node recovery: '//TRIM(em)) + scale = w(1)*l0(1)**2/12.0_wp + WRITE (*, '(A,ES10.3,A,I0,A,ES10.3,A)') ' laid line: max |bend moment| ', MAXVAL(ABS(bend)), & + ' N.m at node ', MAXLOC(ABS(bend)), & + ' (w l0^2/12 = ', scale, ' N.m)' + CALL require(MAXVAL(ABS(bend)) <= 0.02_wp*scale, & + 'laid line: no weight moment at the held ends or along the line') + CALL CD_HermiteCable_Dyn_End(dyn) + END SUBROUTINE check_laid_line_node_recovery + SUBROUTINE check_end_support_energy_gradient() !! The energy line search (energy_minimization) tests the potential against R.dq, so the !! residual must be the gradient of the potential it evaluates. With the end element's