diff --git a/docs/source/user/api_change.rst b/docs/source/user/api_change.rst index 0fdca046fd..2431f9a75c 100644 --- a/docs/source/user/api_change.rst +++ b/docs/source/user/api_change.rst @@ -18,6 +18,8 @@ The generalized support-structure (GS) influence model was added to AeroDyn. Thi OLAF now selects the regularization function for the vortex particles separately from the vortex segments. The new ``RegFunctionPart`` input on line 26 of the OLAF input file controls the particle kernel, and ``RegFunction`` now applies only to the segments. Setting ``RegFunctionPart=1`` reproduces the exponential kernel previously used with ``RegFunction`` greater than 0. Setting ``RegFunctionPart=0`` restores an unregularized particle kernel previously used with ``RegFunction=0``. The newly added default compact particle regularization (``RegFunctionPart=2``) provides significant speed up compared to the exponential kernel (``RegFunctionPart=1``). +MoorDyn now includes slip supports that only constrain motion in the line normal direction, like a hoop the line goes through, and apply only friction in the line tangent direction. Slip supports can be added in the input file under an optional new section called ``Slip Supports``. See :ref:`MD-SlipSupport` for details. + ============================================= ======== ==================== ========================================================================================================================================================================================================================================================================================================== Added in OpenFAST `5.1.0` -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- @@ -39,6 +41,9 @@ AeroDyn \* NumGSMembers 0 AeroDyn \* GSMemberID GSMJointID1 GSMJointID2 GSMDia1 GSMDia2 GSMCd1 GSMCd2 GSMTI1 GSMTI2 GSMDiv AeroDyn \* (-) (-) (-) (m) (m) (-) (-) (-) (-) (m) OLAF 26 RegFunctionPart default RegFunctionPart - Particle regularization function {0: None, 1: Exponential, 2: Compact, "default": 2} [only if VelocityMethod=2,3] (switch) +MoorDyn \* ---------------------- SLIP SUPPORTS ------------------------------------------------ +MoorDyn \* ID Point Line Width K Zeta MuS MuK VTol VS StickDisp +MoorDyn \* (-) (-) (-) (h) (N/m) (-) (-) (-) (m/s) (m/s) (m) ============================================= ======== ==================== ========================================================================================================================================================================================================================================================================================================== The meaning of the HydroDyn ``WaveDisp`` input has changed, although its valid values (0 and 1) are unchanged. diff --git a/docs/source/user/moordyn/index.rst b/docs/source/user/moordyn/index.rst index 9d5d60d15e..89aa53cb3e 100644 --- a/docs/source/user/moordyn/index.rst +++ b/docs/source/user/moordyn/index.rst @@ -15,3 +15,16 @@ and links to publications with the relevant theory. Examples of how to use MoorDynF and MoorDynC can be found here: `MoorDyn Example Uses `_ + + +MoorDyn-F Specific Features +--------------------------- + +The pages below document capabilities that are currently specific to MoorDyn-F, the Fortran version of +MoorDyn in OpenFAST, and are not yet covered by the external MoorDyn documentation. + +.. toctree:: + :maxdepth: 2 + + moordynf_features.rst + slip_support_theory.rst diff --git a/docs/source/user/moordyn/moordynf_features.rst b/docs/source/user/moordyn/moordynf_features.rst new file mode 100644 index 0000000000..ee8b3f9132 --- /dev/null +++ b/docs/source/user/moordyn/moordynf_features.rst @@ -0,0 +1,232 @@ +.. _MD-F-Features: + +MoorDyn-F Specific Features +=========================== + +This page describes input options that are available in MoorDyn-F (the Fortran implementation of MoorDyn +used by OpenFAST) in addition to those documented in the +`MoorDyn input file documentation `_. +All other sections of the MoorDyn input file (line types, points, lines, failures, options, and outputs) are +unchanged and follow the external documentation. + + +.. _MD-SlipSupport: + +Slip Supports +------------- + +A slip support models a line threaded through a hoop, eyelet, or guide that is carried by a MoorDyn Point. +The hoop constrains the line transversely but lets it slide along its length, resisted by friction. +The hoop has no physical size: there is no clearance between the line and the hoop, and no bend radius. +Typical applications include lines running through fairlead guides, mid-line guides on a structure, or a +free weight or buoy riding on a line. + +The line keeps its original discretization while it slides: no nodes are added, removed, or moved along +the line. Instead, the hoop acts on a virtual contact point that moves continuously along the line, and +the contact force is spread to the nearby line nodes with a smooth weighting kernel. As a result, the +contact force varies smoothly as the line nodes slide past. The formulation is given in +:ref:`MD-SlipSupport-Theory`. + +Main capabilities: + +- The hoop can be any Point type: fixed, coupled, free, or attached to a Body. The reaction force from the + line is added to the Point, so it is included in coupled loads passed to the glue code and in the + dynamics of free Points and Bodies. +- Several slip supports may act on the same line, and several lines may pass through the same Point. +- Friction is either a regularized Coulomb model or a LuGre model with true sticking. Both use separate + static and kinetic friction coefficients with a smooth (Stribeck) transition between them. +- If the line slides completely through a hoop, the support releases the line permanently and the + simulation continues. +- Linearization is supported (see :ref:`MD-SlipSupport-Linearization`). + + +.. _MD-SlipSupport-Input: + +Input +~~~~~ + +Slip supports are defined in an optional ``SLIP SUPPORTS`` section of the MoorDyn input file. The section +header is a line containing ``---`` and the text ``SLIP SUPPORTS``, followed by a column-name line, a units +line, and one row per slip support. **The section must come after the** ``LINES`` **section**, because it +refers to Line properties. + +Example (one hoop on Point 2 threading Line 1, with LuGre friction): + +.. code-block:: none + + ---------------------- SLIP SUPPORTS ------------------------------------------------ + ID Point Line Width K Zeta MuS MuK VTol VS StickDisp + (#) (#) (#) (h) (N/m) (-) (-) (-) (m/s) (m/s) (m) + 1 2 1 2 0 0.8 0.6 0.3 0.01 0.1 0.005 + +.. list-table:: Slip support input columns + :header-rows: 1 + :widths: 12 10 78 + + * - Column + - Units + - Description + * - ``ID`` + - (-) + - Slip support number. Must be sequential starting from 1. + * - ``Point`` + - (-) + - ID of the Point carrying the hoop. Any Point type may be used. A free Point used as a hoop + must have a nonzero mass, added mass, or attached lines (see :ref:`MD-SlipSupport-Usage`). + * - ``Line`` + - (-) + - ID of the Line threaded through the hoop. + * - ``Width`` + - (-) + - Contact kernel width, given as a multiple of the segment length :math:`h`. Here :math:`h` is the + longest unstretched segment length among all lines that pass through the same Point, so all lines + through one hoop share the same physical contact width :math:`W = \mathrm{Width}\cdot h`. + Must be at least 1. Values below 2 produce a warning because the contact force may ripple as the + line slides past nodes. **2 is recommended.** :math:`4W` must also be shorter than the + unstretched length of the line. + * - ``K`` + - (N/m) + - Transverse contact stiffness (penalty stiffness) that keeps the line centered in the hoop. Enter 0 + for an automatic value tied to the line node mass and the time step ``dtM`` (recommended). + * - ``Zeta`` + - (-) + - Transverse contact damping ratio, relative to critical damping of the contact stiffness. A value of + about 0.8 is recommended. + * - ``MuS`` + - (-) + - Static friction coefficient. Must satisfy ``MuS`` :math:`\ge` ``MuK``. Enter 0 together with + ``MuK`` = 0 for a frictionless hoop. + * - ``MuK`` + - (-) + - Kinetic (sliding) friction coefficient, :math:`\ge 0`. + * - ``VTol`` + - (m/s) + - Regularization velocity of the Coulomb friction model. Must be positive when ``MuS`` > 0 and + ``StickDisp`` = 0. Ignored when LuGre friction is used. The line creeps through the hoop at a speed of + order ``VTol`` when it should be stuck, so it should be small compared with the sliding speeds of + interest. MoorDyn raises it internally when needed for numerical stability + (see :ref:`MD-SlipSupport-Coulomb`). + * - ``VS`` + - (m/s) + - Stribeck velocity: the slip speed over which friction drops from the static to the kinetic level. + Must be positive when ``MuS`` > ``MuK`` or when LuGre friction is used. + * - ``StickDisp`` + - (m) + - Pre-sliding displacement of the LuGre friction model: the elastic slip the contact can take before + it breaks loose, typically a few millimeters. Enter 0 to use the regularized Coulomb model instead. + +Friction model selection: + +- ``StickDisp`` = 0: regularized Coulomb friction. No additional states. The line can creep slowly through + the hoop while it should be stuck. +- ``StickDisp`` > 0: LuGre friction. One additional continuous state per slip support. The line truly + sticks until the static friction limit is reached. If ``StickDisp`` is too small for the time step, + MoorDyn stiffens the contact less than requested to stay stable and prints a one-time warning. Either + increase ``StickDisp`` or reduce the time step to remove the warning. + + +.. _MD-SlipSupport-Outputs: + +Output Channels +~~~~~~~~~~~~~~~ + +Slip support output channels are requested in the ``OUTPUTS`` section with the prefix ``SLIP`` followed by +the slip support ID and a quantity code, for example ``SLIP1FN``. + +.. list-table:: Slip support output channels + :header-rows: 1 + :widths: 18 10 72 + + * - Channel + - Units + - Description + * - ``SLIPS`` + - (m) + - Contact location: unstretched arc length from line end A to the contact point. + * - ``SLIPVS`` + - (m/s) + - Slip velocity of the line relative to the hoop, along the line. Positive when the line moves toward + end B (``SLIPS`` decreasing). This is a physical velocity of the stretched line, not the rate of + change of ``SLIPS``, which is measured in unstretched arc length: + d(``SLIPS``)/dt :math:`\approx -` ``SLIPVS`` :math:`/(1+\varepsilon)`, with :math:`\varepsilon` + the line strain at the contact point. + * - ``SLIPFN`` + - (N) + - Magnitude of the normal (transverse) contact force applied to the line. Always :math:`\ge 0`. + * - ``SLIPFF`` + - (N) + - Friction force applied to the line, along the line. Positive pushes the line toward end B. It + normally has the opposite sign of ``SLIPVS``. + * - ``SLIPFX``, ``SLIPFY``, ``SLIPFZ`` + - (N) + - Total force exerted by the line on the hoop in global coordinates. This is the reaction added to the + Point. + * - ``SLIPZ`` + - (m) + - LuGre bristle deflection. Zero when Coulomb friction is used. + +After a slip support is released, ``SLIPVS``, ``SLIPFN``, ``SLIPFF``, and ``SLIPFX/FY/FZ`` are +zero, while ``SLIPS`` and ``SLIPZ`` keep their last values. + + +.. _MD-SlipSupport-Usage: + +Usage Notes +~~~~~~~~~~~ + +**Initial line shape.** A line with slip supports is not initialized with the catenary solution. Its +internal nodes are placed at equal distances along a polyline that runs from end A through each slip +support Point, in the order the slip supports of that line are listed, to end B. The initial contact +location therefore follows the ratio of the straight-line distances between the line ends and the hoops. +A hoop does not need to coincide with a node. List the slip supports of a line in order from end A to +end B. + +**Initial condition generation.** During dynamic relaxation (``TmaxIC`` > 0), friction is replaced by a +viscous slip damper so that the line can slide to its frictionless equilibrium. The damping coefficient is +set automatically from the line properties (see :ref:`MD-SlipSupport-IC`). LuGre states are held at +zero and reset after dynamic relaxation. Note that the convergence check of the dynamic relaxation only +monitors fairlead tensions, so a line may still be sliding slowly when it ends; friction then holds the +remaining imbalance. Set ``TmaxIC`` = 0 to start a free hoop (for example, a bead on a line) from its +input position. + +**Release.** A slip support releases its line permanently once the hoop is more than 0.1 segment +length beyond either line end throughout a whole time step. A message is written to the screen and log +file. Supports are not released during dynamic relaxation. + +**Choosing Width.** Larger values give smoother contact forces when the line nodes slide past the support, +but spread the load over a longer length and place the line slightly outside the hoop at bends. +``Width`` = 2 removes node-passing ripple for practical purposes. Because ``Width`` is relative to the +segment length, mesh refinement also narrows the contact and shrinks the offset between the line and the +hoop; to keep a fixed physical contact width while refining, scale ``Width`` with the number of +segments. + +**Bend angle.** The model is intended for moderate bend angles at the hoop. It can work up to a hairpin +turn, but likely with degraded accuracy. To model a line wrapping around a large-radius sheave or bollard, +use several slip supports distributed around its circumference. + +**Free hoop points.** The automatic stiffness ``K`` = 0 is sized from the line node mass. A free hoop Point +that is much lighter than a line node may need an explicit, smaller ``K`` for numerical stability. + +**Active tensioning.** Changes of the last segment length by line control channels are not considered when +checking ``Width``. + + +.. _MD-SlipSupport-Linearization: + +Linearization +~~~~~~~~~~~~~ + +Slip supports may be used with linearization. Each LuGre slip support adds one first-order state (the +bristle deflection) named ``Slip bristle z, m``. The linearized model reflects the numerical model of +the contact: + +- With Coulomb friction, near-zero slip the friction acts as a viscous damper along the line. Because + MoorDyn perturbs velocities by 0.1 m/s, the linearized damping is a secant over that perturbation. It + approaches the tangent value :math:`\mu_s F_N / \tilde v_{tol}` only if the adjusted regularization velocity + :math:`\tilde v_{tol} \ge` ``VTol`` (see :ref:`MD-SlipSupport-Coulomb`) is well above 0.1 m/s; otherwise + the damping is smaller. +- The transverse contact stiffness adds a high-frequency transverse mode. With the automatic stiffness, its + frequency is about :math:`0.3/dt_M` rad/s. +- With LuGre friction, the stuck contact appears as a stiff spring along the line (the bristle stiffness) + with an associated fast first-order mode. The bristle damping is subject to the same 0.1 m/s velocity + perturbation and is underestimated unless ``VS`` is well above 0.1 m/s. diff --git a/docs/source/user/moordyn/slip_support_theory.rst b/docs/source/user/moordyn/slip_support_theory.rst new file mode 100644 index 0000000000..3bfe26b7e0 --- /dev/null +++ b/docs/source/user/moordyn/slip_support_theory.rst @@ -0,0 +1,480 @@ +.. _MD-SlipSupport-Theory: + +Slip Support Theory +=================== + +This page describes the formulation of the slip supports in MoorDyn-F. Input and usage are described in +:ref:`MD-SlipSupport`. + + +Notation +-------- + +A line has nodes :math:`j = 0,\dots,N` with positions :math:`\mathbf{r}_j` and velocities +:math:`\dot{\mathbf{r}}_j`. Segment :math:`k = 1,\dots,N` connects nodes :math:`k-1` and :math:`k` and has +*unstretched* length :math:`l_k`. The unstretched arc-length coordinate of node :math:`j` is + +.. math:: + + s_0 = 0, \qquad s_j = s_{j-1} + l_j, \qquad L = s_N . + +The hoop (slip support) is located at the position :math:`\mathbf{x}_h` of its Point and moves with velocity +:math:`\dot{\mathbf{x}}_h`. The contact location on the line is given by the *unstretched* arc length +:math:`s_c` measured from end A. All directions along the line are oriented from end A toward end B. +:math:`\Delta t` denotes the requested mooring time step ``dtM``. + + +Contact Kernel +-------------- + +The hoop acts on a virtual contact point that is a weighted average of nearby line nodes. The weights come +from a cubic B-spline kernel of fixed width in *unstretched* arc length, + +.. math:: + + B(t) = + \begin{cases} + \left(4 - 6t^2 + 3|t|^3\right)/6, & |t| < 1 \\ + \left(2 - |t|\right)^3/6, & 1 \le |t| < 2 \\ + 0, & |t| \ge 2 + \end{cases} + +with kernel scale :math:`W = \mathrm{Width}\cdot h`, where :math:`h` is the longest segment length among +all lines through the same hoop Point. The unnormalized weight of node :math:`j` is + +.. math:: + + a_j(s_c) = B\!\left(\frac{s_j - s_c}{W}\right)\Delta s_j , + +where :math:`\Delta s_j = (l_j + l_{j+1})/2` is the line length represented by node :math:`j`, with +:math:`\Delta s_0 = l_1` and :math:`\Delta s_N = l_N`. The end nodes use full element lengths instead of +half element lengths due to the mirroring algorithm explained below. + +Near the line ends, the kernel support :math:`[s_c - 2W,\ s_c + 2W]` extends beyond the line. The line is +then extended by point reflection about the end node: every real node :math:`j` whose mirror image falls +inside the support contributes a virtual node, at arc length :math:`-s_j` with position +:math:`2\mathbf{r}_0 - \mathbf{r}_j` beyond end A, or at :math:`2L - s_{N-j}` with position +:math:`2\mathbf{r}_N - \mathbf{r}_{N-j}` beyond end B. Depending on :math:`W/h` and :math:`s_c`, several +virtual nodes may be involved (up to about :math:`2W/h` when the contact is at the end). Because a virtual +position is a linear combination of two real nodes, the weight :math:`a` of a virtual node is folded back +onto the real nodes, as :math:`+2a` on the end node and :math:`-a` on the mirrored node, which leaves the +weight sum unchanged. Because the end node becomes an interior node of the extended line, it represents a +full segment length. The reflection keeps a straight line straight across the end and makes the contact +curve pass exactly through the end node when :math:`s_c` is at the end, instead of being pulled inward as a +truncated kernel would be. The normalized weights are + +.. math:: + + w_j(s_c) = \frac{\tilde a_j(s_c)}{A(s_c)}, \qquad A(s_c) = \sum_j \tilde a_j(s_c) , + +where :math:`\tilde a_j` includes the folded contributions. The contact point, its material velocity, and +the unit tangent are + +.. math:: + + \mathbf{C}(s_c) = \sum_j w_j \mathbf{r}_j, \qquad + \mathbf{V}(s_c) = \sum_j w_j \dot{\mathbf{r}}_j, \qquad + \hat{\mathbf{t}} = \frac{\mathbf{C}'}{\lVert\mathbf{C}'\rVert}, \quad + \mathbf{C}' = \sum_j \frac{\mathrm{d}w_j}{\mathrm{d}s_c}\mathbf{r}_j , + +with :math:`\mathrm{d}w_j/\mathrm{d}s_c = (\tilde a_j' - w_j A')/A`. + +For fixed node positions, :math:`\mathbf{C}(s)` is a smooth curve following the line nodes parameterized by +the *unstretched* arc length :math:`s`, and :math:`\mathbf{C}'` is its derivative with respect to :math:`s` +at the current node positions. Only the weights depend on :math:`s`; the node positions are held fixed. +:math:`\hat{\mathbf{t}}` is therefore the tangent of the smoothed curve at the contact point rather than +the direction of any single segment, and it varies continuously as the contact moves past nodes. The +weight derivatives follow from the kernel, +:math:`a_j' = -B'\!\left((s_j - s_c)/W\right)\Delta s_j/W`, with the same folding at the line ends. If +:math:`\lVert\mathbf{C}'\rVert` vanishes, the direction of the contact segment +:math:`\mathbf{r}_k - \mathbf{r}_{k-1}` is used instead. + +The same weights :math:`w_j` are used both to interpolate the contact kinematics :math:`\mathbf{C}` and +:math:`\mathbf{V}` from the nodes and to spread the contact force :math:`\mathbf{F}` back onto the nodes as +:math:`w_j\mathbf{F}` (see :ref:`MD-SlipSupport-Distribution`). With this choice, the power delivered to the +line is + +.. math:: + + \sum_j w_j\mathbf{F}\cdot\dot{\mathbf{r}}_j = \mathbf{F}\cdot\mathbf{V}, + +so a spring and damper acting on :math:`\mathbf{C}` and :math:`\mathbf{V}` act on the line nodes as a passive +spring and damper: the spring stores exactly its potential energy and the damper only dissipates energy. The +same choice gives the line a single well-defined inertia :math:`m_C` at the contact point (see +:ref:`MD-SlipSupport-Normal`), on which the stiffness, damping, and stability bounds below are based. If the +kinematics were interpolated with weights different from those used to spread the force, the damper could +feed energy into some node motions. The per-node tangent split described in +:ref:`MD-SlipSupport-Distribution` is applied only to the spring and breaks this property for it; the damper +is spread with the weights :math:`w_j` without projection, and the friction velocity is measured with the same +per-node tangents used to apply the friction force, so both remain dissipative. Interpolation and spreading with the same +kernel follow the same principle as in the immersed boundary method (Peskin, 2002). + +Properties of the kernel: + +- With the reflected virtual nodes, :math:`\mathbf{C}(0) = \mathbf{r}_0` and :math:`\mathbf{C}(L) = \mathbf{r}_N`. +- The kernel has an RMS spread of :math:`W/\sqrt{3}` along the line. +- :math:`\mathbf{C}` is a smoothed version of the line, so the line passes slightly outside the contact point + at :math:`\mathbf{C}(s_c)` at a moderate bend, with increased offset for a sharper bend and larger kernel width. + This is one reason why solution accuracy degrades for sharply bent lines. +- On a line with uniform segments away from the ends, the unnormalized weight sum :math:`A` and the mean and + variance of the weight distribution do not depend on the contact position when :math:`W/h` is an integer. + For moderate bends, the node-passing ripple of the contact force is governed mainly by the variation of + :math:`\sum_j w_j^2` as the kernel shifts along the evenly spaced nodes. This variation repeats every segment, + and its relative range over one segment decreases rapidly with :math:`W/h`: + + .. list-table:: + :header-rows: 1 + + * - :math:`W/h` + - 1 + - 1.5 + - 2 + - 3 + - 6 + * - Variation of :math:`\sum_j w_j^2` (max - min)/mean + - 8.3% + - 0.41% + - 0.27% + - 0.043% + - 0.0024% + + +Contact Location +---------------- + +The contact location is defined as the closest point on the curve :math:`\mathbf{C}(s)` to the hoop position +:math:`\mathbf{x}_h`: + +.. math:: + + s_c = \arg\min_{s}\ \lVert \mathbf{C}(s) - \mathbf{x}_h \rVert^2 . + +It is parameterized by segment :math:`k` and a local coordinate :math:`u \in [0,1]` as +:math:`s_c = s_{k-1} + u\,l_k`. For each candidate segment, :math:`u` is initialized by projecting +:math:`\mathbf{x}_h` onto the chord :math:`\mathbf{r}_{k-1}\mathbf{r}_k` and refined by Gauss-Newton +iterations + +.. math:: + + u \leftarrow \operatorname{clamp}_{[0,1]}\!\left(u - \frac{(\mathbf{C} - \mathbf{x}_h)\cdot\mathbf{C}_u} + {\mathbf{C}_u\cdot\mathbf{C}_u}\right), \qquad \mathbf{C}_u = l_k\,\mathbf{C}' . + +The candidate segments are the segment found in the previous evaluation and two segments on either side +of it (all segments on the first evaluation). The candidate with the smallest distance is used. The contact +location is not a state; it is recomputed from the line and hoop positions in every evaluation. + +Searching only neighboring segments keeps the update efficient without checking every segment of the line. +It also prevents the contact location from jumping between evaluations to a distant part of the line that +happens to pass near the hoop, such as the other leg of a sharp bend. + +If the contact moves by more than two segments between evaluations, the search stops at the edge of the +window, so the contact lags behind the hoop and the contact force is applied at the wrong location. The +window follows the contact by up to two segments per evaluation, so the lag clears after a few evaluations. +This is unlikely to be an issue in practice: with an explicit time step limited by the axial wave speed +:math:`c` of the line (:math:`\Delta t \lesssim h/c`), it would require slip speeds comparable to :math:`c`. + + +.. _MD-SlipSupport-Normal: + +Normal Contact Force +-------------------- + +The contact offset perpendicular to the line and the transverse relative velocity are + +.. math:: + + \boldsymbol{\delta} = (\mathbf{C}(s_c) - \mathbf{x}_h) - \left[(\mathbf{C}(s_c) - \mathbf{x}_h)\cdot\hat{\mathbf{t}}(s_c)\right]\hat{\mathbf{t}}(s_c), + \qquad + \mathbf{v}_\perp = (\mathbf{V}(s_c) - \dot{\mathbf{x}}_h) - \left[(\mathbf{V}(s_c) - \dot{\mathbf{x}}_h)\cdot\hat{\mathbf{t}}(s_c)\right]\hat{\mathbf{t}}(s_c). + +The transverse penalty force acting on the contact point is + +.. math:: + + \mathbf{F}_\perp = \mathbf{F}_K + \mathbf{F}_c, \qquad \mathbf{F}_K = -K\boldsymbol{\delta}, \qquad + \mathbf{F}_c = -c\,\mathbf{v}_\perp, \qquad c = 2\zeta\sqrt{K m_C}, + +where :math:`m_C` is the inertia felt at the contact point when a force is shared among the nodes, + +.. math:: + + m_C = \frac{m_{node}}{\sum_j w_j^2}, \qquad + m_{node} = \frac{L}{N}\left(m' + \rho_w C_{an} A_d\right), + +with line mass per unit length :math:`m'`, water density :math:`\rho_w`, transverse added-mass +coefficient :math:`C_{an}`, and :math:`A_d = \pi d^2/4`. + +The expression for :math:`m_C` follows from treating the nodes as independent masses on the short time scale +of the contact. A force :math:`\mathbf{F}` at the contact point is applied to node :math:`j` as +:math:`w_j\mathbf{F}`, which accelerates the node by :math:`w_j\mathbf{F}/m_{node}`. The contact point moves +with the weighted node motion, so its acceleration is +:math:`\sum_j w_j\,(w_j\mathbf{F}/m_{node}) = \mathbf{F}\sum_j w_j^2/m_{node}`, which gives the effective mass +:math:`m_C`. Here, we assumed :math:`m_{node}` is constant along each line, including at the ends, when +estimating :math:`m_C`. Because the weights are nonnegative and sum to one, :math:`\sum_j w_j^2 \le 1`, and +therefore + +.. math:: + + m_C \ge m_{node} . + +Equality holds when the whole weight sits on a single node, such as at a line end. Spreading the force over +more nodes makes the contact point heavier; for a uniform line, +:math:`\sum_j w_j^2 \approx 0.48\,h/W`, so :math:`m_C \approx 4\,m_{node}` for :math:`W/h = 2`. Near the line +ends, the folded contribution :math:`-a` partly cancels the weight of the mirrored node, but the net weight +remains nonnegative because the virtual node is always farther from the contact than its real counterpart. +The bound :math:`m_C \ge m_{node}` therefore also holds near the line ends. + +:math:`m_C` is recomputed from the current weights in every evaluation, so the damping coefficient :math:`c` +and the friction regularization and bounds that depend on :math:`m_C` (see below) follow the contact as it +slides. :math:`K`, in contrast, is fixed at initialization. When ``K`` = 0 is input, the following default +is used: + +.. math:: + + K = m_{node}\left(\frac{0.3}{\Delta t}\right)^2 , + +which keeps the contact frequency :math:`\sqrt{K/m_C}` at or below :math:`0.3/\Delta t`, including at +the line ends where the weight concentrates on the end node. + + +.. _MD-SlipSupport-Distribution: + +Force Distribution to Nodes +--------------------------- + +Each node receives its share of the contact force. The spring and friction parts are split along the node's +own tangent; the damper part is not. The node tangent +:math:`\hat{\mathbf{q}}_j` is the unit bisector of the adjacent segment directions, + +.. math:: + + \hat{\mathbf{q}}_j = \frac{\hat{\mathbf{e}}_j + \hat{\mathbf{e}}_{j+1}}{\lVert\hat{\mathbf{e}}_j + \hat{\mathbf{e}}_{j+1}\rVert}, + \qquad \hat{\mathbf{e}}_k = \frac{\mathbf{r}_k - \mathbf{r}_{k-1}}{\lVert\mathbf{r}_k - \mathbf{r}_{k-1}\rVert}, + +using the single adjacent segment at the end nodes (and :math:`\hat{\mathbf{t}}` if the bisector vanishes). +Near the line ends, the share :math:`+2a` that a virtual node folds onto the end node keeps the tangent of its +mirrored node in these per-node splits, because point reflection preserves the tangent; only the end node's +own share uses :math:`\hat{\mathbf{q}}_0` or :math:`\hat{\mathbf{q}}_N`. +The normal force applied to node :math:`j` is the part of its spring share perpendicular to +:math:`\hat{\mathbf{q}}_j` plus its full damper share, + +.. math:: + + \mathbf{f}^N_j = w_j\left[\mathbf{F}_K - (\mathbf{F}_K\cdot\hat{\mathbf{q}}_j)\,\hat{\mathbf{q}}_j\right] + w_j\mathbf{F}_c, + \qquad + F_N = \left\lVert \sum_j \mathbf{f}^N_j \right\rVert . + +:math:`F_N` is the normal force actually delivered to the line and is used for the friction capacity. With +the friction force :math:`F_f` defined below, the total force on node :math:`j` and the reaction on the hoop +are + +.. math:: + + \mathbf{f}_j = \mathbf{f}^N_j + w_j F_f\,\hat{\mathbf{q}}_j, \qquad + \mathbf{F}_h = -\sum_j \mathbf{f}_j . + +The node forces :math:`\mathbf{f}_j` are added to the net forces of the line nodes; forces on end nodes are +passed to the attached objects as usual. :math:`\mathbf{F}_h` is added to the net force of the hoop Point. + +Splitting the spring force along each node's own tangent ensures that a frictionless hoop at rest applies no +force along the line at any node. Without it, the parts of the kernel acting on nodes away from a sharp bend +would push the line along its length, and the hoop could hold tension without friction. The split makes the +spring force non-conservative (it is not the gradient of a contact potential), but the contact damping likely +far outweighs the resulting energy change. The damper is spread with the same weights used to interpolate +:math:`\mathbf{V}`, so its power delivered to the line and hoop together is +:math:`\mathbf{F}_c\cdot\mathbf{v}_\perp = -c\lVert\mathbf{v}_\perp\rVert^2 \le 0`. Its tangential +components at a bend vanish when the contact is at rest, so they cannot hold tension. Projecting the damper +per node as well would keep it dissipative, but at a sharp bend the projections of both the velocity and the +force would each reduce the damping by about :math:`\cos^2` of the angle between :math:`\hat{\mathbf{q}}_j` and +:math:`\hat{\mathbf{t}}`, leaving the contact underdamped. Friction is dissipative because the slip velocity +:math:`v` below is measured with the same tangents used to apply the friction force, so its power is +:math:`F_f v`, which is :math:`\le 0` for Coulomb friction. + +The friction force is not projected: each node receives its full share :math:`w_j F_f` along its own tangent. +Because the weights sum to one, the tension changes by exactly :math:`F_f` across the support. The projection of the spring force, in +contrast, discards part of each node's share at a bend, so the effective contact stiffness drops +below :math:`K` at sharp bends. The delivered normal force :math:`F_N` is unaffected at +equilibrium, because the penalty offset grows until the line tension is balanced. This loss of stiffness is +the price for avoiding spurious tangential forces. + + +Slip Velocity and Friction Coefficient +-------------------------------------- + +The slip velocity of the line relative to the hoop is + +.. math:: + + v = \sum_j w_j\,\hat{\mathbf{q}}_j\cdot(\dot{\mathbf{r}}_j - \dot{\mathbf{x}}_h) , + +positive when the line material moves toward end B. On a straight line, this is +:math:`(\mathbf{V}(s_c) - \dot{\mathbf{x}}_h)\cdot\hat{\mathbf{t}}`. Because :math:`s_c` is a material coordinate, +:math:`\dot s_c \approx -v/(1+\varepsilon)` with :math:`\varepsilon` the local line strain. + +Friction uses the static friction coefficient :math:`\mu_s` (input ``MuS``), which sets the maximum force +before the line starts sliding, and the kinetic friction coefficient :math:`\mu_k` (input ``MuK``), which +applies during steady sliding. The effective coefficient :math:`\mu(v)` transitions smoothly from +:math:`\mu_s` at zero slip to :math:`\mu_k` at high slip speed over the Stribeck velocity :math:`v_S` +(input ``VS``), + +.. math:: + + \mu(v) = \mu_k + (\mu_s - \mu_k)\,e^{-(v/v_S)^2} . + + +.. _MD-SlipSupport-Coulomb: + +Regularized Coulomb Friction +---------------------------- + +When ``StickDisp`` = 0, + +.. math:: + + F_f = -\mu(v)\,F_N \tanh\!\left(\frac{v}{\tilde v_{tol}}\right), + \qquad + \tilde v_{tol} = \max\!\left(v_{tol},\ \frac{\mu_s F_N \Delta t}{m_C}\right). + +Near zero slip, the friction acts as a viscous damper with coefficient :math:`\mu_s F_N/\tilde v_{tol}`. The +lower bound on :math:`\tilde v_{tol}`, evaluated with the current :math:`F_N` and :math:`m_C`, limits this +damping to :math:`(\mu_s F_N/\tilde v_{tol})\,\Delta t/m_C \le 1`, which is half of the explicit stability +limit of a viscous damper, :math:`c\,\Delta t/m \le 2` for second-order Runge-Kutta integration. Here, +:math:`m_C` is used as the nominal contact inertia. The actual inertia resisting slip differs because friction +acts tangentially, where the added mass differs from the transverse value used in :math:`m_{node}`, and +because of the line's axial stiffness and end conditions. Using half of the stability limit provides a margin +for these uncertainties. + + +LuGre Friction +-------------- + +When ``StickDisp`` :math:`= d > 0`, friction follows the LuGre model (Canudas de Wit et al., 1995) with the +friction normalized by the normal force. The bristle deflection :math:`z` is a continuous state of the +slip support. The bristle stiffness per unit normal force is + +.. math:: + + \sigma_0 = \frac{\mu_s}{d}, \qquad + \tilde\sigma_0 = \min\!\left(\sigma_0,\ \frac{m_C}{F_N}\left(\frac{0.5}{\Delta t}\right)^2\right), + +where the cap, evaluated with the current :math:`F_N` and :math:`m_C`, keeps the nominal bristle frequency +:math:`\sqrt{\tilde\sigma_0 F_N/m_C}` at or below :math:`0.5/\Delta t`. The bristle dynamics and friction +force are + +.. math:: + + \dot z = v - \kappa z, \qquad + \kappa = \min\!\left(\frac{|v|\,\tilde\sigma_0}{\mu(v)},\ \frac{1}{\Delta t}\right), + +.. math:: + + F_f^{*} = -\left(F_N\tilde\sigma_0 z + \sigma_1 F_N\,\dot z\right), \qquad + \sigma_1 F_N = 2\sqrt{\tilde\sigma_0 F_N\, m_C}\; e^{-(v/v_S)^2} , + +with no viscous term (:math:`\sigma_2 = 0`). The bristle damping :math:`\sigma_1` is critical damping of the +bristle stiffness acting on :math:`m_C`. It is faded out as the slip speed exceeds :math:`v_S`, which keeps the +model dissipative during sliding transients. + +In the standard LuGre model, the relaxation rate is :math:`\kappa = |v|\,\sigma_0/\mu(v)` without an upper limit. +It grows with the slip speed and the bristle stiffness, and :math:`\dot z = v - \kappa z` becomes a stiff +equation during fast sliding with a small ``StickDisp``: for example, :math:`v` = 1 m/s and :math:`d` = 1 mm +give :math:`\kappa` of order 1000 1/s. Since the explicit integration of :math:`\dot z = -\kappa z` is +stable only for :math:`\kappa\,\Delta t \le 2`, MoorDyn caps the rate at :math:`1/\Delta t`, half of the +stability limit. For second-order Runge-Kutta integration, this value also gives the fastest decay of +:math:`z` per time step. The cap is not part of the standard model. It is inactive in the sticking regime and at low slip +speeds, where the model is unchanged. When it is active, :math:`z` cannot relax fast enough to stay at the +steady value :math:`\mu(v)/\tilde\sigma_0`; during steady sliding it instead approaches +:math:`|v|\,\Delta t`, which is larger. The corresponding bristle force would exceed :math:`\mu(v)F_N`, so it +is clipped by the friction bound below. The friction force is bounded, + +.. math:: + + F_f = \operatorname{sign}(F_f^{*})\,\min\!\left(|F_f^{*}|,\ F_{max}\right), \qquad + F_{max} = + \begin{cases} + \mu_s F_N, & \kappa < 1/\Delta t \\ + \mu(v)\,F_N, & \kappa = 1/\Delta t + \end{cases} + +The :math:`\mu_s F_N` bound is a safeguard that preserves the frictional lag of the LuGre model at breakaway. +When the relaxation rate :math:`\kappa` is capped (fast sliding with a small ``StickDisp``), :math:`z` can +exceed its steady value and the :math:`\mu(v) F_N` bound restores the kinetic friction level. + +Limiting behavior: + +- Steady sliding: :math:`z \to \mu(v)\operatorname{sign}(v)/\tilde\sigma_0` and + :math:`F_f \to -\mu(v) F_N \operatorname{sign}(v)`. +- Sticking (:math:`|z| \ll \mu/\tilde\sigma_0`): :math:`\dot z \approx v`, so :math:`z` is the elastic slip + and the contact acts as a spring along the line with stiffness :math:`\tilde\sigma_0 F_N = \mu_s F_N/d`. + Breakaway occurs as :math:`|z|` approaches :math:`d`. +- Like the original LuGre model, small oscillations of :math:`v` relax :math:`z` through the :math:`|v|` + term, which can produce slow drift below breakaway when the line vibrates (Dupont et al., 2002). + + +.. _MD-SlipSupport-IC: + +Initial Conditions +------------------ + +**Line layout.** For a line with :math:`n` slip supports, the line end positions and the hoop positions (in +the order the supports are listed) define the polyline vertices + +.. math:: + + \mathbf{P}_0 = \mathbf{r}_0, \quad \mathbf{P}_i = \mathbf{x}_{h,i}\ (i = 1,\dots,n), \quad + \mathbf{P}_{n+1} = \mathbf{r}_N, \qquad + c_0 = 0, \quad c_i = c_{i-1} + \lVert\mathbf{P}_i - \mathbf{P}_{i-1}\rVert , + +where :math:`c_i` is the distance along the polyline from end A to vertex :math:`i`. Internal node :math:`j` is +placed on the polyline at the distance :math:`d_j = (s_j/L)\,c_{n+1}` from end A. With :math:`i` such that +:math:`c_{i-1} \le d_j \le c_i`, + +.. math:: + + \mathbf{r}_j = \mathbf{P}_{i-1} + \frac{d_j - c_{i-1}}{c_i - c_{i-1}}\left(\mathbf{P}_i - \mathbf{P}_{i-1}\right), + +and the initial node velocities are zero. The first force evaluation then locates each contact by searching all +segments of the line. + +**Dynamic relaxation.** During initial condition generation, the friction force is replaced by a viscous +slip damper and the LuGre state is frozen, + +.. math:: + + F_f = -c_{IC}\,v, \qquad \dot z = 0, \qquad + c_{IC} = \min\!\left(2\sqrt{2\,w_{wet}\,M_L},\ \frac{0.5\,m_{node}}{\Delta t}\right), + +where :math:`w_{wet} = |m' - \rho_w A_d|\,g` is the wet weight per unit length and +:math:`M_L = L\,(m' + \rho_w C_{at} A_d)` is the line mass including axial added mass. + +The first term is an ad-hoc damping that slows sliding at the rate :math:`c_{IC}/M_L`, which is +proportional to :math:`\sqrt{g/L}` for lines much heavier than water. The second term limits the damper to +:math:`c_{IC}\,\Delta t/m_{node} \le 0.5`, well within the explicit stability limit. The value of +:math:`c_{IC}` affects only the convergence rate of dynamic relaxation, not the equilibrium it converges to. +:math:`z` is set to zero at the start and at the end of dynamic relaxation. + + +Release +------- + +When the contact point sits at a line end (:math:`k = 1, u = 0` or :math:`k = N, u = 1`), the overshoot of +the hoop beyond the end is + +.. math:: + + o = -(\mathbf{x}_h - \mathbf{C}(s_c))\cdot\hat{\mathbf{t}(s_c)} \ \ \text{(end A)}, \qquad + o = (\mathbf{x}_h - \mathbf{C}(s_c))\cdot\hat{\mathbf{t}(s_c)} \ \ \text{(end B)} . + +The slip support releases the line permanently if :math:`o > 0.1\,l_k` in every force evaluation within an +accepted time step. Requiring this in every evaluation prevents releases triggered by intermediate +Runge-Kutta stages or linearization perturbations. + + +References +---------- + +- C. Canudas de Wit, H. Olsson, K. J. Åström, and P. Lischinsky, "A new model for control of systems with + friction," *IEEE Transactions on Automatic Control*, 40(3):419-425, 1995. +- P. Dupont, V. Hayward, B. Armstrong, and F. Altpeter, "Single state elastoplastic friction models," + *IEEE Transactions on Automatic Control*, 47(5):787-792, 2002. +- C. S. Peskin, "The immersed boundary method," *Acta Numerica*, 11:479-517, 2002. diff --git a/modules/moordyn/CMakeLists.txt b/modules/moordyn/CMakeLists.txt index ef15882385..8d79fd47fa 100644 --- a/modules/moordyn/CMakeLists.txt +++ b/modules/moordyn/CMakeLists.txt @@ -26,6 +26,7 @@ add_library(moordynlib STATIC src/MoorDyn_Misc.f90 src/MoorDyn_Point.f90 src/MoorDyn_Rod.f90 + src/MoorDyn_SlipSupport.f90 src/MoorDyn_Types.f90 ) target_link_libraries(moordynlib seastlib nwtclibs) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index 85674b96a1..78a32d5b86 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -28,6 +28,7 @@ MODULE MoorDyn USE MoorDyn_Rod USE MoorDyn_Body USE MoorDyn_Misc + USE MoorDyn_SlipSupport IMPLICIT NONE @@ -78,6 +79,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er INTEGER(IntKi) :: il ! index INTEGER(IntKi) :: iil ! index INTEGER(IntKi) :: Success ! flag for checking whether line is attached to failure point + REAL(DbKi) :: hSeg ! segment length used to size slip support kernels INTEGER(IntKi) :: I ! Current line number of input file INTEGER(IntKi) :: J ! index INTEGER(IntKi) :: K ! index @@ -443,6 +445,18 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er Line = NextLine(i) END DO + else if (INDEX(Line, "SLIP SUPPORT") > 0) then ! if slip support header + + ! skip following two lines (label line and unit line) + Line = NextLine(i) + Line = NextLine(i) + + Line = NextLine(i) + DO while (INDEX(Line, "---") == 0) ! while we DON'T find another header line + p%nSlips = p%nSlips + 1 + Line = NextLine(i) + END DO + else if (INDEX(Line, "OPTIONS") > 0) then ! if options header @@ -646,6 +660,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ALLOCATE(m%LineList( p%nLines ), STAT = ErrStat2 ); if(AllocateFailed("LineList" )) return ALLOCATE(m%ExtLdList( p%nExtLds ), STAT = ErrStat2 ); if(AllocateFailed("ExtLdList" )) return ALLOCATE(m%FailList( p%nFails ), STAT = ErrStat2 ); if(AllocateFailed("FailList" )) return + ALLOCATE(m%SlipList( p%nSlips ), STAT = ErrStat2 ); if(AllocateFailed("SlipList" )) return ! Allocate associated index arrays (note: some are allocated larger than will be used, for simplicity) @@ -2290,6 +2305,118 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er endif enddo + !------------------------------------------------------------------------------------------- + else if (INDEX(Line, "SLIP SUPPORT") > 0) then ! if slip support header + + IF (wordy > 0) CALL WrScr(" Reading slip supports") + + ! skip following two lines (label line and unit line) + Line = NextLine(i) + Line = NextLine(i) + + DO l = 1,p%nSlips + + Line = NextLine(i) + + ! ID Point Line Width K Zeta MuS MuK VTol VS StickDisp (Width as a multiple of segment length) + READ(Line,*,IOSTAT=ErrStat2) m%SlipList(l)%IdNum, m%SlipList(l)%PointID, m%SlipList(l)%LineID, m%SlipList(l)%width, & + m%SlipList(l)%kPerp, m%SlipList(l)%zeta, m%SlipList(l)%muS, & + m%SlipList(l)%muK, m%SlipList(l)%vTol, m%SlipList(l)%vS, m%SlipList(l)%stickDisp + IF (ErrStat2 /= 0) THEN + CALL SetErrStat( ErrID_Fatal, ' Unable to parse slip support '//trim(Num2LStr(l))//' on row '//trim(Num2LStr(i))//' in input file. Expecting: ID Point Line Width K Zeta MuS MuK VTol VS StickDisp.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + END IF + + IF (m%SlipList(l)%width < 1.0_DbKi) THEN + CALL SetErrStat( ErrID_Fatal, ' Slip support '//trim(Num2LStr(l))//': Width must be at least 1 segment length (2 recommended).', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + ELSE IF (m%SlipList(l)%width < 2.0_DbKi) THEN + CALL SetErrStat( ErrID_Warn, ' Slip support '//trim(Num2LStr(l))//': Width is less than 2 segment lengths; contact forces may ripple as the line slides.', ErrStat, ErrMsg, RoutineName ) + END IF + + IF ( m%SlipList(l)%IdNum .NE. l ) THEN + CALL SetErrStat( ErrID_Fatal, 'Slip support ID numbers must be sequential starting from 1.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + END IF + + IF ((m%SlipList(l)%PointID < 1) .OR. (m%SlipList(l)%PointID > p%nPoints)) THEN + CALL SetErrStat( ErrID_Fatal, ' Slip support '//trim(Num2LStr(l))//': Point number out of bounds.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + END IF + + ! the hoop reaction adds no mass, so a free hoop Point needs mass of its own or attached lines + IF ((m%PointList(m%SlipList(l)%PointID)%typeNum == 0) .AND. & + (m%PointList(m%SlipList(l)%PointID)%pointM <= 0.0_DbKi) .AND. & + (m%PointList(m%SlipList(l)%PointID)%pointV*m%PointList(m%SlipList(l)%PointID)%pointCa <= 0.0_DbKi) .AND. & + (m%PointList(m%SlipList(l)%PointID)%nAttached == 0)) THEN + CALL SetErrStat( ErrID_Fatal, ' Slip support '//trim(Num2LStr(l))//': free Point '// & + trim(Num2LStr(m%SlipList(l)%PointID))//' has no mass and no attached lines.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + END IF + + IF ((m%SlipList(l)%LineID < 1) .OR. (m%SlipList(l)%LineID > p%nLines)) THEN + CALL SetErrStat( ErrID_Fatal, ' Slip support '//trim(Num2LStr(l))//': Line number out of bounds.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + END IF + + IF (.NOT. ALLOCATED(m%LineList(m%SlipList(l)%LineID)%r)) THEN + CALL SetErrStat( ErrID_Fatal, ' The SLIP SUPPORTS section must come after the LINES section.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + END IF + + IF ((m%SlipList(l)%zeta < 0.0_DbKi) .OR. (m%SlipList(l)%muK < 0.0_DbKi) .OR. & + (m%SlipList(l)%muS < m%SlipList(l)%muK) .OR. (m%SlipList(l)%stickDisp < 0.0_DbKi) .OR. & + ((m%SlipList(l)%muS > 0.0_DbKi) .AND. (m%SlipList(l)%stickDisp <= 0.0_DbKi) .AND. (m%SlipList(l)%vTol <= 0.0_DbKi)) .OR. & + ((m%SlipList(l)%muS > m%SlipList(l)%muK) .AND. (m%SlipList(l)%vS <= 0.0_DbKi)) .OR. & + ((m%SlipList(l)%stickDisp > 0.0_DbKi) .AND. ((m%SlipList(l)%muS <= 0.0_DbKi) .OR. (m%SlipList(l)%vS <= 0.0_DbKi)))) THEN + CALL SetErrStat( ErrID_Fatal, ' Slip support '//trim(Num2LStr(l))//': require Zeta >= 0, 0 <= MuK <= MuS, StickDisp >= 0, VTol > 0 when MuS > 0 and StickDisp = 0, '// & + 'VS > 0 when MuS > MuK, and MuS > 0 and VS > 0 when StickDisp > 0.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + END IF + + ! LuGre bristle states are appended after all line states + IF (m%SlipList(l)%stickDisp > 0.0_DbKi) THEN + Nx = Nx + 1 + m%SlipList(l)%stateIs = Nx + END IF + + CALL SlipSupport_Setup(m%SlipList(l), m%LineList(m%SlipList(l)%LineID), p) + + END DO + + ! Width (read as a ratio) is scaled by the longest segment among all lines through the same hoop point + DO l = 1,p%nSlips + hSeg = 0.0_DbKi + DO il = 1,p%nSlips + IF (m%SlipList(il)%PointID == m%SlipList(l)%PointID) hSeg = MAX(hSeg, MAXVAL(m%LineList(m%SlipList(il)%LineID)%l)) + END DO + m%SlipList(l)%width = m%SlipList(l)%width*hSeg + END DO + + DO l = 1,p%nSlips + ! kernel spans +-2*Width; mirroring at the line ends assumes it reaches less than half the line + IF (4.0_DbKi*m%SlipList(l)%width >= m%LineList(m%SlipList(l)%LineID)%UnstrLen) THEN + CALL SetErrStat( ErrID_Fatal, ' Slip support '//trim(Num2LStr(l))//': Width ('//trim(Num2LStr(m%SlipList(l)%width))// & + ' m) must be less than a quarter of the length of Line '//trim(Num2LStr(m%SlipList(l)%LineID))//'. Use more segments or a smaller Width.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + END IF + + if (p%writeLog > 0) then + write(p%UnLog,'(A)') ' Slip support '//trim(Num2LStr(l))//': Point '//trim(Num2LStr(m%SlipList(l)%PointID))// & + ', Line '//trim(Num2LStr(m%SlipList(l)%LineID))//', Width = '//trim(Num2LStr(m%SlipList(l)%width))// & + ' m, K = '//trim(Num2LStr(m%SlipList(l)%kPerp))//' N/m, IC slip damping = '//trim(Num2LStr(m%SlipList(l)%cIC))//' N-s/m' + end if + END DO + !------------------------------------------------------------------------------------------- else if (INDEX(Line, "OUTPUT") > 0) then ! if output header @@ -2815,6 +2942,8 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL CheckError( ErrStat2, ErrMsg2 ) IF (ErrStat >= AbortErrLev) RETURN + IF (p%nSlips > 0) CALL SlipSupport_InitLineShape(m, p, l) + IF (wordy > 2) print *, "Line ", l, " with NumSegs =", N IF (wordy > 2) print *, "its states range from index ", m%LineStateIs1(l), " to ", m%LineStateIsN(l) @@ -2839,6 +2968,10 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er END DO !l = 1, p%NLines + DO l = 1, p%nSlips + IF (m%SlipList(l)%stateIs > 0) x%states(m%SlipList(l)%stateIs) = 0.0_DbKi + END DO + ! -------------------------------------------------------------------- @@ -3071,6 +3204,9 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! Unboost drag coefficient of each line type <<< m%IC_gen = .False. ! turn off IC_gen flag + DO I = 1, p%nSlips + IF (m%SlipList(I)%stateIs > 0) x%states(m%SlipList(I)%stateIs) = 0.0_DbKi + END DO DO I = 1, p%nLines m%LineList(I)%Cdn = m%LineList(I)%Cdn / InputFileDat%CdScaleIC m%LineList(I)%Cdt = m%LineList(I)%Cdt / InputFileDat%CdScaleIC @@ -3095,6 +3231,11 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er other%dummy = 0 xd%dummy = 0 z%dummy = 0 + + ALLOCATE(other%slipKSeg(p%nSlips), STAT = ErrStat2); if(AllocateFailed("other%slipKSeg")) return + ALLOCATE(other%slipStatus(p%nSlips), STAT = ErrStat2); if(AllocateFailed("other%slipStatus")) return + IF (p%nSlips > 0) other%slipKSeg = m%SlipList(:)%kSeg + IF (p%nSlips > 0) other%slipStatus = m%SlipList(:)%status if (InitInp%Linearize .and. ((compVIV) .OR. (compVisco))) then ErrStat2 = ErrID_Fatal @@ -3451,6 +3592,22 @@ subroutine MD_InitVars(Vars, InitInp, u, p, x, z, y, m, InitOut, Linearize, ErrS end do end do + !----------------- + ! first-order states + !----------------- + + ! LuGre bristle deflection of slip supports; perturbation scaled to the pre-sliding range + do l = 1, p%nSlips + if (m%SlipList(l)%stateIs > 0) then + LinStr = 'Slip '//trim(Num2LStr(l)) + call MV_AddVar(Vars%x, LinStr, FieldScalar, DatLoc(MD_x_states), & + iAry=m%SlipList(l)%stateIs, & + Num=1, Flags=VF_DerivOrder1, & + Perturb=0.01_R8Ki*m%SlipList(l)%stickDisp, & + LinNames=[trim(LinStr)//' bristle z, m']) + end if + end do + !------------------------------------------------------------------------- ! Input variables !------------------------------------------------------------------------- @@ -3614,14 +3771,23 @@ SUBROUTINE MD_UpdateStates( t, n, u, t_array, p, x, xd, z, other, m, ErrStat, Er NdtM = ceiling(p%dtCoupling/p%dtM0) ! get number of mooring time steps to do based on desired time step size dtM = p%dtCoupling/REAL(NdtM,DbKi) ! adjust desired time step to satisfy dt with an integer number of time steps + ! restart the contact search and release status from the state at t, in case the glue code has stepped back in time + IF (p%nSlips > 0) THEN + m%SlipList(:)%kSeg = other%slipKSeg + m%SlipList(:)%status = other%slipStatus + END IF !loop through line integration time steps DO I = 1, NdtM ! for (double ts=t; ts<=t+ICdt-dts; ts+=dts) + IF (p%nSlips > 0) CALL SlipSupport_ArmRelease(m, p) + Call MD_Step(t2, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat2, ErrMsg2) IF ( ErrStat2 /= ErrID_None ) THEN CALL CheckError(ErrStat2, ErrMsg2) END IF + + IF (p%nSlips > 0) CALL SlipSupport_CommitRelease(m, p, t2) ! check for NaNs - is this a good place/way to do it? DO J = 1, m%Nx @@ -3639,6 +3805,10 @@ SUBROUTINE MD_UpdateStates( t, n, u, t_array, p, x, xd, z, other, m, ErrStat, Er END DO ! I time steps + IF (p%nSlips > 0) THEN + other%slipKSeg = m%SlipList(:)%kSeg + other%slipStatus = m%SlipList(:)%status + END IF ! destroy dxdt and x2, and u_interp !CALL MD_DestroyContState( dxdt, ErrStat, ErrMsg) @@ -3935,6 +4105,10 @@ SUBROUTINE MD_CalcOutput( t, u, p, x, xd, z, other, y, m, ErrStat, ErrMsg ) ! call CalcContStateDeriv in order to run model and calculate dynamics with provided x and u + IF (p%nSlips > 0) THEN + m%SlipList(:)%kSeg = other%slipKSeg + m%SlipList(:)%status = other%slipStatus + END IF CALL MD_CalcContStateDeriv( t, u, p, x, xd, z, other, m, m%xdTemp, ErrStat, ErrMsg ) ! ! assign net force on fairlead Points to the fairlead force output mesh @@ -4288,6 +4462,16 @@ SUBROUTINE MD_CalcContStateDeriv( t, u, p, x, xd, z, other, m, dxdt, ErrStat, Er CALL Line_SetState(m%LineList(l), x%states(m%LineStateIs1(l):m%LineStateIsN(l)), t, m) END DO + IF (p%nSlips > 0) THEN + DO l = 1, p%nSlips + IF (m%SlipList(l)%stateIs > 0) m%SlipList(l)%z = x%states(m%SlipList(l)%stateIs) + END DO + CALL SlipSupport_CalcForces(m, p) + DO l = 1, p%nSlips + IF (m%SlipList(l)%stateIs > 0) dxdt%states(m%SlipList(l)%stateIs) = m%SlipList(l)%zdot + END DO + END IF + ! calculate dynamics of free objects (will also calculate forces (doRHS()) from any child/dependent objects)... ! calculate line dynamics (and calculate line forces and masses attributed to points) diff --git a/modules/moordyn/src/MoorDyn_IO.f90 b/modules/moordyn/src/MoorDyn_IO.f90 index d853f52748..22160daf84 100644 --- a/modules/moordyn/src/MoorDyn_IO.f90 +++ b/modules/moordyn/src/MoorDyn_IO.f90 @@ -77,6 +77,11 @@ MODULE MoorDyn_IO INTEGER, PARAMETER :: Sub = 26 INTEGER, PARAMETER :: TenA = 27 INTEGER, PARAMETER :: TenB = 28 + INTEGER, PARAMETER :: SlipS = 29 + INTEGER, PARAMETER :: SlipV = 30 + INTEGER, PARAMETER :: SlipFN = 31 + INTEGER, PARAMETER :: SlipFF = 32 + INTEGER, PARAMETER :: SlipZ = 33 ! List of units corresponding to the quantities parameters for QTypes @@ -634,6 +639,10 @@ SUBROUTINE MDIO_ProcessOutList(OutList, p, m, y, InitOut, ErrStat, ErrMsg ) p%OutParam(I)%OType = 4 ! Body object type qVal = let2 ! quantity type string + ELSE IF (Let1(1:1) == 'S') THEN ! Look for S?xxx or Slip?xxx + p%OutParam(I)%OType = 5 ! slip support object type + qVal = let2 + ! should do fairlead option also! ! error @@ -728,6 +737,21 @@ SUBROUTINE MDIO_ProcessOutList(OutList, p, m, y, InitOut, ErrStat, ErrMsg ) ELSE IF (qVal == 'SUB') THEN p%OutParam(I)%QType = Sub p%OutParam(I)%Units = UnitList(Sub) + ELSE IF ((p%OutParam(I)%OType == 5) .AND. (qVal == 'S')) THEN + p%OutParam(I)%QType = SlipS + p%OutParam(I)%Units = UnitList(PosX) + ELSE IF ((p%OutParam(I)%OType == 5) .AND. (qVal == 'VS')) THEN + p%OutParam(I)%QType = SlipV + p%OutParam(I)%Units = UnitList(VelX) + ELSE IF ((p%OutParam(I)%OType == 5) .AND. (qVal == 'FN')) THEN + p%OutParam(I)%QType = SlipFN + p%OutParam(I)%Units = UnitList(Ten) + ELSE IF ((p%OutParam(I)%OType == 5) .AND. (qVal == 'FF')) THEN + p%OutParam(I)%QType = SlipFF + p%OutParam(I)%Units = UnitList(Ten) + ELSE IF ((p%OutParam(I)%OType == 5) .AND. (qVal == 'Z')) THEN + p%OutParam(I)%QType = SlipZ + p%OutParam(I)%Units = UnitList(PosX) ELSE CALL DenoteInvalidOutput(p%OutParam(I)) ! flag as invalid CALL WrScr('Warning: invalid output specifier '//trim(OutListTmp)//'. Quantity type not recognized.') @@ -776,6 +800,12 @@ SUBROUTINE MDIO_ProcessOutList(OutList, p, m, y, InitOut, ErrStat, ErrMsg ) CALL DenoteInvalidOutput(p%OutParam(I)) ! flag as invalid END IF + ELSE IF (p%OutParam(I)%OType==5) THEN ! slip support + IF ((p%OutParam(I)%ObjID < 1) .OR. (p%OutParam(I)%ObjID > p%nSlips)) THEN + CALL WrScr('Warning: output slip support index out of range in requested output '//trim(OutListTmp)//'.') + CALL DenoteInvalidOutput(p%OutParam(I)) ! flag as invalid + END IF + END IF ! is the reverse sign functionality necessary? @@ -1638,6 +1668,30 @@ SUBROUTINE MDIO_WriteOutputs( Time, p, m, y, ErrStat, ErrMsg ) ErrMsg = ' Unsupported output quantity '//TRIM(p%OutParam(I)%Name)//' requested from Body '//TRIM(Num2Lstr(p%OutParam(I)%ObjID))//'.' END SELECT + ELSE IF (p%OutParam(I)%OType == 5) THEN ! if dealing with a slip support output + SELECT CASE (p%OutParam(I)%QType) + CASE (SlipS) + y%WriteOutput(I) = m%SlipList(p%OutParam(I)%ObjID)%s + CASE (SlipV) + y%WriteOutput(I) = m%SlipList(p%OutParam(I)%ObjID)%vSlip + CASE (SlipFN) + y%WriteOutput(I) = m%SlipList(p%OutParam(I)%ObjID)%Fn + CASE (SlipFF) + y%WriteOutput(I) = m%SlipList(p%OutParam(I)%ObjID)%Ff + CASE (SlipZ) + y%WriteOutput(I) = m%SlipList(p%OutParam(I)%ObjID)%z + CASE (FX) + y%WriteOutput(I) = m%SlipList(p%OutParam(I)%ObjID)%Fhoop(1) + CASE (FY) + y%WriteOutput(I) = m%SlipList(p%OutParam(I)%ObjID)%Fhoop(2) + CASE (FZ) + y%WriteOutput(I) = m%SlipList(p%OutParam(I)%ObjID)%Fhoop(3) + CASE DEFAULT + y%WriteOutput(I) = 0.0_ReKi + ErrStat = ErrID_Warn + ErrMsg = ' Unsupported output quantity '//TRIM(p%OutParam(I)%Name)//' requested from slip support '//TRIM(Num2Lstr(p%OutParam(I)%ObjID))//'.' + END SELECT + ELSE ! it must be an invalid output, so write zero y%WriteOutput(I) = 0.0_ReKi diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index a898e70119..e45866fda9 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -279,7 +279,7 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) ! allocate node force vectors ALLOCATE ( Line%W(3, 0:N), Line%Dp(3, 0:N), Line%Dq(3, 0:N), & Line%Ap(3, 0:N), Line%Aq(3, 0:N), Line%B(3, 0:N), Line%Bs(3, 0:N), & - Line%Lf(3, 0:N), Line%Fnet(3, 0:N), STAT = ErrStat ) + Line%Lf(3, 0:N), Line%Fnet(3, 0:N), Line%Fslip(3, 0:N), STAT = ErrStat ) IF ( ErrStat /= ErrID_None ) THEN ErrMsg = ' Error allocating node force arrays.' !CALL CleanUp() @@ -289,6 +289,7 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) ! set gravity and bottom contact forces to zero initially (because the horizontal components should remain at zero) Line%W = 0.0_DbKi Line%B = 0.0_DbKi + Line%Fslip = 0.0_DbKi ! allocate mass and inverse mass matrices for each node (including ends) ALLOCATE ( Line%S(3, 3, 0:N), Line%M(3, 3, 0:N), STAT = ErrStat ) @@ -1882,11 +1883,11 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! total forces IF (I==0) THEN - Line%Fnet(:,I) = Line%T(:,1) + Line%Td(:,1) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) + Line%Fnet(:,I) = Line%T(:,1) + Line%Td(:,1) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) + Line%Fslip(:,I) ELSE IF (I==N) THEN - Line%Fnet(:,I) = -Line%T(:,N) - Line%Td(:,N) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) + Line%Fnet(:,I) = -Line%T(:,N) - Line%Td(:,N) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) + Line%Fslip(:,I) ELSE - Line%Fnet(:,I) = Line%T(:,I+1) - Line%T(:,I) + Line%Td(:,I+1) - Line%Td(:,I) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) + Line%Fnet(:,I) = Line%T(:,I+1) - Line%T(:,I) + Line%Td(:,I+1) - Line%Td(:,I) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) + Line%Fslip(:,I) END IF END DO ! I - done looping through nodes diff --git a/modules/moordyn/src/MoorDyn_Point.f90 b/modules/moordyn/src/MoorDyn_Point.f90 index b64a5efb70..795be83d32 100644 --- a/modules/moordyn/src/MoorDyn_Point.f90 +++ b/modules/moordyn/src/MoorDyn_Point.f90 @@ -286,6 +286,13 @@ SUBROUTINE Point_DoRHS(Point, m, p) Point%Fnet(J) = Point%Fnet(J) - Point%Blin(J) * Point%rd(J) - Point%Bquad(J) * ABS(Point%rd(J)) * Point%rd(J) END DO + ! reaction from any lines threaded through this point via slip supports + DO l = 1, p%nSlips + IF ((m%SlipList(l)%PointID == Point%IdNum) .AND. (m%SlipList(l)%status == 1)) THEN + Point%Fnet = Point%Fnet + m%SlipList(l)%Fhoop + END IF + END DO + ! would this sub ever need to include the m*a inertial term? Is it ever called for coupled points? <<< END SUBROUTINE Point_DoRHS diff --git a/modules/moordyn/src/MoorDyn_Registry.txt b/modules/moordyn/src/MoorDyn_Registry.txt index e394009d73..9f10f125fa 100644 --- a/modules/moordyn/src/MoorDyn_Registry.txt +++ b/modules/moordyn/src/MoorDyn_Registry.txt @@ -303,6 +303,7 @@ typedef ^ ^ DbKi Aq {:}{:} typedef ^ ^ DbKi B {:}{:} - - "node bottom contact force" "[N]" typedef ^ ^ DbKi Bs {:}{:} - - "node force due to bending moments" "[N]" typedef ^ ^ DbKi Lf {:}{:} - - "node crossflow viv lift force. TODO this will have to be allocated for only internal after fix" "[N]" +typedef ^ ^ DbKi Fslip {:}{:} - - "node forces from slip supports" "[N]" typedef ^ ^ DbKi Fnet {:}{:} - - "total force on node" "[N]" typedef ^ ^ DbKi S {:}{:}{:} - - "node inverse mass matrix" "[kg]" typedef ^ ^ DbKi M {:}{:}{:} - - "node mass matrix" "[kg]" @@ -339,6 +340,34 @@ typedef ^ ^ DbKi failTime - typedef ^ ^ DbKi failTen - - - "tension threshold of failure" "N" typedef ^ ^ IntKi failStatus - - - "0 not failed yet, 1 failed, 2 invalid" "-" +# this is the slip support type: a line threaded through a hoop on a Point +typedef ^ MD_SlipSupport IntKi IdNum - - - "integer identifier of this slip support" "-" +typedef ^ ^ IntKi PointID - - - "index of the hoop Point in PointList" "-" +typedef ^ ^ IntKi LineID - - - "index of the threaded Line in LineList" "-" +typedef ^ ^ DbKi width - - - "cubic B-spline contact kernel scale in unstretched arc length (input as a multiple of segment length)" "[m]" +typedef ^ ^ DbKi kPerp - - - "transverse penalty stiffness" "[N/m]" +typedef ^ ^ DbKi zeta - - - "transverse damping ratio" "[-]" +typedef ^ ^ DbKi muS - - - "static friction coefficient" "[-]" +typedef ^ ^ DbKi muK - - - "kinetic friction coefficient" "[-]" +typedef ^ ^ DbKi vTol - - - "friction regularization velocity" "[m/s]" +typedef ^ ^ DbKi vS - - - "Stribeck velocity for static-to-kinetic transition" "[m/s]" +typedef ^ ^ DbKi stickDisp - - - "LuGre pre-sliding displacement at breakaway; 0 = regularized Coulomb" "[m]" +typedef ^ ^ IntKi stateIs - 0 - "index of the LuGre bristle state in x%states (0 if none)" "-" +typedef ^ ^ DbKi z - - - "LuGre bristle deflection" "[m]" +typedef ^ ^ DbKi zdot - - - "LuGre bristle deflection rate" "[m/s]" +typedef ^ ^ IntKi capWarned - 0 - "1 once the bristle stiffness cap warning has been issued" "-" +typedef ^ ^ DbKi mNode - - - "representative line node mass including added mass" "[kg]" +typedef ^ ^ DbKi cIC - - - "viscous slip damping used during IC generation" "[N-s/m]" +typedef ^ ^ IntKi status - - - "1 active, 0 released" "-" +typedef ^ ^ IntKi releasePending - 0 - "1 while every force evaluation in the current time step found the hoop past a line end" "-" +typedef ^ ^ IntKi kSeg - 0 - "active segment index (0 = not yet located)" "-" +typedef ^ ^ DbKi xi - - - "parametric contact coordinate on active segment" "[-]" +typedef ^ ^ DbKi s - - - "unstretched arc length from end A to contact" "[m]" +typedef ^ ^ DbKi vSlip - - - "slip velocity of line relative to hoop along line tangent" "[m/s]" +typedef ^ ^ DbKi Fn - - - "normal (transverse) contact force magnitude" "[N]" +typedef ^ ^ DbKi Ff - - - "tangential friction force on line along tangent" "[N]" +typedef ^ ^ DbKi Fhoop {3} - - "force exerted by the line on the hoop Point" "[N]" + # this is the MDOutParmType - a less literal alternative of the NWTC OutParmType for MoorDyn (to avoid huge lists of possible output channel permutations) typedef ^ MD_OutParmType CHARACTER(10) Name - - - "name of output channel" typedef ^ ^ CHARACTER(10) Units - - - "units string" @@ -370,6 +399,8 @@ typedef ^ ConstraintStateType SiKi dummy - ## ============================== Define Other states here: ===================================================================================================================================== typedef ^ OtherStateType SiKi dummy - - - "Remove this variable if you have other states" - +typedef ^ ^ IntKi slipKSeg {:} - - "contact segment of each slip support at the end of the last state update (search-window center)" - +typedef ^ ^ IntKi slipStatus {:} - - "status of each slip support at the end of the last state update (1: engaged, 0: released)" - ## ============================== Parameters ============================================================================================================================================ typedef ^ ParameterType IntKi nLineTypes - 0 - "number of line types" "" @@ -383,6 +414,7 @@ typedef ^ ^ IntKi nSyropeLineICs - typedef ^ ^ IntKi nExtLds - 0 - "number of external loads or damping" "" typedef ^ ^ IntKi nCtrlChans - 0 - "number of distinct control channels specified for use as inputs" "" typedef ^ ^ IntKi nFails - 0 - "number of failure conditions" "" +typedef ^ ^ IntKi nSlips - 0 - "number of slip supports" "" typedef ^ ^ IntKi nFreeBodies - 0 - "" "" typedef ^ ^ IntKi nFreeRods - 0 - "" "" typedef ^ ^ IntKi nFreePoints - 0 - "" "" @@ -487,6 +519,7 @@ typedef ^ ^ MD_Point PointList {:} typedef ^ ^ MD_Line LineList {:} - - "array of line objects" - typedef ^ ^ MD_ExtLd ExtLdList {:} - - "array of external load objects" - typedef ^ ^ MD_Fail FailList {:} - - "array of line objects" - +typedef ^ ^ MD_SlipSupport SlipList {:} - - "array of slip support objects" - typedef ^ ^ IntKi FreePointIs {:} - - "array of free point indices in PointList vector" "" typedef ^ ^ IntKi CpldPointIs {:}{:} - - "array of coupled/fairlead point indices in PointList vector" "" typedef ^ ^ IntKi FreeRodIs {:} - - "array of free rod indices in RodList vector" "" diff --git a/modules/moordyn/src/MoorDyn_SlipSupport.f90 b/modules/moordyn/src/MoorDyn_SlipSupport.f90 new file mode 100644 index 0000000000..dbaca5d970 --- /dev/null +++ b/modules/moordyn/src/MoorDyn_SlipSupport.f90 @@ -0,0 +1,645 @@ +!********************************************************************************************************************************** +! LICENSING +! Copyright (C) 2026 National Laboratory of the Rockies +! +! This file is part of MoorDyn. +! +! Licensed under the Apache License, Version 2.0 (the "License"); +! you may not use this file except in compliance with the License. +! You may obtain a copy of the License at +! +! http://www.apache.org/licenses/LICENSE-2.0 +! +! Unless required by applicable law or agreed to in writing, software +! distributed under the License is distributed on an "AS IS" BASIS, +! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +! See the License for the specific language governing permissions and +! limitations under the License. +! +!********************************************************************************************************************************** +! Slip supports: a line threaded through a hoop on a Point. The hoop acts on a virtual contact point on the +! line; contact forces are spread to nearby nodes with a cubic B-spline kernel of fixed width in arc length. +MODULE MoorDyn_SlipSupport + + USE MoorDyn_Types + USE NWTC_Library + + IMPLICIT NONE + + PRIVATE + + PUBLIC :: SlipSupport_Setup + PUBLIC :: SlipSupport_InitLineShape + PUBLIC :: SlipSupport_CalcForces + PUBLIC :: SlipSupport_ArmRelease + PUBLIC :: SlipSupport_CommitRelease + +CONTAINS + + !-------------------------------------------------------------- + SUBROUTINE SlipSupport_Setup(Slip, Line, p) + + TYPE(MD_SlipSupport), INTENT(INOUT) :: Slip + TYPE(MD_Line), INTENT(IN ) :: Line + TYPE(MD_ParameterType), INTENT(IN ) :: p + + REAL(DbKi) :: area, w, wWet, mLine + + area = 0.25_DbKi*Pi*Line%d**2 + w = Line%rho*area + Slip%mNode = Line%UnstrLen/REAL(Line%N, DbKi) * (w + p%rhoW*Line%Can*area) + + IF (Slip%kPerp <= 0.0_DbKi) Slip%kPerp = Slip%mNode*(0.3_DbKi/p%dtM0)**2 + + ! ad-hoc slip damping, rate ~sqrt(g/L) for heavy lines, capped for explicit stability + wWet = ABS(w - p%rhoW*area)*p%g + mLine = Line%UnstrLen*(w + p%rhoW*Line%Cat*area) + Slip%cIC = MIN(2.0_DbKi*SQRT(2.0_DbKi*wWet*mLine), 0.5_DbKi*Slip%mNode/p%dtM0) + + Slip%status = 1 + Slip%releasePending = 0 + Slip%z = 0.0_DbKi + Slip%zdot = 0.0_DbKi + Slip%capWarned = 0 + Slip%kSeg = 0 + Slip%xi = 0.0_DbKi + Slip%s = 0.0_DbKi + Slip%vSlip = 0.0_DbKi + Slip%Fn = 0.0_DbKi + Slip%Ff = 0.0_DbKi + Slip%Fhoop = 0.0_DbKi + + END SUBROUTINE SlipSupport_Setup + !-------------------------------------------------------------- + + + !-------------------------------------------------------------- + ! Lay out a line's internal nodes along the polyline A -> hoop(s) -> B, with hoops in input order. + SUBROUTINE SlipSupport_InitLineShape(m, p, iLine) + + TYPE(MD_MiscVarType), INTENT(INOUT) :: m + TYPE(MD_ParameterType), INTENT(IN ) :: p + INTEGER(IntKi), INTENT(IN ) :: iLine + + INTEGER(IntKi) :: is, nV, I, J, N + REAL(DbKi), ALLOCATABLE :: rV(:,:), cV(:) + REAL(DbKi) :: sNode, dTarget, frac + + N = m%LineList(iLine)%N + + nV = 0 + DO is = 1, p%nSlips + IF (m%SlipList(is)%LineID == iLine) nV = nV + 1 + END DO + IF (nV == 0) RETURN + + ALLOCATE(rV(3, 0:nV+1), cV(0:nV+1)) + + rV(:,0) = m%LineList(iLine)%r(:,0) + J = 0 + DO is = 1, p%nSlips + IF (m%SlipList(is)%LineID == iLine) THEN + J = J + 1 + rV(:,J) = m%PointList(m%SlipList(is)%PointID)%r + END IF + END DO + rV(:,nV+1) = m%LineList(iLine)%r(:,N) + + cV(0) = 0.0_DbKi + DO J = 1, nV+1 + cV(J) = cV(J-1) + SQRT(SUM((rV(:,J) - rV(:,J-1))**2)) + END DO + + IF (cV(nV+1) <= 0.0_DbKi) THEN + DEALLOCATE(rV, cV) + RETURN + END IF + + ! node arc length is mapped proportionally onto the polyline chord length + sNode = 0.0_DbKi + DO I = 1, N-1 + sNode = sNode + m%LineList(iLine)%l(I) + dTarget = sNode/m%LineList(iLine)%UnstrLen*cV(nV+1) + DO J = 1, nV+1 + IF ((dTarget <= cV(J)) .OR. (J == nV+1)) EXIT + END DO + frac = (dTarget - cV(J-1))/MAX(cV(J) - cV(J-1), 1.0E-12_DbKi) + m%LineList(iLine)%r(:,I) = rV(:,J-1) + frac*(rV(:,J) - rV(:,J-1)) + m%LineList(iLine)%rd(:,I) = 0.0_DbKi + END DO + + J = 0 + DO is = 1, p%nSlips + IF (m%SlipList(is)%LineID == iLine) THEN + J = J + 1 + m%SlipList(is)%s = cV(J)/cV(nV+1)*m%LineList(iLine)%UnstrLen + ! seed the contact search at the hoop's segment so the first search stays local + sNode = 0.0_DbKi + DO I = 1, N + sNode = sNode + m%LineList(iLine)%l(I) + IF (sNode >= m%SlipList(is)%s) EXIT + END DO + m%SlipList(is)%kSeg = MIN(I, N) + END IF + END DO + + ! viscoelastic (non-Syrope) lines: reset the static stretch to the new geometry, as in Line_Initialize + IF ((m%LineList(iLine)%ElasticMod > 1) .AND. (m%LineList(iLine)%ElasticMod < 4)) THEN + DO I = 1, N + m%LineList(iLine)%dl_1(I) = SQRT(SUM((m%LineList(iLine)%r(:,I) - m%LineList(iLine)%r(:,I-1))**2)) & + - m%LineList(iLine)%l(I) + END DO + END IF + + DEALLOCATE(rV, cV) + + END SUBROUTINE SlipSupport_InitLineShape + !-------------------------------------------------------------- + + + !-------------------------------------------------------------- + ! Compute contact forces of all slip supports, filling Line%Fslip and Slip%Fhoop. + ! Must be called after line and point kinematics are set and before Line_GetStateDeriv. + SUBROUTINE SlipSupport_CalcForces(m, p) + + TYPE(MD_MiscVarType), INTENT(INOUT) :: m + TYPE(MD_ParameterType), INTENT(IN ) :: p + + INTEGER(IntKi) :: is, iL, N, k, kLo, kHi, kBest, j, jj, nW, jEnd + REAL(DbKi) :: xh(3), vh(3), u, d2, d2Best, uBest + REAL(DbKi) :: tC(3), tCmag, rStar(3), delta(3), vRelj(3), vRelPerp(3) + REAL(DbKi) :: Fperp(3), Fd(3), F(3), Fsum(3), mC, cPerp, vTolEff, muEff, overshoot + REAL(DbKi) :: sig0, kCap, kappa, sig1Fn + LOGICAL :: rateCapped + CHARACTER(200) :: Msg + INTEGER(IntKi), ALLOCATABLE :: idx(:) + REAL(DbKi), ALLOCATABLE :: w(:), dw(:), wg(:), wN(:), sNode(:), qN(:,:), Fnode(:,:) + + DO is = 1, p%nSlips + m%LineList(m%SlipList(is)%LineID)%Fslip = 0.0_DbKi + END DO + + DO is = 1, p%nSlips + + m%SlipList(is)%zdot = 0.0_DbKi + IF (m%SlipList(is)%status /= 1) CYCLE + + iL = m%SlipList(is)%LineID + N = m%LineList(iL)%N + xh = m%PointList(m%SlipList(is)%PointID)%r + vh = m%PointList(m%SlipList(is)%PointID)%rd + + IF (ALLOCATED(sNode)) THEN + IF (SIZE(sNode) /= N+1) DEALLOCATE(idx, w, dw, wg, wN, sNode, qN, Fnode) + END IF + IF (.NOT. ALLOCATED(sNode)) ALLOCATE(idx(N+1), w(N+1), dw(N+1), wg(N+1), wN(N+1), sNode(0:N), qN(3,N+1), Fnode(3,N+1)) + sNode(0) = 0.0_DbKi + DO j = 1, N + sNode(j) = sNode(j-1) + m%LineList(iL)%l(j) + END DO + + ! search neighbors of the cached segment; full search when not yet located + IF (m%SlipList(is)%kSeg < 1) THEN + kLo = 1 + kHi = N + ELSE + kLo = MAX(1, m%SlipList(is)%kSeg - 2) + kHi = MIN(N, m%SlipList(is)%kSeg + 2) + END IF + + d2Best = HUGE(1.0_DbKi) + kBest = kLo + uBest = 0.0_DbKi + DO k = kLo, kHi + CALL ProjectOnSegment(m%LineList(iL), m%SlipList(is)%width, sNode, k, xh, u, d2, idx, w, dw) + IF (d2 < d2Best) THEN + d2Best = d2 + kBest = k + uBest = u + END IF + END DO + + k = kBest + u = uBest + m%SlipList(is)%kSeg = k + m%SlipList(is)%xi = u + + CALL KernelWeights(m%LineList(iL), m%SlipList(is)%width, sNode, k, u, nW, idx, w, dw, wg, jEnd) + + rStar = 0.0_DbKi + tC = 0.0_DbKi + DO j = 1, nW + rStar = rStar + w(j) *m%LineList(iL)%r(:,idx(j)) + tC = tC + dw(j)*m%LineList(iL)%r(:,idx(j)) + END DO + + tCmag = SQRT(DOT_PRODUCT(tC, tC)) + IF (tCmag > 1.0E-6_DbKi) THEN + tC = tC/tCmag + ELSE + tC = UnitOrZero(m%LineList(iL)%r(:,k) - m%LineList(iL)%r(:,k-1)) + IF (DOT_PRODUCT(tC, tC) <= 0.0_DbKi) CYCLE + END IF + + ! a release is only committed by SlipSupport_CommitRelease after an accepted time step + overshoot = 0.0_DbKi + IF ((k == 1) .AND. (u <= 0.0_DbKi)) overshoot = -DOT_PRODUCT(xh - rStar, tC) + IF ((k == N) .AND. (u >= 1.0_DbKi)) overshoot = DOT_PRODUCT(xh - rStar, tC) + IF (overshoot <= 0.1_DbKi*m%LineList(iL)%l(k)) m%SlipList(is)%releasePending = 0 + + delta = rStar - xh + delta = delta - DOT_PRODUCT(delta, tC)*tC + + ! the spring force is split per node along each node's own tangent, so a frictionless hoop never + ! pushes nodes along the line (avoids tension trapping at sharp bends); friction uses the same + ! per-node tangents for its velocity and force so that it can only dissipate energy. + ! The end node's folded share wg(j) of a mirrored node carries that node's tangent (point reflection). + wN(1:nW) = w(1:nW) + IF (jEnd > 0) wN(jEnd) = w(jEnd) - SUM(wg(1:nW)) + + vRelPerp = 0.0_DbKi + m%SlipList(is)%vSlip = 0.0_DbKi + DO j = 1, nW + qN(:,j) = NodeTangent(m%LineList(iL), idx(j)) + IF (DOT_PRODUCT(qN(:,j), qN(:,j)) <= 0.0_DbKi) qN(:,j) = tC + vRelj = m%LineList(iL)%rd(:,idx(j)) - vh + vRelPerp = vRelPerp + w(j)*vRelj + m%SlipList(is)%vSlip = m%SlipList(is)%vSlip + wN(j)*DOT_PRODUCT(vRelj, qN(:,j)) + END DO + IF (jEnd > 0) THEN + vRelj = m%LineList(iL)%rd(:,idx(jEnd)) - vh + DO j = 1, nW + IF (wg(j) == 0.0_DbKi) CYCLE + m%SlipList(is)%vSlip = m%SlipList(is)%vSlip + wg(j)*DOT_PRODUCT(vRelj, qN(:,j)) + END DO + END IF + vRelPerp = vRelPerp - DOT_PRODUCT(vRelPerp, tC)*tC + + ! inertia felt at the contact point when a force is shared among nodes of mass mNode + mC = m%SlipList(is)%mNode/SUM(w(1:nW)**2) + cPerp = 2.0_DbKi*m%SlipList(is)%zeta*SQRT(m%SlipList(is)%kPerp*mC) + Fperp = -m%SlipList(is)%kPerp*delta + ! the damper acts on the interpolated velocity and is spread with the same weights, so it stays + ! dissipative without the per-node projections that weaken it at sharp bends + Fd = -cPerp*vRelPerp + + Fsum = 0.0_DbKi + DO j = 1, nW + Fnode(:,j) = wN(j)*(Fperp - DOT_PRODUCT(Fperp, qN(:,j))*qN(:,j)) + w(j)*Fd + Fsum = Fsum + Fnode(:,j) + END DO + IF (jEnd > 0) THEN + DO j = 1, nW + IF (wg(j) == 0.0_DbKi) CYCLE + F = wg(j)*(Fperp - DOT_PRODUCT(Fperp, qN(:,j))*qN(:,j)) + Fnode(:,jEnd) = Fnode(:,jEnd) + F + Fsum = Fsum + F + END DO + END IF + + ! normal load actually delivered to the line sets the friction capacity + m%SlipList(is)%Fn = SQRT(DOT_PRODUCT(Fsum, Fsum)) + + IF (m%IC_gen) THEN + m%SlipList(is)%Ff = -m%SlipList(is)%cIC*m%SlipList(is)%vSlip + ELSE IF (m%SlipList(is)%muS > 0.0_DbKi) THEN + IF (m%SlipList(is)%muS > m%SlipList(is)%muK) THEN + muEff = m%SlipList(is)%muK + (m%SlipList(is)%muS - m%SlipList(is)%muK)*EXP(-(m%SlipList(is)%vSlip/m%SlipList(is)%vS)**2) + ELSE + muEff = m%SlipList(is)%muK + END IF + + IF (m%SlipList(is)%stateIs > 0) THEN + ! LuGre with friction normalized by Fn; bristle stiffness capped for explicit stability + sig0 = m%SlipList(is)%muS/m%SlipList(is)%stickDisp + IF (m%SlipList(is)%Fn > 0.0_DbKi) THEN + kCap = mC*(0.5_DbKi/p%dtM0)**2/m%SlipList(is)%Fn + IF (sig0 > kCap) THEN + sig0 = kCap + IF (m%SlipList(is)%capWarned == 0) THEN + m%SlipList(is)%capWarned = 1 + Msg = 'Slip support '//TRIM(Num2LStr(m%SlipList(is)%IdNum))//': StickDisp too small for dtM; '// & + 'bristle stiffness capped, so sticking is softer than specified.' + CALL WrScr(' '//TRIM(Msg)) + IF (p%writeLog > 0) WRITE(p%UnLog,'(A)') TRIM(Msg) + END IF + END IF + END IF + + ! relaxation rate |v|/z_ss with z_ss = muEff/sig0, capped at 1/dtM (fastest RK2 decay, half the stability limit) + kappa = 1.0_DbKi/p%dtM0 + rateCapped = .TRUE. + IF (ABS(m%SlipList(is)%vSlip)*sig0 < kappa*muEff) THEN + kappa = ABS(m%SlipList(is)%vSlip)*sig0/muEff + rateCapped = .FALSE. + END IF + m%SlipList(is)%zdot = m%SlipList(is)%vSlip - kappa*m%SlipList(is)%z + + ! critically damped bristles, faded out once sliding to keep LuGre dissipative + sig1Fn = 2.0_DbKi*SQRT(sig0*m%SlipList(is)%Fn*mC)*EXP(-(m%SlipList(is)%vSlip/m%SlipList(is)%vS)**2) + + m%SlipList(is)%Ff = -(m%SlipList(is)%Fn*sig0*m%SlipList(is)%z + sig1Fn*m%SlipList(is)%zdot) + ! muS bound keeps LuGre's frictional lag; muEff bound only where the capped rate lets z overshoot z_ss + IF (rateCapped) THEN + m%SlipList(is)%Ff = SIGN(MIN(ABS(m%SlipList(is)%Ff), muEff*m%SlipList(is)%Fn), m%SlipList(is)%Ff) + ELSE + m%SlipList(is)%Ff = SIGN(MIN(ABS(m%SlipList(is)%Ff), m%SlipList(is)%muS*m%SlipList(is)%Fn), m%SlipList(is)%Ff) + END IF + ELSE + ! floor vTol so the peak friction slope muS*Fn/vTol stays within explicit stability + vTolEff = MAX(m%SlipList(is)%vTol, m%SlipList(is)%muS*m%SlipList(is)%Fn*p%dtM0/mC) + m%SlipList(is)%Ff = -muEff*m%SlipList(is)%Fn*TANH(m%SlipList(is)%vSlip/vTolEff) + END IF + ELSE + m%SlipList(is)%Ff = 0.0_DbKi + END IF + + Fsum = 0.0_DbKi + DO j = 1, nW + F = Fnode(:,j) + wN(j)*m%SlipList(is)%Ff*qN(:,j) + IF (j == jEnd) THEN + DO jj = 1, nW + IF (wg(jj) /= 0.0_DbKi) F = F + wg(jj)*m%SlipList(is)%Ff*qN(:,jj) + END DO + END IF + m%LineList(iL)%Fslip(:,idx(j)) = m%LineList(iL)%Fslip(:,idx(j)) + F + Fsum = Fsum + F + END DO + m%SlipList(is)%Fhoop = -Fsum + + m%SlipList(is)%s = sNode(k-1) + u*m%LineList(iL)%l(k) + + END DO + + END SUBROUTINE SlipSupport_CalcForces + !-------------------------------------------------------------- + + + !-------------------------------------------------------------- + ! Closest contact parameter u in [0,1] on segment k of the contact curve, and squared distance to xh. + SUBROUTINE ProjectOnSegment(Line, width, sNode, k, xh, u, d2, idx, w, dw) + + TYPE(MD_Line), INTENT(IN ) :: Line + REAL(DbKi), INTENT(IN ) :: width + REAL(DbKi), INTENT(IN ) :: sNode(0:) + INTEGER(IntKi), INTENT(IN ) :: k + REAL(DbKi), INTENT(IN ) :: xh(3) + REAL(DbKi), INTENT( OUT) :: u + REAL(DbKi), INTENT( OUT) :: d2 + INTEGER(IntKi), INTENT(INOUT) :: idx(:) + REAL(DbKi), INTENT(INOUT) :: w(:) + REAL(DbKi), INTENT(INOUT) :: dw(:) + + INTEGER(IntKi) :: it, j, nW + REAL(DbKi) :: a(3), lvec(3), L2, C(3), Cp(3), h, uNew + + a = Line%r(:,k-1) + lvec = Line%r(:,k) - a + L2 = DOT_PRODUCT(lvec, lvec) + IF (L2 > 0.0_DbKi) THEN + u = MIN(1.0_DbKi, MAX(0.0_DbKi, DOT_PRODUCT(xh - a, lvec)/L2)) + ELSE + u = 0.0_DbKi + END IF + + ! Gauss-Newton on the smoothed contact curve, started from the chord projection + DO it = 1, 10 + CALL KernelWeights(Line, width, sNode, k, u, nW, idx, w, dw) + C = 0.0_DbKi + Cp = 0.0_DbKi + DO j = 1, nW + C = C + w(j) *Line%r(:,idx(j)) + Cp = Cp + dw(j)*Line%r(:,idx(j)) + END DO + h = DOT_PRODUCT(Cp, Cp) + IF (h <= 0.0_DbKi) EXIT + uNew = MIN(1.0_DbKi, MAX(0.0_DbKi, u - DOT_PRODUCT(C - xh, Cp)/h)) + IF (ABS(uNew - u) < 1.0E-10_DbKi) THEN + u = uNew + EXIT + END IF + u = uNew + END DO + + CALL KernelWeights(Line, width, sNode, k, u, nW, idx, w, dw) + C = 0.0_DbKi + DO j = 1, nW + C = C + w(j)*Line%r(:,idx(j)) + END DO + d2 = SUM((C - xh)**2) + + END SUBROUTINE ProjectOnSegment + !-------------------------------------------------------------- + + + !-------------------------------------------------------------- + ! Weights from a cubic B-spline kernel of fixed scale 'width' in unstretched arc length, normalized to + ! sum to one. Nodes are mirrored about each line end (r -> 2*r_end - r) and folded back, so the contact + ! curve passes through the end nodes. dw is d/du on segment k. Optional wg(j) is the share of the folded + ! end node (entry jEnd, 0 if none) that comes from the mirror image of node idx(j). + SUBROUTINE KernelWeights(Line, width, sNode, k, u, nW, idx, w, dw, wg, jEnd) + + TYPE(MD_Line), INTENT(IN ) :: Line + REAL(DbKi), INTENT(IN ) :: width + REAL(DbKi), INTENT(IN ) :: sNode(0:) + INTEGER(IntKi), INTENT(IN ) :: k + REAL(DbKi), INTENT(IN ) :: u + INTEGER(IntKi), INTENT( OUT) :: nW + INTEGER(IntKi), INTENT( OUT) :: idx(:) + REAL(DbKi), INTENT( OUT) :: w(:) + REAL(DbKi), INTENT( OUT) :: dw(:) + REAL(DbKi), INTENT( OUT), OPTIONAL :: wg(:) + INTEGER(IntKi), INTENT( OUT), OPTIONAL :: jEnd + + INTEGER(IntKi) :: N, j, jLo, jHi + REAL(DbKi) :: s, R, Ltot, Asum, dAsum, aK, daK + + N = Line%N + Ltot = sNode(N) + s = sNode(k-1) + u*Line%l(k) + R = 2.0_DbKi*width + + jLo = k-1 + DO WHILE (jLo > 0) + IF (sNode(jLo-1) <= s - R) EXIT + jLo = jLo - 1 + END DO + jHi = k + DO WHILE (jHi < N) + IF (sNode(jHi+1) >= s + R) EXIT + jHi = jHi + 1 + END DO + + nW = jHi - jLo + 1 + DO j = jLo, jHi + idx(j-jLo+1) = j + END DO + w(1:nW) = 0.0_DbKi + dw(1:nW) = 0.0_DbKi + Asum = 0.0_DbKi + dAsum = 0.0_DbKi + IF (PRESENT(wg)) wg(1:nW) = 0.0_DbKi + IF (PRESENT(jEnd)) jEnd = 0 + + DO j = jLo, jHi + CALL Sample(sNode(j), j) + CALL Add(j, aK, daK) + END DO + + ! the kernel reaches less than half the line, so at most one end is mirrored + IF (s - R < 0.0_DbKi) THEN + IF (PRESENT(jEnd)) jEnd = 1 - jLo + DO j = 1, N + IF (-sNode(j) <= s - R) EXIT + CALL Sample(-sNode(j), j) + CALL Add(0, 2.0_DbKi*aK, 2.0_DbKi*daK) + CALL Add(j, -aK, -daK) + IF (PRESENT(wg)) wg(j-jLo+1) = wg(j-jLo+1) + 2.0_DbKi*aK + END DO + END IF + + IF (s + R > Ltot) THEN + IF (PRESENT(jEnd)) jEnd = N - jLo + 1 + DO j = 1, N + IF (2.0_DbKi*Ltot - sNode(N-j) >= s + R) EXIT + CALL Sample(2.0_DbKi*Ltot - sNode(N-j), N-j) + CALL Add(N, 2.0_DbKi*aK, 2.0_DbKi*daK) + CALL Add(N-j, -aK, -daK) + IF (PRESENT(wg)) wg(N-j-jLo+1) = wg(N-j-jLo+1) + 2.0_DbKi*aK + END DO + END IF + + ! normalize; Width >= longest segment guarantees Asum > 0 + w(1:nW) = w(1:nW)/Asum + dw(1:nW) = (dw(1:nW) - w(1:nW)*dAsum)/Asum * Line%l(k) + IF (PRESENT(wg)) wg(1:nW) = wg(1:nW)/Asum + + CONTAINS + + ! raw kernel value (a) and its d/ds (da) for a node at arc position sj with tributary length of node jt + SUBROUTINE Sample(sj, jt) + REAL(DbKi), INTENT(IN) :: sj + INTEGER(IntKi), INTENT(IN) :: jt + REAL(DbKi) :: t, x, Kx, dKx, ds + + t = (sj - s)/width + x = ABS(t) + IF (x < 1.0_DbKi) THEN + Kx = (4.0_DbKi - 6.0_DbKi*x**2 + 3.0_DbKi*x**3)/6.0_DbKi + dKx = (-12.0_DbKi*x + 9.0_DbKi*x**2)/6.0_DbKi + ELSE IF (x < 2.0_DbKi) THEN + Kx = (2.0_DbKi - x)**3/6.0_DbKi + dKx = -(2.0_DbKi - x)**2/2.0_DbKi + ELSE + Kx = 0.0_DbKi + dKx = 0.0_DbKi + END IF + + ! tributary length; end nodes count as interior because of mirroring + IF (jt == 0) THEN + ds = Line%l(1) + ELSE IF (jt == N) THEN + ds = Line%l(N) + ELSE + ds = 0.5_DbKi*(Line%l(jt) + Line%l(jt+1)) + END IF + + aK = Kx*ds + daK = -SIGN(1.0_DbKi, t)*dKx*ds/width + END SUBROUTINE Sample + + SUBROUTINE Add(jn, ai, dai) + INTEGER(IntKi), INTENT(IN) :: jn + REAL(DbKi), INTENT(IN) :: ai, dai + w( jn-jLo+1) = w( jn-jLo+1) + ai + dw(jn-jLo+1) = dw(jn-jLo+1) + dai + Asum = Asum + ai + dAsum = dAsum + dai + END SUBROUTINE Add + + END SUBROUTINE KernelWeights + !-------------------------------------------------------------- + + + !-------------------------------------------------------------- + ! Arm release checks before a time step; CalcForces clears the flag at any evaluation that is not past a line end. + SUBROUTINE SlipSupport_ArmRelease(m, p) + + TYPE(MD_MiscVarType), INTENT(INOUT) :: m + TYPE(MD_ParameterType), INTENT(IN ) :: p + + INTEGER(IntKi) :: is + + DO is = 1, p%nSlips + m%SlipList(is)%releasePending = m%SlipList(is)%status + END DO + + END SUBROUTINE SlipSupport_ArmRelease + !-------------------------------------------------------------- + + + !-------------------------------------------------------------- + ! Permanently release supports whose hoop was past a line end in every force evaluation of the last time step. + SUBROUTINE SlipSupport_CommitRelease(m, p, t) + + TYPE(MD_MiscVarType), INTENT(INOUT) :: m + TYPE(MD_ParameterType), INTENT(IN ) :: p + REAL(DbKi), INTENT(IN ) :: t + + INTEGER(IntKi) :: is + CHARACTER(200) :: Msg + + DO is = 1, p%nSlips + IF ((m%SlipList(is)%status == 1) .AND. (m%SlipList(is)%releasePending == 1)) THEN + m%SlipList(is)%status = 0 + m%SlipList(is)%releasePending = 0 + m%SlipList(is)%vSlip = 0.0_DbKi + m%SlipList(is)%Fn = 0.0_DbKi + m%SlipList(is)%Ff = 0.0_DbKi + m%SlipList(is)%Fhoop = 0.0_DbKi + Msg = 'Slip support '//TRIM(Num2LStr(m%SlipList(is)%IdNum))//' released Line '// & + TRIM(Num2LStr(m%SlipList(is)%LineID))//' at t = '//TRIM(Num2LStr(t))//' s.' + CALL WrScr(' '//TRIM(Msg)) + IF (p%writeLog > 0) WRITE(p%UnLog,'(A)') TRIM(Msg) + END IF + END DO + + END SUBROUTINE SlipSupport_CommitRelease + !-------------------------------------------------------------- + + + !-------------------------------------------------------------- + ! unit bisector of the adjacent segment directions at node i (end nodes use their single segment) + FUNCTION NodeTangent(Line, i) RESULT(q) + + TYPE(MD_Line), INTENT(IN) :: Line + INTEGER(IntKi), INTENT(IN) :: i + REAL(DbKi) :: q(3) + + q = 0.0_DbKi + IF (i > 0) q = q + UnitOrZero(Line%r(:,i) - Line%r(:,i-1)) + IF (i < Line%N) q = q + UnitOrZero(Line%r(:,i+1) - Line%r(:,i)) + q = UnitOrZero(q) + + END FUNCTION NodeTangent + !-------------------------------------------------------------- + + + !-------------------------------------------------------------- + FUNCTION UnitOrZero(v) RESULT(u) + + REAL(DbKi), INTENT(IN) :: v(3) + REAL(DbKi) :: u(3) + REAL(DbKi) :: vMag + + vMag = SQRT(DOT_PRODUCT(v, v)) + IF (vMag > 0.0_DbKi) THEN + u = v/vMag + ELSE + u = 0.0_DbKi + END IF + + END FUNCTION UnitOrZero + !-------------------------------------------------------------- + +END MODULE MoorDyn_SlipSupport diff --git a/modules/moordyn/src/MoorDyn_Types.f90 b/modules/moordyn/src/MoorDyn_Types.f90 index 2f6bfd87bf..c5f4cf1400 100644 --- a/modules/moordyn/src/MoorDyn_Types.f90 +++ b/modules/moordyn/src/MoorDyn_Types.f90 @@ -324,6 +324,7 @@ MODULE MoorDyn_Types REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: B !< node bottom contact force [[N]] REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: Bs !< node force due to bending moments [[N]] REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: Lf !< node crossflow viv lift force. TODO this will have to be allocated for only internal after fix [[N]] + REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: Fslip !< node forces from slip supports [[N]] REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: Fnet !< total force on node [[N]] REAL(DbKi) , DIMENSION(:,:,:), ALLOCATABLE :: S !< node inverse mass matrix [[kg]] REAL(DbKi) , DIMENSION(:,:,:), ALLOCATABLE :: M !< node mass matrix [[kg]] @@ -363,6 +364,36 @@ MODULE MoorDyn_Types INTEGER(IntKi) :: failStatus = 0_IntKi !< 0 not failed yet, 1 failed, 2 invalid [-] END TYPE MD_Fail ! ======================= +! ========= MD_SlipSupport ======= + TYPE, PUBLIC :: MD_SlipSupport + INTEGER(IntKi) :: IdNum = 0_IntKi !< integer identifier of this slip support [-] + INTEGER(IntKi) :: PointID = 0_IntKi !< index of the hoop Point in PointList [-] + INTEGER(IntKi) :: LineID = 0_IntKi !< index of the threaded Line in LineList [-] + REAL(DbKi) :: width = 0.0_R8Ki !< cubic B-spline contact kernel scale in unstretched arc length (input as a multiple of segment length) [[m]] + REAL(DbKi) :: kPerp = 0.0_R8Ki !< transverse penalty stiffness [[N/m]] + REAL(DbKi) :: zeta = 0.0_R8Ki !< transverse damping ratio [[-]] + REAL(DbKi) :: muS = 0.0_R8Ki !< static friction coefficient [[-]] + REAL(DbKi) :: muK = 0.0_R8Ki !< kinetic friction coefficient [[-]] + REAL(DbKi) :: vTol = 0.0_R8Ki !< friction regularization velocity [[m/s]] + REAL(DbKi) :: vS = 0.0_R8Ki !< Stribeck velocity for static-to-kinetic transition [[m/s]] + REAL(DbKi) :: stickDisp = 0.0_R8Ki !< LuGre pre-sliding displacement at breakaway; 0 = regularized Coulomb [[m]] + INTEGER(IntKi) :: stateIs = 0 !< index of the LuGre bristle state in x%states (0 if none) [-] + REAL(DbKi) :: z = 0.0_R8Ki !< LuGre bristle deflection [[m]] + REAL(DbKi) :: zdot = 0.0_R8Ki !< LuGre bristle deflection rate [[m/s]] + INTEGER(IntKi) :: capWarned = 0 !< 1 once the bristle stiffness cap warning has been issued [-] + REAL(DbKi) :: mNode = 0.0_R8Ki !< representative line node mass including added mass [[kg]] + REAL(DbKi) :: cIC = 0.0_R8Ki !< viscous slip damping used during IC generation [[N-s/m]] + INTEGER(IntKi) :: status = 0_IntKi !< 1 active, 0 released [-] + INTEGER(IntKi) :: releasePending = 0 !< 1 while every force evaluation in the current time step found the hoop past a line end [-] + INTEGER(IntKi) :: kSeg = 0 !< active segment index (0 = not yet located) [-] + REAL(DbKi) :: xi = 0.0_R8Ki !< parametric contact coordinate on active segment [[-]] + REAL(DbKi) :: s = 0.0_R8Ki !< unstretched arc length from end A to contact [[m]] + REAL(DbKi) :: vSlip = 0.0_R8Ki !< slip velocity of line relative to hoop along line tangent [[m/s]] + REAL(DbKi) :: Fn = 0.0_R8Ki !< normal (transverse) contact force magnitude [[N]] + REAL(DbKi) :: Ff = 0.0_R8Ki !< tangential friction force on line along tangent [[N]] + REAL(DbKi) , DIMENSION(1:3) :: Fhoop = 0.0_R8Ki !< force exerted by the line on the hoop Point [[N]] + END TYPE MD_SlipSupport +! ======================= ! ========= MD_OutParmType ======= TYPE, PUBLIC :: MD_OutParmType CHARACTER(10) :: Name !< name of output channel [-] @@ -405,6 +436,8 @@ MODULE MoorDyn_Types ! ========= MD_OtherStateType ======= TYPE, PUBLIC :: MD_OtherStateType REAL(SiKi) :: dummy = 0.0_R4Ki !< Remove this variable if you have other states [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: slipKSeg !< contact segment of each slip support at the end of the last state update (search-window center) [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: slipStatus !< status of each slip support at the end of the last state update (1: engaged, 0: released) [-] END TYPE MD_OtherStateType ! ======================= ! ========= MD_ParameterType ======= @@ -420,6 +453,7 @@ MODULE MoorDyn_Types INTEGER(IntKi) :: nExtLds = 0 !< number of external loads or damping [] INTEGER(IntKi) :: nCtrlChans = 0 !< number of distinct control channels specified for use as inputs [] INTEGER(IntKi) :: nFails = 0 !< number of failure conditions [] + INTEGER(IntKi) :: nSlips = 0 !< number of slip supports [] INTEGER(IntKi) :: nFreeBodies = 0 !< [] INTEGER(IntKi) :: nFreeRods = 0 !< [] INTEGER(IntKi) :: nFreePoints = 0 !< [] @@ -518,6 +552,7 @@ MODULE MoorDyn_Types TYPE(MD_Line) , DIMENSION(:), ALLOCATABLE :: LineList !< array of line objects [-] TYPE(MD_ExtLd) , DIMENSION(:), ALLOCATABLE :: ExtLdList !< array of external load objects [-] TYPE(MD_Fail) , DIMENSION(:), ALLOCATABLE :: FailList !< array of line objects [-] + TYPE(MD_SlipSupport) , DIMENSION(:), ALLOCATABLE :: SlipList !< array of slip support objects [-] INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: FreePointIs !< array of free point indices in PointList vector [] INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: CpldPointIs !< array of coupled/fairlead point indices in PointList vector [] INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: FreeRodIs !< array of free rod indices in RodList vector [] @@ -2164,6 +2199,18 @@ subroutine MD_CopyLine(SrcLineData, DstLineData, CtrlCode, ErrStat, ErrMsg) end if DstLineData%Lf = SrcLineData%Lf end if + if (allocated(SrcLineData%Fslip)) then + LB(1:2) = lbound(SrcLineData%Fslip) + UB(1:2) = ubound(SrcLineData%Fslip) + if (.not. allocated(DstLineData%Fslip)) then + allocate(DstLineData%Fslip(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstLineData%Fslip.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstLineData%Fslip = SrcLineData%Fslip + end if if (allocated(SrcLineData%Fnet)) then LB(1:2) = lbound(SrcLineData%Fnet) UB(1:2) = ubound(SrcLineData%Fnet) @@ -2360,6 +2407,9 @@ subroutine MD_DestroyLine(LineData, ErrStat, ErrMsg) if (allocated(LineData%Lf)) then deallocate(LineData%Lf) end if + if (allocated(LineData%Fslip)) then + deallocate(LineData%Fslip) + end if if (allocated(LineData%Fnet)) then deallocate(LineData%Fnet) end if @@ -2463,6 +2513,7 @@ subroutine MD_PackLine(RF, Indata) call RegPackAlloc(RF, InData%B) call RegPackAlloc(RF, InData%Bs) call RegPackAlloc(RF, InData%Lf) + call RegPackAlloc(RF, InData%Fslip) call RegPackAlloc(RF, InData%Fnet) call RegPackAlloc(RF, InData%S) call RegPackAlloc(RF, InData%M) @@ -2561,6 +2612,7 @@ subroutine MD_UnPackLine(RF, OutData) call RegUnpackAlloc(RF, OutData%B); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%Bs); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%Lf); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Fslip); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%Fnet); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%S); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%M); if (RegCheckErr(RF, RoutineName)) return @@ -2690,6 +2742,119 @@ subroutine MD_UnPackFail(RF, OutData) call RegUnpack(RF, OutData%failStatus); if (RegCheckErr(RF, RoutineName)) return end subroutine +subroutine MD_CopySlipSupport(SrcSlipSupportData, DstSlipSupportData, CtrlCode, ErrStat, ErrMsg) + type(MD_SlipSupport), intent(in) :: SrcSlipSupportData + type(MD_SlipSupport), intent(inout) :: DstSlipSupportData + integer(IntKi), intent(in ) :: CtrlCode + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'MD_CopySlipSupport' + ErrStat = ErrID_None + ErrMsg = '' + DstSlipSupportData%IdNum = SrcSlipSupportData%IdNum + DstSlipSupportData%PointID = SrcSlipSupportData%PointID + DstSlipSupportData%LineID = SrcSlipSupportData%LineID + DstSlipSupportData%width = SrcSlipSupportData%width + DstSlipSupportData%kPerp = SrcSlipSupportData%kPerp + DstSlipSupportData%zeta = SrcSlipSupportData%zeta + DstSlipSupportData%muS = SrcSlipSupportData%muS + DstSlipSupportData%muK = SrcSlipSupportData%muK + DstSlipSupportData%vTol = SrcSlipSupportData%vTol + DstSlipSupportData%vS = SrcSlipSupportData%vS + DstSlipSupportData%stickDisp = SrcSlipSupportData%stickDisp + DstSlipSupportData%stateIs = SrcSlipSupportData%stateIs + DstSlipSupportData%z = SrcSlipSupportData%z + DstSlipSupportData%zdot = SrcSlipSupportData%zdot + DstSlipSupportData%capWarned = SrcSlipSupportData%capWarned + DstSlipSupportData%mNode = SrcSlipSupportData%mNode + DstSlipSupportData%cIC = SrcSlipSupportData%cIC + DstSlipSupportData%status = SrcSlipSupportData%status + DstSlipSupportData%releasePending = SrcSlipSupportData%releasePending + DstSlipSupportData%kSeg = SrcSlipSupportData%kSeg + DstSlipSupportData%xi = SrcSlipSupportData%xi + DstSlipSupportData%s = SrcSlipSupportData%s + DstSlipSupportData%vSlip = SrcSlipSupportData%vSlip + DstSlipSupportData%Fn = SrcSlipSupportData%Fn + DstSlipSupportData%Ff = SrcSlipSupportData%Ff + DstSlipSupportData%Fhoop = SrcSlipSupportData%Fhoop +end subroutine + +subroutine MD_DestroySlipSupport(SlipSupportData, ErrStat, ErrMsg) + type(MD_SlipSupport), intent(inout) :: SlipSupportData + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'MD_DestroySlipSupport' + ErrStat = ErrID_None + ErrMsg = '' +end subroutine + +subroutine MD_PackSlipSupport(RF, Indata) + type(RegFile), intent(inout) :: RF + type(MD_SlipSupport), intent(in) :: InData + character(*), parameter :: RoutineName = 'MD_PackSlipSupport' + if (RF%ErrStat >= AbortErrLev) return + call RegPack(RF, InData%IdNum) + call RegPack(RF, InData%PointID) + call RegPack(RF, InData%LineID) + call RegPack(RF, InData%width) + call RegPack(RF, InData%kPerp) + call RegPack(RF, InData%zeta) + call RegPack(RF, InData%muS) + call RegPack(RF, InData%muK) + call RegPack(RF, InData%vTol) + call RegPack(RF, InData%vS) + call RegPack(RF, InData%stickDisp) + call RegPack(RF, InData%stateIs) + call RegPack(RF, InData%z) + call RegPack(RF, InData%zdot) + call RegPack(RF, InData%capWarned) + call RegPack(RF, InData%mNode) + call RegPack(RF, InData%cIC) + call RegPack(RF, InData%status) + call RegPack(RF, InData%releasePending) + call RegPack(RF, InData%kSeg) + call RegPack(RF, InData%xi) + call RegPack(RF, InData%s) + call RegPack(RF, InData%vSlip) + call RegPack(RF, InData%Fn) + call RegPack(RF, InData%Ff) + call RegPack(RF, InData%Fhoop) + if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine MD_UnPackSlipSupport(RF, OutData) + type(RegFile), intent(inout) :: RF + type(MD_SlipSupport), intent(inout) :: OutData + character(*), parameter :: RoutineName = 'MD_UnPackSlipSupport' + if (RF%ErrStat /= ErrID_None) return + call RegUnpack(RF, OutData%IdNum); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%PointID); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%LineID); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%width); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%kPerp); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%zeta); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%muS); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%muK); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%vTol); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%vS); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%stickDisp); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%stateIs); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%z); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%zdot); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%capWarned); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%mNode); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%cIC); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%status); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%releasePending); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%kSeg); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%xi); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%s); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%vSlip); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Fn); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Ff); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Fhoop); if (RegCheckErr(RF, RoutineName)) return +end subroutine + subroutine MD_CopyOutParmType(SrcOutParmTypeData, DstOutParmTypeData, CtrlCode, ErrStat, ErrMsg) type(MD_OutParmType), intent(in) :: SrcOutParmTypeData type(MD_OutParmType), intent(inout) :: DstOutParmTypeData @@ -3047,10 +3212,36 @@ subroutine MD_CopyOtherState(SrcOtherStateData, DstOtherStateData, CtrlCode, Err integer(IntKi), intent(in ) :: CtrlCode integer(IntKi), intent( out) :: ErrStat character(*), intent( out) :: ErrMsg + integer(B4Ki) :: LB(1), UB(1) + integer(IntKi) :: ErrStat2 character(*), parameter :: RoutineName = 'MD_CopyOtherState' ErrStat = ErrID_None ErrMsg = '' DstOtherStateData%dummy = SrcOtherStateData%dummy + if (allocated(SrcOtherStateData%slipKSeg)) then + LB(1:1) = lbound(SrcOtherStateData%slipKSeg) + UB(1:1) = ubound(SrcOtherStateData%slipKSeg) + if (.not. allocated(DstOtherStateData%slipKSeg)) then + allocate(DstOtherStateData%slipKSeg(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstOtherStateData%slipKSeg.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstOtherStateData%slipKSeg = SrcOtherStateData%slipKSeg + end if + if (allocated(SrcOtherStateData%slipStatus)) then + LB(1:1) = lbound(SrcOtherStateData%slipStatus) + UB(1:1) = ubound(SrcOtherStateData%slipStatus) + if (.not. allocated(DstOtherStateData%slipStatus)) then + allocate(DstOtherStateData%slipStatus(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstOtherStateData%slipStatus.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstOtherStateData%slipStatus = SrcOtherStateData%slipStatus + end if end subroutine subroutine MD_DestroyOtherState(OtherStateData, ErrStat, ErrMsg) @@ -3060,6 +3251,12 @@ subroutine MD_DestroyOtherState(OtherStateData, ErrStat, ErrMsg) character(*), parameter :: RoutineName = 'MD_DestroyOtherState' ErrStat = ErrID_None ErrMsg = '' + if (allocated(OtherStateData%slipKSeg)) then + deallocate(OtherStateData%slipKSeg) + end if + if (allocated(OtherStateData%slipStatus)) then + deallocate(OtherStateData%slipStatus) + end if end subroutine subroutine MD_PackOtherState(RF, Indata) @@ -3068,6 +3265,8 @@ subroutine MD_PackOtherState(RF, Indata) character(*), parameter :: RoutineName = 'MD_PackOtherState' if (RF%ErrStat >= AbortErrLev) return call RegPack(RF, InData%dummy) + call RegPackAlloc(RF, InData%slipKSeg) + call RegPackAlloc(RF, InData%slipStatus) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -3075,8 +3274,13 @@ subroutine MD_UnPackOtherState(RF, OutData) type(RegFile), intent(inout) :: RF type(MD_OtherStateType), intent(inout) :: OutData character(*), parameter :: RoutineName = 'MD_UnPackOtherState' + integer(B4Ki) :: LB(1), UB(1) + integer(IntKi) :: stat + logical :: IsAllocAssoc if (RF%ErrStat /= ErrID_None) return call RegUnpack(RF, OutData%dummy); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%slipKSeg); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%slipStatus); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine MD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) @@ -3103,6 +3307,7 @@ subroutine MD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%nExtLds = SrcParamData%nExtLds DstParamData%nCtrlChans = SrcParamData%nCtrlChans DstParamData%nFails = SrcParamData%nFails + DstParamData%nSlips = SrcParamData%nSlips DstParamData%nFreeBodies = SrcParamData%nFreeBodies DstParamData%nFreeRods = SrcParamData%nFreeRods DstParamData%nFreePoints = SrcParamData%nFreePoints @@ -3501,6 +3706,7 @@ subroutine MD_PackParam(RF, Indata) call RegPack(RF, InData%nExtLds) call RegPack(RF, InData%nCtrlChans) call RegPack(RF, InData%nFails) + call RegPack(RF, InData%nSlips) call RegPack(RF, InData%nFreeBodies) call RegPack(RF, InData%nFreeRods) call RegPack(RF, InData%nFreePoints) @@ -3611,6 +3817,7 @@ subroutine MD_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%nExtLds); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nCtrlChans); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nFails); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%nSlips); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nFreeBodies); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nFreeRods); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nFreePoints); if (RegCheckErr(RF, RoutineName)) return @@ -4307,6 +4514,22 @@ subroutine MD_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) if (ErrStat >= AbortErrLev) return end do end if + if (allocated(SrcMiscData%SlipList)) then + LB(1:1) = lbound(SrcMiscData%SlipList) + UB(1:1) = ubound(SrcMiscData%SlipList) + if (.not. allocated(DstMiscData%SlipList)) then + allocate(DstMiscData%SlipList(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstMiscData%SlipList.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + do i1 = LB(1), UB(1) + call MD_CopySlipSupport(SrcMiscData%SlipList(i1), DstMiscData%SlipList(i1), CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + end do + end if if (allocated(SrcMiscData%FreePointIs)) then LB(1:1) = lbound(SrcMiscData%FreePointIs) UB(1:1) = ubound(SrcMiscData%FreePointIs) @@ -4651,6 +4874,15 @@ subroutine MD_DestroyMisc(MiscData, ErrStat, ErrMsg) end do deallocate(MiscData%FailList) end if + if (allocated(MiscData%SlipList)) then + LB(1:1) = lbound(MiscData%SlipList) + UB(1:1) = ubound(MiscData%SlipList) + do i1 = LB(1), UB(1) + call MD_DestroySlipSupport(MiscData%SlipList(i1), ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + end do + deallocate(MiscData%SlipList) + end if if (allocated(MiscData%FreePointIs)) then deallocate(MiscData%FreePointIs) end if @@ -4803,6 +5035,15 @@ subroutine MD_PackMisc(RF, Indata) call MD_PackFail(RF, InData%FailList(i1)) end do end if + call RegPack(RF, allocated(InData%SlipList)) + if (allocated(InData%SlipList)) then + call RegPackBounds(RF, 1, lbound(InData%SlipList), ubound(InData%SlipList)) + LB(1:1) = lbound(InData%SlipList) + UB(1:1) = ubound(InData%SlipList) + do i1 = LB(1), UB(1) + call MD_PackSlipSupport(RF, InData%SlipList(i1)) + end do + end if call RegPackAlloc(RF, InData%FreePointIs) call RegPackAlloc(RF, InData%CpldPointIs) call RegPackAlloc(RF, InData%FreeRodIs) @@ -4955,6 +5196,19 @@ subroutine MD_UnPackMisc(RF, OutData) call MD_UnpackFail(RF, OutData%FailList(i1)) ! FailList end do end if + if (allocated(OutData%SlipList)) deallocate(OutData%SlipList) + call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return + if (IsAllocAssoc) then + call RegUnpackBounds(RF, 1, LB, UB); if (RegCheckErr(RF, RoutineName)) return + allocate(OutData%SlipList(LB(1):UB(1)),stat=stat) + if (stat /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating OutData%SlipList.', RF%ErrStat, RF%ErrMsg, RoutineName) + return + end if + do i1 = LB(1), UB(1) + call MD_UnpackSlipSupport(RF, OutData%SlipList(i1)) ! SlipList + end do + end if call RegUnpackAlloc(RF, OutData%FreePointIs); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%CpldPointIs); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%FreeRodIs); if (RegCheckErr(RF, RoutineName)) return diff --git a/openfast_io/openfast_io/FAST_reader.py b/openfast_io/openfast_io/FAST_reader.py index 23b629f087..fae0a8ecea 100644 --- a/openfast_io/openfast_io/FAST_reader.py +++ b/openfast_io/openfast_io/FAST_reader.py @@ -3532,7 +3532,7 @@ def read_MoorDyn(self, moordyn_file): f.readline() # Lines - self.fst_vt['MoorDyn']['Line_ID'] = [] + self.fst_vt['MoorDyn']['Line_ID'] = [] self.fst_vt['MoorDyn']['LineType'] = [] self.fst_vt['MoorDyn']['AttachA'] = [] self.fst_vt['MoorDyn']['AttachB'] = [] @@ -3551,6 +3551,38 @@ def read_MoorDyn(self, moordyn_file): data_line = readline_filterComments(f).split() data_line = ''.join(data_line) # re-join for reading next section uniformly + elif '---' in data_line and 'slipsupport' in data_line: + f.readline() + f.readline() + + # Lines + self.fst_vt['MoorDyn']['SlipSupport_ID'] = [] + self.fst_vt['MoorDyn']['SlipSupport_Point'] = [] + self.fst_vt['MoorDyn']['SlipSupport_Line'] = [] + self.fst_vt['MoorDyn']['SlipSupport_Width'] = [] + self.fst_vt['MoorDyn']['SlipSupport_K'] = [] + self.fst_vt['MoorDyn']['SlipSupport_Zeta'] = [] + self.fst_vt['MoorDyn']['SlipSupport_MuS'] = [] + self.fst_vt['MoorDyn']['SlipSupport_MuK'] = [] + self.fst_vt['MoorDyn']['SlipSupport_VTol'] = [] + self.fst_vt['MoorDyn']['SlipSupport_VS'] = [] + self.fst_vt['MoorDyn']['SlipSupport_StickDisp'] = [] + data_line = readline_filterComments(f).split() + while data_line[0] and data_line[0][:3] != '---': # OpenFAST searches for ---, so we'll do the same + self.fst_vt['MoorDyn']['SlipSupport_ID'].append(int(data_line[0])) + self.fst_vt['MoorDyn']['SlipSupport_Point'].append(int(data_line[1])) + self.fst_vt['MoorDyn']['SlipSupport_Line'].append(int(data_line[2])) + self.fst_vt['MoorDyn']['SlipSupport_Width'].append(float(data_line[3])) + self.fst_vt['MoorDyn']['SlipSupport_K'].append(float(data_line[4])) + self.fst_vt['MoorDyn']['SlipSupport_Zeta'].append(float(data_line[5])) + self.fst_vt['MoorDyn']['SlipSupport_MuS'].append(float(data_line[6])) + self.fst_vt['MoorDyn']['SlipSupport_MuK'].append(float(data_line[7])) + self.fst_vt['MoorDyn']['SlipSupport_VTol'].append(float(data_line[8])) + self.fst_vt['MoorDyn']['SlipSupport_VS'].append(float(data_line[9])) + self.fst_vt['MoorDyn']['SlipSupport_StickDisp'].append(float(data_line[10])) + data_line = readline_filterComments(f).split() + data_line = ''.join(data_line) # re-join for reading next section uniformly + elif 'failure' in data_line.lower(): f.readline() f.readline() diff --git a/openfast_io/openfast_io/FAST_writer.py b/openfast_io/openfast_io/FAST_writer.py index cda77e046a..b4aefbf14e 100644 --- a/openfast_io/openfast_io/FAST_writer.py +++ b/openfast_io/openfast_io/FAST_writer.py @@ -2753,6 +2753,25 @@ def write_MoorDyn(self): ln.append('{:^11}'.format(self.fst_vt['MoorDyn']['Outputs'][i])) f.write(" ".join(ln) + '\n') + if 'SlipSupport_ID' in self.fst_vt['MoorDyn'] and self.fst_vt['MoorDyn']['SlipSupport_ID']: # there are slip support inputs + f.write('---------------------- SLIP SUPPORTS ----------------------\n') + f.write(" ".join(['{:<11s}'.format(i) for i in ['ID', 'Point', 'Line', 'Width', 'K', 'Zeta', 'MuS', 'MuK', 'VTol', 'VS', 'StickDisp']])+'\n') + f.write(" ".join(['{:<11s}'.format(i) for i in ['(-)', '(-)', '(-)', '(h)', '(N/m)', '(-)', '(-)', '(-)', '(m/s)', '(m/s)', '(m)']])+'\n') + for i in range(len(self.fst_vt['MoorDyn']['SlipSupport_ID'])): + ln = [] + ln.append('{:<7d}'.format(self.fst_vt['MoorDyn']['SlipSupport_ID'][i])) + ln.append('{:^11d}'.format(self.fst_vt['MoorDyn']['SlipSupport_Point'][i])) + ln.append('{:^11d}'.format(self.fst_vt['MoorDyn']['SlipSupport_Line'][i])) + ln.append('{:^11.4f}'.format(self.fst_vt['MoorDyn']['SlipSupport_Width'][i])) + ln.append('{:^11.4e}'.format(self.fst_vt['MoorDyn']['SlipSupport_K'][i])) + ln.append('{:^11.4f}'.format(self.fst_vt['MoorDyn']['SlipSupport_Zeta'][i])) + ln.append('{:^11.4f}'.format(self.fst_vt['MoorDyn']['SlipSupport_MuS'][i])) + ln.append('{:^11.4f}'.format(self.fst_vt['MoorDyn']['SlipSupport_MuK'][i])) + ln.append('{:^11.4e}'.format(self.fst_vt['MoorDyn']['SlipSupport_VTol'][i])) + ln.append('{:^11.4e}'.format(self.fst_vt['MoorDyn']['SlipSupport_VS'][i])) + ln.append('{:^11.4e}'.format(self.fst_vt['MoorDyn']['SlipSupport_StickDisp'][i])) + f.write(" ".join(ln) + '\n') + if 'Failure_ID' in self.fst_vt['MoorDyn'] and self.fst_vt['MoorDyn']['Failure_ID']: # there are failure inputs f.write('---------------------- FAILURE ----------------------\n') f.write('FailureID Point Line(s) FailTime FailTen\n') diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index 68671d506e..9348eafd67 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -384,6 +384,7 @@ of_regression("IEA22MW_ModalDamping" "openfast;beamdyn;servody of_regression("IEA22MW_ModalDampingLoose" "openfast;beamdyn;servodyn") of_regression("OC6_phaseII" "openfast;soildyn;subdyn;hydrodyn;offshore;stc") of_regression("MinimalExample" "openfast;elastodyn") +of_regression("Floating_Cylinder_Guide_Ropes" "openfast;hydrodyn;moordyn;offshore") of_aeromap_regression("5MW_Land_AeroMap" "aeromap;elastodyn;aerodyn") @@ -582,6 +583,7 @@ py_md_regression("py_md_5MW_OC4Semi" "moordyn;python") #md_regression("md_Single_Line_Quasi_Static_Test" "moordyn") md_regression("md_viscoelastic" "moordyn") md_regression("md_syrope" "moordyn") +md_regression("md_SlipSupport_LuGre" "moordyn") # OpenFAST IO Library regression tests py_openfast_io_library_pytest("openfast_io_library" "openfast_io;python") diff --git a/reg_tests/r-test b/reg_tests/r-test index f735abb001..77895d4f91 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit f735abb0011ded84cfac4ec414b72949bda539a5 +Subproject commit 77895d4f91ad511a6aa0c01807759473fa76d92b diff --git a/vs-build/modules/MoorDyn.vfproj b/vs-build/modules/MoorDyn.vfproj index d2e951ba95..73736f560a 100644 --- a/vs-build/modules/MoorDyn.vfproj +++ b/vs-build/modules/MoorDyn.vfproj @@ -118,6 +118,7 @@ +