Generalized Quasi-Static Mooring System Modeling with Analytic Jacobians

Abstract

This paper presents a generalized and efficient method for quasi-static analysis of mooring systems, including complex scenarios such as when shared mooring lines interconnect multiple floating wind or wave energy devices. While quasi-static mooring models are well established, most published formulations are focused on specific applications, and no publicly available implementations provide efficient handling of large mooring system networks. The present formulation addresses these gaps by: (1) formulating solutions for edge cases not typically supported by quasi-static models; (2) creating a fully generalized model structure such that any combination of mooring lines, point masses, and floating bodies can be assembled; and (3) deriving analytic expressions for the system Jacobians (stiffness matrices) so that systems with many degrees of freedom can be solved efficiently. These techniques form the theory basis of MoorPy, an open-source mooring analysis library. The model is demonstrated on nine scenarios of increasing complexity with features of interest for offshore renewable energy applications. When compared with steady-state results from a lumped-mass dynamic model, the results show that the quasi-static formulation accurately calculates profiles and tensions and that its analytic approach provides more efficient and reliable computation of system stiffness matrices than finite-differencing methods. These results verify the accuracy of the MoorPy model.

1. Introduction

This paper presents a generalized and efficient formulation for quasi-static modeling of mooring systems. Quasi-static models are useful for a variety of uses in mooring system design and analysis because of their computational efficiency and simplicity. There is considerable literature available on quasi-static mooring models. However, many formulations are limited to specific mooring configurations, unable to model unusual mooring line profiles, or limited in their ability to model interconnected arrays of floating structures.

New mooring system designs for floating wind turbines and wave energy converters have raised the need for more generalized quasi-static models. Spar platforms for floating wind turbines typically use bridle mooring configurations where each mooring line splits into two before attaching to the platform [1]. Shared mooring systems have been proposed that interconnect floating wind turbine platforms with mooring lines [2,3]. Shared mooring systems have also been considered for wave energy converters, for which complex mooring configurations with weights and floats are also common [4,5]. Efficient and versatile mooring analysis methods for such scenarios are needed to support the advancement of these offshore renewable energy technologies.

Responding to this need, the present quasi-static model formulation is general in that it can apply to any mooring system arrangement, including multiple interconnected floating bodies. It is also highly efficient thanks to the use of analytic Jacobians (stiffness matrices) for every part of the system. This formulation, which has not been presented previously, is the foundation of the open-source quasi-static mooring model MoorPy [6].

1.1. Motivation for Quasi-Static Mooring Modeling

Quasi-static mooring models solve the static equilibrium of a mooring system while neglecting inertial and hydrodynamic phenomena. The term quasi-static reflects the fact that such models can be used in time domain simulations where a floating body is in motion; they simply compute the instantaneous mooring system positions and reaction forces by making the approximation that the mooring system is in static equilibrium at each instant in time. As such, quasi-static mooring models can be incorporated into coupled dynamics models of floating systems where wave-, wind-, and current-induced loads are present on other parts of the system, while including only static tension and stiffness-related restoring behavior for the mooring system.

Quasi-static mooring models have long been a staple of mooring system design, as they are simpler and computationally faster than dynamic models [7]. The quasi-static approximation was historically used in mooring design for offshore oil and gas structures [8], and has also been used for research and development of floating offshore renewable energy devices. However, floating wind turbines and wave energy converters generally experience larger environmental loads relative to their size than offshore oil and gas structures, meaning that neglecting mooring system dynamic effects becomes less justifiable. For wave energy converters, Davidson and Ringwood [5] have provided a thorough review of mooring modeling options and the trade-off between speed and accuracy. For floating wind turbines, Hall et al. [9] compared quasi-static and dynamic mooring models within a coupled aero-hydro-servo-elastic simulation of three floating wind turbine designs (one from each stability class). The results showed that while platform motions had little sensitivity to the mooring model, the mooring line tensions could be significantly affected, particularly for catenary mooring systems undergoing larger motions. Further comparisons of quasi-static and dynamic mooring models for semisubmersible floating wind turbines [10,11] reinforce the conclusion that while quasi-static models provide a good approximation of the restoring forces that affect floating platform motions, they can significantly underpredict mooring tension amplitudes.

Despite their accuracy limitations, quasi-static models have unsurpassed computational efficiency that is ideal for initial concept evaluation and iteration and any other uses where fast analysis is important. As a result, quasi-static models are still widely used and remain an active research topic in the literature.

1.2. Quasi-Static Mooring Model Developments and Limitations

Quasi-static models work by calculating the instantaneous state of a mooring line, assuming that it is under static equilibrium and free of hydrodynamic loads. In other words, the forces involved are limited to net weight and axial stiffness. With these simplifications, each homogeneous mooring line section will have the shape of a catenary, for which a set of nonlinear analytic equations is well established (e.g., [12,13], and as presented in the next section). The catenary equations are often adjusted to account for a portion of the line resting along the seabed, sometimes including the effects of seabed friction or a sloped seabed as well. Regardless, the solution process generally involves numerically solving a system of two or three nonlinear equations for each mooring line section using an iterative process.

Many implementations have been made for quasi-static mooring models, some of which are notable for supporting more complicated situations. Huang and Vassalos [14] provided a semi-analytic formulation for fully suspended cables that works in three dimensions and can include current drag forces. Chai et al. [15] developed a generalized quasi-static model with seabed contact, including seabed slope and the ability to incorporate multiple segments in order to handle discontinuities in the seabed surface. These are examples of more general and advanced quasi-static implementations; however, they are not publicly available. On the other hand, Jonkman [16] implemented a straightforward two-dimensional formulation that includes elasticity and friction on a level seabed. This solution was included in the FAST open-source floating wind turbine simulator. Various commercial codes, such as OrcaFlex, are also known to include quasi-static mooring models.

One preexisting example of a general open-source quasi-static mooring system model is the Mooring Analysis Program (MAP) [17]. It follows a general multisegmented approach, meaning that it supports any arrangement of mooring line sections connected together. MAP is part of the OpenFAST floating wind turbine simulator, and is also available as a Python package that facilitates use in other applications. Shortcomings of MAP include that it cannot solve certain line profiles, it does not support more than one floating body (meaning that it cannot simulate coupled floating systems), and its approach for computing the mooring stiffness matrix of a body uses finite differences.

It is more efficient and accurate to compute a mooring system’s stiffness matrix analytically. Al Solihat and Nahon formulated analytic stiffness matrices for several mooring system topologies [18]. Amaral, Pesce, and Franzini presented a more general formulation for analytic mooring system stiffness matrices that rely on the force–displacement relations of individual mooring lines [19,20]. Neither of these papers provided a fully generalized approach or a publicly available implementation.

Considering the many methods described in the literature and the comparatively few publicly available models, the range of quasi-static mooring models available to researchers has been very limited, and is missing an efficient means of analyzing novel scenarios involving multiple floating bodies.

1.3. Objectives of This Work

The quasi-static mooring model formulation presented in this paper was developed to address the limitations of existing available models. Specifically, it aims to provide the following:

• Robust modeling of individual mooring line sections with support for less-common situations, such as when lines on the seabed become full slack or rest entirely on the seabed.
• A fully general mooring system formulation allowing any number of floating bodies to be specified.
• Analytic computation of system Jacobians for greater efficiency and robustness.

This paper presents the quasi-static model formulation that was developed to meet these objectives. The formulation was implemented in the Python-based MoorPy mooring analysis library. Section 2 describes the quasi-static formulation for individual mooring line sections in a range of different situations; Section 3 describes the formulation for connecting line sections with point or rigid body objects to form a complete mooring system, along with how the system Jacobians can be computed analytically and how to solve for the equilibrium state; and Section 4 presents results for a set of test scenarios of increasing complexity, comparing the results against equilibrium solutions from MoorDyn to verify the approach.

2. Quasi-Static Line Section Model

The first step of a quasi-static model of a mooring system is to model each mooring line section. Here, “mooring line section” refers to a length of mooring line over which the mechanical properties are uniform and there are no applied loads aside from weight, buoyancy, and seabed contact. Loads or imposed motions are only applied at the ends of the line section. Attachments or changes in properties along a mooring line are handled by connecting multiple sections together using the capabilities described in Section 3.

2.1. Modeling a Basic Mooring Line Section

As is typical of most quasi-static mooring models, a mooring line section is defined by its unstretched length L and two cross-sectional properties: wet weight per unit length w and axial stiffness $EA$. The wet weight per unit length (or net weight considering buoyancy) can be calculated as

$$w=({\rho }_{line}-{\rho }_{water})g\frac{\pi d^2}{4},$$

(1)

where $\rho$ is density of the mooring line material or the water, g is the gravitational constant, and d is the section’s volume-equivalent diameter (the diameter of a circle that corresponds to the displaced volume per unit length of the mooring line [11]). The elasticity, represented by the cross-sectional stiffness coefficient $EA$, can be considered as a product of the elasticity modulus and cross-sectional area; however, as this neglects elastic contributions from the line material’s internal structure, $EA$ is usually treated as a single constant for a given mooring line type.

In the formulation presented here, hydrodynamics and seabed friction are neglected and the seabed is assumed to be flat and horizontal. Although seabed friction is often included in quasi-static models, it is excluded here because it adds causality (or path-dependence) to the quasi-static solution of more complex mooring systems, which is at odds with the generalized approach developed in this paper.

With these simplifications, the only forces acting on the suspended portion of a mooring line are tension and weight. The shape and tension distribution of a mooring line can be derived by considering a small element of unstretched length $ds$ along a free-hanging catenary mooring line (Figure 1). The tension at one end of the element is denoted as T, resulting in an axial strain of $T/EA$ that stretches the element. The element has weight $wds$, which must be equal to the difference in the vertical component of tension at each end of the element.

Figure 1 shows some conventions used in this work. The line ends are labeled A and B, of which A is generally the lower end from which measurements are made. The vertical tension components at each end are denoted $V_A$ and $V_B$, and are considered positive in the directions shown in the figure. The horizontal tension component H is assumed to be constant over the section length, as there are no other horizontal forces applied over the line section length under the present modeling assumptions. The span between the mooring line ends, or the horizontal and vertical distances from end A to end B, are denoted by $x_B$ and $z_B$, respectively. The unstretched arc length measured from end A is denoted by s.

In combination with the definition of uniform material properties, the existence of solely vertical loads along the mooring line section allows for a quasi-static model to be constructed on the basis of several key principles:

• The line segment profile falls along a vertical plane, and as such can be analyzed in only two dimensions.
• Because only vertical loads are applied along the segment length, the horizontal component of the tension along the segment will be uniform.
• Because the horizontal tension component H is uniform, the local incline angle $\phi$ and total tension magnitude T are interrelated at any point on the line by

  $$Tcos\phi =H.$$

  (2)
• For suspended portions of the line segment, the change in vertical tension component is equal to the wet weight per unit length of the line segment.

These principles result in the characteristic catenary curve that describes the profile of a suspended cable hanging under its own weight. Solving the differential relations illustrated in Figure 1 over the length of a suspended line gives rise to the well-known catenary equations (see e.g., [16]), which describe the profile of a suspended cable hanging under its own weight. The following equations describe the horizontal and vertical coordinates along the profile as a function of the unstretched arc length s, as measured from end A of the line:

$$x\left(s\right)=\frac{H}{w}\left\{ln\left[\frac{V_A+ws}{H}+\sqrt{1+{\left(\frac{V_A+ws}{H}\right)}^2}\right]-ln\left[\frac{V_A}{H}+\sqrt{1+{\left(\frac{V_A}{H}\right)}^2}\right]\right\}+\frac{Hs}{EA}$$

(3)

$$z\left(s\right)=\frac{H}{w}\left[\sqrt{1+{\left(\frac{V_A+ws}{H}\right)}^2}-\sqrt{1+{\left(\frac{V_A}{H}\right)}^2}\right]+\frac{1}{EA}\left(V_{A}s+\frac{w{s}^2}{2}\right).$$

(4)

The first term in each equation comes from the effect of the weight of the line, while the second term in each equation is due to the elasticity of the line. A more compact version of (3) is sometimes written using the identity ${sinh}^{-1}a=ln\left(a+\sqrt{1+a^2}\right)$.

The term $V_A+ws$ is the vertical tension component, which comes from a balance of forces in the vertical direction:

$$V\left(s\right)=V_A+ws.$$

(5)

The vertical tension component at end B is $V_B=V_A+wL$. As noted earlier, the horizontal tension component remains constant. From (5) and (2), the tension magnitude is then

$$\tau \left(s\right)=\sqrt{H^2+{(V_A+ws)}^2}.$$

(6)

When applying Equations (3)–(6) to model the shape and tension of a mooring line segment, the end positions of the segment are usually known, while the tension components H and $V_A$ are typically unknown. Equations (3) and (4) need to be solved numerically in order to determine these unknowns. This can be accomplished by setting $s=L$, which results in (3) and (4) being equal to the horizontal and vertical spans between the segment end points ($x_B=x\left(L\right)$ and $z_B=z\left(L\right)$, respectively). As described by Jonkman [16], the equations can then be numerically solved with an iterative Newton–Raphson method to find the tension components (Figure 2). Denoting the known end point spacings by $x_b^0$ and $z_b^0$, the iteration is

$$\left\{\begin{array}{c}{H}_{n+1}\\ V_{B,n+1}\end{array}\right\}=\left\{\begin{array}{c}{H}_n\\ V_{B,n}\end{array}\right\}-{\mathbf{J}}_n^{-1}\left\{\begin{array}{c}{x}_{b,n}-x_b^0\\ z_{b,n}-z_b^0\end{array}\right\},$$

(7)

where n is the iteration number and ${\mathbf{J}}_n$ is the Jacobian of the positions with respect to the tension components evaluated with the tension component values at iteration n:

$${\mathbf{J}}_n=\left[\begin{array}{cc}\frac{\partial x_B}{\partial H}(H_n,V_{B,n})& \frac{\partial x_B}{\partial V_B}(H_n,V_{B,n})\\ \frac{\partial z_B}{\partial H}(H_n,V_{B,n})& \frac{\partial z_B}{\partial V_B}(H_n,V_{B,n})\end{array}\right].$$

(8)

The entries for the inverse of this matrix are the partial derivatives of (3) and (4), which are straightforward to express analytically, as has been done by [16] and others. After H and $V_B$ are solved for, the distributed line profile and tensions can be calculated directly from (3)–(6).

Thus far, the presented equations apply only to a mooring line with a pure catenary shape. To support other cases, the following subsections detail the quasi-static modeling of a mooring line section with different degrees of spacing, seabed contact, and tension. These cases apply the catenary approach introduced above along with simpler models for when mooring lines are slack or lying along the seabed. In combination, these methods fully describe the mooring section profile and tension distribution across the range of likely scenarios that could be encountered.

2.2. Fully Suspended Sections

The simplest case is a mooring line section that does not contact the seabed at any part of its length. As such, the segment has a catenary shape over its entire length and (3)–(6) apply to its entirety. Figure 3 shows three scenarios that fit this case.

Substituting $s=L$ and $V_B=V_A+wL$ into (3) and (4) provides the equations for the catenary horizontal and vertical dimensions as a function of the horizontal and vertical tension at end B:

$$\begin{array}{cc}\hfill x_B(H,V_B)& =\frac{H}{w}\left\{ln\left[\frac{V_B}{H}+\sqrt{1+{\left(\frac{V_B}{H}\right)}^2}\right]-ln\left[\frac{V_B-wL}{H}+\sqrt{1+{\left(\frac{V_B-wL}{H}\right)}^2}\right]\right\}+\frac{HL}{EA}\hfill \\ \hfill & =\frac{H}{w}\left\{{sinh}^{-1}\left[\frac{V_B}{H}\right]-{sinh}^{-1}\left[\frac{V_B-wL}{H}\right]\right\}+\frac{HL}{EA},\hfill \end{array}$$

(9)

$$z_B(H,V_B)=\frac{H}{w}\left[\sqrt{1+{\left(\frac{V_B}{H}\right)}^2}-\sqrt{1+{\left(\frac{V_B-wL}{H}\right)}^2}\right]+\frac{1}{EA}\left(V_{B}L-\frac{w{L}^2}{2}\right).$$

(10)

Using these equations and their partial derivatives, the tension components at end B can be solved according to (7).

2.3. Sections Fully on the Seabed

A very simple case is when a mooring line segment lies entirely along the seabed. With the friction on the seabed neglected, the tension through the line segment is uniform:

$$\tau \left(s\right)=\mathrm{Maximum}\left(0,\left(\frac{x_B}{L}-1\right)EA\right).$$

(11)

The section profile’s z-coordinates are uniform at the seabed depth, while the x-coordinates are evenly stretched:

$$x\left(s\right)=\frac{x_B}{L}s.$$

(12)

As a result, the end tensions are $V_B=0$ and $H=\tau$, for which no iterative solution is required.

2.4. Sections Partly on the Seabed

A common scenario is when a mooring line section lies partially along the seabed (Figure 4). In this case, the two previous cases are combined using a piecewise approach. The portion of the line resting on the seabed is denoted as $L_b$ and the suspended portion is denoted as $L_s$. These unstretched lengths sum to provide the total unstretched line length L.

As with the fully suspended case, determining the line profile requires solving nonlinear profile shape equations for the end tension components. Because the end of the line that rests on the seabed will have no vertical tensions, the suspended length of the line can be expressed from (5) as

$$L_s=\frac{V_B}{w}.$$

(13)

If the vertical force at the fairlead is greater than the total weight of the line ($L_s>L$ or $V_B>Lw$), then the line will be fully suspended and have no contact with the seabed, matching the case discussed in Section 2.2. Otherwise (i.e., when $V_B<wL$), the vertical force at end B is less than the weight of the line, meaning that a length of line will rest on the seabed:

$$L_b=L-L_S=L-\frac{V_B}{w}.$$

(14)

As long as this resting length of line is less than the horizontal span of the profile ($L_b<x_B$), the line will have some horizontal tension and can be solved with the catenary equations; otherwise, the next case in Section 2.5 applies.

The equations governing the profile and tension distribution of a line that is partially on the seabed consist of a combination of those for a suspended line and a line on the seabed. For the portion on the seabed, the equations are identical to those in Section 2.3. For the suspended portion, the equations are based on those of Section 2.1 with $V_A=0$, and most instances of s are replaced with $s-L_b$ to account for the portion of line on the seabed. These substitutions result in the following equations for the line profile:

$$x\left(s\right)=\left\{\begin{array}{cc}s\left(1+\frac{H}{EA}\right)\hfill & 0\le s\le L_b\hfill \\ L_b+\frac{H}{w}ln\left[\frac{w(s-L_b)}{H}+\sqrt{1+{\left(\frac{w(s-L_b)}{H}\right)}^2}\right]+\frac{Hs}{EA}\hfill & L_b<s\le L\hfill \end{array}\right.$$

(15)

$$z\left(s\right)=\left\{\begin{array}{cc}0\hfill & 0\le s\le L_b\hfill \\ \frac{H}{w}ln\left[\sqrt{1+{\left(\frac{w(s-L_b)}{H}\right)}^2}-1\right]+\frac{w{(s-L_b)}^2}{2EA}\hfill & L_b<s\le L\hfill \end{array}\right..$$

(16)

The tensions along the line are

$$\tau \left(s\right)=\left\{\begin{array}{cc}\mathrm{Maximum}(0,H)\hfill & 0\le s\le L_b\hfill \\ \sqrt{H^2+{\left(w(s-L_b)\right)}^2}\hfill & L_b<s\le L.\hfill \end{array}\right.$$

(17)

Similar to Section 2.2, the tension components H and $V_B$ can found by a numerical solution of the equations for the known end coordinates ($x_B$ and $z_B$). These equations are obtained by substituting L for s in (15) and (16):

$$x_B(H,V_B)=L_b+\frac{H}{w}ln\left[\frac{V_B}{H}+\sqrt{1+{\left(\frac{V_B}{H}\right)}^2}\right]+\frac{HL}{EA},$$

(18)

$$z_B(H,V_B)=\frac{H}{w}\left[\sqrt{1+{\left(\frac{V_B}{H}\right)}^2}-1\right]+\frac{V_B^2}{2EAw}.$$

(19)

As in in Section 2.2, these equations can be solved using (7).

2.5. Sections Partly on the Seabed and Slack

If the line section from the previous case is sufficiently slack that the length lying on the seabed equals or exceeds the horizontal spacing of the endpoints ($L_b\ge x_B$), then there will be no horizontal tension in the line. In this case, the suspended portion of the line will hang vertically from end B (Figure 5) and the catenary equations do not apply.

The tension where the vertically hanging portion contacts the seabed will be zero, and the tension will increase up the length of the suspended portion according to

$$\tau \left(s\right)=w(s-L_b).$$

(20)

The strain along the vertical line can then be expressed as

$$\frac{dz}{ds}=1+\frac{T\left(s\right)}{EA}=1+\frac{w(s-L_b)}{EA},$$

(21)

which can be integrated to provide the z-coordinate of a given point along the unstretched length of the suspended portion:

$$z\left(s\right)=(s-L_b)+\frac{w{(s-L_b)}^2}{2EA},s\ge L_b.$$

(22)

Setting $s=L$ provides an equality that relates the height of end B off the seabed ($z_B$) to the suspended length $L_h=L-L_b$. Solving this using a quadratic formula provides the following equation for the suspended length of the line:

$$L_h=L-L_b=\frac{EA}{w}\left(\sqrt{1+\frac{2w{z}_B}{EA}}-1\right).$$

(23)

In a typical application, this equation can be used immediately to determine how much of the line is suspended versus how much is laying on the seabed, at which point (22) and (20) can be used determine the positions and tensions of the suspended portion. The tension at end B is

$$V_B=w{L}_b=EA\left(\sqrt{1+\frac{2w{z}_B}{EA}}-1\right).$$

(24)

No iterative solution process is required in this case. Along the horizontal portion of the line that sits on the seabed, the tensions are zero and the node positions are assumed to be evenly spaced.

2.6. Vertical Sections

If there is no horizontal separation between the line section end points, then the catenary equations do not apply; instead, the section behavior needs to be modeled purely in the vertical direction. A key calculation for this case is the stretched length of an elastic line section hanging freely from one end. Similar to (22), this stretched length, or height, of a line section of length L is

$$h=L+\frac{1}{2EA}{L}^2.$$

(25)

A vertical line attached at both ends can be either taut or slack (Figure 6). Applying (25) to the full unstretched length provides $h_L$, a measure of the vertical extent of the line section if it had zero tension at the lower end. If this distance is greater than the actual distance between the end points ($h_L>z_B$), then the line is slack; otherwise, it is taut.

If the line section is slack, then it will have two portions that hang vertically from the end points and come together at the lowest point of the profile, which is the point of zero tension. Labeling these two portions A and B, by applying (25) to each of them and then considering that $L_A+L_B=L$ and $h_B-h_A=z_B$, the unstretched length of one portion can be solved for as follows:

$$L_B=\frac{z_B+L+\frac{W{L}^2}{2EA}}{2+\frac{WL}{EA}}.$$

(26)

Applying the same for the other side, the tension on each end point can then be calculated from $V_A=w{L}_A$ and $V_B=w{L}_B$. Equations (22) and (20) can be adapted to compute the vertical location and tension at any point along either portion’s unstretched length.

If the segment is taut ($h_L<z_B$), then a uniform additional strain term ${ϵ}_u$ can be calculated; this represents how much the segment is stretched beyond a free-hanging state due to the end point locations:

$${ϵ}_u=\frac{L_B-h_L}{L}.$$

(27)

This strain is superimposed with the strain implied by (25) to provide the total strain throughout the line. The corresponding tension is

$$\tau \left(s\right)={ϵ}_{u}EA+ws.$$

(28)

2.7. Sections on the Seabed with Both Ends Suspended

A special case arises if a line segment running between two points that are off the seabed hangs in such a way that it contacts the seabed. This situation does not lend itself to any of the solving techniques already discussed, as it will have two catenary portions, one at each side. Thus, determining the suspended lengths of these catenaries and the forces at the line ends is more difficult than previous cases, as the system of equations is larger. This scenario is solved using an iterative approach that performs two separate catenary solutions for each side of the line.

This special case is detected when the lowest point of a suspended line section is found to be below the seabed. The lowest point is located where the vertical component of the line tension crosses zero ($s_{low}=L-V_B/w$), which is substituted into the equation for the z-coordinate (10) to obtain

$$z_{low}=z\left(s=L-\frac{V_B}{w}\right)=\left[1-\sqrt{1+\frac{{(V_B-wL)}^2}{H^2}}\right]\frac{H}{w}-\frac{{(V_B-wL)}^2}{2EAw}.$$

(29)

When this theoretical low point $z_{low}$ is below the seabed, the mooring line is split into two portions at this location along its unstretched length $s_{low}$ (Figure 7). Each portion undergoes its own catenary solution, following the approach described in Section 2.4, while an additional outer-loop iteration is applied to couple the two line portions. In order for the two line portion solutions to represent the full U-shaped line profile, the horizontal tension component in each portion must be equal. The outer-loop iteration achieves this by adjusting the x-location of the split point on the seabed. A Newton–Raphson iteration method is used to make adjustments until the line tensions are equal:

$$x_{split,n+1}=x_{split,n}+\frac{H_B-H_A}{k_{x,A}+k_{x,B}}$$

(30)

where $k_x$ is the stiffness in the x-direction from either portion. Figure 8 illustrates this process.

If the line is fully slack ($H=0$), then the above process is not required; instead, the equations of Section 2.5 can be applied to both hanging ends of the line.

In either case, the solutions of the line profile coordinates and tensions from each side of the U-shaped line are combined, along with appropriate coordinate transformations, to provide the final coordinate and tension results for the U-shaped line.

2.8. Generalization to Other Cases

The above methods cover the range of anticipated scenarios for a uniform mooring line section with axial elasticity over an even seabed with no friction or external applied loads. However, they assume that the section is negatively buoyant and that end B is higher than end A. A few simple adjustments can be made to handle cases in which these assumptions are reversed.

In the case of a positively buoyant line, the fully suspended line calculations of Section 2.2 are applied with the buoyancy reversed. Then, the results in terms of vertical forces and z-coordinates are flipped. Similarly, in cases where a mooring line is specified with end A higher than end B, the scenario can be transformed to match the existing assumptions by swapping the end coordinates and reversing the sign of $z_B$.

All of the equations and methods described above have been incorporated into MoorPy in its catenary solver module. From there, they form the basis for how MoorPy computes the profile and forces of any given mooring line section.

In practice, any mooring line section connecting two points in three dimensions can be looked at from a two-dimensional view that labels the higher end as end B and the lower end as end A. The specifics of switching between two-dimensional and three-dimensional views and how multiple mooring line sections can be connected are discussed in the next section.

3. Mooring System Assembly and Solution

A general quasi-static mooring system model requires the ability to connect individual mooring line sections together in various arrangements in order to represent different types of mooring systems. Additional methods are then required to solve for the equilibrium state of this interconnected system of mooring line sections.

3.1. Mooring System Object Hierarchy

MoorPy represents a mooring system using a framework of three distinct object types, namely, lines, points, and bodies, which are also used in MoorDyn:

• A line is a single mooring line section with uniform distributed properties, as described in Section 2.
• A point is an entity with three translational degrees of freedom (DOFs); it can have weight and buoyancy properties, and serves as the attachment mechanism at the ends of lines. Any number of line ends can be attached to a point, while exactly one point must be assigned to every line end.
• A body is a representation of a rigid body that can both translate and rotate. It can have weight, buoyancy, and hydrostatic properties such that it can represent a floating platform. Points can be attached to a body, which then allows any lines attached to those points to impart forces and moments on the body’s six DOFs.

Following the same convention used in MoorDyn, the bodies and points in MoorPy can be fixed, coupled, or free. Fixed objects cannot move, although they can be tied to another object and move with it. Coupled objects are expected to be manipulated by an external program, for example to represent the motions of a floating platform. Free objects have degrees of freedom that are solved for within the model.

MoorPy uses a hierarchy in its objects when assigning positions or computing mooring system reaction forces (Figure 9). Bodies are at the highest level, followed by points, then lines. Positions are assigned in a top-down manner: first, the body positions are set, which then determines the positions of any attached points, which then determines the positions of any attached line ends. Points not attached to a body have their positions set separately. These are constant for a fixed point (such as an anchor), specified externally for a coupled point, or solved for within the model for a free point. When all point positions are known, these define the line end coordinates, allowing the line profiles and forces to be computed. Then, the forces are computed in a bottom-up manner: the line end forces are summed to compute the net forces on the points, then the body forces and moments are summed from those points attached to each body.

Computing the system stiffness matrix analytically is a key part of this formulation’s efficiency; it avoids the computational expense of finite difference calculations for every degree of freedom, providing a more robust and accurate calculation. When using an analytic approach, the mooring system stiffness is computed from the bottom up in the same way as for the forces in Figure 9. First, the stiffness matrix of each mooring line is computed based on the Jacobians that have already been used to solving each mooring line’s profile. Next, the stiffness matrix of each point is computed based on those of the attached lines. Finally, the stiffness matrix of each body is computed based on those of the attached points. How this process works is covered in the next subsections.

3.2. Line Forces and Stiffnesses

Section 2 has described the solution method for computing the profile and distributed tension of any mooring line section in two dimensions. The same formulation can be extended to three dimensions by a rotational transformation about the Z-axis (Figure 10). The mooring line’s heading $\theta$ is defined here as the direction from end A to end B as measured as a rotation about the Z-axis relative to the direction of the global X-axis. In the following, the local horizontal and vertical coordinates in the catenary plane discussed in Section 2 are denoted with a prime as $x^{\prime }$ and $y^{\prime }$ to distinguish them from the global coordinates.

A heading rotation matrix $\mathbf{R}$ can be defined as the direction cosines matrix that would rotate a line’s orientation from along the global X-axis to along its true heading (along $x^{\prime }$):

$$\mathbf{R}=\left[\begin{array}{ccc}cos\left(\theta \right)& -sin\left(\theta \right)& 0\\ sin\left(\theta \right)& cos\left(\theta \right)& 0\\ 0& 0& 1\end{array}\right].$$

(31)

Based on a line section’s end coordinates in three-dimensional space (${\mathbf{r}}_A$ and ${\mathbf{r}}_B$), the horizontal and vertical distances for the two-dimensional catenary calculations of Section 2 are

$$x_B^{\prime }=\sqrt{{(r_{Bx}-r_{Ax})}^2+{(r_{By}-r_{Ay})}^2}$$

(32)

and

$$z_B^{\prime }=r_{Bz}-r_{Az},$$

(33)

where nomenclature such as $r_{Az}$ denotes the z-component of the ${\mathbf{r}}_A$ position vector. After the mooring line section profiles are solved in two dimensions ($x^{\prime }\left(s\right)$ and $z^{\prime }\left(s\right)$), they can be converted into three-dimensional profiles ($\mathbf{r}\left(s\right)$) with

$$\mathbf{r}\left(s\right)=\left[\begin{array}{c}{x}^{\prime }\left(s\right)cos\left(\theta \right)+r_{Ax}\\ x^{\prime }\left(s\right)sin\left(\theta \right)+r_{Ay}\\ z^{\prime }\left(s\right)+r_{Az}\end{array}\right].$$

(34)

Similarly, the tension force vectors at the line ends in three dimensions are

$${\mathbf{f}}_A=\left[\begin{array}{c}Hcos\left(\theta \right)\\ Hsin\left(\theta \right)\\ V_A\end{array}\right]$$

(35)

and

$${\mathbf{f}}_B=\left[\begin{array}{c}-Hcos\left(\theta \right)\\ -Hsin\left(\theta \right)\\ -V_B\end{array}\right].$$

(36)

The mooring line section’s stiffness matrix needs to be known in order to compute the overall system stiffness matrix. For mooring line sections that are suspended at both ends, each end has a stiffness matrix. The matrix can be computed using a finite difference approach by perturbing each line end position and measuring how the reaction forces change. However, a more efficient approach is to use the Jacobian computed analytically from the equations in Section 2. These Jacobians are already used for solving each line section profile, as in (7). Inverting the matrix provides the corresponding two-dimensional stiffness matrix for each mooring line end:

$${\mathbf{K}}_A^{\left(2D\right)}={\mathbf{J}}_A^{-1}=\left[\begin{array}{cc}\frac{\partial H}{\partial x_A^{\prime }}& \frac{\partial H}{\partial z_A^{\prime }}\\ \frac{\partial V_A}{\partial x_A^{\prime }}& \frac{\partial V_A}{\partial z_A^{\prime }}\end{array}\right]$$

(37)

and

$${\mathbf{K}}_B^{\left(2D\right)}={\mathbf{J}}_B^{-1}=\left[\begin{array}{cc}\frac{\partial H}{\partial x_B^{\prime }}& \frac{\partial H}{\partial z_B^{\prime }}\\ \frac{\partial V_B}{\partial x_B^{\prime }}& \frac{\partial V_B}{\partial z_B^{\prime }}\end{array}\right].$$

(38)

It should be noted that the off-diagonal terms in the Jacobians often play a large role, meaning that $\frac{\partial H}{\partial x_B^{\prime }}$ can be very different from ${\left(\frac{\partial x_B^{\prime }}{\partial H}\right)}^{-1}$.

While the above matrices only consider motion in the plane of the catenary, a mooring line will also apply a reaction force to out-of-plane displacements. As discussed in [21], the out-of-plane (or transverse) stiffness term $k_t$ can be computed from the horizontal tension and span as

$$k_t=\frac{H}{x_B^{\prime }}.$$

(39)

This fraction is undefined if the mooring line is vertical, in which case the same fraction can instead be made with vertical components. The full $3\times 3$ matrix for a mooring line end is

$${\mathbf{K}}_B^{\prime }=\left[\begin{array}{ccc}\frac{\partial H}{\partial x_B^{\prime }}& 0& \frac{\partial H}{\partial z_B^{\prime }}\\ 0& k_t& 0\\ \frac{\partial V_B}{\partial x_B^{\prime }}& 0& \frac{\partial V_B}{\partial z_B^{\prime }}\end{array}\right].$$

(40)

As this matrix is relative to the line’s local heading, the next step is to transform it into the global reference frame using the heading rotation matrix from (31):

$${\mathbf{K}}_B=\mathbf{R}{\mathbf{K}}_B^{\prime }{\mathbf{R}}^{\top }$$

(41)

where ${\mathbf{R}}^{\top }$ denotes the transpose of the rotation matrix. The result is the stiffness matrix in global coordinates, which can be used in equations such as $\Delta {\mathbf{f}}_B={\mathbf{K}}_B\Delta {\mathbf{r}}_B$.

In MoorPy, the approach of (40) and (41) is used to make three stiffness matrices for every line, that is, for end A, end B, and the coupling between ends A and B:

$$\begin{array}{cc}\hfill {\mathbf{K}}_A& =\frac{\partial {\mathbf{f}}_A}{\partial {\mathbf{r}}_A}\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill {\mathbf{K}}_B& =\frac{\partial {\mathbf{f}}_B}{\partial {\mathbf{r}}_B}\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill {\mathbf{K}}_{BA}& =\frac{\partial {\mathbf{f}}_B}{\partial {\mathbf{r}}_A}.\hfill \end{array}$$

(44)

When a line is fully suspended, ${\mathbf{K}}_A={\mathbf{K}}_B$ and ${\mathbf{K}}_{AB}=-{\mathbf{K}}_A$; however, when a line contacts the seabed, the matrices differ, as the reaction force is affected by a change in how much line weight is supported by the seabed. These differences are captured when deriving the Jacobians for each case in Section 2. This means that the formulation presented thus far can provide the three-dimensional profiles, tensions, and stiffness properties of any mooring line section within the stated assumptions.

3.3. Point Forces and Stiffnesses

When evaluating multiple mooring line sections, point objects provide the means of connection. Points signify the end points for line objects, allowing line ends to be fixed, moved, or connected together. Points can also have their own mass and volume properties, meaning that they can add weight or buoyancy forces to the mooring system. The forces on a point are then the sum of its own weight and buoyancy, plus the summed tension vectors from any attached line ends ${\mathbf{f}}_{linei}$, plus any specified external applied force ${\mathbf{f}}_{ext}$:

$${\mathbf{f}}_{point}=\left[\begin{array}{c}0\\ 0\\ (\rho v_{point}-m_{point})g\end{array}\right]+\sum _i{\mathbf{f}}_{linei}+{\mathbf{f}}_{ext},$$

(45)

where $\rho$ is the water density, g is the acceleration due to gravity, v is the volume, and m is the mass.

The stiffness of a point is the sum of the stiffnesses of any attached mooring line ends, from (41):

$${\mathbf{K}}_{point}=\sum _i{\mathbf{K}}_{linei}.$$

(46)

The summation approach means that points can support any number of attached mooring lines.

3.4. Body Forces and Stiffnesses

Rotational moments and stiffness effects are provided by body objects. Similar to points, bodies serve to sum up forces and stiffnesses from what is attached to them. However, bodies work with attached point objects, which can be at offset locations to allow for moments; in addition, bodies can rotate.

A body’s orientation is described by a rotation matrix $\mathbf{R}\left(\mathsf{\theta }\right)$, where $\mathsf{\theta }$ is a vector of the roll, pitch, and yaw angles. MoorPy uses a z-y-x intrinsic Tait–Bryan rotation sequence, meaning that a body’s orientation comes from first applying a yaw rotation, then a pitch rotation (about the rotated pitch axis), then a roll rotation (about the twice-rotated roll axis):

$$\mathbf{R}\left(\mathsf{\theta }\right)=\left[\begin{array}{ccc}cos{\theta }_{1}cos{\theta }_2& cos{\theta }_{1}sin{\theta }_{2}sin{\theta }_3-cos{\theta }_{3}sin{\theta }_1& sin{\theta }_{1}sin{\theta }_3+cos{\theta }_{1}cos{\theta }_{3}sin{\theta }_2\\ cos{\theta }_{2}s1& cos{\theta }_{1}cos{\theta }_3+sin{\theta }_{1}sin{\theta }_{2}sin{\theta }_3& cos{\theta }_{3}s1sin{\theta }_2-cos{\theta }_{1}sin{\theta }_3\\ -sin{\theta }_2& cos{\theta }_{2}sin{\theta }_3& cos{\theta }_{2}cos{\theta }_3\end{array}\right].$$

(47)

The three moments corresponding to the rotation directions of $\mathsf{\theta }$ are denoted by the vector $\mathsf{\mu }$. A combined force and moment vector is denoted as ${\mathbf{f}}^{\left(6\right)}$:

$${\mathbf{f}}^{\left(6\right)}=\left\{\begin{array}{c}\mathbf{f}\\ \mathsf{\mu }\end{array}\right\}.$$

(48)

A 6-DOF pose vector combining positions and rotations is denoted as ${\mathbf{r}}^{\left(6\right)}$:

$${\mathbf{r}}^{\left(6\right)}=\left\{\begin{array}{c}\mathbf{r}\\ \mathsf{\theta }\end{array}\right\}.$$

(49)

The translation, force, and moment components in these expressions are always assumed to be aligned with the global reference frame ($x,y,z$), while the orientation angles follow the Tait–Bryan rotation sequence mentioned earlier.

The relative position of an attached point i relative to the body’s reference frame is denoted as ${\mathbf{b}}_i^{\prime }$. When rotated into the global reference frame, This relative position vector is

$${\mathbf{b}}_i=\mathbf{R}\left(\mathsf{\theta }\right){\mathbf{b}}_i^{\prime }.$$

(50)

The net force on a point will also apply a moment on the body of ${\mathsf{\mu }}_{pointi}={\mathbf{b}}_i\times {\mathbf{f}}_{pointi}$, resulting in a size-6 load vector of

$${\mathbf{f}}_{pointi}^{\left(6\right)}=\left\{\begin{array}{c}{\mathbf{f}}_{pointi}\\ {\mathsf{\mu }}_{pointi}\end{array}\right\}=\left\{\begin{array}{c}{\mathbf{f}}_{pointi}\\ {\mathbf{b}}_i\times {\mathbf{f}}_{pointi}\end{array}\right\},$$

(51)

which is the cross-product of the relative position of the point from the body reference point and the total force on the point (which may include any number of attached lines from Equation (45)).

The total force on a body comes from the net forces on the attached points (which are in turn attached to lines) along with any intrinsic body forces. Bodies can have mass, buoyancy, and hydrostatic properties, defined by a mass m centered at a position ${\mathbf{r}}_{CG}^{\prime }$ in the body reference frame, a volume v assumed to act at a metacenter ${\mathbf{r}}_M^{\prime }$, and a waterplane area A with a centroid that is assumed to lie along the z-axis of the body. These hydrostatic properties are a simple implementation included for convenience, and are not the focus of the model. Combining the factors listed above results in the following expression for the total force and moment on a body:

$${\mathbf{f}}_{body}^{\left(6\right)}=\left\{\begin{array}{c}\left\{\begin{array}{c}0\\ 0\\ -mg\end{array}\right\}+\left\{\begin{array}{c}0\\ 0\\ \rho gv\end{array}\right\}+\left\{\begin{array}{c}0\\ 0\\ \rho g{A}_{WP}\Delta z\end{array}\right\}\\ \left(\mathbf{R}\left(\mathsf{\theta }\right){\mathbf{r}}_{CG}^{\prime }\right)\times \left\{\begin{array}{c}0\\ 0\\ -mg\end{array}\right\}+\left(\mathbf{R}\left(\mathsf{\theta }\right){\mathbf{r}}_M^{\prime }\right)\times \left\{\begin{array}{c}0\\ 0\\ \rho gv\end{array}\right\}\end{array}\right\}+\sum _i\left\{\begin{array}{c}{\mathbf{f}}_{pointi}\\ \mathbf{R}\left(\mathsf{\theta }\right){\mathbf{b}}_i^{\prime }\times {\mathbf{f}}_{pointi}\end{array}\right\}+{\mathbf{f}}_{ext}^{\left(6\right)},$$

(52)

where $\Delta z$ is the body’s vertical displacement and ${\mathbf{f}}_{ext}^{\left(6\right)}$ is any externally applied force and moment (e.g., the thrust force and overturning moment from a wind turbine).

The body’s stiffness matrix is more complicated than the other stiffness matrices, as it must account for rotations and moments. Its derivation makes use of the concept of the alternator matrix, denoted as $\mathbf{H}$, which is defined as follows:

$$\mathbf{H}\left(\mathbf{b}\right)=\left[\begin{array}{ccc}0& b_3& -b_2\\ -b_3& 0& b_1\\ b_2& -b_1& 0\end{array}\right].$$

(53)

The alternator matrix has two roles. First, it provides the partial derivative of the position of a point attached to a body with respect to the rotation angles of the body:

$$\frac{\partial \mathbf{b}}{\partial \mathsf{\theta }}=\mathbf{H}\left(\mathbf{b}\right).$$

(54)

Second, it provides a matrix that can be multiplied by a vector to perform the same operation as the cross product, for example,

$$\mathbf{b}\times \mathbf{f}=-\mathbf{H}\left(\mathbf{b}\right)\mathbf{f}.$$

(55)

The body’s $6\times 6$ stiffness matrix is based on a summation of the body’s hydrostatic stiffness along with the stiffness contribution from each point attached to the body and any mooring lines attached to those points:

$${\mathbf{K}}_{body}^{\left(6\right)}=\left[\begin{array}{cc}\left[\begin{array}{ccc}0& & \\ & 0& \\ & & \rho g{A}_{WP}\end{array}\right]& \mathbf{0}\\ \mathbf{0}& \mathbf{H}\left(\begin{array}{c}0\\ 0\\ -mg\end{array}\right)\mathbf{H}\left({\mathbf{r}}_{\mathbf{CG}}\right)+\mathbf{H}\left(\begin{array}{c}0\\ 0\\ \rho gV\end{array}\right)\mathbf{H}\left({\mathbf{r}}_{\mathbf{M}}\right)\end{array}\right]+\sum _i{\mathbf{K}}_{pointi}^{\left(6\right)}$$

(56)

where the $6\times 6$ stiffness contribution from each attached point is

$${\mathbf{K}}_{pointi}^{\left(6\right)}=-\frac{\partial {\mathbf{f}}_{pointi}^{\left(6\right)}}{\partial {\mathbf{r}}_{body}^{\left(6\right)}}=\left[\begin{array}{cc}-\frac{\partial {\mathbf{f}}_{pointi}}{\partial {\mathbf{r}}_{body}}& -\frac{\partial {\mathbf{f}}_{pointi}}{\partial {\mathsf{\theta }}_{body}}\\ -\frac{\partial {\mathsf{\mu }}_{pointi}}{\partial {\mathbf{r}}_{body}}& -\frac{\partial {\mathsf{\mu }}_{pointi}}{\partial {\mathsf{\theta }}_{body}}\end{array}\right].$$

(57)

This represents the stiffness contribution of point i as four $3\times 3$ submatrices. The upper left-hand submatrix is the same as the stiffness matrix of the attached point (or line end (40) if only one line is attached) because it is the Jacobian of the force exerted by the point/line on the body with respect to the attachment point translation:

$$-\frac{\partial {\mathbf{f}}_{pointi}}{\partial {\mathbf{r}}_{body}}={\mathbf{K}}_{pointi}.$$

(58)

The upper right-hand submatrix represents the change in the mooring line force with respect to the body rotation, which depends on the relative location of the point on the body (${\mathbf{b}}_i={\mathbf{r}}_{pointi}-{\mathbf{r}}_{body}$):

$$-\frac{\partial {\mathbf{f}}_{pointi}}{\partial {\mathsf{\theta }}_{body}}=-\frac{\partial {\mathbf{f}}_{pointi}}{\partial {\mathbf{r}}_{pointi}}\frac{\partial {\mathbf{r}}_{pointi}}{\partial {\mathsf{\theta }}_{body}}={\mathbf{K}}_{pointi}\mathbf{H}\left({\mathbf{b}}_i\right).$$

(59)

The bottom left-hand submatrix represents how the moment from the mooring line changes with body translation:

$$-\frac{\partial {\mathsf{\mu }}_{pointi}}{\partial {\mathbf{r}}_{body}}=-{\mathbf{b}}_i\times \frac{\partial {\mathbf{f}}_{pointi}}{\partial {\mathbf{r}}_{body}}=-\mathbf{H}\left({\mathbf{b}}_i\right){\mathbf{K}}_{pointi}=\mathbf{H}{\left({\mathbf{b}}_i\right)}^{\top }{\mathbf{K}}_{pointi}.$$

(60)

The bottom right-hand submatrix represents the reaction moments to body rotations. It is the least straightforward to calculate, as it involves a partial derivative of the cross-product from (51):

$$-\frac{\partial {\mathsf{\mu }}_{pointi}}{\partial {\mathsf{\theta }}_{body}}=-\frac{\partial }{\partial {\mathsf{\theta }}_{body}}({\mathbf{b}}_i\times {\mathbf{f}}_{pointi}).$$

(61)

Both sides of the cross-product are functions of $\mathsf{\theta }$. The expression can be evaluated using higher-order tensor math; however, a more intuitive explanation comes from considering one rotation at a time. Replacing $\mathsf{\theta }$ with ${\theta }_j$ allows straightforward use of the product rule when differentiating the cross-product:

$$\frac{\partial \mathsf{\mu }}{\partial {\theta }_j}=\frac{\partial }{\partial {\theta }_j}(\mathbf{b}\times \mathbf{f})=\frac{\partial \mathbf{b}}{\partial {\theta }_j}\times \mathbf{f}+\mathbf{b}\times \frac{\partial \mathbf{f}}{\partial {\theta }_j}$$

(62)

where the subscripts denoting point i have been dropped for brevity, as the entire equation pertains to point i. This equation produces column j of the $\partial \mathsf{\mu }/\partial \mathsf{\theta }$ matrix, showing how this part of the stiffness is composed of two terms: one related to the change in moment arm of the line tension force, and one related to the change in line tension force.

The first term in (62) represents the effect of the moment arm changing when the body rotates while the line tension vector remains constant. It can be rewritten as

$$\frac{\partial \mathbf{b}}{\partial {\theta }_j}\times \mathbf{f}=\mathbf{H}{\left(\mathbf{b}\right)}_j\times \mathbf{f}=-\mathbf{f}\times \mathbf{H}{\left(\mathbf{b}\right)}_j=\mathbf{H}\left(\mathbf{f}\right)\mathbf{H}{\left(\mathbf{b}\right)}_j,$$

(63)

where $\mathbf{H}{\left(\mathbf{b}\right)}_j$ represents a column of the $\mathbf{H}\left(\mathbf{b}\right)$ matrix corresponding to ${\theta }_j$. The second term in (62) represents the change in tension due to rotation, and it can be rewritten as

$$\mathbf{b}\times \frac{\partial \mathbf{f}}{\partial {\theta }_j}=-\mathbf{H}\left(\mathbf{b}\right)\frac{\partial \mathbf{f}}{\partial {\theta }_j}=-\mathbf{H}\left(\mathbf{b}\right)\frac{\partial \mathbf{f}}{\partial \mathbf{b}}\frac{\partial \mathbf{b}}{\partial {\theta }_j}=\mathbf{H}\left(\mathbf{b}\right)\mathbf{K}\mathbf{H}{\left(\mathbf{b}\right)}_j,$$

(64)

where $\mathbf{K}$ is the stiffness matrix of the point (46), recognizing that $\partial \mathbf{r}=\partial \mathbf{b}$. Adding (63) and (64) provides the final expression for (62) as $\mathbf{H}\left(\mathbf{f}\right)\mathbf{H}{\left(\mathbf{b}\right)}_j+\mathbf{H}\left(\mathbf{b}\right)\mathbf{K}\mathbf{H}{\left(\mathbf{b}\right)}_j$.

The extension to all three rotation directions to produce the $3\times 3$ submatrix of (61) simply requires removing the j index and using the full $\mathbf{H}\left(\mathbf{b}\right)$ matrix. Restoring the subscripts of point i, we have

$$-\frac{\partial {\mathsf{\mu }}_{pointi}}{\partial {\mathsf{\theta }}_{body}}=-\mathbf{H}\left({\mathbf{f}}_{pointi}\right)\mathbf{H}\left({\mathbf{b}}_i\right)-\mathbf{H}\left({\mathbf{b}}_i\right){\mathbf{K}}_{pointi}\mathbf{H}\left({\mathbf{b}}_i\right).$$

(65)

Now that expressions for the four submatrices of (57) are established, the completed matrix for the stiffness contribution from a point on a body is

$${\mathbf{K}}_{pointi}^{\left(6\right)}=\sum _i\left[\begin{array}{cc}{\mathbf{K}}_{pointi}& {\mathbf{K}}_{pointi}\mathbf{H}\left({\mathbf{b}}_i\right)\\ \mathbf{H}{\left({\mathbf{b}}_i\right)}^{\top }{\mathbf{K}}_{pointi}& -\mathbf{H}\left({\mathbf{f}}_{pointi}\right)\mathbf{H}\left({\mathbf{b}}_i\right)-\mathbf{H}\left({\mathbf{b}}_i\right){\mathbf{K}}_{pointi}\mathbf{H}\left({\mathbf{b}}_i\right)\end{array}\right].$$

(66)

Because ${\mathbf{K}}_{line}$ (which makes up ${\mathbf{K}}_{point}$) is symmetric and $\mathbf{H}$ is anti-symmetric, all of the matrix except for the bottom right-hand quadrant must be symmetric. The asymmetry in the bottom right-hand quadrant is explained in Section 4.4 of the results.

At this point, it is worth highlighting a sometimes-overlooked characteristic of the body stiffness matrix. The term $\mathbf{H}\left({\mathbf{f}}_{pointi}\right)\mathbf{H}\left({\mathbf{b}}_i\right)$, which relates to the change in moment arm of a constant tension force (63), is not an inherent stiffness characteristic to which transformations can be applied. It is dependent on external tension vectors relative to the reference point location, and is unrelated to any other terms in the stiffness matrix. As a result, a body’s mooring system stiffness matrix, once calculated, cannot be transformed to be about a different reference point. Instead, the stiffness matrix needs to be recomputed whenever a new reference point is used.

The final stiffness matrix for the body, incorporating all of the above, is

$${\mathbf{K}}_{body}^{\left(6\right)}=\left[\begin{array}{cc}\left[\begin{array}{ccc}0& & \\ & 0& \\ & & \rho g{A}_{WP}\end{array}\right]& \mathbf{0}\\ \mathbf{0}& \mathbf{H}\left(\begin{array}{c}0\\ 0\\ -mg\end{array}\right)\mathbf{H}\left({\mathbf{r}}_{\mathbf{CG}}\right)+\mathbf{H}\left(\begin{array}{c}0\\ 0\\ \rho gV\end{array}\right)\mathbf{H}\left({\mathbf{r}}_{\mathbf{M}}\right)\end{array}\right]+\sum _i{\left[\begin{array}{cc}\mathbf{K}& \mathbf{K}\mathbf{H}\left(\mathbf{b}\right)\\ \mathbf{H}{\left(\mathbf{b}\right)}^{\top }\mathbf{K}& -\mathbf{H}\left(\mathbf{f}\right)\mathbf{H}\left(\mathbf{b}\right)-\mathbf{H}\left(\mathbf{b}\right)\mathbf{K}\mathbf{H}\left(\mathbf{b}\right)\end{array}\right]}_{pointi}.$$

(67)

3.5. Cross-Coupling Stiffness Terms

Forming a full mooring system stiffness matrix for a system that has more than one free object requires combining the individual point or body stiffness matrices with cross-coupling terms that represent how displacements of one object will generate forces on another object. These cross-coupling terms only need to be calculated when a mooring line runs between one free object and another free object.

For any mooring line, a cross-coupling matrix ${\mathbf{K}}_{BA}$ is defined to represent the reaction forces on end B of the line due to motions of end A. To consider the opposite direction (forces on end A due to motions of end B), the matrix is simply transposed:

$${\mathbf{K}}_{BA}={\mathbf{K}}_{AB}^{\top }.$$

(68)

These $3\times 3$ matrices are used directly when representing the coupled stiffness effects between two points connected by a mooring line.

When a mooring line runs between two bodies, there is a $6\times 6$ cross-coupling matrix between the bodies, of the form

$${\mathbf{K}}_{1,2}=\left[\begin{array}{cc}{\mathbf{K}}_{AB}& {\mathbf{K}}_{AB}\mathbf{H}\left({\mathbf{r}}_2\right)\\ \mathbf{H}{\left({\mathbf{r}}_1\right)}^{\top }{\mathbf{K}}_{AB}& \mathbf{H}\left({\mathbf{r}}_2\right){\mathbf{K}}_{AB}\mathbf{H}{\left({\mathbf{r}}_1\right)}^{\top }\end{array}\right].$$

(69)

This matrix represents the loads on body 1 due to motions of body 2 caused by a mooring line between the bodies when end A is attached to body 1 and end B is attached to body 2. This matrix is similar to (66), except that it omits the term that represents a body rotation causing a change in the moment arm. For coupling between bodies, this term would represent the rotation of one body causing a change in the moment arm of the other body, which is a higher-order effect that is negligible except when the bodies are very close together.

When a mooring line runs between a body and a point, we have a $6\times 3$ cross-coupling matrix

$${\mathbf{K}}_{1,2}=\left[\begin{array}{c}{\mathbf{K}}_{AB}\\ \mathbf{H}{\left({\mathbf{r}}_1\right)}^{\top }{\mathbf{K}}_{AB}\end{array}\right].$$

(70)

Following (68), when dealing with mooring lines in the reverse orientation, ${\mathbf{K}}_{AB}$ should be transposed. Thus, to represent the coupling effects of body 1 on body 2 (rather than the other way around), the entire ${\mathbf{K}}_{1,2}$ matrix should be transposed. At this point, all the components of the full system stiffness matrix have now been defined.

3.6. System Stiffness Matrix and Equilibrium Solution

Beyond the individual line segment computations discussed in Section 2, the main challenge in creating a quasi-static mooring system model is solving for the full system’s equilibrium state, meaning the equilibrium values of all the free object’s degrees of freedom. If $\mathbf{r}$ is the vector of system DOFs and $\mathbf{f}$ is the vector of net forces in each DOF, then solving for equilibrium means determining the vector $\mathbf{r}$ that will satisfy

$$\mathbf{f}\left(\mathbf{r}\right)=\mathbf{0}.$$

(71)

Because the individual mooring line force–displacement relations are nonlinear, solving (71) requires an iterative solution process over the individual mooring lines’ iterative solutions. Similar to many other models, MoorPy uses a Newton–Raphson approach. Iteration is as follows:

$${\mathbf{r}}_{n+1}={\mathbf{r}}_n-\mathbf{K}{\left({\mathbf{r}}_n\right)}^{-1}\mathbf{f}\left({\mathbf{r}}_n\right)$$

(72)

where n is the iteration number and $\mathbf{K}$ is the full system stiffness matrix at the current position. This process involves computing the individual line profiles and the full system stiffness matrix at each iteration (Figure 11).

MoorPy is general in that it supports solving the equilibrium for multiple bodies and any combination of attachments in the mooring system. The vector $\mathbf{r}$ includes any enabled DOFs. DOFs can be disabled to improve efficiency; for example, DOFs in the y-direction can be disabled to perform a two-dimensional analysis.

The simplest way to calculate the system stiffness matrix is via finite differencing. Each system DOF j can be perturbed in turn to compute column j of the stiffness matrix, denoted as $\mathbf{K}{\left(\mathbf{r}\right)}_j$. The following is an example using a forward difference:

$$\mathbf{K}{\left(\mathbf{r}\right)}_j=\frac{\mathbf{f}\left({\mathbf{r}}^{j+}\right)-\mathbf{f}\left(\mathbf{r}\right)}{d{r}_j}$$

(73)

where $\mathbf{r}$ is the unperturbed DOF vector, ${\mathbf{r}}^{j+}$ is the DOF vector with DOF j perturbed, and $d{r}_j$ is the size of the perturbation.

A more efficient approach is to use the analytic stiffness matrix calculations detailed in the previous subsections to assemble the full system stiffness matrix analytically. The full system stiffness matrix consists of the stiffness matrix of each object along with the cross-coupling matrices between each pair of objects (appearing twice as transposes of each other in the upper and lower halves of the matrix), as shown in Figure 12.

When considering the computational expense of these calculations, a given calculation of the system forces $\mathbf{f}\left(\mathbf{r}\right)$ can be considered roughly proportional in expense to the number of DOFs in the system m. From (73), populating the stiffness matrix using forward differencing requires $m+1$ evaluations, resulting in a total computational expense of order $m^2$ for constructing the system stiffness matrix.

In comparison, the analytic approach requires the system forces and corresponding Jacobians of each object to be evaluated only once, meaning that the system stiffness matrix computation is of order m when using the analytic approach.

Considering that the stiffness matrix must be inverted, there is a clear nonlinear increase in computational complexity with system size. However, the analytic approach for the system Jacobians is one order lower in terms of computational expense than the finite difference approach. This is a significant difference, particularly for more complex mooring systems.

3.7. Stiffness of Coupled DOFs with Free DOFs

There is often a need to find the stiffness matrix for certain DOFs (such as those attached to a floating structure) while leaving other DOFs (such as those of the connection points along the mooring lines) free to equilibrate. MoorPy allows coupled and free DOFs to be distinguished for this purpose.

There are two ways to find the stiffness of only coupled DOFs (Figure 13). The simplest approach is to use a finite difference method, as mentioned in the previous subsection, where each coupled DOF is perturbed while the rest of the mooring system is allowed to settle into equilibrium at each perturbation.

An analytic approach is also possible, and can be explained by considering the meaning of a stiffness matrix and its inverse. In the equation for a stiffness matrix ($\Delta \mathbf{f}=\mathbf{K}\Delta \mathbf{r}$), any zero entry in the $\Delta \mathbf{r}$ vector means that the displacement in that DOF is zero; therefore, the corresponding column in $\mathbf{K}$ has no bearing on the reaction forces. In the inverse of the equation ($\Delta \mathbf{x}={\mathbf{K}}^{-1}\Delta \mathbf{f}$), any zero entry in $\Delta \mathbf{f}$ means that there is zero net force in that degree of freedom; thus, the corresponding column of ${\mathbf{K}}^{-1}$ has no bearing on the displacements. Free DOFs in equilibrium will have zero net force, and we do not track their motions when computing the stiffness matrix of coupled DOFs. Therefore, the rows and columns corresponding to free DOFs in $\Delta \mathbf{r}={\mathbf{K}}^{-1}\mathbf{f}$ will have no bearing on the stiffness matrix of coupled DOFs. This leads to the following process for computing the stiffness matrix for coupled DOFs:

• Perform a system equilibrium solution for only the free DOFs.
• Compute the inverse of the full system stiffness matrix (including both free and coupled DOFs).
• Delete the rows and columns of free DOFs from the inverse of the matrix, leaving a matrix for only the coupled DOFs.
• Invert the reduced inverse stiffness matrix to obtain the stiffness matrix for coupled DOFs with other DOFs in equilibrium.

The above process is for cases in which free DOFs are assumed to equilibrate. Alternatively, if the free DOFs are assumed to be frozen, meaning that they do not shift to stay in equilibrium, then the coupled stiffness matrix can be obtained by simply removing the rows and columns of the free DOFs from the full stiffness matrix. While this does not require any matrix inversion, it is not usually a very useful matrix.

The above approach for finding the stiffness matrix of coupled DOFs while the other DOFs are equilibrated complements the general analytic approach for computing stiffness matrices while avoiding the computational expense and numerical approximation of finite difference approaches. This completes the core set of methods used in MoorPy for quasi-static mooring system analysis with analytic Jacobians, providing both generality and efficiency.

4. Demonstration and Verification

This section presents nine cases that demonstrate and verify the methods presented above, beginning with the simplest scenarios representing individual mooring lines and then progressively increasing in scenario complexity. The quasi-static results were computed directly in MoorPy software based on the formulations described in the previous sections. For verification, these quasi-static results are compared with static equilibrium results from MoorDyn, a lumped-mass mooring dynamics model that has been validated and verified in previous studies [11,22]. MoorDyn’s general formulation allows it to simulate both simple and highly intricate mooring configurations with equal effectiveness, making it well suited to testing the validity of MoorPy’s implementation of different types of mooring line profiles and interconnections. The implementation of individual mooring lines in Section 2 was originally verified against results from other quasi-static models; however, it is the more complicated capabilities of MoorPy that are novel and challenging. Verifying those capabilities is the focus here, for which a general model such as MoorDyn allows comparisons to be made for every case.

4.1. Overview of the Verification Cases

Figure 14 illustrates the nine verification cases. Case 1 tests basic catenary behavior. Case 2 tests the combination of a catenary section and a taut section connected together. Cases 3–6 test the handling of scenarios with clump weights and floats, buoyant sections (as occur in dynamic power cables), slack mooring lines, and midline seabed contact, all of which are situations that may arise in unconventional mooring systems (e.g., for certain wave energy converters). Cases 7–9 test the mooring system stiffness effects on floating bodies, including bridle configurations and shared mooring connections, which are applicable for floating wind turbines and wave energy converters.

Table 1. Mooring system design parameters for the cases.

	Case	1	2	3	4	5	6	7 (a, b)	8	9
Depth		300	300	300	300	300	300	50	100	100
Anchor	$r_x$ (m)	−800	−800	−800	−800	−400	−400	100	200	−141.4
	$r_y$ (m)	0	0	0	0	0	0	0	0	−141.4
	$r_z$ (m)	−300	−300	−300	−300	−300	0	−50	−100	−100
Section 1	Type	Reg.	Reg.	Reg.	Reg.	Reg.	Reg.	Reg.	Rope	Reg.
	L (m)	900	500	400	360	890	660	112	160	210
Point 1	$r_x$ (m)		−400	−400	−300		0		40
	$r_y$ (m)		0	0	0		0		0
	$r_z$ (m)		−100	−100	−200		−100		−50
	m (kg)		0	0	0		0		0
	v (m3)		0	200	0		200		0
Section 2	Type		Rope	Reg.	Buoy.		Reg.		Rope
	L (m)		350	250	240		330		45.36
Point 2	$r_x$ (m)			−200	−100		0
	$r_y$ (m)			0	0		200
	$r_z$ (m)			−100	−200		−200
	m (kg)			100,000	0		0
	v (m3)			0	0		0
Section 3	Type			Reg.	Reg.		Reg.
	L (m)			250	360		330
Fairlead	$r_x$ (m)	0	0	0	0	0	0	0, 5	3.94	−14.14
	$r_y$ (m)	0	0	0	0	0	400	0, 3	±6.82	−14.14
	$r_z$ (m)	0	0	0	0	0	0	−10	−21	−20

provides the dimensions and quantitative properties of the test cases. To reduce complexity, only three types of cross-sectional properties were used in these cases:

• “Regular” (Reg.) has high weight and stiffness to represent chains that are commonly used for mooring systems, as well as a dynamic power cable in Case 4.
• “Rope” is nearly neutrally buoyant and has lower stiffness to represent fiber rope mooring line sections used in semi-taut and taut mooring configurations.
• “Buoyant” (Buoy.) has the same weight as Regular, but with much greater diameter to represent the buoyancy section of a dynamic power cable.

The diameter, linear density, and elasticity of these section types are provided in

Table 2. Mooring line section types and properties.

Name	Regular	Rope	Buoyant
Diameter (m)	0.2	0.15	0.7
Linear density (kg/m)	500	25	500
Stiffness, EA (MN)	2000	30	2000

. More qualitative descriptions of each test case are provided when the test case results are introduced in the following subsections. The input files for the test cases are available in the open-source MoorPy repository.

4.2. Numerical Modeling Settings

In the results that follow, MoorPy is assigned a convergence tolerance of 1 mm for its equilibrium solution process. Meanwhile, each MoorDyn result comes from a time domain simulation of at least 600 s with elevated damping coefficients to accelerate convergence to static equilibrium. Unless otherwise noted, each mooring line segment in MoorDyn is subdivided into at least 20 internal segments.

Cases 7–9 include results of stiffness matrices to provide further insight and verification on how MoorPy solves system equilibrium. When computing these stiffness matrices using the finite difference approach, the perturbation sizes in MoorPy are 0.02 m and 0.001 radians, whie those in MoorDyn are 0.2 m and 0.01 radians. Smaller values are not used in MoorDyn, as the lumped mass discretization is vulnerable to small bumps when nodes contact or leave the seabed, an effect that can be magnified when using very small step sizes. MoorDyn is assigned 100 s to settle into equilibrium after each perturbation, and the equilibrium result is averaged from the final 20 s of that period. The accuracy obtained with this method is similar to that of the level of discretization used in the MoorDyn models. The comparisons generally use the MoorDyn results to ensure that the MoorPy results are close to the correct value (i.e., that there are no bugs), after which a comparison of the MoorPy finite difference and analytical results can provide a more precise measure of the accuracy of the results, as MoorPy is not limited by a lumped-mass discretization.

4.3. Cases 1–6: Individual Mooring Line Assemblies

The first six cases consider individual mooring lines, i.e., an assembly of mooring line sections running between an anchor point and an attachment point, which would be on a floating platform. Such an assembly involves line objects connected and terminated by point objects. If the fairlead end (the end that would be attached to a floating platform) is assumed to be stationary, then the system has $3n$ degrees of freedom, where n is the number of connecting points along the mooring line assembly (i.e., there would be $n+1$ sections). These cases focus on translational effects only.

Case 1 is a basic catenary mooring line with a portion laying along the seabed. This is the most generic mooring line configuration and the easiest case to model. Figure 15 shows a comparison with the equilibrated MoorDyn results for this case. The continuous quasi-static results for the line section’s profile, distributed horizontal and vertical positions, and tension along its arc length are shown in green. The results from MoorDyn are plotted as black points, corresponding to the MoorDyn model’s discretization of the line section; the profile and positions are plotted with points at the node locations, while the tensions are computed at each segment between nodes, making them offset from the position measurements. The plots show very close agreement in all results, confirming MoorPy’s profile and tension prediction accuracy for a single catenary mooring line section.

Case 2 is a two-segment mooring line that corresponds to a catenary chain at the seabed and a rope section running up to the floating platform. This approximates semi-taut mooring configurations commonly used for floating wind turbines. The same type of results as shown for the last case are plotted in Figure 16. The comparison is again very close, confirming the accuracy of MoorPy with catenary and taut sections, including an equilibrium solution of their connection point.

Case 3 comprises three chain segments connected by a float and a clump weight; this combination is sometimes used in wave energy converters to provide compliance, as it permits a high degree of device motion. This case demonstrates finding equilibrium with multiple catenary segments and the presences of point forces. The results are plotted in Figure 17. As seen in the last subplot, the presence of the clump weight and float cause large jumps in tension when transitioning from one section to the next. The agreement between MoorPy and MoorDyn remains very close, demonstrating that the point weight and buoyancy forces are modeled correctly.

Case 4 represents a dynamic power cable in the common “lazy wave” configuration, with a catenary section rising from the seabed, a buoyant section (representing the length-averaged sectional properties of a dynamic cable section fitted with evenly spaced buoyancy modules) shaped in an inverse catenary, and then another non-buoyant catenary section rising up to the platform. The results are shown in Figure 18. The last subplot shows how the peak tensions occur at the transition points between sections, where the abrupt change in net weight/buoyancy causes a sharp reversal in the tension slope. The results agree closely between models, demonstrating that these phenomena are modeled well.

Case 5 is a fully slack mooring line, with the suspended portion hanging vertically and the portion on the seabed having no tension. This represents an extreme case that is often not supported in catenary models but that is useful to include when simulating unconventional designs or mooring line failures in nonredundant mooring systems, as these situations that can be of high relevance for wave energy or floating wind turbine systems. The results are shown in Figure 19. While the profiles are identical in most cases, there is a large difference in how MoorPy and MoorDyn position the nodes lengthwise along the seabed when the tension is zero, meaning that there is not enough force to move the nodes. In MoorDyn, those positions are causal (dependent on the initial conditions), while in MoorPy they are assumed to be evenly spaced. This difference does not effect any important quantities for this level of modeling fidelity, while the profile shape and tension are matched well between models.

Case 6 features a mooring line assembly with two different portions making contact with the seabed and a float halfway between them keeping the midpoint suspended in the middle of the water column. There is also a massless point partway along the second portion that lies on the seabed, which tests the ability for different seabed contact scenarios. This case tests more unusual scenarios that could occur, for example if there is a failure in a shared mooring system. Additionally, this mooring system is assigned a 45 degree heading, meaning that it has both x- and y-axis components to test support for different headings. In both MoorPy and MoorDyn, the system is initialized with the buoyant midpoint located at 0,0 (far from its equilibrium location); thus, solving for the equilibrium of the floating point involves motion in all three dimensions. The results are shown in Figure 20. Despite the complexity of this case, all the results agree closely, demonstrating MoorPy’s ability to handle each of these issues.

To provide some quantification of the agreement levels between MoorPy and MoorDyn,

Table 3. Root mean square differences between MoorPy and MoorDyn results in Cases 1–6.

Case	1	2	3	4	5	6
Profile positions (m)	0.062	0.020	0.023	0.062	43.4	2.67
Distributed tension (kN)	1.57	0.09	2.82	0.77	4.02	18.6

presents the root mean square (RMS) differences in the node positions and tensions predicted by the two models. These RMS values are all quite small relative to the distance and tension magnitudes of the case. The close agreement suggests that the MoorPy results are accurate, and that even closer agreement could be obtained by using finer discretization levels in MoorDyn. Taken together, these results provide confidence in the accuracy of the quasi-static modeling method for individual section profiles, point forces, and translational degrees of freedom.

4.4. Case 7: Mooring Line Attached to a Body

The rotational effects that arise from mooring lines attached at offset positions on a body are complicated (Section 3.4), and call for careful verification. To do this, Case 7 features a single catenary mooring line where the fairlead is attached to a body at an offset position. Two different positions are used: (a) an offset in z, and (b) an offset in all three dimensions. This provides visibility on the complex 6-DOF stiffness effects that arise with just one direction of offset before verifying the stiffness calculations when the fairlead is offset in three dimensions. The unstretched mooring line length differs between the two cases, providing roughly the same tension despite the two different fairlead locations. Both cases are shown in Figure 21, illustrating the fairlead location relative to the body reference frame.

Figure 22 shows the stiffness matrices calculated for the two fairlead positions using three different methods: finite differencing of the MoorDyn results, finite differencing of the MoorPy results, and MoorPy’s analytic stiffness calculation method. These stiffness matrix tables use darker and bolder text to highlight the most significant matrix values; near-zero values can generally be ignored. The background shading provides a reminder of the three different types of stiffness coefficients, namely, the translational stiffness, rotational stiffness, and translational–rotational coupling stiffness terms.

Comparing the results of the two variations in fairlead position illustrates which stiffness terms are dominant from an offset in a single direction and those that are added by other offset directions. In Case 7a, which lies in the x–z plane, the stiffness matrix has significant values in all the entries that do not involve out-of-plane DOFs (sway, roll, and yaw). In addition, there are some significant terms involving these out-of-plane DOFs, such as the roll stiffness and the coupling term between sway and roll. In Case 7b, with its three-dimensional fairlead position, none of the 36 entries in the stiffness matrix are negligible. This illustrates the complexity of three-dimensional mooring system stiffness effects.

In both variations of Case 7, the z-coordinate has the largest effect when it comes to rotations and moments. This can be used to explain the nature of asymmetry in the lower right-hand quadrant of the matrix. A rotation about the z-axis will have little effect on the moment about x, while a rotation about x will have a large effect on the moment about z, as the x-rotation (combined with the z-offset of the fairlead) causes a large change in the fairlead’s y-coordinate. Accordingly, $|K_{46}|<|K_{64}|$. This demonstrates why stiffness matrix asymmetry occurs, something that was not understood in the discussion of [19], where the method produced symmetric stiffness matrices.

To provide a visual comparison of the results from the three stiffness calculation methods, Figure 23 plots selected stiffness coefficients predicted by MoorPy normalized by the results from MoorDyn. In general, these results show small differences from MoorDyn, on the order of 2%, especially for cross-coupling stiffness terms. The close alignment of the circles and crosses shows that the MoorPy analytic and finite difference results agree closely. Closer agreement is possible by using smaller perturbation sizes, further confirming that both MoorPy methods provide accurate results. The MoorDyn results could be increased in accuracy by using a finer lumped-mass discretization, which would bring all the results further into alignment. The decent agreement between MoorPy and MoorDyn for all terms suggests that the MoorPy formulations are correct, and the close agreement between the MoorPy finite difference and analytic calculations further supports the accuracy of the stiffness calculations.

4.5. Case 8: Floating Wind Turbine Bridle Mooring System

Case 8 considers the stiffness effects of a full mooring system on a floating platform. The combined stiffness of multiple mooring lines attached to a floating platform is simply a summation of the stiffness provided by each line; thus, this case adds an additional complexity in the form of bridle connections. As typical for spar–buoy floating wind turbine support structures, this mooring system has three main mooring lines, each dividing into two lines before attaching to the platform, as shown in Figure 24. This bridle configuration increases yaw stiffness and results in a more complex mooring system behavior. The width of the bridle attachments is such that the adjacent mooring lines share the same attachment points, spaced at 120 degrees and a radius of 7.88 m around the body centerline. Their depth is 21 m.

This system has fifteen degrees of freedom: six for the body (which can be considered coupled DOFs) and three for each of the three bridle points (which can be considered free DOFs). While the full $15\times 15$ stiffness matrix is not shown for the sake of space, the correct stiffness behavior can be tested by comparing the results for the $6\times 6$ coupled stiffness matrix of the floating body, as accurately computing this rigid-body stiffness depends on accounting for the couplings with the three bridle points, and through them the entire system of nine mooring line sections.

Figure 25 presents the coupled body stiffness matrix computed using four different methods:

• MoorDyn finite differencing, where each body DOF is purturbed and then the rest of the system is allowed to settle into equilibrium.
• MoorPy analytic with free DOFs frozen, where only body DOFs are varied and the free DOFs of the bridle points are assumed to not move.
• MoorPy finite differencing, where each body DOF is perturbed while the bridle points are equilibrated.
• MoorPy analytic, where the full system stiffness matrix is used to directly find the coupled stiffness using the analytic technique described in Section 3.7.

The method with the bridle point DOFs frozen is valuable for demonstrating the significance of the couplings between all the DOFs of the mooring system. The predicted body stiffness coefficients are drastically higher than from the other methods, showing how critical it is to properly account for the coupling effects with free DOFs in the system. The other three methods all provide results that are in good agreement with each other, verifying the coupled stiffness calculation methods described in Section 3.7. All the dominant coefficients from these three other methods agree to within 1%, though some smaller terms have greater disagreement. The heave–roll coupling stiffness terms deviate by a large factor; however, these terms are three orders of magnitude smaller than the other terms in their quadrant, indicating that they are unimportant. In fact, for a mooring system that is supposed to be axisymmetric, the presence of these terms is simply a result of loss of precision from entering the positions numerically in the input file.

Figure 26 provides a visual comparison of the important stiffness terms calculated by each of the methods, normalizing by the MoorDyn results. The figure provides two different y-axis magnifications to first show the massive difference when the free DOFs are frozen and then to zoom in on the small differences between the three other methods. As with the previous case, differences on the order of 1% are considered to be due to numerical approximations from finite differencing and the MoorDyn discretization, both of which could be improved at the cost of computation time. The presence of some zero-valued normalized results correspond to DOF pairs for which the MoorPy methods accurately predicted zero stiffness while the numerical approximations in the MoorDyn method caused a nonzero stiffness result. Overall, these results show the effectiveness of the coupled stiffness calculation methods of Section 3.7, and demonstrate that the analytic method provides the best results at much less computational expense.

4.6. Case 9: Two Floating Bodies and a Shared Mooring Line

The remaining capability to be demonstrated is the case where multiple floating bodies are coupled by mooring lines. A good test is provided by the simplest shared mooring system arrangement: two floating bodies, each with two anchored mooring lines, and one shared mooring line running between the bodies (Figure 27).

In this example, the bodies are spaced 200 m apart (center-to-center) along the x-axis. The mooring lines are attached to the bodies at radii of 20 m and depths of 20 m. The water depth is 100 m; the anchors are positioned at a radius of 200 m from the platforms’ centerlines and spread at an angle of 90 degrees. The rationale for this geometry comes from [21]. As a best practice, the shared line length is carefully selected so that the system is in equilibrium at the specified spacing, rather than allowing the platforms to drift into unintended equilibrium positions. The anchored mooring lines have unstretched lengths of 210 m, while the shared line has an unstretched length of 168.2 m.

The degrees of freedom for this case are the six rigid-body degrees of freedom for each of the two bodies. As such, the correct behavior of the coupling between bodies can be checked by looking at the full $12\times 12$ system stiffness matrix. Figure 28 presents the system stiffness matrices calculated using the same three methods as in Case 7. There is no difference in a stiffness calculation method whether dealing with one body or two. Finite difference methods simply pass through and perturb each degree of freedom one-by-one, and the reaction forces and moments are subjected to finite differencing in the same way to form the stiffness matrix.

Comparing the results shows decent agreement between all three methods and very good agreement between the two quasi-static methods. Figure 29 plots the normalized results of key entries in the stiffness matrix to compare between methods. The upper pane shows entries for the stiffness on one of the bodies, while the lower pane shows entries of the cross-coupling terms between the two bodies. The level of agreement is similar in the upper and lower panes, indicating that the cross-coupling stiffness calculations (including the analytic method described in Section 3.5) works as well as the other calculations explored in previous cases. The MoorPy results differ from the MoorDyn results by up to 6% in the worst case (for the small coupling term between heave and pitch); however, the agreement between finite difference and analytic MoorPy results is with 0.3%, and often much closer. As in previous cases, these trends indicate that the quasi-static calculation methods are correct, and that agreement with the MoorDyn results could be improved asymptotically by using finer discretization in the MoorDyn modeling.

An important general observation from the stiffness matrices in Figure 28 is that the presence of a single shared mooring line between two bodies creates very many cross-coupling stiffness terms. As is visible from the highlighting in the tables, the off-diagonal $6\times 6$ matrices indicating the coupling terms between bodies have sixteen non-negligible entries. That these entries are captured equally well by the MoorPy finite difference and analytic methods demonstrates that the most difficult parts of the analytic stiffness method have been successfully formulated and implemented.

Recalling from Section 3.6 that the finite difference method requires about m times more computations than the analytic method (where m is the number of degrees of freedom), the ability of the analytic method to capture all stiffness terms is very useful, allowing for twelve times faster computation in the present example and even greater speed gains for more complicated mooring systems.

5. Conclusions

In this paper, a general approach for quasi-static mooring system modeling has been formulated, implemented, and verified. The well-established catenary equations for the profile and tension of a mooring line section have been augmented to include additional scenarios, such as when a mooring line section is fully slack, fully along the seabed, or has two suspended ends with seabed contact in the middle. These additions allow for robust modeling capability supporting the many scenarios that can arise with offshore renewable energy systems.

To represent any mooring system topology, the formulation includes point objects to connect the ends of mooring line sections as well as body objects that group point objects at relative positions to account for rotations and moments on rigid floating bodies. These objects are formulated generally, meaning that any combination of attachments can be modeled.

Lastly, the formulation includes an analytic approach for computing the stiffness matrix of an entire mooring system, beginning with the stiffness of each mooring line section and building up to the combined translational and rotational stiffness terms of multiple floating bodies in a shared mooring system. The derivation of this approach provides mathematical clarity on the various ways in which mooring lines provide stiffness effects, including based on how their tension vectors change when a body moves. One particular note is that a body’s mooring system stiffness matrix, once calculated, cannot be transformed to be about a different reference point. In addition, the rotational portion of a stiffness matrix can be asymmetric. The analytic stiffness approach is more computationally efficient than finite differencing approaches by a factor of the number of degrees of freedom in the system, making it significantly more practical for analyzing interconnected arrays of floating wind turbines or wave energy converters.

We have presented nine test cases to demonstrate and verify the quasi-static modeling formulation, using equilibrated MoorDyn results as a basis for comparison. The general trend across all results is that the quasi-static results computed by MoorPy agree with the results from MoorDyn, verifying the quasi-static formulation. More specific conclusions are as follows:

• Mooring line profiles and tensions agree very closely between MoorPy and MoorDyn, within 1% in most configurations, verifying the quasi-static line section model (Section 2) and the handling of point objects and equilibrium solution (Section 3.2, Section 3.3 and Section 3.6).
• Mooring system stiffness matrices typically agree to around 1–2% between MoorPy and MoorDyn, and agree very closely between the MoorPy finite difference and analytic methods, verifying the quasi-static stiffness formulations for bodies and coupled systems (Section 3.4, Section 3.5 and Section 3.6), including cases where certain DOFs are allowed to equilibrate (Section 3.7).
• Small differences in the results between MoorPy and MoorDyn are attributable to MoorDyn’s lumped-mass discretization and the numerical approximations inherent in computing finite differences from equilibrium results of a time domain model. All of these factors could be asymptotically improved at the cost of computation time.
• The many off-diagonal stiffness terms arising in the tests with rigid bodies are captured similarly in the various methods, indicating that no non-negligible terms are missing from the analytic stiffness computation method, even when there are couplings between bodies or when certain DOFs are allowed to equilibrate.

Apart from the general verification provided by the results, perhaps the most important observation is that even relatively simple mooring configurations can have many off-diagonal stiffness terms when there are rigid bodies in a three-dimensional scenario. The presented quasi-static formulation can resolve these terms without the errors that arise when using a finite difference approach and with significantly greater computational efficiency.

The presented formulation is publicly available in the open-source MoorPy library, along with many other functions not covered in this paper that can aid in the efficient analysis of mooring systems.