# Numerical Stability and Accuracy of Temporally Coupled Multi-Physics Modules in Wind Turbine CAE Tools∗ Amir Gasmi,† Michael A. Sprague,‡ Jason M. Jonkman,§ and Wesley B. Jones¶ National Renewable Energy Laboratory, Golden, Colorado, 80401, USA In this paper we examine the stability and accuracy of numerical algorithms for coupling time-dependent multi-physics modules relevant to computer-aided engineering (CAE) of wind turbines. This work is motivated by an in-progress major revision of FAST, the National Renewable Energy Laboratory’s (NREL’s) premier aero-elastic CAE simulation tool. We employ two simple examples as test systems, while algorithm descriptions are kept general. Coupled-system governing equations are framed in monolithic and partitioned representations as differential-algebraic equations. Explicit and implicit loose partition coupling is examined. In explicit coupling, partitions are advanced in time from known information. In implicit coupling, there is dependence on other-partition data at the next time step; coupling is accomplished through a predictor-corrector (PC) approach. Numerical time integration of coupled ordinary-differential equations (ODEs) is accomplished with one of three, fourth-order fixed-time-increment methods: Runge-Kutta (RK), Adams-Bashforth (AB), and Adams-Bashforth-Moulton (ABM). Through numerical experiments it is shown that explicit coupling can be dramatically less stable and less accurate than simulations performed with the monolithic system. However, PC implicit coupling restored stability and fourth-order accuracy for ABM; only second-order accuracy was achieved with RK integration. For systems without constraints, explicit time integration with AB and explicit loose coupling exhibited desired accuracy and stability. # I. Introduction In this paper we examine the numerical stability and accuracy of several temporal coupling methods for multi-physics modules relevant to computer-aided engineering (CAE) of wind turbines. The wind turbine system is composed of many physics that have significantly different characteristic time and length scales. Figure 1 shows an illustration of a floating wind turbine with some possible physical interactions. The wind turbine industry relies heavily on CAE tools for analyzing wind turbine performance, loading, and stability. Over the past two decades, the U.S. Department of Energy has sponsored the National Renewable Energy Laboratory’s (NREL’s) development of CAE tools for wind turbine analysis. NREL’s premier suite of tools is FAST,1 which is a modular assembly of advanced CAE tools. FAST is an open-source, professionalgrade software package. As shown in Figure 2, FAST encompasses modules for aerodynamics (AeroDyn $2,3$ ), platform hydrodynamics (HydroDyn $4,5$ ) for offshore systems, control and electrical systems (servo), and structural dynamics. The modules are coupled to allow for nonlinear analysis of aero-hydro-servo-elastic interactions in the time domain. The FAST tool enables the analysis of a range of wind turbine configurations, including two- or three-blade horizontal-axis rotors, pitch or stall regulation, rigid or teetering hub, upwind or downwind rotor, and lattice or tubular towers. The wind turbine can be modeled on land or offshore on fixed-bottom or floating substructures. ![](images/9ebb259cf9beeeee4c9999aafec5594b229b6bb795f1b96c75a096099a4d7eef.jpg) Figure 1. Illustration showing possible physical interactions experienced by an offshore floating wind turbine. ![](images/637bd13efa7ff97e36d5a388dbd16cdd71d37d7f76a64e217f34b5a3f6d55b63.jpg) Figure 2. Schematic illustrating primary simulation modules of NREL’s wind turbine CAE tool FAST. This work is motivated by a major restructuring of the FAST tool suite, which is described in a companion paper.6 FAST-restructuring goals include (1) improving the ability to read, implement, and maintain source code; (2) increasing module sharing and shared-code development across the wind community; (3) improving numerical performance and robustness; and (4) greatly enhancing flexibility and expandability to enable further developments of functionality without the need to recode established modules. It is envisioned that the new modularization framework will transform FAST into a powerful, robust, and flexible wind turbine modeling tool with a large number of developers and a range of modeling fidelities across the aerodynamic, hydrodynamic, servo-dynamic, and structural-dynamic components. We note that the proposed modularization framework has some features similar to those of other model-coupling frameworks (See, e.g., Refs. $_{7-10}$ ). In general, when modeling wind turbine multi-physics, each physics is represented as a system of nonlinear, time-dependent equations. These equations may be some combination of partial-differential equations (PDEs), ordinary-differential equations (ODEs), or, more generally, differential-algebraic equations (DAEs). Here, we restrict our scope to situations where spatial differential operators have been discretized and only temporal differential operators and/or algebraic constraints remain. This is known as the method-of-lines approach in PDE analysis. In the established version of FAST, the wind turbine system is partitioned $^{11}$ into its natural domains, and each partition is treated with a particular software module, e.g., aerodynamics are modeled with the package AeroDyn. The time-dependent ODEs or DAEs associated with each module are time integrated with various integrators, e.g., AeroDyn is integrated with a multi-step predictor-corrector algorithm (fourth-order Adams-Bashforth-Moulton). The modules interact in FAST in an explicit loosely coupled configuration, where partition solutions are updated from time $t$ to time $t+\Delta t$ using coupling information known at time $t$ , where $\Delta t$ is the time increment of the numerical integrator. Our taxonomy of models for a time-dependent multi-physics system is summarized schematically in Figure 3. On the left is a multi-physics system that can be modeled in either a monolithic or a partitioned approach. In a monolithic approach, where physics are naturally and inherently coupled, the system is represented by a single system of differential or differential-algebraic equations, and time integration is accomplished with a single solver; all quantities are advanced simultaneously. In a partitioned approach, the system is decomposed into an assembly of subsystems, each being represented by a partition that is subsequently coupled to other subsystem partitions through partition input-output relationships. We refer readers to Felippa et al.11 for a comparison of monolithic and partitioned approaches. A partitioned system can be assembled into a tightly coupled system that is appropriate for time integration by a single integrator (likely as a system of DAEs). Alternatively, partitions can be loosely coupled, where they are time integrated in a conventional serial staggered procedure $^{12}$ (i.e., sequential integration). Loose coupling is appealing because it allows for flexibility in modeling each subsystem. This flexibility may allow for different time integrators, different time increments, and different spatial discretization. Further, it allows for the use of existing software modules such as dedicated fluid-dynamics and structural-dynamics software. However, linearization of loosely coupled partitioned systems can be problematic.6 Finally, loosely coupled partitions can be coupled either explicitly or implicitly (right side of Figure 3). Under explicit loose coupling, as described above, advancement of all partition solutions from time $t$ to $t+\Delta t$ is accomplished based on shared information at $t$ . While appealing in its simplicity, loose explicit coupling may suffer numerical instabilities. Further, explicit coupling presents an inherent challenge for modules that employ a DAE or ODE solver based on implicit numerical time integration, or where a module employs a high-order single step method (like Runge-Kutta), which performs right-hand-side (RHS) evaluations in the semi-closed interval $(t,t+\Delta t]$ . Alternatively, implicit partition coupling is formulated such that one or more partitions require other-module data in $(t,t+\Delta t]$ before the solution can be advanced; the partitions are thus tied together in a linear or nonlinear system that must be solved for time advancement. However, approximate solution to the implicitly coupled systems can be pursued through a predictor-corrector (PC) approach, where each partition is advanced from $t$ to $t+\Delta t$ one or more times using predicted/corrected values from other partitions. A PC approach can maintain the modularity and flexibility of loose coupling. ![](images/2503c29c1ba0d9f11e9586438c62c3cfcd84ddffde95cee5af83be93109f4cbc.jpg) Figure 3. Schematic illustrating our “taxonomy” for models that describe a multi-physics systems. In this paper, we examine time integration of monolithic systems and equivalent partitioned systems with explicit and implicit loose coupling. Our work is guided by a motivating question: How can we couple multi-physics modules of a partitioned system, such as those in FAST, in a manner that is accurate and numerically stable? Our application examples include a simple two-mass damped oscillator, and a damped oscillator coupled to a quasi-static nonlinear cable. We examine various partition-coupling strategies where the underlying fourth-order fixed-step time integrator is either Runge-Kutta, Adams-Bashforth, or AdamsBashforth-Moulton. # II. Model Formulation In this section we present a general framework for describing time-dependent monolithic and partitioned systems that is adopted by the new modular framework in FAST, $^6$ and we provide specifics of two simple example applications that will be used as test cases to explore the stability and accuracy of the proposed coupling methods described in Section III. # II.A. General Formulation II.A.1. Monolithic System Representation In a monolithic representation of a continuous-time-dependent physical system, there is a single system of inherently tightly coupled DAEs. Using a generic state-space representation, we can write the monolithicsystem DAEs as follows: $$ \begin{array}{r l}&{\dot{\mathbf{x}}=\mathbf{X}\left(t,\mathbf{x}(t),\mathbf{z}(t)\right),}\\ &{\mathbf{0}=\mathbf{Z}\left(t,\mathbf{x}(t),\mathbf{z}(t)\right),}\\ &{\mathbf{y}=\mathbf{Y}\left(t,\mathbf{x}(t),\mathbf{z}(t)\right),}\end{array} $$ where an overdot denotes time differentiation, $\mathbf{X}$ , $\mathbf{Z}$ , and $\mathbf{Y}$ are multi-variable vector functions, corresponding to the state, constraint, and output equations, respectively. As with the formulation described in the companion paper (which describes the new modular framework for FAST $^6$ ), we restrict our formulation to semi-explicit DAEs of index 1. While not included here, the companion paper $^6$ includes discrete-time states. The associated vector dependent variables, $\mathbf{x}$ and $\mathbf{z}$ , are the state and constraint, respectively; y is the output-vector variable. While in some systems it is natural to setup the monolithic equations in the form represented by Eqs. (1)–(3), often the complete multi-physics system involves independent spatial domains and the interaction is achieved through shared-boundary interfaces, which makes it well suited for partitioning. In the next subsection, we present the generalized framework for such multi-physics-system representation. # II.A.2. Partitioned System Representation Let us consider that the entire multi-physics system is subdivided into $N$ subsystems; each of these subsystems can be represented in a general DAE form: $$ {\dot{\bf x}}_{i}={\bf X}_{i}\left(t,{\bf x}_{i}(t),{\bf z}_{i}(t),{\bf u}_{i}(t)\right), $$ $$ {\bf0}={\bf Z}_{i}\left(t,{\bf x}_{i}(t),{\bf z}_{i}(t),{\bf u}_{i}(t)\right), $$ $$ {\bf y}_{i}(t)={\bf Y}_{i}\left(t,{\bf x}_{i}(t),{\bf z}_{i}(t),{\bf u}_{i}(t)\right), $$ where the subscript $i$ corresponds to the $i^{t h}$ subsystem, and $\mathbf{u}_{i}$ is a vector of inputs derived from outputs (and, in general, inputs) of coupled partitions. Here, the input to partition $i$ is determined through the additional implicit input-output relationship $$ {\bf0}={\bf U}_{i}\left({\bf u}_{1}(t)\,,\dots,{\bf u}_{N}(t)\,,{\bf y}_{1}(t)\,,\dots\,,{\bf y}_{N}(t)\right). $$ Equations (4)–(7) for $i\in\{1,\ldots,N\}$ constitute what we consider a partitioned-system representation. # II.B. Application Examples In this subsection we demonstrate how the above representations can be applied to two application examples. These examples will serve as our test systems for illustrating the essence of various coupling schemes. We consider first a linear two degrees-of-freedom (DOF) damped oscillator shown by Figure 4. This is a common system for examining numerical coupling algorithms (see, e.g., Ref. $^{13}$ ). A monolithic-system representation is shown in Figure 4(a), and our partitioned-system representation is shown in Figure 4(b). The partitioning was accomplished by introducing a virtual interaction node between the System 1 mass, $m_{1}$ , and spring-damper coupling, whose displacement and velocity are denoted $q_{I}(t)$ and ${\dot{q}}_{I}(t)$ , respectively. System 1 output (related to System 2 input) is displacement $q_{1}(t)$ and velocity ${\dot{q}}_{1}(t)$ of mass 1, while System 1 input (System 2 output) is the coupling force $f_{c}$ . ![](images/191021bb7e45e7de6d35c0e729a556e7c0f6d35a68ae6eb58e8da48e065ec707.jpg) Figure 4. Schematic of the two-mass damped oscillator: (a) monolithic system, and (b) partitioned system. The governing equation for displacement $q_{1}(t)$ of the point mass in System 1 is given by $$ m_{1}\ddot{q}_{1}+c_{1}\dot{q}_{1}+k_{1}q_{1}=f_{c}+f_{1}, $$ where $f_{c}$ is the coupling force due to the relative motion between Systems 1 and 2; $c_{1}$ and $k_{1}$ are the damping and stiffness coefficients, respectively, and $f_{1}$ is a prescribed external force. Following the general formulation of the partitioned-system representation, we recast the System 1 governing equations as follows, $$ \begin{array}{r l}&{\dot{\mathbf{x}}_{1}(t)=\mathbf{X}_{1}\left(t,\mathbf{x}_{1},\mathbf{u}_{1}\right),}\\ &{\mathbf{y}_{1}(t)=\mathbf{Y}_{1}\left(t,\mathbf{x}_{1}\right),\mathbf{u}_{1}\right),}\end{array} $$ where $\mathbf{x}_{1}=\left[\begin{array}{l l}{q_{1}}&{\dot{q}_{1}}\end{array}\right]^{T}$ is the state vector (a $T$ -superscript denotes a transpose), $\mathbf{u}_{1}=\left[\begin{array}{l}{f_{c}}\end{array}\right]$ is the input vector, and $\mathbf{y}_{1}=\left[\begin{array}{l l}{q_{1}}&{\dot{q}_{1}}\end{array}\right]^{T}$ is the output vector. We note that in the above formulation, we have neither a constraint equation nor an associated constraint-state vector. For this linear system, we can write the right-hand sides in common state-space form as $$ \begin{array}{r l}&{{{\bf{X}}_{1}}={{\bf{A}}_{1}}{\bf{x}}_{1}+{{\bf{B}}_{1}}{\bf{u}}_{1}+{{\bf{f}}_{1}},}\\ &{{{\bf{Y}}_{1}}={{\bf{C}}_{1}}{\bf{x}}_{1}+{{\bf{D}}_{1}}{\bf{u}}_{1},}\end{array} $$ where $\mathbf{f}_{1}=\left[\begin{array}{l l}{0}&{f_{1}/m_{1}}\end{array}\right]^{T}$ , and the state matrix $\mathbf{A}_{1}$ , input matrix $\mathbf{B}_{1}$ , output matrix $\mathbf{C}_{1}$ , and feed-through matrix $\mathbf{D}_{1}$ , are given by $$ \mathbf{A}_{1}=\left[\begin{array}{c c}{0}&{1}\\ {-\frac{k_{1}}{m_{1}}}&{-\frac{c_{1}}{m_{1}}}\end{array}\right],\qquad\mathbf{B}_{1}=\left[\begin{array}{c}{0}\\ {\frac{1}{m_{1}}}\end{array}\right],\qquad\mathbf{C}_{1}=\left[\begin{array}{c c}{1}&{0}\\ {0}&{1}\end{array}\right],\qquad\mathbf{D}_{1}=\left[\begin{array}{c}{0}\\ {0}\end{array}\right], $$ respectively. Similarly, the governing equation for displacement $q_{2}(t)$ of the point mass $m_{2}$ in System 2 is given by $$ m_{2}\ddot{q}_{2}+\left(c_{c}+c_{2}\right)\dot{q}_{2}+\left(k_{c}+k_{2}\right)q_{2}=c_{c}\dot{q}_{I}+k_{c}q_{I}+f_{2}, $$ for which a general state-space form may be written $$ \begin{array}{r}{\dot{\mathbf{x}}_{2}=\mathbf{X}_{2}\left(t,\mathbf{x}_{2},\mathbf{u}_{2}\right),}\\ {\mathbf{y}_{2}=\mathbf{Y}_{2}\left(t,\mathbf{x}_{2},\mathbf{u}_{2}\right).}\end{array} $$ Here, $c_{2}$ and $k_{2}$ are damping and stiffness coefficients, respectively, $c_{c}$ and $k_{c}$ are damping- and stiffnesscoupling coefficients, respectively, $f_{2}$ is a prescribed force, $\mathbf{x}_{2}=\left[\begin{array}{l l}{q_{2}}&{\dot{q}_{2}}\end{array}\right]^{T}\mathbf{,\;u}_{2}=\left[\begin{array}{l l}{q_{I}}&{\dot{q}_{I}}\end{array}\right]^{T}$ , and $$ \begin{array}{r l}&{{\bf{X}}_{2}={\bf{A}}_{2}{\bf{x}}_{2}+{\bf{B}}_{2}{\bf{u}}_{2}+{\bf{f}}_{2},}\\ &{{\bf{Y}}_{2}={\bf{C}}_{2}{\bf{x}}_{2}+{\bf{D}}_{2}{\bf{u}}_{2},}\end{array} $$ where $\mathbf{f}_{2}=\left[\begin{array}{l l}{0}&{f_{2}/m_{2}}\end{array}\right]^{T}$ , $$ \mathbf{A}_{2}=\left[\begin{array}{c c}{0}&{1}\\ {-\frac{k_{c}+k_{2}}{m_{2}}}&{-\frac{c_{c}+c_{2}}{m_{2}}}\end{array}\right],\quad\mathbf{B}_{2}=\left[\begin{array}{c c}{0}&{0}\\ {\frac{k_{c}}{m_{2}}}&{\frac{c_{c}}{m_{2}}}\end{array}\right],\quad\mathbf{C}_{2}=\left[\begin{array}{c c}{k_{c}}&{c_{c}}\end{array}\right],\quad\mathbf{D}_{2}=\left[\begin{array}{c c}{-k_{c}}&{-\sigma_{2}}\\ {-k_{c}}&{-\sigma_{2}}\end{array}\right], $$ Hence, $\mathbf{y}_{2}=\left[\begin{array}{c}{c}{c_{c}\left(\dot{q}_{2}-\dot{q}_{I}\right)+k_{c}\left(q_{2}-q_{I}\right)}\end{array}\right]=\left[\begin{array}{c}{f_{c}}\end{array}\right]$ . In this form, the input-output relationships governing the data dependence between Systems $^{1}$ and $2$ are ${\bf U}_{1}\left({\bf u}_{1},{\bf y}_{2}\right)={\bf0}$ and $\mathbf{U}_{2}\left(\mathbf{u}_{2},\mathbf{y}_{1}\right)=\mathbf{0}$ , respectively, where $$ {\bf U}_{1}={\bf u}_{1}-{\bf y}_{2}\,,\qquad{\bf U}_{2}={\bf u}_{2}-{\bf y}_{1}. $$ Alternatively, in the monolithic approach, the governing equations (8) and (14) can be written as a single state-space representation, e.g., $$ \begin{array}{r l}&{\dot{\mathbf{x}}_{1-2}=\mathbf{X}_{1-2}\left(t,\mathbf{x}_{1-2}(t)\right),}\\ &{\mathbf{y}_{1-2}=\mathbf{x}_{1-2}\,,}\end{array} $$ where $$ \mathbf{X}_{1-2}=\mathbf{A}\mathbf{x}_{1-2}+\mathbf{f}_{1-2}, $$ $$ \mathbf{x}_{1-2}=\left[\begin{array}{l l l l}{q_{1}}&{\dot{q}_{1}}&{q_{2}}&{\dot{q}_{2}}\end{array}\right]^{T},\;\mathbf{f}_{1-2}=\left[\begin{array}{l l l l}{0}&{f_{1}/m_{1}}&{0}&{f_{2}/m_{2}}\end{array}\right]^{T},\;\mathrm{and} $$ $$ \mathbf{A}=\left[\begin{array}{c c c c}{0}&{1}&{0}&{0}\\ {-\frac{k_{c}+k_{1}}{m_{1}}}&{-\frac{c_{c}+c_{1}}{m_{1}}}&{\frac{k_{c}}{m_{1}}}&{\frac{c_{c}}{m_{1}}}\\ {0}&{0}&{0}&{1}\\ {\frac{k_{c}}{m_{2}}}&{\frac{c_{c}}{m_{2}}}&{-\frac{k_{c}+k_{2}}{m_{2}}}&{-\frac{c_{c}+c_{2}}{m_{2}}}\end{array}\right]. $$ We note that, here, we do not have constraints or inputs. II.B.2. Linear 1-DOF Damped Oscillator Connected with a Nonlinear Quasi-Static Cable We consider here the case where System 2 in the coupled system, as described above, is replaced with System 3 (shown in Figure 5), which is a quasi-static catenary cable where the right end is fixed and the left end is coupled to the displacement of the mass in System 1. The nonlinear relationship between cable horizontal force $H(t)$ and left-end displacement $q_{3}(t)$ is given by Jonkman and Buhl $^{14}$ as $$ 0=q_{3}-D+\frac{H L}{A E}+\frac{H}{w}\left(\ln\left(\sqrt{\frac{L^{2}w^{2}}{4H^{2}}+1}+\frac{L w}{2H}\right)-\ln\left(\sqrt{\frac{L^{2}w^{2}}{4H^{2}}+1}-\frac{L w}{2H}\right)\right), $$ where $D$ is the distance between the left and right cable ends when $q_{3}=0$ , $L$ is the unstretched length of the cable, $w$ is the weight per unit length, $E$ is the modulus of elasticity, and $A$ is the cross-sectional area. In our framework, this partition may be represented as having zero state variables, $\mathbf{x}_{3}=\varnothing$ , one constraint variable, $\mathbf{z}_{3}=\left[\begin{array}{l}{H}\end{array}\right]$ , with input $\mathbf{u}_{3}={\left[\begin{array}{l}{q_{3}}\end{array}\right]}$ , and output $\mathbf{y}_{3}=\mathbf{Y}_{3}=\left[\begin{array}{l}{R}\end{array}\right]$ . Under the notation described above, the constraint equation ${\bf Z}_{3}\left(t,{\bf z}_{3},{\bf u}_{3}\right)$ is given by the right side of Eq. (25). For System 1 coupled to System 3, the input-output relationships are ${\bf U}_{1}\left({\bf u}_{1},{\bf y}_{3}\right)={\bf0}$ and $\bf U_{3}\left(u_{3},y_{1}\right)=\bf0$ , where $$ \mathbf{U}_{1}=\mathbf{u}_{1}-\mathbf{y}_{3},\quad\mathbf{U}_{3}=\mathbf{u}_{3}-\left[\begin{array}{l l}{1}&{0}\end{array}\right]\mathbf{y}_{1}. $$ ![](images/1603526a8cfa553ea725bb6fd4277ac128aed0126d88232ca3ca9ff185fdd037.jpg) Figure 5. Schematic of System 3, a nonlinear quasi-static cable. The monolithic representation of System 1 coupled to System 3 with $q_{3}=q_{1}$ is given by $$ \begin{array}{r}{\dot{{\mathbf x}}_{1}={\mathbf X}_{1-3}\left(t,{\mathbf x}_{1}(t),{\mathbf z}_{3}(t)\right)\,,}\\ {{\mathbf0}={\mathbf Z}_{1-3}\left(t,{\mathbf x}_{1}(t),{\mathbf z}_{3}(t)\right)\,,}\\ {{\mathbf y}_{1-3}={\mathbf Y}_{1-3}\,,\qquad\qquad\qquad\mathbf{\mathbf{y}}_{1}.}\end{array} $$ where $$ \begin{array}{r l}&{\mathbf{X}_{1-3}=\mathbf{A}_{1}\mathbf{x}_{1}+\left[\begin{array}{l l}{0}&{1/m_{1}}\end{array}\right]^{T}\mathbf{z}_{3}+\mathbf{f}_{1}\,,}\\ &{\mathbf{Z}_{1-3}=\left[\begin{array}{l l}{q_{1}-D+\frac{H L}{A E}+\frac{H}{w}\left(\ln\left(\sqrt{\frac{L^{2}w^{2}}{4H^{2}}+1}+\frac{L w}{2H}\right)-\ln\left(\sqrt{\frac{L^{2}w^{2}}{4H^{2}}+1}-\frac{L w}{2H}\right)\right)}\end{array}\right]\,,}\\ &{\mathbf{Y}_{1-3}=\left[\begin{array}{l l l}{q_{1}}&{\dot{q}_{1}}&{H}\end{array}\right]\,.}\end{array} $$ # III. Numerical Coupling Methods We are interested in solving initial-value problems associated with the monolithic-system representation given by Eqs. (1)–(3) and the partitioned-system representation given by Eqs. (4)–(6). While it is straight forward to time-advance the solution to the monolithic DAEs, time-advancing the partitioned-system solution requires care, as it entails the manipulation and transfer of data between partitions at each time step. In the following subsections we present in detail the different coupling strategies employed in our numerical investigation. For clarity, we refer to the mathematical model for a given subsystem as a partition, while the numerical-algorithm implementation for time-advancement of a partition solution is called a module. # III.A. Monolithic System The monolithic-system DAEs are well suited for implementation in a general DAE solver. Alternatively, in the absence of constraints, any number of first-order-ODE solvers can be employed, e.g, Runge-Kutta or Adams methods. Here, we denote the time advance from discrete-time station $t^{n}=t^{0}+n\Delta t$ to $t^{n+1}$ as $$ \mathbf{x}^{n}\,,{\dot{\mathbf{x}}}^{n}\,,\dots\,,{\dot{\mathbf{x}}}^{n-m+1}\,,\mathbf{z}^{n}\,,\mathbf{X}(t)\,,\mathbf{Z}(t)\,\xrightarrow{\mathrm{ADV}}\mathbf{x}^{n+1},\mathbf{z}^{n+1}, $$ where $t^{0}$ is the initial time and a $n$ superscript denotes the discrete-time station, where $n\in\{0,\dots,n_{m a x}\}$ . The notation on the left hand side (LHS) explicitly indicates the data accessible to the underlying integrator, and where $\mathbf{X}(t)$ and $\mathbf{Z}(t)$ express access to the continuous-time forms of state and constraint right-hand sides, respectively; ADV denotes the execution of the underlying ${\boldsymbol{r}}m$ -step DAE- or ODE-solver algorithm over a single time step. In our formulation, we do not restrict how a particular module will handle its constraint-variable update. For example, a module may simply time lag its constraints. In that case, ${\mathbf z}^{n}$ is first calculated to be consistent with the states at $t^{\pi}$ , and continuous states are then updated; $\mathbf z^{n+1}$ is calculated independently and is seen as a guess for the constraint at the next time step. Alternatively, in an implicit DAE solver, $\mathbf{z}^{n+1}$ and $\mathbf{x}^{n+1}$ are determined concurrently. For more discussion on DAE solvers, we refer interested readers to Brenan et al.15 We consider here a multi-physics system that has been modeled as $N$ coupled partitions following the representation shown in Eqs. (4)–(7). To enhance modularity, we assume that data may be communicated between modules only at discrete-time stations. Further, we restrict our analysis to lock-step time integration; all partitions are time advanced with the same time increment $\Delta t$ . Two data-transfer schemes are considered: explicit and implicit partition coupling, where the latter system is solved and advanced in a predictorcorrector approach. # III.B.1. Explicit Loose Explicit coupling and time-advancement of partitions is straight forward, as each partition is advanced based on information known at $t^{\prime\prime}$ . Our algorithm is summarized as follows for advancement of partition solutions from $t^{\prime\prime}$ to $t^{n+1}$ using an ${\boldsymbol{r}}n$ -step integrator. For each partition $i$ , where $i\in\{1,\ldots,N\}$ , we assume that we have previous-solution information required by the underlying time-integrator, as well as output data, and input data, i.e., $$ \mathbf{x}_{i}^{n}\,,{\dot{\mathbf{x}}}_{i}^{n}\,,\dots\,,{\dot{\mathbf{x}}}_{i}^{n-m+1}\,,\mathbf{z}_{i}^{n}\,,\mathbf{y}_{i}^{n}\,,\mathbf{u}_{i}^{n}. $$ In general, these input and output values (associated with $t=t^{n}$ ) are determined as a solution to the system $$ \begin{array}{r l}&{\mathbf{0}=\mathbf{U}_{1}\left(\mathbf{u}_{1}^{n},\dots,\mathbf{u}_{N}^{n},\mathbf{y}_{1}^{n},\dots,\mathbf{y}_{N}^{n}\right)\,,}\\ &{\phantom{\quad}\dots\,}\\ &{\mathbf{0}=\mathbf{U}_{N}\left(\mathbf{u}_{1}^{n},\dots,\mathbf{u}_{N}^{n},\mathbf{y}_{1}^{n},\dots,\mathbf{y}_{N}^{n}\right)\,,}\\ &{\mathbf{y}_{1}^{n}=\mathbf{Y}_{1}\left(t^{n},\mathbf{x}_{1}^{n},\mathbf{z}_{1}^{n},\mathbf{u}_{1}^{n}\right)\,,}\\ &{\phantom{\quad}\dots\,}\\ &{\mathbf{y}_{N}^{n}=\mathbf{Y}_{N}\left(t^{n},\mathbf{x}_{N}^{n},\mathbf{z}_{N}^{n},\mathbf{u}_{N}^{n}\right)\,,}\end{array} $$ which, in general, is nonlinear, and the state time derivative at $t^{\prime\iota}$ is calculated as $$ \begin{array}{r}{\dot{\mathbf{x}}_{i}^{n}=\mathbf{X}_{i}\left(t^{n},\mathbf{x}_{i}^{n},\mathbf{z}_{i}^{n},\mathbf{u}_{i}^{n}\right),}\end{array} $$ for each $i\in\{1,\ldots,N\}$ . We denote time update under explicit coupling as follows: $$ \mathbf{x}_{i}^{n}\,,\dot{\mathbf{x}}_{i}^{n}\,,\hdots\,,\dot{\mathbf{x}}_{i}^{n-m+1}\,,\mathbf{z}_{i}^{n}\,,\mathbf{u}_{i}^{n}\,,\mathbf{X}_{i}(t)\,,\mathbf{Z}_{i}(t)\xrightarrow{\mathrm{ADV}}\mathbf{x}_{i}^{n+1}\,,\mathbf{z}_{i}^{n+1}\,, $$ where, again, the notation on the left hand side (LHS) explicitly indicates the data accessible to the underlying integrator, and where $\mathbf{X}_{i}(t)$ and $\mathbf{Z}_{i}(t)$ express access to the continuous-time forms of state and constraint right-hand sides, respectively. # III.B.2. Implicit Loose Time advancement of modules in our implicit loose coupling is based on a predictor-corrector (PC) approach. The approach described here deals with the communication of data between partitions, and is independent of the underlying ODE or DAE solver employed by each partition. This approach is related to one described elsewhere, $11,13$ and can be described as follows for a two-partition system: (1) predict the partition 1 state and constraint at $t^{n+1}$ using known information; (2) substitute those predicted values into partition 2 through the coupling terms, and calculate a predicted value of the partition 2 state and constraint at $t^{n+1}$ ; (3) correct the partition 1 state by using the predicted partition 2 data. Steps (2) and (3) may then be repeated in order to improve the solution to the associated partitioned system. This approach is associated with the Gauss-Seidel iterative-solution method,16 in that Steps (2) and (3) are applied to each partition with the most updated predicted values. If the PC procedure was modified such that each partition was updated with data known from the previous implicit-solve iteration, the method would be associated with the Jacobi iterative-solution method, which offers slower convergence rates but is more suitable to distributed-memory parallel computation. The step-by-step procedure for advancing the solution from $t^{n}$ to $t^{n+1}$ is described below. We denote this algorithm $\mathrm{PC}(j_{m a x})$ , where $j$ is the correction counter and $j_{m a x}\ge1$ is the user-defined number of corrector steps to be taken over each time step. Starting known information: We assume that the following data are known at time $t^{\pi}$ : $$ \mathbf{x}_{i}^{n}\,,\dot{\mathbf{x}}_{i}^{n}\,,\hdots\,,\dot{\mathbf{x}}_{i}^{n-m+1}\,,\mathbf{z}_{i}^{n}\,,\mathbf{y}_{i}^{n}\,,\mathbf{y}_{i}^{n-1}\,,\mathbf{u}_{i}^{n}\,,\mathbf{u}_{i}^{n-1}\,, $$ where the state-derivative RHS is calculated as $$ \begin{array}{r}{\dot{\mathbf{x}}_{i}^{n}=\mathbf{X}_{i}\left(t^{n},\mathbf{x}_{i}^{n},\mathbf{z}_{i}^{n},\mathbf{u}_{i}^{n}\right),}\end{array} $$ for each $i\in\{1,\ldots,N\}$ . Step 1 (Predict): Let $j=0$ , and predict the constraint, state, and output of partition $^{1}$ at $t^{n+1}$ . (Step 1.a) Predict the input and output at $t^{n+1}$ of all partitions through linear extrapolation (over constant $\Delta t$ ): $$ \begin{array}{r}{\mathbf{u}_{i}^{n+1(j)}=2\mathbf{u}_{i}^{n}-\mathbf{u}_{i}^{n-1},}\\ {\mathbf{y}_{i}^{n+1(j)}=2\mathbf{y}_{i}^{n}-\mathbf{y}_{i}^{n-1},}\end{array} $$ for each $i\in\{1,\ldots,N\}$ , which is an $O(\Delta t^{2})$ accurate prediction. Here, a superscript in parentheses indicates the correction iteration. (Step 1.b) Solve the input-output equation for the predicted input at $t^{n+1}$ for partition 1, i.e., solve $$ {\bf0}={\bf U}_{1}\left({\bf u}_{1}^{n+1(j+1)}\,,{\bf u}_{2}^{n+1(j)}\,,\dots\,,{\bf u}_{N}^{n+1(j)}\,,{\bf y}_{1}^{n+1(j)},\dots,{\bf y}_{N}^{n+1(j)}\right), $$ for u $\mathbf u_{1}^{n+1(j+1)}$ (Step 1.c) Let $\mathbf{u}_{1}^{\alpha}\;=\;\mathbf{L}\left(\alpha_{1},\mathbf{u}_{1}^{n},\mathbf{u}_{1}^{n+1(j+1)}\right)$ be a linearly interpolated input value based on known and predicted inputs, where $$ \mathbf{L}\left(\alpha,\ \mathbf{u}_{1}^{n},\ \mathbf{u}_{1}^{n+1}\right)=(1-\alpha)\mathbf{u}_{1}^{n}+\alpha\mathbf{u}_{1}^{n+1}. $$ The value of $\alpha_{1}$ is tied to the underlying ODE/DAE integrator. When the underlying ODE or DAE integrator requires other-partition data in the interval $t^{n}\;<\;t\;\leq\;t^{n+1}$ , $\alpha_{1}~>~0$ . For example, for implicit AdamsMoulton update, or predictor-corrector Adams-Bashforth-Moulton update, $\alpha_{1}\,=\,1$ . However, for explicit Adams-Bashforth, $\alpha_{1}=0$ (Step 1.d) Advance the solution to yield predicted state and constraint values, i.e., $$ {\bf x}_{1}^{n}\,,\dot{\bf x}_{1}^{n}\,,\dots\,,\dot{\bf x}_{1}^{n-m+1}\,,{\bf z}_{1}^{n}\,,{\bf u}_{1}^{\alpha}\,,{\bf X}_{1}(t)\,,{\bf Z}_{1}(t)\,\xrightarrow{\mathrm{ADV}}{\bf x}_{1}^{n+1(j+1)}\,,{\bf z}_{1}^{n+1(j+1)}\,, $$ while $\mathbf{u}_{1}^{n+1(j+1)}$ is used to calculated the output, i.e., $$ \mathbf{y}_{1}^{n+1(j+1)}=\mathbf{Y}_{1}\left(t^{n+1},\mathbf{x}_{1}^{n+1(j+1)},\mathbf{z}_{1}^{n+1(j+1)},\mathbf{u}_{1}^{n+1(j+1)}\right). $$ # Step 2 (Substitute): (Step 2.a) Solve the input-output equation for the input for partition 2, i.e., solve $$ \begin{array}{r}{0=\mathbf{U}_{2}\left(\mathbf{u}_{1}^{n+1(j+1)},\mathbf{u}_{2}^{n+1(j+1)},\mathbf{u}_{3}^{n+1(j)},\dots,\mathbf{u}_{N}^{n+1(j)},\mathbf{y}_{1}^{n+1(j+1)},\mathbf{y}_{2}^{n+1(j)},\mathbf{y}_{3}^{n+1(j)},\dots,\mathbf{y}_{N}^{n+1(j)}\right),}\end{array} $$ to yield u2n+1(j+1). (Step 2.b) Advance the solution with input $\mathbf{u}_{2}^{\alpha}=\mathbf{L}\left(\alpha_{2},\mathbf{u}_{2}^{n},\mathbf{u}_{2}^{n+1(j+1)}\right)$ to yield state, constraint, and output values, i.e., $$ {\bf x}_{2}^{n}\,,\dot{\bf x}_{2}^{n}\,,\ldots\,,\dot{\bf x}_{2}^{n-m+1}\,,{\bf z}_{2}^{n}\,,{\bf u}_{2}^{\alpha}\,,{\bf X}_{2}(t)\,,{\bf Z}_{2}(t)\,\xrightarrow{\mathrm{ADV}}{\bf x}_{2}^{n+1(j+1)}\,,{\bf z}_{2}^{n+1(j+1)}\,, $$ $$ \mathbf{y}_{2}^{n+1(j+1)}=\mathbf{Y}_{2}\left(t^{n+1},\mathbf{x}_{2}^{n+1(j+1)},\mathbf{z}_{2}^{n+1(j+1)},\mathbf{u}_{2}^{n+1(j+1)}\right). $$ Again, the value of $\alpha_{2}$ is related to the underlying ODE/DAE integrator. (Step 2.c) In the same manner as in Step 2.b, advance the solution of the remaining $N\,-\,2$ partitions $(\mathbf{z}_{i}^{n+1(j+1)}$ , $\mathbf{x}_{i}^{n+1(j+1)}$ for $i\in\{3,\ldots,N\})$ ), always calculating the input with the most updated outputs from other partitions, i.e., the input equation for partition $i\in\{3,\ldots,N\}$ is $$ \begin{array}{r}{\mathbf{0}=\mathbf{U}_{i}\left(\mathbf{u}_{1}^{n+1(j+1)},\allowbreak\cdot\cdot,\mathbf{u}_{i}^{n+1(j+1)},\mathbf{u}_{i+1}^{n+1(j)},\allowbreak\cdot\cdot,\mathbf{u}_{N}^{n+1(j)},\mathbf{y}_{1}^{n+1(j+1)},\allowbreak\cdot\cdot,\mathbf{y}_{i-1}^{n+1(j+1)},\mathbf{y}_{i}^{n+1(j)},\allowbreak\cdot\cdot,\mathbf{u}_{N}^{n+1(j)},\allowbreak\cdot\cdot\cdot,\mathbf{u}_{N}^{n+1(j+1)},\allowbreak\cdot\cdot\cdot,\mathbf{u}_{N}^{n+1(j+1)},\allowbreak\cdot\cdot\cdot,\mathbf{u}_{N}^{n+1(j+1)},\allowbreak\cdot\cdot\cdot,\mathbf{u}_{N}^{n+1(j)},\allowbreak\cdot\cdot\cdot\times,\mathbf{u}_{N}^{n+1(j)},}\end{array} $$ which is solved to yield uin+1(j+1). Step 3 (Correct): Correct the solution for partition 1. (Step 3.a) Solve the input-output equation for the input for partition 1, i.e., solve $$ \mathbf{0}=\mathbf{U}_{1}\left(\mathbf{u}_{1}^{n+1(j+2)},\mathbf{u}_{2}^{n+1(j+1)},\ldots,\mathbf{u}_{N}^{n+1(j+1)},\mathbf{y}_{1}^{n+1(j+1)},\ldots,\mathbf{y}_{N}^{n+1(j+1)}\right), $$ to yield u1n+1(j+2). (Step 3.b) Advance the solution with input $\mathbf{u}_{1}^{\alpha}=\mathbf{L}\left(\alpha_{1},\mathbf{u}_{1}^{n},\mathbf{u}_{1}^{n+1(j+2)}\right)$ to yield state, constraint, and output values, i.e., $$ {\bf x}_{1}^{n}\,,\dot{\bf x}_{1}^{n}\,,\dots\,,\dot{\bf x}_{1}^{n-m+1}\,,{\bf z}_{1}^{n}\,,{\bf u}_{1}^{\alpha}\,,{\bf X}_{1}(t)\,,{\bf Z}_{1}(t)\,\xrightarrow{\mathrm{ADV}}{\bf x}_{1}^{n+1(j+2)}\,,{\bf z}_{1}^{n+1(j+2)}\,, $$ $$ \mathbf{y}_{1}^{n+1(j+2)}=\mathbf{Y}_{1}\left(t^{n+1},\mathbf{x}_{1}^{n+1(j+2)},\mathbf{z}_{1}^{n+1(j+2)},\mathbf{u}_{1}^{n+1(j+2)}\right). $$ Let $j=j+1$ . If $j
System 1System 2
m1C1k1m2C2k2Cckc
1.00.11.01.00.11.00.011.0
![](images/24e59ff5dab881e45a25b389400718028240d135e1f1143ed1d7206a8ff091f5.jpg) Figure 6. Benchmark displacement histories $q_{1}(t)$ and $q_{2}(t)$ for the system parameters listed in Table 1. # IV.A.1. Stability Analysis Figure 7 shows numerically determined critical time increments $\Delta t^{c r i t}$ as a function of coupling-spring stiffness ( $k_{c}$ ) for the three coupling schemes with the three time integrators. Here, $\Delta t^{c r i t}$ is the largest time increment for which bounded numerical solutions were produced. Also included are the critical time increments associated with uncoupled integration of System 2, $\Delta t_{2}^{c r i t}$ , which was determined with the displacement of the interface node held fixed. For all cases shown, $\Delta t_{2}^{c r i t}<\Delta t_{1}^{c r i t}$ , where $\Delta t_{1}^{c r i t}$ is the critical time increment for the uncoupled System 1. In regard to the numerical-stability data shown in Figure 7, we make the following observations: # • $\Delta t^{m a x}=\mathrm{min}(\Delta t_{1}^{c r i t},\Delta t_{2}^{c r i t})$ is an upper bound on stable time increments for all cases examined. • For integration of the monolithic system, RK4 was the most stable, while ABM4 was significantly more stable than AB4. However, RK4 requires four RHS evaluations per time step, ABM4 requires two, while AB4 requires only one. From an accuracy vs. computational-cost perspective, the three integrators are comparable performers for the monolithic system. • Explicit coupling is prone to numerical instability for RK4 and ABM4 integration schemes, and can require $\Delta t\,\ll\,\Delta t^{m a x}$ for stable solutions. This is, perhaps, not surprising. RK4 and ABM4 both have an implicit dependence on input (from the other partition) over $t^{n}\,<\,t\,\leq\,t^{n+1}$ , but the input is held constant in the implementation here, which introduces additional stiffness. Explicitly coupled AB4, however, exhibits stability characteristics indistinguishable from those of the AB4 monolithic treatment. • For RK4 and ABM4 integration, PC coupling is significantly more stable than explicit coupling. The stability of PC(2) is significantly better than that of PC(1). For some RK4 and ABM4 cases, PC(2) exhibited critical time increments that exceeded those of the monolithic approach (but were still bounded by $\Delta t^{m a x}$ ). • While not shown, PC(3) solutions exhibited stability characteristics virtually indistinguishable from those of PC(2). ![](images/40d8b47233f9d7a7ec00b26176d22ede5b30feb534516abf8c187c7ad79a6c71.jpg) Figure 7. Numerically determined critical time increments as a function of coupling-spring stiffness for (a) RK4, (b) ABM4, and (c) AB4 time integration and with the three coupling schemes. Only monolithic and explicit coupling is shown for AB4, for which PC-implicit and explicit are equivalent. $\Delta t_{2}^{c r i t}$ corresponds to the critical increment associated with uncoupled integration of System 2. # IV.A.2. Accuracy Analysis We examine here the accuracy of the three numerical coupling methods for the three numerical integrators, and with the baseline parameters listed in Table 1. Figure 8 shows normalized root-mean-square (RMS) error of the numerical solutions for the displacement of System 1 over the time interval $0\le t\le30$ . Normalized RMS error for $n_{m a x}$ numerical response values $q_{1}^{n}$ , where $q_{1}^{n}\approx q_{1}(t^{n})$ , was calculated as $$ \varepsilon\left({{q}_{1}}\right)=\sqrt{\frac{\sum_{k=0}^{{n_{m a x}}}\left[{{q}_{1}^{k}}-{{q}_{b}}(t^{k})\right]^{2}}{\sum_{k=0}^{{n_{m a x}}}\left[{{q}_{b}}(t^{k})\right]^{2}}}, $$ where $q_{b}(t)$ is a known benchmark solution; $q_{b}(t)$ is either known analytically or is interpolated from a highly resolved numerical solution. In regard to the numerical-accuracy data shown in Figure 8, we make the following observations: • Integration of the monolithic systems exhibits expected accuracy (convergence rates) for the three fourth-order time integrators. RK4 and ABM4 exhibit similar accuracy that is better than that of AB4. However, RK4 requires four ODE-RHS evaluations per time step, ABM4 requires two, while AB4 requires only one. From this perspective, the three time integrators applied to the monolithic system exhibit similar computational cost for a given accuracy. • Explicit coupling with RK4 and ABM4 exhibits only first-order accuracy, while fourth-order accuracy is maintained with AB4 integration. RK4 and ABM4 integrators are clearly hindered by their implicit dependence on other-partition data. • For RK4 integration, PC(1) and PC(2) provide second-order accuracy. Numerical experiments showed that RK4 with PC( $j_{m a x}\ge1$ ) is no better than second-order accurate. This limitation in accuracy is attributed to the fact that input data is held constant over the time step, even though a high-order multi-step method like RK4 performs RHS evaluations within the time step. • For ABM4, PC(1) exhibits second-order accuracy, while PC(2) exhibits the desired fourth-order accuracy of the underlying time integrator. • For this simple partitioned test system, AB4 integration with explicit coupling is clearly the best performer. While ABM4 with PC(2) exhibits slightly better accuracy for a given $\Delta t$ , it requires five “steps” with ABM4 (one prediction, two substitutions, and two correctors) per time step. Alternatively, AB4 with explicit coupling requires only two “steps” (one per module) per time step. However, this strong performance is largely tied to the simplicity of this non-stiff ODE system; AB4 and explicit coupling may well be a poor choice for stiff systems and/or DAE systems. # IV.B. A Linear 1-DOF Damped Oscillator Connected with a Nonlinear Quasi-Static Cable We examine here numerical-accuracy results for System 1, a one-DOF damped oscillator, coupled to System 3, a nonlinear quasi-static cable, as described in Section II.B.2. Dimensionless property values for both systems are listed in Table 2. Benchmark solutions for the coupled System 1 and System 3 were generated with IDA, from the Sundials $^{17}$ software suite. IDA employs variable-step variable-order backwards differentiation formulas. Cable displacement and force benchmark histories are shown in Figure 9(a) and (b), respectively, for $0\,\leq\,t\,\leq\,20$ and for the following initial conditions: $q_{1}(0)\,=\,1$ , ${\dot{q}}_{1}(0)\,=\,0$ , $q_{3}(0)\,=\,1$ . For comparison, Figure 9(a) includes $q_{1}(t)$ calculated in the absence of force from the attached cable; the cable force clearly has a significant impact on the response of $q_{1}(t)$ . Table 2. System 1 and 3 dimensionless parameters associated with the test simulations.
System 1System 3
m1C1k1DLAEM
1.00.13.01.53.01.0500002.9033
In our numerical implementation, System 1 is treated in the same manner as in Section IV.A, where one of several time-integrators can be employed. The nonlinear constraint equation in System 3 is solved with a standard Newton-Raphson (NR) algorithm with a tolerance of $10^{-10}$ . The System 3 module is purely a root-finding algorithm; the output will correspond to the same time as that of the input. Under PC-coupling, partition 1 (System 3) employs $\alpha_{1}=1$ . Due to the limited performance exhibited by RK4 in our coupling schemes, we focus here on ABM4 and AB4 integrators coupled to the NR constraint solve. Figure 10 shows RMS errors of $q_{1}(t)$ histories produced with ABM4 and AB4 integration of System 1 and with explicit, PC(1), and PC(2) coupling to System 3. Explicit coupling yields only first-order convergence, as the force from the cable (a constraint variable) is time lagged by one step, which introduced $O(\Delta t)$ error. This first-order error can be reduced to $O(\Delta t^{2})$ if the cable force constraint is calculated based on a linearly extrapolated input as in Eq. (40) in PC coupling. However, PC(1) coupling provides second-order accuracy for ABM4-NR and fourth-order accuracy for AB4-NR. AB4-NR achieves desired accuracy with only PC(1) because System 1 is advanced with the consistent value of the cable force (at $t^{n}$ ), and the System 3 cable force is then corrected to be consistent with displacement at $t^{n+1}$ . For this coupled system, ABM4-NR required PC(2) coupling for fourth-order accuracy. This is because, in the first cycle (embodied in PC(1)), ABM4 is advanced with a $O(\Delta t^{2})$ prediction of the constraint; a second correction cycle is required to obtain desired accuracy. Here, solutions converge to error of $O(10^{-6})$ , rather than machine precision, due to the discrete nature of the benchmark solution and the interpolation of that solution for comparison. ![](images/dff1cce8fca48fd5ced67a98438fa70c0711fae59db6a09cf6c2447e0ddd01c7.jpg) Figure 8. Normalized RMS error of System 1 displacement histories over $0\,\leq\,t\,\leq\,30$ as a function of time increment for (a) RK4, (b) ABM4, and (c) AB4 time integration and the three coupling schemes. Only monolithic and explicit coupling is shown for AB4, for which PC-implicit and explicit are equivalent. Dashed lines show reference convergence rates. (c) AB4 ![](images/f36378482ad06867bf3a47d47e4bdc07ce3188024dc5a85eb97aa2b0ab8dfa76.jpg) Figure 9. Numerically calculated benchmark response histories for (a) System 1 displacement, and (b) System 3 cable force. System 1 displacements are shown with and without the effects of the attached cable. ![](images/00f9a2d8bdbe9870a0c2175fc0e033b26798ccaebf3fec098ea121b2ed7b1b26.jpg) Figure 10. Normalized RMS errors of $q_{1}(t)$ histories calculated for $0\le t\le20$ with (a) ABM4-NR and (b) AB4-NR and with various coupling schemes; PC(1) and $\mathbf{PC}(\mathbf{2})$ results are indistinguishable for AB4-NR. Dashed lines show reference convergence rates. # V. Conclusions This work was motivated by a major restructuring of the FAST wind turbine CAE tool suite; details of the restructuring are described in a companion paper.6 Our goal was to guide FAST-module developers and to provide them with stable, accurate, and efficient numerical coupling strategies for multi-physics multi-module simulations. Using two simple example systems, we described possible partition choices and methods for coupling those partitions, under the restriction that data is communicated between partitions only at full time steps. We examined two loose-coupling schemes, explicit and implicit. Loose coupling implies that modules may be updated sequentially, rather than concurrently, which provides flexibility in modularization and parallelization. In explicit coupling, partitions are advanced in time using only known information. In implicit coupling, one or more partitions require data from other partitions at the next step. Here, implicit coupling is accomplished in a loose-coupling manner through a predictor-corrector approach. Time integration was examined with three fixed-step fourth-order schemes: Runge-Kutta (explicit singlestep method), Adams-Bashforth-Moulton (predictor-corrector multi-step method), and Adams-Bashforth (explicit multi-step method). Through numerical experiments, it was shown that, if the underlying ODE integrator has an implicit dependence on other-partition data (as in RK and ABM), or if there is an inherent implicit dependence due to constraint equations (as in our System 3 coupled to System 1), explicit coupling can be dramatically less stable and less accurate than simulations performed with the monolithic system. However, predictor-corrector implicit coupling can restore desired accuracy and stability. For coupled partitions without constraints, explicit time integration with AB and explicit loose coupling exhibited desired accuracy and stability. Based on this work we make the following recommendations. First, for general multi-partition systems, we recommend implicit coupling through a predictor-corrector approach, which is significantly more robust than explicit coupling. Second, high-order single step methods such as Runge-Kutta should be avoided under the restriction that partitions may only communicate at whole time steps. In future work, we will extend our analysis to include time-step subcycling,13,18,19 where one or more partitions is updated with a time increment that is a fraction of the largest time increment used by other systems. Further, we will examine accuracy and stability associated with higher-order extrapolation of input and output data in our PC algorithm; e.g., Eqs. (40) and (41), respectively. Linear extrapolation was employed in this work; however, quadratic extrapolation may provide significant improvement in accuracy with limited stability degradation.20 Regardless of the extrapolation accuracy, desired accuracy of the underlying multi-step time integrators appears attainable if a sufficient number of correction steps are applied. In our examples, fourth-order accuracy with multi-step methods was attained with two correction iterations (and only one in some cases). We will also compare the loose-coupling schemes described in this paper with tight implicit coupling, where the partitioned systems are assembled into a system suitable for concurrent time update with a single DAE solver, and where constraints are naturally introduced from the input-output relations. Future work may also include the examination of implicit numerical integrators, and models with discrete-time states. # Acknowledgments The authors acknowledge useful discussions with Bonnie Jonkman and John Michalakes. Thanks also to Dr. Ray Grout and Dr. Peter Graf for reviewing the paper prior to submission. This work was performed at NREL in support of the U.S. Department of Energy under contract number DE-AC36-08GO28308 and under a project funded through topic area 1.1 of the Funding Opportunity Announcement (FOA) number DE-FOA-0000415. # References $^{1}$ Jonkman, J. M. and Buhl, M. L., “FAST User’s Guide,” Tech. Rep. NREL/EL-500-38230, National Renewable Energy Laboratory, Golden, CO, August 2005. $^2$ Laino, D. J. and Hansen, A. C., “Users Guide to the Wind Turbine Dynamics Aerodynamics Computer Software AeroDyn,” December 2002, Windward Engineering LLC. Prepared for the National Renewable Energy Laboratory under Subcontract No. TCX-9-29209-01. $^3$ Moriarty, P. J. and Hansen, A. C., “AeroDyn Theory Manual,” Tech. Rep. NREL/EL-500-36881, National Renewable Energy Laboratory, Golden, CO, December 2005. $^4$ Jonkman, J. M., Dynamics Modeling and Loads Analysis of an Offshore Floating Wind Turbine, Ph.D. thesis, Department of Aerospace Engineering Sciences, University of Colorado, Boulder, CO, 2007, also published in tech. report NREL/TP-500- 41958, National Renewable Energy Laboratory, Golden, CO. $^{5}$ Jonkman, J. M., “Dynamics of Offshore Floating Wind Turbines–Model Development and Verification,” Wind Energy, Vol. 12, 2009, pp. 459–492, also published in tech. report NREL/JA-500-45311 National Renewable Energy Laboratory, Golden, CO. $^6$ Jonkman, J. M., “The New Modularization Framework for the FAST Wind Turbine CAE Tool,” 2012, To appear in the proceedings of the 32nd ASME Wind Energy Symposium, 51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition. 7Larson, J. W., Jacob, R. L., Foster, I., and Guo, J., “The Model Coupling Toolkit,” Proceedings of the 2001 International Conference on Computational Science, 2001. $^8$ Bettencourt, M. T., “Distributed Model Coupling Framework,” Proceedings of the 11th IEEE Symposium on High Performance Distributed Computing, Edinburgh, UK, 2002, pp. 284–290. $^{9}$ Elwasif, W. R., Bernholdt, D. E., Shet, A. G., Foley, S. S., Bramley, R., Batchelor, D. B., and Berry, L. A., “The design and implementation of the SWIM Integrated Plasma Simulator,” Proceedings of the 18th Euromicro Conference on Parallel, Distributed and Network-based Processing, 2010. $^{10}$ Pawlowski, R., Bartlett, R., Belcourt, N., Hooper, R., and Schmidt, R., “A Theory Manual for Multi-Physics Code Coupling in LIME,” Tech. Rep. SAND2011-2195, Sandial National Laboratories, Albuquerque, New Mexico and Livermore, California, March 2011. $\mathrm{~11}$ Felippa, C. A., Park, K. C., and Farhat, C., “Partitioned analysis of coupled mechanical systems,” Computer Methods in Applied Mechanics and Engineering, Vol. 190, 2001, pp. 3247–3270. ${\mathrm{12}}$ Farhat, C., van der Zee, K. G., and Geuzaine, P., “Provably second-order time-accurate loosely-coupled solution algorithms for transient nonlinear computational aeroelasticity,” Computer Methods in Applied Mechanics and Engineering, Vol. 195, 2006, pp. 1973–2001. $^{13}$ Belytschko, T., Yen, H. J., and Mullen, R., “Mixed Methods for Time Integration,” Computer Methods in Applied Mechanics and Engineering, Vol. 17, 1979, pp. 259–275. $^{14}$ Jonkman, J. M. and Buhl, M. L., “Development and Verification of a Fully Coupled Simulator for Offshore Wind Turbines,” Proceedings of the 45th AIAA Aerospace Sciences Meeting, Reno, NV, 2007, also published in tech. report NREL/CP500-40979 National Renewable Energy Laboratory, Golden, CO. $^{15}$ Brenan, K. E., Campbell, S. L., and Petzold, L. R., Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, Society of Industrial and Applied Mathematics, 1996. $^{16}$ Matthies, H. G., Nekamp, R., and Steindorf, J., “Algorithms for strong coupling procedures,” Computer Methods in Applied Mechanics and Engineering, Vol. 195, 2006, pp. 2028–2049. 17Hindmarsh, A. C., Brown, P. N., Grant, K. E., Lee, S. L., Serban, R., Shumaker, D. E., and Woodward., C. S., “SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers,” ACM Transactions on Mathematical Software, Vol. 31, 2005, pp. 363–396. $^{18}$ Belytschko, T., “Partitioned and Adaptive Algorithms for Explicit Time Integration,” Nonlinear Finite Element Analysis in Structural Mechanics, edited by W. Wunderlich, E. Stein, and K. J. Bath, Springer-Verlag, 1981, pp. 572–584. $^{19}$ Hulbert, G. M. and Hughes, T. J. R., “Numerical Evaluation and Comparison of Subcycling Algorithms for Structural Dynamics,” Tech. Rep. DNA-TR-88-8, Defense Nuclear Agency, Washington, D.C., 1988. 20Elliott, A. S., “A Highly Efficient, General-Purpose Approach for Co-Simulation with ADAMS $^\mathrm{\textregistered}$ ,” Proceedings of the 15th European ADAMS Users’ Conference, 2000.