vault backup: 2025-07-22 09:09:33

书籍/力学书籍/OpenFast/Numerical Stability and Accuracy of Coupled Multi-Physics/auto/Numerical Stability and Accuracy of Coupled Multi-Physics.md

@@ -0,0 +1,453 @@
+# 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<j_{m a x}$ , repeat Steps 2 and 3. If $j=j_{m a x}$ , proceed to Step 4.  
+
+Step 4: Save all the final variables,  
+
+$$
+\mathbf{x}_{i}^{n+1}=\mathbf{x}_{i}^{n+1(j_{m a x})},\quad\mathbf{z}_{i}^{n+1}=\mathbf{z}_{i}^{n+1(j_{m a x})},\quad\mathbf{y}_{i}^{n+1}=\mathbf{y}_{i}^{n+1(j_{m a x})},\quad\mathbf{u}_{i}^{n+1}=\mathbf{u}_{i}^{n+1(j_{m a x})},
+$$  
+
+for each $i\in\{2,\dots,N\}$ , and  
+
+$$
+\begin{array}{r}{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\cdot\mathbf{x}^{n+1}=\mathbf{x}_{1}^{n+1(j_{m a x}+1)},\quad\mathbf{z}_{1}^{n+1}=\mathbf{z}_{i}^{n+1(j_{m a x}+1)},\quad\mathbf{y}_{1}^{n+1}=\mathbf{y}_{1}^{n+1(j_{m a x}+1)},\quad\mathbf{u}_{1}^{n+1}=\mathbf{u}_{1}^{n+1(j_{m a x}+1)},}\end{array}
+$$  
+
+which completes solution advancement to time $t^{n+1}$  
+
+$\mathrm{PC}(j_{m a x})$ , as described above, requires $(j_{m a x}N+1)$ ADV evaluations per time step. Hence, for a given $\Delta t$ , the PC approach requires more operations than explicit coupling. However, this additional cost per time step is negated if the PC-approach is significantly more numerically stable and/or more accurate than explicit coupling. The above algorithm could be modified to stop iterating after convergence to within some tolerance, e.g., iterations continue until $\left|\mathbf{x}_{i}^{n+1(j+1)}-\mathbf{x}_{i}^{n+1(j)}\right|\;<\;\epsilon$ and $\left|\mathbf{z}_{i}^{n+1(j+1)}-\mathbf{z}_{i}^{n+1(j)}\right|\ <\ \epsilon$ for each $i\in\{1,\ldots,N\}$ , where $\epsilon$ is a user-specified￿ tolerance.  
+
+# IV. Application Example Results and Discussions  
+
+In this section, we examine the numerical performance of the three coupling approaches, described in Section III, using the two application examples, described in Section II.B. We examine the numerical stability and accuracy of three standard, fixed-time-increment methods: (1) fourth-order, explicit Runge-Kutta (RK4) which is a single-step (though “multi-stage”) method, (2) fourth-order Adams-Bashforth-Moulton (ABM4), which is an explicit multi-step method with an implicit dependence on partition-interaction data, and (3) fourth-order explicit Adams-Bashforth (AB4). ABM4 also serves as a proxy for a multi-step implicit time integrator as it has the same implicit dependence on other-partition data.  
+
+We note that the RK4 method is self-starting, while ABM4 and AB4 are not. In our simulations with ABM4 and AB4, the first four time steps are initialized with benchmark solutions. In practice, ABM4 and AB4 could be initialized with RK4 integrated with sufficiently small time increments. In our PC implementation of RK4, numerical experiments showed that $\alpha=0.5$ provides the best performance. RK4 performs four RHS evaluations over a time step; providing a midpoint average of the input seems prudent. For ABM4, PC coupling is accomplished with $\alpha=1$ , while $\alpha=0$ is employed for AB4. As such, AB4 PC simulations for System 1 coupled to System 2 are equivalent to those under explicit coupling. However, in the PC simulations where System 1 is coupled to System 3, the constraint is updated with $\alpha=1$ . Importantly, PC-coupling (as described in Section III.B.2) was implemented with System 2 (or 3) as the first partition, and with System 1 as the second partition. Hence, System 2 (or 3) was predicted (Step 1), System 1 was substituted (Step 2), and then System 2 (or 3) was corrected (Step 3).  
+
+# IV.A. Linear 2-DOF Damped Oscillator System  
+
+We implemented the three numerical-coupling algorithms described in Section III (monolithic, explicit loose, and PC implicit loose) for the linear 2-DOF damped oscillator system shown in Figure 4. Throughout this subsection, solutions were calculated for $0\le t\le30$ in the absence of external forcing ( $f_{1}=f_{2}=0$ ), with (dimensionless) initial conditions given by  
+
+$$
+q_{1}(0)=1,\quad\dot{q}_{1}(0)=0,\quad q_{2}(0)=0,\quad\dot{q}_{2}(0)=0,
+$$  
+
+and with baseline system parameters listed in Table 1. For accuracy assessment, a closed-form exact solution was calculated; Figure 6 shows the exact benchmark displacement histories.  
+
+Table 1. System 1 and 2 baseline parameters associated with the test simulations.   
+
+
+<html><body><table><tr><td colspan="3">System 1</td><td colspan="5">System 2</td></tr><tr><td>m1</td><td>C1</td><td>k1</td><td>m2</td><td>C2</td><td>k2</td><td>Cc</td><td>kc</td></tr><tr><td>1.0</td><td>0.1</td><td>1.0</td><td>1.0</td><td>0.1</td><td>1.0</td><td>0.01</td><td>1.0</td></tr></table></body></html>  
+
+![](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.   
+
+
+<html><body><table><tr><td colspan="3">System 1</td><td colspan="5">System 3</td></tr><tr><td>m1</td><td>C1</td><td>k1</td><td>D</td><td>L</td><td>A</td><td>E</td><td>M</td></tr><tr><td>1.0</td><td>0.1</td><td>3.0</td><td>1.5</td><td>3.0</td><td>1.0</td><td>50000</td><td>2.9033</td></tr></table></body></html>  
+
+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.  

书籍/力学书籍/OpenFast/The New Modularization Framework for the FAST/auto/The New Modularization Framework for the FAST.md

@@ -0,0 +1,518 @@
+# The New Modularization Framework for the FAST Wind Turbine CAE Tool  
+
+Preprint  
+
+J. Jonkman  
+
+To be presented at the 51st AIAA Aerospace Sciences Meeting, including the New Horizons Forum and Aerospace Exposition Dallas, Texas   
+January 7-10, 2013  
+
+# NOTICE  
+
+The submitted manuscript has been offered by an employee of the Alliance for Sustainable Energy, LLC (Alliance), a contractor of the US Government under Contract No. DE-AC36-08GO28308. Accordingly, the US Government and Alliance retain a nonexclusive royalty-free license to publish or reproduce the published form of this contribution, or allow others to do so, for US Government purposes.  
+
+This report was prepared as an account of work sponsored by an agency of the United States government. Neither the United States government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights.  Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or any agency thereof.  The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States government or any agency thereof.  
+
+Available electronically at http://www.osti.gov/bridge  
+
+Available for a processing fee to U.S. Department of Energy   
+and its contractors, in paper, from: U.S. Department of Energy Office of Scientific and Technical Information P.O. Box 62 Oak Ridge, TN 37831-0062 phone:  865.576.8401 fax: 865.576.5728 email:  mailto:reports@adonis.osti.gov  
+
+Available for sale to the public, in paper, from:  
+
+U.S. Department of Commerce   
+National Technical Information Service   
+5285 Port Royal Road   
+Springfield, VA 22161   
+phone:  800.553.6847   
+fax:  703.605.6900   
+email: orders@ntis.fedworld.gov   
+online ordering:  http://www.ntis.gov/help/ordermethods.aspx  
+
+Cover Photos: (left to right) PIX 16416, PIX 17423, PIX 16560, PIX 17613, PIX 17436, PIX 17721  
+
+Printed on paper containing at least $50\%$ wastepaper, including $10\%$ post consumer waste.  
+
+# The New Modularization Framework for the FAST Wind Turbine CAE Tool  
+
+Jason M. Jonkman National Renewable Energy Laboratory (NREL), Golden, Colorado, 80401  
+
+NREL recently has put considerable effort into improving the overall modularity of its FAST wind turbine aero-hydro-servo-elastic tool to (1) improve the ability to read, implement, and maintain source code; (2) increase module sharing and shared code development  across  the  wind  community;  (3)  improve  numerical  performance  and robustness;  and  (4)  greatly  enhance  flexibility  and  expandability  to  enable  further developments of functionality without the need to recode established modules.  The new FAST modularization framework supports module-independent inputs, outputs, states, and parameters; states in continuous-time, discrete-time, and constraint form; loose and tight coupling; independent time and spatial discretizations; time marching, operating-point determination, and linearization; data encapsulation; dynamic allocation; and save/retrieve capability.  This paper explains the features of the new FAST modularization framework, as well as the concepts and mathematical background needed to understand and apply it correctly.  It is envisioned that the new modularization framework will transform FAST into a powerful, flexible, and robust wind turbine modeling tool with a large number of developers and a range of modeling fidelities across the aerodynamic, hydrodynamic, servodynamic, and structural-dynamic components.  
+
+# Nomenclature  
+
+A $=$ continuous-state matrix for linearized continuous-state equations   
+AADC $=$ continuous-state matrix for linearized discrete-state equations   
+Ad $=$ discrete-state matrix for linearized discrete-state equations   
+ADAC $=$ discrete-state matrix for linearized continuous-state equations   
+$B$ $=$ input matrix for linearized continuous-state equations   
+BADC $=$ input matrix for linearized discrete-state equations   
+$C$ $=$ continuous-state matrix for linearized output equations   
+CDAC $=$ discrete-state matrix for linearized output equations   
+$D$ $=$ input-transmission matrix for linearized output equations   
+$G$ $=$ the matrix formed by $\frac{\partial U}{\partial u}-\frac{\partial U}{\partial z}\bigg[\frac{\partial Z}{\partial z}\bigg]^{-I}\frac{\partial Z}{\partial u}$   
+$H(~)$ $=$ Heaviside-step (unit-step) function   
+$n$ $=$ discrete-step counter   
+$N$ $=$ total number of modules coupled together   
+$p(t)$ $=$ parameters   
+$t$ $=$ time   
+$u(t)$ $=$ inputs   
+$U(\,)$ $=$ input-output transformation functions   
+$x(t)$ $=$ continuous states   
+${\dot{x}}\left(t\right)$ $=$ first time derivative of the continuous states   
+$X(\,)$ $=$ continuous-state functions   
+$x^{d}[n]$ $=$ discrete states   
+$X^{d}(\mathbf{\xi})$ $=$ discrete-state functions   
+$y(t)$ $=$ outputs   
+$Y(\,)$ $=$ output functions   
+$z(t)$ $=$ constraint (algebraic) states   
+$Z(\,)$ $=$ constraint-state (algebraic) functions   
+$\varDelta t$ $=$ discrete time step (increment)  
+
+# I. Introduction  
+
+T HE wind industry relies extensively on computer-aided-engineering (CAE) tools for wind turbine performance, loads, and stability analyses.  Limitations—and consequent inaccuracies—in the tools slow the advancement of wind power.  Accurate tools are required for the wind industry to develop more innovative, optimized, reliable, and cost-effective wind technology.  Overcoming current modeling limitations increases in importance as turbines scale up to larger sizes, incorporate novel architectures and load-control technologies, and are installed on offshore support platforms.  
+
+Over the past two decades, the U.S. Department of Energy (DOE) has sponsored NREL’s development of CAE tools for wind turbine analysis.  The tools are developed as free, publicly available, open-source, professional-grade products as a resource for the wind industry.  The tools are used by thousands of wind turbine designers, manufacturers, consultants, certifiers, researchers, educators, and students throughout the world.  The open-source approach facilitates the tools’ credibility and adaptability by the wind industry.  The tools are modular, well documented, and supported by NREL through workshops and an on-line forum.  They have been verified through model-to-model comparisons and validated with test measurements.  
+
+Analyzing wind-energy structures requires CAE tools that model the important physical phenomena (illustrated in Figure 1) and system couplings, including the environmental excitation (wind, waves, and current) and fullsystem dynamic response (rotor, drivetrain, nacelle, support structure, and controller).  
+
+omputers because the certificationdriven design process is iterative and must consider a vast set of environmental conditions  and operational scenarios. Computationally intensive solutions are generally unsuitable for these applications because the long run times make it impossible to  consider  all  of  the  required cases. So,  CAE  tools  cannot solely  be  massively  discretized high-performance computing (HPC) solutions of the fundamental laws of physics (e.g., computational fluid dynamics (CFD)  solutions  of  the  NavierStokes equations). Instead, NREL’s core CAE tool, FAST,1,2 is based on advanced engineering models—derived from fundamental laws, but with appropriate simplifications and assumptions,  and  supplemented where applicable with computational solutions and test data, as illustrated in Figure 2.  
+
+![](images/2c5c06548259b5732cafbb396ac5abbe751727f5326b37835cf9b9f6ef4825f1.jpg)  
+Figure 1.  Physical phenomena affecting a floating wind turbine system.  
+
+As shown in Figure 3 and Figure 4, FAST joins a rotor aerodynamics module (AeroDyn ), a platform hydrodynamics module (HydroDyn ) for offshore systems, a control and electrical system (servo) dynamics module, and a structural (elastic) dynamics module to enable coupled nonlinear aerohydro-servo-elastic analysis in the time domain.  The FAST tool enables the analysis of a range of wind turbine configurations, including two- or threeblade horizontal-axis rotor, pitch or stall regulation, rigid or teetering hub, upwind or downwind rotor, and lattice or tubular tower.  The wind turbine can be modeled on land or offshore on fixed-bottom or floating substructures.  
+
+![](images/02d1fb86013489bfa381fa216ec1ea11414cbf20fc8a92cb5ef4275598d08f7e.jpg)  
+Figure 2.  CAE Tools derived from  theory,  computational solutions, and test data.  
+
+AeroDyn uses wind-inflow data and solves for the rotor-wake effects and blade-element  aerodynamic  loads,  including  dynamic  stall.    HydroDyn simulates the regular or irregular incident waves and currents and solves for the hydrostatic, radiation, diffraction, and viscous loads on the offshore substructure.    The  control  and  electrical  system  module  simulates  the controller logic, sensors, and actuators of the blade-pitch, generator-torque, nacelle-yaw, and other control devices, as well as the generator and powerconverter components of the electrical drive.  The structural-dynamics module applies the control and electrical system reactions, applies the aerodynamic and hydrodynamic loads, adds gravitational loads, and simulates the elasticity of the rotor, drivetrain, and support structure.  The modular interface and coupler enables interactions between all modules.  All of the modules are  
+
+available fo free on the NREL website7 and include compiled executables, source  code,  and  sample  input  data.    User  and  theory documentation for each module are available.1-6  
+
+One of the benefits of a modular framework is that it allows modules to be interchanged; this feature is important for  benchmarking,  research,  and  industrial  applications because the required model fidelity is dependent on the application.  For example, the fidelity of FAST’s structural module is currently limited by its degrees of freedom (DOF). However, the AeroDyn and HydroDyn modules have also been interfaced to the nearly unlimited DOF multibodydynamics tool MSC.ADAMS—a commercially available and general-purpose tool from MSC Software Corporation— through  the  ADAMS-to-AeroDyn  (A2AD)  interface  to enable higher-fidelity structural analysis of wind systems with proper aerodynamic and hydrodynamic loading.  A range of  
+
+![](images/4390c8817fa65963abbea64a95ba9f81092fefe958202221f09ce329d73e5cf6.jpg)  
+Figure 3.  NREL’s modular CAE tool, FAST.  
+
+theories of varying model fidelity is also plausible for  other  physics. In rotor-wake   
+aerodynamics, for example, available theories of  different fidelity include  (from lowest to highest fidelity) blade  
+element/momentum   
+(BEM) theory, generalized dynamic wake (GDW) theory, vortex-wake methods, and CFD. Available theories in hydrodynamics include  
+
+![](images/359e99c2c9103f07dca06fec0914b6e72225e83d540677e6b55947bea34d1979.jpg)  
+Figure 4.  Coupled aero-hydro-servo-elastic interaction.  
+
+Morison’s  equation  with  linear waves;  Morison’s  equation  with higher-order  waves;  potential-flow theory  with linear  radiation, diffraction, and hydrostatics; secondorder  potential-flow  theory;  and CFD.  Many of these theories have not been implemented within FAST, but would be of great value.  There are  also  many  CAE  tools  where these individual theories have been implemented (e.g., the Aerodynamic Wind  turbine  Simulation  Module (AWSM)9 from the Energy research Centre of the Netherlands (ECN) for the vortex-wake method), but most have not yet been interfaced with  
+
+![](images/11d14e77d0fae682992bd1d1eeb70512fb049c1d58d05245a4d67d859c8ff5d4.jpg)  
+Figure 5.  Conceptualization of module interchange.  
+
+FAST.  The concept of interchanging modules is illustrated in Figure 5, although in practice all modules are interfaced through a common driver program (i.e., modular interface and coupler or “glue code”) as illustrated in Figure 3.  
+
+While the FAST CAE tool has always been modular, NREL recently has put considerable effort into improving its overall modularity to accomplish the following benefits: (1) improved ability to read, implement, and maintain source code; (2) increased module sharing and shared code development across the wind community; (3) improved numerical performance and robustness; and (4) greatly enhanced flexiblity 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, flexible, and robust wind turbine CAE tool with a large number of developers and a range of modeling fidelities across the aerodynamic, hydrodynamic, servodynamic, and structural-dynamic components.  
+
+This paper explains the features of the new FAST modularization framework, as well as the concepts and mathematical background needed to understand and apply it correctly.  Specific problems with earlier versions of FAST that the new modularization framework addresses are summarized in Table 1 and are discussed where  
+
+Table 1.  Summary of problems and solutions addressed by the new FAST modularization framework.   
+
+
+<html><body><table><tr><td>Problem</td><td>Solution</td></tr><tr><td>Limited range of modelingfidelity</td><td>Framework allowing modules to be exchanged Developmentofnewmodulesofhigherfidelity</td></tr><tr><td>Solutiondrivenbystructuralsolver</td><td>Separatemoduleinterfaceandcoupler</td></tr><tr><td>·Inabilitytoisolate a given model</td><td>Modulesthatcanbecalledbyseparatedriverprograms orinterfacedtogethertoformacoupledsolution</td></tr><tr><td>Dependentspatialdiscretizationsand timestepsacrossmodules</td><td>Library of spatial elements and mesh-to-mesh mapping ·Datatransferwithinterpolation/extrapolationintime</td></tr><tr><td>Inabilitytolinearizeallsystem equations</td><td>Tightcouplingwithoptionsforoperating-point determinationandlinearization</td></tr><tr><td>Focusonsingleturbine</td><td>Dynamicallocationofmodulesforwind-plant simulation</td></tr><tr><td>"Spaghetticode"due touncleardata transferandglobaldata</td><td>Modularizationwithdataencapsulation</td></tr><tr><td>Limitednumberofdevelopersdueto codesize&complexity</td><td>Modularizationofcodeintoseparatecomponents ·Programmer'shandbookexplainingcodedevelopment requirementsandbestpractices</td></tr><tr><td>Potentiallypoornumericalaccuracy andstability</td><td>Multiple coupling schemes and integration/solver options</td></tr></table></body></html>  
+
+appropriate in this  paper. While  not covered  in  this paper, a programmer’s handbook10 explaining the code development requirements and best practices has been written to support shared code development across the wind community. The  handbook includes a Fortran-based‡ source-code template for developing modules within the framework, including the necessary data structures and interface  procedures. Please note  that  the  current  FAST source  code  at  the  time  of publication (v7.02.00d-bjj) does not yet conform to this framework;  effort  is  being made to convert to this new framework.  
+
+# II.  Features of the New Modularization Framework  
+
+At its core, the new FAST modularization framework is a means by which various mathematical systems (modeling the important physical phenomena as illustrated in Figure 1) are implemented in distinct modules and interconnected to solve for the global, coupled, dynamic response of a wind turbine system.  The new FAST modularization  framework  supports  module-independent  inputs,  outputs,  states,  and  parameters;  states  in continuous-time, discrete-time, and constraint form; loose and tight coupling; independent time and spatial discretizations; time marching, operating-point determination, and linearization; data encapsulation; dynamic allocation; and save/retrieve capability.  While not covered in this paper, modularization also establishes a basis for mixed-language programming, multicore processing, co-simulation across a network, hiding the details of individual model components to protect intellectual property (IP), and developing hardware-in-the-loop (HIL) simulation.  
+
+# A. Inputs, Outputs, States, and Parameters  
+
+Mathematical models of dynamic systems can be developed in terms of inputs, outputs, states, and parameters. Inputs (identified in this paper by $u$ ) are a set of values supplied to a system (i.e., module) and that, together with the states, are needed to calculate future states and/or the system’s output.  Outputs $(y)$ are a set of values calculated by and returned from a system and dependent on the system’s states, inputs, and/or parameters through output equations (with functions Y).  States are a set of internal values of a system influenced by inputs and/or time and used to calculate future state values and/or the system’s output.  There are three types of states.  Continuous states $(x)$ are states that are differentiable in time and characterized by continuous-time differential equations (with functions $X$ ).  Discrete states (with functions $x^{d}$ ) are states that only have a value at discrete steps in time and are characterized by discrete-time difference (recurrence) equations (with functions $X^{d}$ ).  Constraint states $(z)$ are states that are not differentiated or discrete (i.e., constraint states are algebraic variables) and are characterized by algebraic constraint equations (i.e., equations without time derivatives) (with functions $Z$ ).  Parameters $(p)$ are a set of internal system values, independent of the states and inputs, that can be fully defined at initialization (possibly with timedependence that can be fully prescribed at initialization) and characterize a system’s state equations (differential, difference, and/or constraint) and output equations.  
+
+As examples in the wind turbine application, structural displacements and velocities are typically represented as continuous states (structural accelerations are not states themselves, but are time derivatives of the velocity states), control-system logic and dynamic-stall formulations are often implemented with discrete states, and structural joints impose constraint states; BEM wake models and other quasi-static or quasi-steady formulations also involve constraint states.  For a rotor-aerodynamics module, examples of inputs and outputs include blade-element motion (position, orientation, and translational and rotational velocity) and blade-element loads (forces and moments), respectively.  Examples of parameters include structural definitions such as blade and tower length, mass, and stiffness and time-only-dependent environmental conditions such as wind inflow and incident waves not influenced by the structural response.  These examples, along with a few others, are summarized in Table 2.  
+
+Table 2.  Examples of inputs, outputs, states, and parameters in the wind turbine application.   
+
+
+<html><body><table><tr><td>Variable</td><td>Aerodynamics</td><td>Hydrodynamics</td><td>Controller</td><td>Structural Dynamics</td></tr><tr><td>Inputs</td><td>Turbine displacements · Turbine velocities</td><td>Substructuredisplacements ·Substructurevelocities</td><td>·Structural accelerations ·Reaction loads</td><td>Aerodynamicloads Hydrodynamicloads Controllercommands</td></tr><tr><td>· Outputs</td><td>·Aerodynamicloads</td><td>·Hydrodynamicloads</td><td>·Controllercommands</td><td>Displacements Velocities Accelerations ReactionLoads</td></tr><tr><td>Continuousstates</td><td>·Induction in GDW</td><td>·State-space-basedradiation "memory"</td><td>·Analog control signals</td><td>Displacements Velocities</td></tr><tr><td>·Discrete states</td><td>Beddoes-Leishman dynamic-stall states</td><td>Numerical-convolution- based radiation"memory"</td><td>·Digital control signals</td><td></td></tr><tr><td>Constraintstates</td><td>InductioninBEM</td><td></td><td></td><td>Constraintloads atjoints Quasi-staticmooringsystem</td></tr><tr><td>·Parameters</td><td>·Turbine geometry ·Staticairfoildata Undisturbedwindinflow</td><td>·Substructuregeometry Hydrodynamiccoefficients Undisturbedincidentwaves</td><td>·Controller gains Controllerlimits</td><td>Geometry Mass/inertia Stiffness coefficients Damping coefficients</td></tr></table></body></html>  
+
+The developer can choose a module’s specific inputs, outputs, states, and parameters in the new FAST modularization framework.  When choosing inputs for a module as part of the module development process, the only restriction is that a module’s inputs must be algebraically derivable from the available outputs of the modules coupled together (including, perhaps, from the module under development—a recursive formulation)—see Section II.D.  
+
+# B. Loose and Tight Coupling  
+
+Before its modularization was improved, FAST applied a loosely coupled time-integration scheme, where data (inputs and outputs) are exchanged between the modules at each coupling step, but where each module tracks its own states and integrates its own equations with its own solver.  Figure 6 illustrates the difference between loose  
+
+![](images/b9e565aa1d1b8e82b442d65847e3a90dc1e32482bd0e743bb11b97b0c1ab9708.jpg)  
+
+and tight-coupling schemes.  In a tightly coupled timeintegration scheme, each module sets up its own equations, but the states are tracked and integrated by a solver common to all of the modules.  The new FAST modularization allows for both loose and tight coupling.  
+
+Loose  coupling  is  convenient  for  introducing legacy code and can be quite computationally efficient because the choice of solver and time steps can be fitfor-purpose to the module.  But loose coupling can lead to numerical errors (e.g., drift) or numerical stability problems in the coupled solution in some cases. 11,§ These numerical problems can sometimes  
+
+be resolved through predictor-corrector-based loose coupling schemes.12  In the new FAST modularization framework, a given loosely coupled module must have a fixed coupling step, and the continuous-time and discretetime states within a given module must share this time step.  
+
+Tight coupling has numerical advantages over loose coupling if the common solver is appropriately suited for the problem; however, the modules have to be developed in a form amenable to tight coupling, which is less common among legacy tools and profoundly impacts how new modules must be developed.  Computational performance can be lost in tight coupling because the same (perhaps variable) time step must be applied to all continuous-time states of all interconnected modules.  In the new FAST modularization framework, separate modules can have different discrete-time steps, but all of the discrete-time states within a given module must share the same discrete time step.  In addition, in the new FAST modularization framework, tight coupling also permits operating-point determination and linearization—see Section II.F.  
+
+An initial assessment of the numerical stability, numerical accuracy, and computational performance of various coupling schemes is provided in a companion paper.  
+
+The loose or tight coupling of individual modules is achieved through the modular interface and coupler described in Section II.G.  Because of the modularization, it is also possible to isolate the dynamics of an individual  
+
+module in an uncoupled way.  A module in the new FAST modularization framework does not run by itself, but it is called by a separate driver program.  The uncoupled solution of a module intended for loose and tight coupling is illustrated in Figure 7.  
+
+![](images/bb1377abdc38eda71de02c266722ad221901b33d80cc23eef613d242d98e9777.jpg)  
+
+# C. Module Form  
+
+The mathematical formulation of a module permitted within the new FAST modularization framework is very general, making the overall framework extremely powerful and flexible enough to be considered for almost any system.  A general (need-not-be-linear) state-space formulation is considered with any combination of continuoustime-state, discrete-time-state, constraint-state, and output equations.  No assumption is made about the theory from which the state-space formulation was derived.  A system described by partial differential equations (PDEs) in space and time can be written in a general state-space formulation once the spatial dimensions have been discretized.  If no states are present, the system is often referred to as a “feed-forward only system.”  If both continuous and discrete states are present, the system is often referred to as a “sampled system” or “hybrid system.” If there are no constraint states, the continuous-time state equations form ordinary differential equations (ODEs), which have wellunderstood numerical solutions. With constraint states, the system is characterized by differential algebraic equations (DAEs) —that is, differential equations combined with algebraic constraint equations— which are much harder to solve.13  
+
+The tight-coupling feature of the new FAST modularization framework permits systems of the form of a hybrid semi-explicit DAE of index 1.  This form—as described in more detail below—has only the following limitations: (1) the continuous-time state derivatives and discrete-time state updates must be written as an explicit function of the states, inputs, and parameters and (2) the constraints must be of index 1.  These are the same limitations imposed within even the most advanced solvers available in the popular MATLAB/Simulink commercial computing package. 14 The most general form allowed by the framework in tight coupling is represented mathematically as  
+
+$$
+\begin{array}{r c l}{\dot{x}=X\Big(x,x^{d},z,u,t\Big)\,,}&\\ {x^{d}\big[n+I\big]=X^{d}\left(x\big|_{t=n\mathscr{A}}\,,x^{d}\left[n\right],z\big|_{t=n\mathscr{A}}\,,u\big|_{t=n\mathscr{A}}\,,t\big|_{t=n\mathscr{A}}\right),}&\\ {{\theta}=Z\Big(x,x^{d},z,u,t\Big)\,\,\mathrm{with}\,\left|\displaystyle\frac{\partial Z}{\partial z}\right|\ne\theta\,,\mathrm{and}}&\\ {y=Y\Big(x,x^{d},z,u,t\Big)\,.}&\end{array}
+$$  
+
+The continuous states, $x.$ , are defined by the explicit first-order ODEs of Eq. (1a) with the continuous-state functions, $X,$ on the right-hand side (RHS).  The continuous states are time dependent, so, $x(t)$ is implied by $x$ where $t$ is time and $\dot{x}$ is the first time derivative of $x$ .  The continuous-state equations are written in first-order form without loss of generality.   The continuous-state functions—as well as the other functions of Eq. (1)—can depend on the continuous states, $x,$ discrete states, $x^{d}$ , constraint states, $z$ , inputs, $\boldsymbol{u}$ , and time, $t$ (and, of course, are characterized by the parameters, $p$ , not shown).  Just as $x(t)$ is implied by $x$ , the time-dependency of $x^{d},z,u.$ , and $p$ is also implied.  (The direct time dependency of the continuous-state functions—as well as the other functions of Eq. (1)—results from the parameters having direct continuous-time-dependence, $p(t)$ , fully prescribed at initialization.) While $t$ is a scalar, it should be understood that $x,\,x_{\mathrm{~,~}}^{d}z,\,u$ , and $p$ may be one-dimensional arrays (vectors) of variables, each of different size. $X$ is a vector of the same size as $x\textbf{S O}$ that the number of equations matches the number of states.  
+
+The discrete states, $x^{d}$ , are defined by the explicit difference (recurrence) equations of Eq. (1b) from step $n$ to $n+l$ with the discrete-state functions, $\boldsymbol{X}^{d}.$ , on the RHS.  In this paper, square brackets $[]$ in functions denote discrete \*\*Second-order (or higher) form can easily be reduced to first-order form.  For example, if $q$ and $\dot{\boldsymbol{q}}$ are continuous states of the second-order system described by $\stackrel{..}{q}=Q\Big(q,\dot{q},x^{d},z,u,t\Big)$ where $\ddot{q}$ is the second time derivative of $q$ and $\boldsymbol{Q}$ are the second-order continuous-state equations, first-order form can be realized by using $x={\Biggl\{}{q\atop{\dot{q}}}{\Biggr\}}$ , with  
+
+$$
+X=\left\lbrace\begin{array}{c}{\dot{q}}\\ {Q\Bigl(q,\dot{q},x^{d},z,u,t\Bigr)\Bigr\rbrace.}\end{array}\right.
+$$  
+
+In the existing structural module of FAST, the equations of motion take on the generalized form of Newton’s second law, $M\big(\boldsymbol{q},\boldsymbol{u},t\big)\ddot{\boldsymbol{q}}=F\big(\boldsymbol{q},\dot{\boldsymbol{q}},\boldsymbol{u},t\big)$ , where $M$ is the generalized mass matrix and $F$ is the generalized force vector (no discrete or constraint states are included, but $M$ depends on the displacements and control inputs in general); in this case, $Q\big(q,\dot{q},u,t\big)\!=\!M^{-I}\big(q,u,t\big)F\big(q,\dot{q},u,t\big)$ .  
+
+time whereas round brackets () in functions denote continuous time; the brackets are dropped when implied.  The discrete state and discrete-state functions are evaluated only at discrete steps in time based on the fixed discrete time step (interval), $\varDelta t$ , which is greater than zero and the same for all discrete states of a given module.  While $n$ and $\varDelta t$ are scalars, $X^{d}$ is a vector of the same size as $x^{d}$ so that the number of equations matches the number of states. Because the discrete-state functions can depend on continuous-time variables, an analog-to-digital conversion (ADC) is required whereby the continuous-time variables are sampled (evaluated) at the discrete time step, as represented by $\Big|_{t=n\varDelta t}$ in Eq. (1b).  Likewise, the zero-order hold (ZOH) method of digital-to-analog conversion (DAC)  is  used  when  continuous-time  functions  depend  on  discrete  states. Mathematically, $x^{d}\left(t\right)=\sum_{n=-\infty}^{\infty}x^{d}\left[n\right]\!\left(H\left(t-n{\cal{A}}t\right)-H\left(t-\left(n+I\right){\cal{A}}t\right)\right)$ is  used  in  place  of $x^{d}[n]$ in  all  continuous-time functions, where $x^{d}(t)$ are the discrete states expressed in continuous time and are piecewise constant and $H$ is the Heaviside-step (unit-step) function; although the summation does not have to be implemented in practice, effectively, if $x^{d}$ exists, $x^{d}[n]$ is applied over $n\varDelta t<=t<\big(n+l\big)\varDelta t$ in all continuous-time functions.  Even if discrete states are present in a module, the constraint states, $z$ , inputs, $u$ , outputs, $y$ , and parameters, $p$ , are always expressed in continuous time.  (Because they are characterized by continuous-time parameters, $p$ , the discrete-state functions, $X^{d}$ , are expressed in continuous time even though they are evaluated at discrete time steps.)  
+
+The constraint (algebraic) states, $z.$ , are defined implicitly in continuous-time form by Eq. (1c) with the constraint-state (algebraic) functions, $Z.$ , on the RHS. $Z$ is a vector of the same size as $z$ so that the number of equations matches the number of states.  Equation (1c) can be understood in the following context: given the continuous states, $x$ , discrete states, $x^{d}$ , inputs, $u$ , and time, $t$ (and, of course, the parameters, $p$ , not shown), $Z$ can be used to solve for constraint states, $z,$ , at time $t$ for use in the other functions of Eq. (1).  In tight coupling, the constraints must be of index 1, which means that the constraint-state function must be invertible such that the constraint states could be written as an explicit function of the other states, inputs, and parameters, guaranteeing the local (but not global) existence and uniqueness of a solution.  Although the inverse of $Z$ with respect to $z$ does not have to be formulated in practice, it must exist and be bounded in a neighborhood around a solution. Mathematically, this requires that $\left|{\frac{\partial Z}{\partial z}}\right|\neq O$ , where $|\ |$ is used to represent the determinant of the Jacobian matrix $\frac{\partial Z}{\partial z}$ The requirement $\left|{\frac{\partial Z}{\partial z}}\right|\neq O$ also means that the matrix inverse of the Jacobian, $\left[\frac{\partial Z}{\partial z}\right]^{-I}$ exists and is bounded in a neighborhood around a solution, which will be used later.  For a system whose DOF result in a naturally higher constraint index (e.g., constraints imposed through Lagrange multipliers in multi-body dynamics often result in index-3 constraints), an index reduction technique must be applied to reduce the index to 1 before the system can be implemented within a tightly coupled module in the new FAST modularization framework.  In some systems, it is possible that alternate DOF will naturally result in a lower constraint index.  
+
+The outputs, $y$ , are defined explicitly in continuous time by Eq. (1d) with the output functions, $Y,$ on the RHS. Variables $y$ and $Y$ may be one-dimensional arrays (vectors), but of the same size so that the number of equations matches the number of states; the time-dependency of $y$ is also implied.  The output equations permit direct feedthrough of input.  The output functions need not depend explicitly on the first time derivative of the continuous states, $\dot{x}$ , because $\dot{x}$ itself depends on the same variables Y does.  
+
+Unlike in tight coupling, a loosely coupled module in the new FAST modularization framework is not limited to a semi-explicit DAE of index 1.  An even more general state-space formulation is available in loose coupling (e.g., an index-3 constraint is allowed), but the numerical solution (including time integration) of the loosely coupled equations is implemented by the module developer within the module and the overall solvability, numerical stability, and convergence of the coupled solution is not guaranteed (these can be verified in tight coupling—see Section II.D).  
+
+The semi-explicit DAE of index 1 formulation available in Eq. (1) for tight coupling is a subset of the more general formulation available in loose coupling. As such, it is  possible (likely desirable) in the FAST modularization framework to develop a given module for both loose and tight coupling when the system is expressible in the semi-explicit DAE of index 1 formulation of Eq. (1).  
+
+# D. Input-Output Transformations and Coupled Solution  
+
+The module interface and coupler described in Section II.G interconnects all of the individual modules together and drives the overall coupled solution forward.  One key role of the module interface and coupler is to derive the inputs to the individual modules from the available outputs of the modules coupled together.  
+
+Given $N$ total number of modules coupled together, the inputs, $\boldsymbol{u}$ , from each individual module can be combined into a single one-dimensional array (vector) across all coupled modules; likewise for the outputs, $y$ ,  
+
+$$
+\boldsymbol{u}=\left\{\begin{array}{c}{\boldsymbol{u}^{(I)}}\\ {\boldsymbol{u}^{(2)}}\\ {\vdots}\\ {\boldsymbol{u}^{(N)}}\end{array}\right\}\mathrm{~and~}\,\boldsymbol{y}=\left\{\begin{array}{c}{\boldsymbol{y}^{(I)}}\\ {\boldsymbol{y}^{(2)}}\\ {\vdots}\\ {\boldsymbol{y}^{(N)}}\end{array}\right\}.
+$$  
+
+In Eq. (2), superscript (i) is used to identify the $i^{\mathrm{th}}$ module and $u$ without a superscript denotes the combined vector of inputs across all coupled modules; the same convention is used for the outputs, $y$ .  In this paper, it should be clear from the context whether a variable refers to an individual module or the specific variable combined from all modules into a single vector.  Also in this paper, curly brackets $\{\}$ in arrays denote one-dimensional vectors whereas square brackets $[]$ in arrays denote two-dimensional matrices; the brackets are dropped when implied.  
+
+The input-output transformation equations used by the new FAST modularization framework in both loose and tight coupling are represented mathematically in their most general form as ( $\tilde{u}$ is defined following Eq. (4) below)  
+
+$$
+O=U\left(u,y,t\right)\,{\mathrm{with}}\,\left|{\frac{\partial U}{\partial{\tilde{u}}}}\right|\neq O\,.
+$$  
+
+The input-output transformation equations of Eq. (3) are algebraic (that is, without derivatives) and expressed implicitly in continuous time with the input-output transformation functions, $U.$ , on the RHS. $U$ is a vector of the same size as $\boldsymbol{u}$ so that the number of equations matches the number of inputs across all modules.  If two or more of the modules coupled together share some of the same exact inputs, the size of $u$ and $U$ can be reduced by that same number.  Equation (3) can be understood in the following context: given the outputs across all coupled modules, $y,$ and time, $t$ , $U$ can be used to solve for the inputs across all coupled modules, $u$ , at time $t,$ .  Of course, if modules have direct feedthrough of input to output, the solve can present a problem, as discussed below.  The input-output transformation functions must be defined specifically for each collection of modules coupled together.  Because no assumption is made about the theory from which each module was derived, it is possible to mix methodologies if the inputs and outputs are compatible.  Equation (3) implies that the input of any module is permitted to depend on any module’s output, including, perhaps, a module’s own output—a recursive formulation.  To clarify this point, it is useful to identify the input-output transformation equations for each individual module as shown in Eq. (4) below. (Equation (4) is equivalent to Eq. (3), but with input-output transformation for each module identified.)  It should be clear by Eq. (4) that while there are separate input-output transformation equations for each module, each transformation can depend on the inputs and outputs across all modules, and so cannot be solved independently from the other modules’ input-output transformations in general.  The input-output transformation functions may also depend on time-dependent parameters, implied by the dependency on time, $t$ , in Eqs. (3) and (4).  While the general form of the input-output transformation equations expressed in Eqs. (3) and (4) is powerful, it is fairly common for the transformation to be quite trivial.††  
+
+††For example, modules are often developed so that the input of one module equals the output of another, as implied by Table 2.  For the case of a two module system each with its input equaling the output of the other, the inputoutput transformation equations are $\left\{\!\!\begin{array}{c}{\!\!\displaystyle0\!\!\right\}\!=\!\!\left\{U^{(l)}\!\left(u,y,t\right)\!\right\}\!\!=\!\!\left\{u^{(l)}-y^{(2)}\right\}\!\!\!\!}\\ {\!\!\!\displaystyle\left\{U^{(2)}\!\left(u,y,t\right)\!\right\}\!\!\!}\end{array}\!\!\!\right.=\!\!\left\{u^{(l)}-y^{(l)}\right\}$ with ${\Bigg|}{\frac{\partial U}{\partial{\tilde{u}}}}{\Bigg|}={\Bigg[}{\begin{array}{l l}{I}&{\theta}\\ {\theta}&{I}\end{array}}{\Bigg]}=I\neq0$ and ${\frac{\partial U}{\partial y}}\!=\!\!{\left[\begin{array}{l l}{0}&{-I}\\ {-I}&{\;\,\theta}\end{array}\right]},$ , where $I$ is the identify matrix.  
+
+$$
+\left\{\begin{array}{l}{\displaystyle0}\\ {\displaystyle0}\\ {\displaystyle\vdots}\\ {\displaystyle0}\end{array}\right\}=\left\{\begin{array}{c}{\displaystyle U^{(I)}\left(u,y,t\right)}\\ {\displaystyle U^{(2)}\left(u,y,t\right)}\\ {\vdots}\\ {\displaystyle U^{(N)}\left(u,y,t\right)}\end{array}\right\}\,\mathrm{with}\,\left|\frac{\partial U}{\partial\tilde{u}}\right|\ne O
+$$  
+
+Similar to the restriction on constraint index in tight coupling, the input-output transformation functions must be invertible such that the inputs could $b e$ written as an explicit function of the outputs, guaranteeing the local (but not global) existence and uniqueness of a solution.  This restriction on the input-output transformation functions exists for both loose and tight coupling in the FAST modularization framework.  Although the function inverse does not have to be formulated in practice, it must exist and be bounded in a neighborhood around a solution. Mathematically, this requires that $\left|\frac{\partial U}{\partial\tilde{u}}\right|\neq0$ , where $\tilde{u}$ is a dummy variable representing $\boldsymbol{u}$ that is needed to clarify that the partial derivative of $U$ is with respect to the $\boldsymbol{u}$ explicitly identified in Eq. (3), or $\begin{array}{r}{\theta\!=\!U\!\left(\tilde{u},y,t\right)}\end{array}$ .  That is, the partial derivative of $U$ in Eq. (3) does not involve the chain rule with the $u$ that is included in the calculation of $y$ by Eq. (1d).  Making use of the chain rule, however, reveals that  
+
+$$
+\frac{\partial U}{\partial u}\!=\!\frac{\partial U}{\partial\widetilde{u}}\!+\!\frac{\partial U}{\partial y}\frac{\partial Y}{\partial u}\mathrm{~and~}\frac{\partial U}{\partial z}\!=\!\frac{\partial U}{\partial y}\frac{\partial Y}{\partial z},
+$$  
+
+which are used below.  
+
+Because inputs are derived from outputs and the new FAST modularization framework permits the output of individual modules to depend directly on their input (i.e., direct feedthrough of input to output), the input-output transformations themselves form algebraic constraint equations.  This condition is equivalent to what is referred to as “algebraic loops” in MATLAB/Simulink.  In some cases, it is possible that there is no global solution, no local solution, or no solution at all.  
+
+In a tightly coupled system, the existence of a solution can be checked, but this is not possible in a loosely coupled system because of  the  possibility  of  the  more  general  state-space  formulation available with loose coupling.  To clarify this point, again consider $N$ total number of modules coupled together, as illustrated in Figure 8, which follows the organization of Figure 6 (without integration shown), but which uses the nomenclature for the state, output, and input functions and inputs and outputs.  Just as the inputs and outputs from each individual module were combined across all modules in Eq. (2), the states, state functions, and output functions can be combined into vectors across all modules as well  
+
+![](images/c987be072f1ca829944c6a2788b8a3f80ca19d396c69c8c431da09fb655bc2ab.jpg)  
+Figure 8.  Coupling of $N$ modules.  
+
+$$
+x=\left\{\begin{array}{c}{\left[{x}^{(I)}\right]}\\ {\left[{x}^{(2)}\right]}\\ {\vdots}\\ {\left[{x}^{(N)}\right]}\end{array}\right\},\;x^{d}=\left\{\begin{array}{c}{{x}^{d(I)}}\\ {{x}^{d(2)}}\\ {\vdots}\\ {{x}^{d(N)}}\end{array}\right\},\;z=\left\{\begin{array}{c}{{z}^{(I)}}\\ {{z}^{(2)}}\\ {\vdots}\\ {{z}^{(N)}}\end{array}\right\},
+$$  
+
+$$
+\begin{array}{r}{X=\left\{X^{\left(I\right)}\right\},\;X^{d}=\left\{X^{d\left(I\right)}\atop{\begin{array}{c}{{X^{\left({{\boldsymbol{\it1}}}\right)}}}\\ {{X^{\left({{\boldsymbol{2}}}\right)}}}\\ {{\vdots}}\\ {{X^{\left(N\right)}}}\end{array}}\right\},\;X^{d}=\left\{X^{d\left(I\right)}\atop{\begin{array}{c}{{X^{d\left({{\boldsymbol{2}}}\right)}}}\\ {{X^{d\left({{\boldsymbol{2}}}\right)}}}\\ {{\vdots}}\\ {{X^{d\left(N\right)}}}\end{array}}\right\},\;Z=\left\{\begin{array}{c}{{Z^{\left(I\right)}}}\\ {{Z^{\left({{\boldsymbol{2}}}\right)}}}\\ {{\vdots}}\\ {{Z^{\left(N\right)}}}\end{array}\right\},\;\mathrm{and}\;\;Y=\left\{Y^{\left({{\boldsymbol{1}}}\right)}\atop{\begin{array}{c}{{Y^{\left({{\boldsymbol{2}}}\right)}}}\\ {{Y^{\left({{\boldsymbol{2}}}\right)}}}\\ {{\vdots}}\\ {{Y^{\left(N\right)}}}\end{array}}\right\}.}\end{array}
+$$  
+
+Substituting the last of Eq. (6b) and Eq. (1d) into Eq. (3), results in $\scriptstyle O\,=\,U\left(u,Y\left(x,x^{d},z,u,t\right),t\right)$ , meaning that—like the other functions of Eq. (1)—the input-output transformation functions, $U,$ , could be written as functions of the continuous states, $x$ , discrete states, $x^{d}$ , constraint states, $z$ , inputs, $u$ , combined across all modules and time, $t$ —that is, $O\!=\!U\!\left(x,x^{d},z,u,t\right)$ .  Taking this form of $U$ with Eqs. (2) and (6) together with Eq. (1) and grouping $Z$ and $U$ reveals that in tight coupling, the coupled solution of all modules forms a global hybrid semi-explicit DAE  
+
+$$
+\begin{array}{c}{\displaystyle\dot{x}=X\left(x,x^{d},z,u,t\right),}\\ {\displaystyle x^{d}\left[n+I\right]=X^{d}\left(x\Big|_{t=n\mathcal{A}},x^{d}\left[n\right],z\Big|_{t=n\mathcal{A}},u\Big|_{t=n\mathcal{A}},t\Big|_{t=n\mathcal{A}}\right),}\\ {\displaystyle\left\{\theta\right\}=\left\{Z\!\left(x,x^{d},z,u,t\right)\right\}_{\displaystyle\left[\begin{array}{c c}{\displaystyle\hat{\omega}Z}&{\displaystyle\hat{\omega}Z}\\ {\displaystyle\hat{\omega}U}&{\displaystyle\hat{\omega}U}\end{array}\right]\!\right\}\!=\!0,\mathrm{~and~}}\\ {\displaystyle y=Y\!\left(x,x^{d},z,u,t\right).}\end{array}
+$$  
+
+Although they appear nearly identical, Eq. (7) should not be confused with Eq. (1).  Equation (7) applies to the coupled solution of all modules in tight coupling whereas Eq. (1) applies only to an individual module in tight coupling.  Importantly, even if none of the individual modules themselves have internal constraint states (meaning that $z$ and $Z$ are absent from Eqs. (1) and (7)), the coupled solution of all modules still forms a DAE in tight coupling.  The inputs across all modules, $u$ , act as additional constraint states defined by the constraint-state (algebraic) equations (which are actually  the  input-output  transformation  equations)  in  Eq.  (7c)  with  the  input-output transformation functions, $U_{s}$ , on the RHS.  Effectively, the coupling illustrated in Figure 8 has been reduced in tight coupling to the system illustrated in Figure 9.  Please note that if the discrete time step (interval), $\varDelta t$ , differs between individual modules, separate discrete-state and discrete-state-function groupings are needed, even though they have been shown grouped together in Eq. (7).  
+
+![](images/65f87af27203e4d0882fe88ccb490d29a5f21a6e5615d12e7c78d0e227e8144c.jpg)  
+Figure 9.  Effective tightly coupled  
+
+To ensure that the global semi-explicit DAE of the coupled solution of all modules in  system. tight coupling has an index of 1, the condition identified by the determinant in Eq. (7c) is needed in addition to the conditions required of the determinants in Eqs. (1c) and (3).  Through the properties of determinants of block matrices, it can be shown that  
+
+$$
+\left[\frac{\partial Z}{\partial z}\quad\frac{\partial Z}{\partial u}\right]\!\!=\!\!\left|\frac{\partial Z}{\partial z}\right|\!\!\left|\frac{\partial U}{\partial u}\!-\!\frac{\partial U}{\partial z}\!\left[\frac{\partial Z}{\partial z}\right]^{\!\!-1}\!\frac{\partial Z}{\partial u}\right|,
+$$  
+
+which is nonzero by Eq. (7c).  The first determinant on the right of Eq. (8) is nonzero by Eq. (1c) (expanded across all coupled modules), which means that the second determinant on the right of Eq. (8) must also be nonzero.  The  
+
+matrix in the second determinant on the right of Eq. (8) will be identified with symbol $G$ from now on.  Inserting Eq. (5) into $G$ , it can thus be shown that the condition identified by the determinant in Eq. (7c) is equivalent to  
+
+$$
+\begin{array}{r}{G=\frac{\partial U}{\partial\tilde{u}}+\frac{\partial U}{\partial y}\left[\begin{array}{c c c c c}{\left[\displaystyle{\frac{\partial Y^{(i)}}{\partial u^{(i)}}-\frac{\partial Y^{(i)}}{\partial z^{(i)}}\left[\frac{\partial Z^{(i)}}{\partial z^{(i)}}\right]^{-1}\frac{\partial Z^{(i)}}{\partial u^{(i)}}}\right]}&{0}&{\cdots}&{0}\\ {0}&{\left[\displaystyle{\frac{\partial Y^{(2)}}{\partial u^{(i)}}-\frac{\partial Y^{(2)}}{\partial z^{(i)}}\left[\frac{\partial Z^{(i)}}{\partial{\hat{z}}^{(i)}}\right]^{-1}\frac{\partial Z^{(i)}}{\partial u^{(i)}}}\right]}&{0}&{\cdots}&{0}\\ {\vdots}&{\vdots}&{\ddots}&{\ddots}&{\vdots}\\ {0}&{0}&{0}&{\left[\displaystyle{\frac{\partial Y^{(N)}}{\partial u^{(N)}}-\frac{\partial Y^{(N)}}{\partial z^{(N)}}\left[\frac{\partial Z^{(N)}}{\partial{\hat{z}}^{(N)}}\right]^{-1}\frac{\partial Z^{(N)}}{\partial u^{(N)}}}\right]}\end{array}\right]\mathrm{wit}}\end{array}
+$$  
+
+which means that the matrix inverse $G^{-I}$ exists as will be used later.  
+
+The Jacobian matrices $\frac{\partial U}{\partial\tilde{u}}$ and $\frac{\partial U}{\partial y}$ from Eq. (9) are written out in terms of their relationships between individual modules in Eq. (10) below.  In general, these Jacobian matrices could be fully populated, but in practice they are likely quite sparse.  
+
+$$
+\frac{\partial U}{\partial\tilde{u}}=\left[\begin{array}{l l l l}{\displaystyle\frac{\partial U^{(l)}}{\partial\tilde{u}^{(l)}}}&{\displaystyle\frac{\partial U^{(l)}}{\partial\tilde{u}^{(2)}}}&{\cdots}&{\displaystyle\frac{\partial U^{(l)}}{\partial\tilde{u}^{(N)}}}\\ {\displaystyle\frac{\partial U^{(2)}}{\partial\tilde{u}^{(l)}}}&{\displaystyle\frac{\partial U^{(2)}}{\partial\tilde{u}^{(2)}}}&{\cdots}&{\displaystyle\frac{\partial U^{(2)}}{\partial\tilde{u}^{(N)}}}\\ {\vdots}&{\ddots}&\\ {\displaystyle\frac{\partial U^{(N)}}{\partial\tilde{u}^{(l)}}}&{\displaystyle\frac{\partial U^{(N)}}{\partial\tilde{u}^{(2)}}}&{\displaystyle\frac{\partial U^{(N)}}{\partial\tilde{u}^{(N)}}\right]{\mathrm{and~}}\frac{\partial U}{\partial y}=\left[\begin{array}{l l l l}{\displaystyle\frac{\partial U^{(l)}}{\partial y^{(l)}}}&{\displaystyle\frac{\partial U^{(l)}}{\partial y^{(2)}}}&{\cdots}&{\displaystyle\frac{\partial U^{(N)}}{\partial y^{(N)}}}\\ {\displaystyle\frac{\partial U^{(N)}}{\partial y^{(l)}}}&{\displaystyle\frac{\partial U^{(2)}}{\partial y^{(2)}}}&{\displaystyle\frac{\partial U^{(2)}}{\partial y^{(N)}}}\\ {\vdots}&{\ddots}&\\ {\displaystyle\frac{\partial U^{(N)}}{\partial y^{(l)}}}&{\displaystyle\frac{\partial U^{(N)}}{\partial y^{(N)}}}&{\displaystyle\frac{\partial U^{(N)}}{\partial\tilde{u}^{(N)}}.}\end{array}\right]
+$$  
+
+The condition identified by Eq. (9) determines the existence of a solution in tight coupling, and can be checked by the module interface and coupler, provided that the Jacobian matrices $\frac{\partial Y}{\partial u},\,\frac{\partial Y}{\partial z},\,\frac{\partial Z}{\partial z}$ , and an $\frac{\partial Z}{\partial u}$ are known for each individual module and $\frac{\partial U}{\partial\tilde{u}}$ and $\frac{\partial U}{\partial y}$ are known.  The important role that direct feedthrough of input to output plays in determining whether a solution exists should be clear by Eq. (9). $^{\ddagger\ddagger}$ Of course, it is best to check the  
+
+$$
+\begin{array}{r l}&{\stackrel{\mathrm{tesuple}}{=}\left[\begin{array}{c c c}{\mathrm{the}}&{\mathrm{two}}&{\mathrm{module}}&{\mathrm{example}}\\ {I}&{-\left[\frac{\hat{\sigma}Y^{(2)}}{\hat{\sigma}u^{(2)}}-\frac{\hat{\sigma}Y^{(2)}}{\hat{\sigma}z^{(2)}}\Bigg[\frac{\hat{\sigma}Z^{(2)}}{\hat{\sigma}z^{(2)}}\Bigg]^{-I}\frac{\partial Z^{(2)}}{\hat{\sigma}u^{(2)}}\right]\right]_{\boldsymbol{\cdot}}}\\ {-\left[\frac{\hat{\sigma}Y^{(I)}}{\hat{\sigma}u^{(I)}}-\frac{\hat{\sigma}Y^{(I)}}{\hat{\sigma}z^{(I)}}\Bigg[\frac{\hat{\sigma}Z^{(I)}}{\hat{\sigma}z^{(I)}}\Bigg]^{-I}\frac{\hat{\sigma}Z^{(I)}}{\hat{\sigma}u^{(I)}}\right]}&{I}\end{array}\right]\stackrel{\mathrm{~,~}}{=}\left[\begin{array}{c c c}{\mathrm{~of~}}&{\mathrm{~of~}}\\ {I}&{-\left[\frac{\hat{\sigma}Y^{(2)}}{\hat{\sigma}u^{(2)}}-\frac{\hat{\sigma}Y^{(2)}}{\hat{\sigma}z^{(2)}}\Bigg[\frac{\hat{\sigma}Z^{(I)}}{\hat{\sigma}}\Bigg]^{-I}\frac{\partial Z^{(I)}}{\partial t^{(2)}}-\frac{\hat{\sigma}Y^{(2)}}{\hat{\sigma}z^{(2)}}\Bigg[\frac{\hat{\sigma}Z^{(2)}}{\hat{\sigma}z^{(2)}}\Bigg]^{-I}\frac{\partial Z^{(2)}}{\hat{\sigma}u^{(2)}}\right]\right]\,\boldsymbol{\cdot}}\\ {-\left[\frac{\hat{\sigma}Y^{(I)}}{\hat{\sigma}u^{(I)}}-\frac{\hat{\sigma}Y^{(I)}}{\hat{\sigma}z^{(I)}}\Bigg[\frac{\hat{\sigma}Z^{(I)}}{\hat{\sigma}z^{(I)}}\Bigg]^{-I}\frac{\partial Z^{(I)}}{\hat{\sigma}u^{(I)}}\right]\left[\frac{\hat{\sigma}Y^{(2)}}{\hat{\sigma}u^{(2)}}-\frac{\hat{\
+$$  
+
+condition identified by Eq. (9) in the module-development process if possible, rather than waiting for the solution process.  To minimize the potential that a solution does not exist, it is best to develop modules such that the output does not have direct feedthrough from at least one of the modules coupled together.§§  When a solution does exist, the tight coupling solvers being developed for the new FAST modularization framework can solve the coupled DAE of Eq. (7) robustly (again, limited to index 1).  
+
+In loose coupling, the more general formulation permitted prohibits the use of Eq. (9).  For example, in an index3 formulation, the determinant of the Jacobian matrix $\frac{\partial Z}{\partial z}$ equals zero, meaning that the matrix inverse of the Jacobian, $\left[\frac{\partial Z}{\partial z}\right]^{-I}$ , used in Eq. (9) cannot be calculated.  In reality, there is no reason to form the combined system arrays of Eqs. (6) through (10) in loose coupling.  Instead, the new FAST modularization framework in loose coupling uses a root-finding algorithm to solve Eq. (3), but it is not possible to check for the existence of a solution in the process.  The numerical problems related to loose coupling can sometimes be resolved through predictorcorrector-based loose coupling schemes,  as discussed in Section II.B.  
+
+Descriptions of the root-finding algorithms being developed for the loose-coupling feature of FAST and of the index-1 DAE time-integration schemes being developed for the tight-coupling feature of FAST are outside the scope of this paper.  
+
+Some CAE tools introduce a time delay between the input and output to avoid the complications related to the input-output transformations (e.g., to use values of the output from the prior time step to derive the inputs at the current time step).  This is a common solution approach, but this approach may introduce numerical errors that may adversely affect the accuracy and stability of the coupled solution and is not the approach taken in the new FAST modularization framework.  
+
+# E. Independent Time and Spatial Discretizations  
+
+In the new FAST modularization framework, a system’s (i.e., module’s) inputs and outputs are allowed to (but need not) reside on a discretized spatial boundary characterizing the outer extent of the system; likewise, the states (continuous, discrete, and/or constraint) and parameters are allowed to (but need not) reside within a system’s discretized domain.  Before its modularization was improved, FAST required that the spatial discretization of interface boundaries in the aerodynamic, hydrodynamic, and structural modules be identical.  In the new modularization framework, independent spatial discretizations are allowed.  Allowing each module to use its own appropriate discretization will greatly improve the computational efficiency and provide more flexibility.  Using too coarse  a  discretization  reduces  solution  accuracy,  and  using  too  fine  a  resolution  reduces  computational performance.  Finer discretizations are needed in areas of significant property or response gradient, such as mass and stiffness variations for a structural model or the exponential decay of hydrodynamic loads with depth for a hydrodynamic model.  Figure 10 illustrates the mapping of independent structural and hydrodynamic discretizations.  
+
+A library of spatial elements, operations on those elements, and functions to map between meshes of different discretizations has been developed based on the isoparametric formulations popular in finite-element analysis (FEA).15  The mesh library allows for varying spatial dimension in motions and loads, including point (lumped, e.g., rigid bodies and concentrated loads), line (one-dimensional, e.g., beams and forces per unit length), surface (twodimensional, e.g., shells and pressure forces), and volume (three-dimensional, e.g., solids and body forces) discretizations.  The mapping allows for the discretizations to conform to boundaries moving due to, for example, structural deflection or variations of the fluid surface.  
+
+![](images/38abe392fcb44b2c49996716c6ebc1df78afa95d08ad7ef4d5513280f7804ef4.jpg)  
+Figure  10.    Mapping  independent  structural  and hydrodynamic discretizations.  
+
+When module inputs, $\boldsymbol{u}$ , and/or outputs, $y$ , reside on a discretized spatial boundary characterizing the outer extent of the system, the mesh library will be used by the modular interface and coupler in the formulation  of  the  input-output  transformation functions, $U_{s}$ , in Eq. (3).  The mathematical details of the mesh library are outside the scope of this paper.  
+
+Likewise,  the  new  FAST  modularization framework allows for distinct time steps between individual modules. As identified in Sections II.B and II.C, there are three types of time steps in the FAST modularization framework: (1) the discrete time step (interval), Δt, (2) the (fixed) coupling step of a loosely coupled  module,  and  (3)  the  common  (perhaps variable) time step used to integrate the continuoustime states of all interconnected modules in tight coupling.  In loose coupling, time-step types (1) and  
+
+(2) must be equal within a given loosely coupled module although separate modules may have different steps.  In tight coupling, separate modules may have different discrete time steps.  
+
+When coupled modules have different time steps, the modular interface and coupler will interpolate and extrapolate the module inputs and outputs in time.  This is done to ensure that the input-output transformation functions, $U,$ , are solved at a given time (as discussed in Section II.D), enabling modules to be called at appropriate times.  The mathematical details of this time-based interpolation and extrapolation are outside the scope of this paper.  
+
+# F. Time Marching, Operating-Point Determination, and Linearization  
+
+The primary purpose of FAST is to perform time-domain simulations of the aero-hydro-servo-elastic response of wind-energy systems.  Mathematically, the coupled system equations form an initial value problem (IVP) whereby the response of the system can be found in time if the parameters of all modules are known for all time, $p(t)$ , and initial values (i.e., initial conditions (ICs)) are given for the states of all modules.  To clarify, ICs need to be provided only for the continuous states, $x(O)$ , and discrete states, $x^{d}[O]$ .  The initial values of constraint states, inputs, and outputs can be derived from these ICs and the parameters, but the solution is aided by initial guesses for the constraint states, $z^{G u e s s}(O)$ , and inputs, $u^{G u e s s}(O)$ .  The new FAST modularization framework supports this timemarching calculation with both loose and tight coupling.  
+
+FAST’s new tight-coupling feature, however, permits two additional types of calculations.  The first of the new calculations is operating-point (OP, or fixed-point) determination.  Several types of OPs can be found, including static equilibrium (constant displacement), steady state (constant velocity), and periodic steady state (periodic variation in response).  These OPs can be found with or without trim of inputs to achieve a desired condition.  Time marching can be performed from given ICs or from an OP.  The mathematical details of the OP calculation are outside the scope of this paper.  The second of the new calculations is linearization of the underlying nonlinear system equations, which is valuable for full-system modal analysis (e.g., determining natural frequencies, damping, and mode shapes), linear-system-based controls design (e.g., developing linear state-space representations of a wind turbine plant), and linear-system-based stability analysis.  Linearization can be performed about an OP defined by given ICs, a given time in the time-marching process, or an OP found through the OP calculation discussed above. Before its modularization was improved, only the structural module of FAST could be linearized (without states in other modules).  With the new formulation, the OP and linearization calculations can take place across the entire coupled system (including aerodynamics, hydrodynamics, servo dynamics, and structural dynamics)—see Section III for the mathematical details.  OP and linearization calculations are not available in loose coupling because the loosely coupled state-space formulation does not have the limitations of tight coupling, and more general formulations may not be suitable for OP and linearization calculations.  
+
+# G. Module Interface and Coupler  
+
+Before its modularization was improved, the structural module of FAST also functioned to couple the aerodynamic, hydrodynamic, and control and electrical system dynamics modules together.  This functionality meant that the coupled solution was dictated by the solution of the structural module and made it difficult to make modifications to the structural module.  The new FAST modularization framework introduces the module interface and coupler—also known as the “glue” code—that is distinct from individual modules.  The module interface and coupler’s role is to interconnect all of the individual modules, algebraically derive inputs from outputs as discussed in Section II.D (including mapping between different spatial discretizations and time interpolation and extrapolation as discussed in Section II.E), and drive the overall coupled solution forward.  In tight coupling, the module interface and coupler has the added tasks of integrating the coupled system equations using one of its own solvers and driving the OP and linearization calculations when those options are selected.  
+
+All modules are intended to interface directly to the module interface and coupler as shown in Figure 3, Figure 6, and Figure 8.  A given module can, itself, be further modularized into separate modules that are interfaced directly with each other if and only if they behave collectively and interface with the module interface and coupler as an individual module would.  When these conditions cannot be met or are too cumbersome to implement, the separate modules should be interfaced directly to the module interface and coupler.  There is no limit to the number of modules, $N.$ , that the module interface and coupler can interconnect.  
+
+# H. Data Encapsulation and Dynamic Allocation  
+
+The new FAST modularization framework also supports data encapsulation.  Specifically, the new framework requires that there be no global variables, so, no actual data (whether inputs, outputs, states, or parameters) are stored inside  the modules.  Instead, the data are stored in the driver  (main) program—in this case the module interface and coupler—using  data structures defined in the modules, and the data are passed to/from  the modules through subroutine arguments.  Access to the data is restricted as much as possible in the Fortran programming language.  
+
+Data  encapsulation has two main consequences. First, as previously mentioned, the inputs must be algebraically  derivable  from the available outputs of the modules coupled together; a module  cannot  access the states or parameters of another  module  unless  they have been copied as outputs. Data transfer between modules  is  now  clear,  an improvement  from  previous versions  of  FAST  where global data led to complicated interactions between modules—what is sometimes  
+
+referred  to  as  “spaghetti  code”—that  made  reading  and maintaining of the source code difficult.  Second, dynamic allocation is possible, such that multiple instances of a module can exist simultaneously.  The ability to dynamically allocate modules is of tremendous value.  Dynamic allocation will assist NREL’s coupling of FAST with the OpenFOAM CFD tool for modeling  multiple  turbines  in  a  wind  farm,  including  the modeling  of  wake  and  array  effects  and  their  aeroelastic interaction,  as  illustrated  in  Figure  11. (Before  its modularization was improved, FAST could only be used to model the dynamics of a single turbine.)  Dynamic allocation is also currently being used to implement a new nonlinear beam FEA-based structural model focused only on the dynamics of a single blade but with the interaction of all (two, three, or more) blades of a rotor included as part of the coupled solution; the concept is illustrated in Figure 12.  
+
+![](images/03d0f2b805ed4c070672f1e5bf4deef2cd01c0cff0b2cb6b133578aa826d77a4.jpg)  
+Figure 11.  Coupled simulation between FAST and OpenFOAM.16  
+
+![](images/06504bdf30e30ee49d99b055fdbe292580ff88d94bc657e555db5461c2eff1a2.jpg)  
+Figure 12.  From one blade to an entire rotor.  
+
+# I.  Save/Retrieve Capability  
+
+The new FAST modularization framework also supports the ability to save the data across all modules at a given instance during the course of a simulation.  The data can be written to a file and retrieved later to continue the simulation starting from that point (i.e., restart).  This save/retrieve feature is especially useful in the solution of computationally expensive systems to minimize the performance impact of HPC disruptions and/or the rerunning of common segments of the simulation.  
+
+# III. Linearization  
+
+As discussed in Section II.F, the tight coupling functionality of the new FAST modularization framework supports OP determination and linearization of individual modules as well as of the overall coupled system; the linearization functionality is not available in loose coupling.  This section summarizes the mathematical details of the linearization functionality for tight coupling, which is important to its understanding and proper application.  The linearization mathematics is also illustrative of the coupling properties inherent in coupled system modeling within the new FAST modularization framework whether the underlying modules are fundamentally linear in nature or not.  
+
+# A. Linearization of a Module  
+
+A linear representation of a nonlinear system model is valid only for small deviations (perturbations) from an OP.  As discussed in Section II.F, an OP in the new FAST modularization framework can be defined by given ICs, a given time in the time-marching process, or an OP found through the OP determination calculation.  Assuming OP values are given for the continuous states, $x\big|_{o p}$ , discrete states, $\left.x^{d}\right|_{o p}$ , inputs, $u\Big|_{o p}$ , and time, $t{\Big|}_{o p}$ , of a tightly coupled module, Eqs. (1a), (1c), and (1d) can be used to calculate the OP values of the first time derivative of the continuous state, $\dot{x}\big|_{o p}$ , constraint states, $z\Big|_{o p}$ , and outputs, ${y\Big|}_{o p}$ .  In this paper, $\big|_{o p}$ denotes an OP value of a variable or the evaluation of a function about the OP.  Except for the OP time, $t{\Big|}_{o p}$ , each of these variables can be perturbed (represented by $\varDelta$ ) about their respective OP values  
+
+$$
+\begin{array}{r l}&{\left.\vphantom{\int}{\d x}=x\right|_{o p}+\varDelta x\,,\,\left.x^{d}=x^{d}\,\right|_{o p}+\varDelta x^{d}\,,\,u=u\big|_{o p}+\varDelta u\,,}\\ &{\left.\dot{x}=\dot{x}\right|_{o p}+\varDelta\dot{x}\,,\,\left.z=z\right|_{o p}+\varDelta z\,,\mathrm{and}\,\,y=y\big|_{o p}+\varDelta y\,.}\end{array}
+$$  
+
+Substituting the expressions of Eq. (11) into Eq. (1), expanding as a Taylor-series approximation, and keeping only the linear terms (i.e., neglecting products of perturbations), it can be shown that the most general linearized form of Eq. (1) is represented mathematically as  
+
+$$
+\begin{array}{c}{{\varDelta\dot{x}=A\varDelta x+A^{D d C}\varDelta x^{d}+B A u\,,}}\\ {{\varDelta x^{d}\left[n+l\right]=A^{A D C}\left.\varDelta x\right|_{t=n\varDelta t}+A^{d}\varDelta x^{d}\left[n\right]+B^{A D C}\left.\varDelta u\right|_{t=n\varDelta t},\mathrm{and}}}\\ {{\varDelta y=C\varDelta x+C^{D A C}\varDelta x^{d}+D\varDelta u\,,}}\end{array}
+$$  
+
+where  
+
+$$
+\begin{array}{r}{A\!=\!\left[\!\!\left[\frac{\partial{\cal X}}{\partial x}\!-\!\frac{\partial{\cal X}}{\partial z}\!\left[\frac{\partial{\cal Z}}{\partial z}\right]^{-I}\!\frac{\partial{\cal Z}}{\partial x}\!\right]\!\right]_{o p},}\\ {A^{D d^{C}}\!=\!\left[\!\!\left[\frac{\partial{\cal X}}{\partial x^{d}}\!-\!\frac{\partial{\cal X}}{\partial z}\!\left[\frac{\partial{\cal Z}}{\partial z}\right]^{-I}\!\frac{\partial{\cal Z}}{\partial x^{d}}\!\right]\!\right]_{o p},}\end{array}
+$$  
+
+$$
+B=\left[{\frac{\partial X}{\partial u}}\!-\!{\frac{\partial X}{\partial z}}\!\left[{\frac{\partial Z}{\partial z}}\right]^{-I}\!{\frac{\partial Z}{\partial u}}\right]_{o p},
+$$  
+
+$$
+A^{A D C}=\left[\frac{\partial X^{d}}{\partial x}\!-\!\frac{\partial X^{d}}{\partial z}\!\left[\frac{\partial Z}{\partial z}\right]^{\!-\!I}\frac{\partial Z}{\partial x}\right]_{\!\!\!o p},
+$$  
+
+$$
+\begin{array}{r l}&{A^{\prime}=\left[\frac{\partial X^{\prime}}{\partial X^{\prime}}-\frac{\partial X^{\prime}}{\partial\zeta}\left[\widehat{\mathcal{Q}}\right]^{\prime}\frac{\partial Z^{\prime}}{\partial X^{\prime}}\right]_{\varphi_{\varphi}},}\\ &{B^{\prime\prime}=\left[\frac{\partial X^{\prime}}{\partial u}-\frac{\partial X^{\prime}}{\partial\zeta^{2}}\left[\widehat{\mathcal{Q}}\right]^{\prime}\frac{\partial Z^{\prime}}{\partial u}\right]_{\varphi_{\varphi}},}\\ &{\ C=\left[\widehat{\mathcal{Q}}\right]_{\infty}^{\mathcal{N}}-\widehat{\mathcal{Q}}^{\Gamma}\left[\widehat{\mathcal{Q}}\right]^{\prime}\frac{1}{\widehat{\mathcal{Q}}}\mathscr{L}\right]_{\varphi_{\varphi}},}\\ &{C^{\partial\epsilon}=\left[\frac{\partial Y^{\prime}}{\partial\epsilon^{2}}-\frac{\partial Y^{\prime}}{\partial Z^{\prime}}\left[\widehat{\mathcal{Q}}\right]^{\prime}\frac{\partial Z^{\prime}}{\partial\widehat{\mathcal{Q}}^{\epsilon}}\right]_{\varphi_{\varphi}},}\\ &{D=\left[\frac{\partial Y^{\prime}}{\partial\epsilon}-\frac{\partial Y^{\prime}}{\partial Z^{\prime}}\left[\widehat{\mathcal{Q}}\right]^{\prime}\frac{\partial Z^{\prime}}{\partial\widehat{\mathcal{Q}}^{\epsilon}}\right]_{\varphi_{\varphi}}.}\end{array}
+$$  
+
+The perturbations of continuous states, $\varDelta x$ , are defined by the linear continuous-state equations of Eq. (12a) expressed as explicit first-order ODEs, where the RHS is the linearized form of the nonlinear continuous-state functions, $X,$ of Eq. (1a). $A$ is the continuous-state matrix, $A^{D A C}$ is the discrete-state matrix, and $B$ is the input matrix of the linearized continuous-state equations. $A$ is a square matrix with the number of rows and columns equal to the number of continuous states.  The number of rows in $A^{D A C}$ equals the number of continuous states, and the number of columns equals the number of discrete states.  The number of rows in $B$ equals the number of continuous states, and the number of columns equals the number of inputs.  These matrices, as well as the other matrices of Eq. (12), are time invariant and depend on the Jacobians of the functions in Eq. (1) as shown in Eq. (13) and discussed below. Like the continuous states themselves, the perturbations of continuous states are time dependent, so, $\varDelta x(t)$ is implied by $\varDelta x$ where $t$ is time and $\varDelta\dot{x}$ is the first time derivative of $_{\varDelta x}$ .  And like the continuous-state equations themselves, the linear continuous-state equations are written in first-order form without loss of generality.  
+
+\*\*\*In the second-order example of footnote \*\*, the linearized form of a second-order system without discrete and constraint states, $\ddot{q}=Q\big(q,\dot{q},u,t\big)$ , is $\varDelta\ddot{q}=\frac{\partial Q}{\partial q}\Bigg\vert_{o p}\,\varDelta q+\frac{\partial Q}{\partial\dot{q}}\Bigg\vert_{o p}\,\varDelta\dot{q}+\frac{\partial Q}{\partial u}\Bigg\vert_{o p}\;\varDelta u$ .  The linearized first-order form of this system can be realized by using $\varDelta x={\binom{A q}{A{\dot{q}}}}$ , with $\varDelta\dot{x}=\left\{\varDelta\dot{q}\right\},\;\,A=\left[\dot{\sigma}\dot{Q}\right]_{o p}\;\;\left.\frac{\partial Q}{\partial\dot{q}}\right|_{o p},$ and  
+
+$$
+\boldsymbol{B}=\left[\begin{array}{c}{\boldsymbol{\theta}}\\ {\hat{\displaystyle{\sigma}}\boldsymbol{Q}}\\ {\hat{\displaystyle{\sigma}}\boldsymbol{u}}\end{array}\right].
+$$  
+
+The perturbations of discrete states, $\varDelta x^{d}$ , are defined by the linear discrete-state equations of (12b) expressed as explicit difference (recurrence) equations from step $n$ to $n+l$ , where the RHS is the linearized form of the nonlinear discrete-state functions, $X^{d}.$ , of Eq. (1b).  Matrices $A^{A D C}$ , $A^{d}.$ , and $B^{A D C}$ are the continuous-state, discretestate, and input matrices of the linearized discrete-state equations, respectively.  While $A^{d}$ is a square matrix with the number of rows and columns equal to the number of discrete states, the number of rows in $A^{A D C}$ and $B^{A D C}$ equals the number of discrete states, the number of columns in $A^{A D C}$ equals the number of continuous states, and the number of columns in $B^{A D C}$ equals the number of inputs.  Like the discrete states and discrete-state functions themselves, the perturbations of discrete states and linearized discrete-state equations are evaluated only at discrete steps in time based on the fixed discrete time step (interval), $\varDelta t$ , which is greater than zero and the same for all discrete states of a given module.  The ADC and DAC discussed in Section II.C also apply to the linearized system.  (The superscripts ADC and DAC in the linearized system matrices indicate analog-to-digital and digital-to-analog conversions, respectively.)  The OP time, $t{\Big|}_{o p}$ , need not be an integer multiple of the discrete time step, $\varDelta t.$ .  Even if discrete states are present in a module, the perturbations of inputs, $\varDelta u$ , and outputs, $\varDelta y$ , are always expressed in continuous time.  
+
+The perturbations of outputs, $\varDelta y$ , are defined by the linear output equations of (12c) expressed explicitly in continuous time, where the RHS is the linearized form of the nonlinear output functions, $Y,$ , of Eq. (1d).  Matrices $C$ ,  
+
+For the generalized form of Newton’s second law used by the existing structural module of FAST, $Q\big(q,\dot{q},u,t\big)\!=\!M^{-I}\big(q,u,t\big)F\big(q,\dot{q},u,t\big)$ ,  it  can  be  shown  that ${\frac{\partial Q}{\partial q}}\biggr|_{o p}=\left[M^{-I}\left[{\frac{\partial F}{\partial q}}{-{\frac{\partial M}{\partial q}}{M^{-I}}F}\right]\right]_{o p},$ $\left.\frac{\partial\mathcal{Q}}{\partial\dot{q}}\right|_{o p}=\left[M^{-I}\frac{\partial F}{\partial\dot{q}}\right]\biggr|_{o p}\,,\,\mathrm{and}\;\frac{\partial\mathcal{Q}}{\partial u}\biggr|_{o p}=\left[M^{-I}\left[\frac{\partial F}{\partial u}\!-\!\frac{\partial M}{\partial u}M^{-I}F\right]\right]_{o p}$ .  The standard generalized linear form of Newton’s second law is typically written as $M\varDelta\ddot{q}+C\varDelta\dot{q}+K\varDelta q=F^{u}\varDelta u$ , where $M$ is the generalized linear mass matrix, $C$ is the generalized linear damping matrix (not to be confused with $C$ from Eq. (13g)), $K$ is the generalized linear stiffness matrix, and $F^{u}$ is the generalized linear input forcing matrix.  Relating these matrices to the  J ${\begin{array}{l}{\displaystyle{\mathrm{acobians~~of~}}\ Q,\ {\mathrm{~it~}}\ {\mathrm{is~~~clear~}}\ {\mathrm{that}}\ {\mathrm{~}}M=M{\Big|}_{o p},\ \ {}C=-{\frac{\partial F}{\partial{\dot{q}}}}{\Bigg|}_{o p},\ \ {}K=-{\left[{\frac{\partial F}{\partial q}}-{\frac{\partial M}{\partial q}}M^{-{\prime}}F\right]}}\\ {\displaystyle\left[{\frac{\partial F}{\partial u}}-{\frac{\partial M}{\partial u}}M^{-{\prime}}F\right]_{o p}.\ {\mathrm{~It~may~be~surprising~that~}}K{\mathrm{~doesn"t~equal~}}-{\frac{\partial F}{\partial q}}{\Bigg|}_{o p}\ {\mathrm{and~}}F^{\prime}\,{\mathrm{doesn"t~equal~}}}\end{array}}$ and l $\left.\frac{\partial F}{\partial u}\right|_{o p}$ which are common expressions.  Noticing that $\left\{M^{-I}F\right\}\biggr|_{o p}=\ddot{q}\bigr|_{o p}$ , it is seen that—because the mass matrix depends on the displacements and control inputs (whereby $\frac{\partial M}{\partial q}\Bigg|_{o p}\neq0$ and $\frac{\hat{\sigma}M}{\hat{\sigma}u}\Bigg|_{o p}\neq0$ in general)—the stiffness and input forcing matrices are impacted by an OP that is not the static-equilibrium or steady-state condition (whereby $\left.\ddot{q}\right|_{o p}\neq0\,.$ ).  While the linear model is still valid for the OP that is not the static-equilibrium or steady-state condition about which the model was linearized, it is of less practical use than when $K=-\frac{\partial F}{\partial q}\Biggr|_{o p}$ and $F^{u}=\frac{\partial F}{\partial u}\bigg|_{o p}$ .  As such, it is usually important for the OP to be a static-equilibrium or steady-state condition (whereby $\left.\ddot{q}\right|_{o p}=0\,,\;K=-\frac{\partial F}{\partial q}\right|_{o\mu}$ , and Fu= $F^{u}=\frac{\partial F}{\partial u}\Bigg|_{o p}).$  
+
+$C^{D A C}$ , and $D$ are the continuous-state, discrete-state, and input-transmission matrices of the linearized output equations, respectively.  The number of rows in C, $C^{D A C}$ , and $D$ equals the number of outputs, and the number of columns in $\bar{C_{*}}\;\bar{C}^{D A C}.$ , and $D$ equals the number of continuous states, discrete states, and inputs, respectively.  The input-transmission matrix, $D$ , is present if the module has direct feedthrough of input to output.  
+
+All matrices of Eq. (12) depend on the Jacobians of the functions in Eq. (1), as shown in Eq. (13).  To be able to form the linearized system, a tightly coupled module must be able to calculate and return the OP values of the 16 $\begin{array}{l}{\mathrm{obian~matrices},\thinspace\thinspace\displaystyle\frac{\partial X}{\partial x}\Bigg|_{\sigma_{p}},\thinspace\thinspace\frac{\partial X}{\partial x^{d}}\Bigg|_{\sigma_{p}},\thinspace\frac{\partial X}{\partial z}\Bigg|_{\sigma_{p}},\thinspace\frac{\partial X}{\partial u}\Bigg|_{\sigma_{p}},\thinspace\frac{\partial X^{d}}{\partial x}\Bigg|_{\sigma_{p}},\thinspace\frac{\partial X^{d}}{\partial x^{d}}\Bigg|_{\sigma_{p}},\thinspace\frac{\partial X^{d}}{\partial z}\Bigg|_{\sigma_{p}},\thinspace\frac{\partial X^{d}}{\partial z}\Bigg|_{\sigma_{p}},\thinspace\frac{\partial X^{d}}{\partial u}\Bigg|_{\sigma_{p}},\thinspace\frac{\partial Z}{\partial x}\Bigg|_{\sigma_{p}},\thinspace\frac{\partial Z}{\partial x}\Bigg|_{\sigma_{p}},\thinspace\frac{\partial Z}{\partial x}\Bigg|_{\sigma_{p}},\thinspace}\\ {\thinspace\frac{\gamma}{\partial\sigma},\thinspace\frac{\partial Z}{\partial u}\Bigg|_{\sigma_{p}},\thinspace\frac{\partial Y}{\partial x}\Bigg|_{\sigma_{p}},\thinspace\frac{\partial Y}{\partial x^{d}}\Bigg|_{\sigma_{p}},\thinspace\frac{\partial Y}{\partial z}\Bigg|_{\sigma_{p}},\thinspace\mathrm{and}\thinspace\frac{\partial Y}{\partial u}\Bigg|_{\sigma_{p}}.\thinspace\thinspace\mathrm{These~Jacobian~matrices\mathrm{-}~a n d~t h u s~a l l~m a t r i c e s},}\end{array}$ ∂Z Jac ∂x ∂Z trices of Eq. ∂ (12)—can be functions of the OP values of the continuous states, $x\big|_{o p}$ , discrete states, $\left.x^{d}\right|_{o p}$ , constraint states, $z\Big|_{o p}$ , inputs, $u\Big|_{o p}$ , and time, $t\Big|_{o p}$ , (and, of course, are characterized by the parameters evaluated at the OP time, $p\Big(t\vert_{o p}\Big)\,)$ , but are themselves time invariant.  (If the nonlinear system of Eq. (1) was characterized with timedependent parameters, $p(t)$ , the parameters are effectively replaced with time-invariant parameters, $p\big(\ t\big|_{o p}\big)$ , in the linearized system of Eq. (12).)  Module developers can choose to compute these Jacobian matrices in the module either numerically (e.g., through a central-difference perturbation technique) or analytically.  The analytical approach is preferred when practical because it is more accurate than numerical approaches, whose accuracy is dictated by the perturbation size of the numerical technique.  
+
+The constraint-state (algebraic) equations have been eliminated from the linearized system of Eq. (12) because once linearized, the constraint-state equations can be easily solved for the perturbations of constraint states, $\varDelta z$ .  The perturbations of constraint states are then substituted into the remaining equations of Eq. (12), essentially eliminating $\varDelta z$ as a separate variable.  By eliminating the constraint-state equations, each of the matrices of Eq. (12) depend on Jacobians that are taken with respect to the constraint states and of the Jacobians of the constraint-state functions as shown in Eq. (13).  The existence of the matrix inverse of the Jacobian of the constraint-state functions with respect to constraint states, $\left[\frac{\partial Z}{\partial z}\right]^{-I}$ , resulted from the limitation in tight coupling that the constraints must be of index 1, as discussed in Section II.C.  
+
+Equation (12) is readily identifiable as a general state-space representation of a linear time-invariant (LTI) system with a combination of continuous and discrete states (i.e., a “sampled” or “hybrid” LTI system).  A linearized system with continuous states but without discrete states reduces to a standard continuous LTI state-space model characterized by matrices $A,\,B,\,C_{5}$ , and $D$ .  A linearized system with discrete states but without continuous states reduces to a discrete LTI state-space model, but with continuous inputs and outputs, characterized by matrices $A^{d},B^{A D C}$ , $C^{D A C}$ , and $D$ .  
+
+It is usually important for the OP to be a static-equilibrium or steady-state condition.  An OP that is not the static-equilibrium or steady-state condition may have an undesirable effect on the linear system matrices, as discussed in footnote  .  But, if the system implemented in a module was naturally linear to begin with, the linearization process will simply result in the same linear system regardless of the OP.  
+
+# B. Linearization of the Overall Coupled System  
+
+Once an OP has been determined across all coupled modules and each individual module has been linearized as discussed in Section III.A, linearization of the overall coupled system is possible.  Given $N$ total number of linearized modules coupled together, the perturbations of inputs and outputs can be combined into vectors across all modules, just as the inputs and outputs from each individual module were combined across all modules in Eq. (2),  
+
+$$
+\begin{array}{r}{\varDelta u=\left\{\begin{array}{c}{\left.A u^{(I)}\right.}\\ {\left.A u^{(2)}\right.}\\ {\left.\vdots\quad\right\}}\\ {\left.A u^{(N)}\right]}\end{array}\right.\mathrm{and}\;\;\varDelta y=\left\{\begin{array}{c}{\left.A y^{(I)}\right.}\\ {\left.A y^{(2)}\right.}\\ {\left.\vdots\quad\right\}}\\ {\left.A y^{(N)}\right\}}\end{array}\right\}.}\end{array}
+$$  
+
+Similar to how Eq. (1) was linearized to form Eq. (12), the input-output transformation functions, $U.$ of Eq. (3) can be linearized  
+
+$$
+{\cal O}=\frac{\hat{\sigma}U}{\hat{\sigma}\tilde{u}}\Bigg\vert_{o p}\,{\cal A}u+\frac{\hat{\sigma}U}{\hat{\sigma}y}\Bigg\vert_{o p}\,{\cal A}y\,\,\,\mathrm{with}\,\left\vert\frac{\hat{\sigma}U}{\hat{\sigma}\tilde{u}}\right\vert_{o p}\ne0\,,
+$$  
+
+where  the  Jacobian  matrices $\frac{\partial U}{\partial\tilde{u}}$ and $\frac{\partial U}{\partial y}$ written out in Eq. (10)—have been evaluated at the OP and are time invariant. Equation (15) can be used to solve for the perturbations  of  inputs  given the perturbations of outputs across all coupled modules.  
+
+It is convenient to see the effect of externally provided inputs on the linear coupled system model.  For example, the linear model of the coupled system could be used as a linear plant model from which to develop an advanced linear state-spacebased  controller, where  the  linear  plant model includes the coupled dynamics of wind  turbine  aero-elastics  plus  the dynamics of a baseline controller; it would be useful in this case to see the effect of additional  control  inputs  on  top  of  the baseline control signals in the linear system response.  To accommodate this feature, the input perturbations derived from Eq.  
+
+![](images/8dda2db6d9c0d4e61f2f474aef5907faa96e089f02eea2048d6e82b9ad6af3b6.jpg)  
+Figure 13.  Coupling of $N$ linearized modules.  
+
+(15) are augmented with additional (but not quantified) input perturbations, $\varDelta u^{+}$ , before being sent to each module. The concept is illustrated in Figure 13, which follows the organization of Figure 8 but uses the nomenclature of the matrices  from  the  linearized  state,  output,  and  input  equations  and  perturbations  of  inputs  and  outputs. Mathematically,  
+
+$$
+\varDelta u=\varDelta\tilde{u}+\varDelta u^{+}\,,
+$$  
+
+where $\varDelta\tilde{u}$ is a dummy variable representing the input perturbations derived from Eq. (15), $\varDelta u^{+}$ are the additional input perturbations, and $\varDelta u$ are the actual input perturbations sent to each module.  All of the input perturbations of Eq. (16) are vectors combining the perturbations across all coupled modules.  The perturbations of states can also be combined into vectors across all modules  
+
+$$
+\begin{array}{r}{\varDelta u^{+}=\left\{\begin{array}{c}{\displaystyle\mathcal{A}u^{+(l)}}\\ {\displaystyle\mathcal{A}u^{+(2)}}\\ {\vdots}\\ {\displaystyle\mathcal{A}u^{+(N)}}\end{array}\right\},\;\varDelta x=\left\{\begin{array}{c}{\displaystyle\varDelta x^{(l)}}\\ {\displaystyle\varDelta x^{(2)}}\\ {\vdots}\\ {\displaystyle\varDelta x^{(N)}}\end{array}\right\},\;\mathrm{and}\;\varDelta x^{d}=\left\{\begin{array}{c}{\displaystyle\varDelta x^{d(l)}}\\ {\displaystyle\varDelta x^{d(2)}}\\ {\vdots}\\ {\displaystyle\varDelta x^{d(N)}}\end{array}\right\}.}\end{array}
+$$  
+
+Combining Eqs. (14) through (17) with the linearized module equations of Eq. (12) and using the module coupling identified in Figure 13 yields the linearized system model of the complete coupled system  
+
+$$
+\begin{array}{c}{{\left.\varDelta\dot{x}=A\varDelta x+A^{D d C}\varDelta x^{d}+B\varDelta u^{+}\right.,}}\\ {{\left.\qquad\quad\varDelta x^{d}\left[n+l\right]=A^{A D C}\left.\varDelta x\right|_{t=n\varDelta t}+A^{d}\varDelta x^{d}\left[n\right]+B^{A D C}\left.\varDelta u^{+}\right|_{t=n\varDelta t},\mathrm{and}}}\\ {{\varDelta y=C\varDelta x+C^{D A C}\varDelta x^{d}+D\varDelta u^{+},}}\end{array}
+$$  
+
+# where  
+
+$$
+\begin{array}{r l}{\left[\begin{array}{l l l l}{\rho(\rho-\phi_{1}-\phi_{2})}&&&\\ &{\rho+\phi_{3}}&&\\ &{\rho+\phi_{4}}&&\\ &{\rho+\phi_{5}}&&\\ &{\rho+\phi_{6}}&&\\ &&{\rho+\phi_{7}}&&\\ &&{\rho+\phi_{8}}&&\end{array}\right]\left[\begin{array}{l l l l}{\rho(\rho-\phi_{1}-\phi_{2})}&&&{\rho}\\ &{\phi_{1}-\rho(\rho)}&&{\rho}\\ {\rho}&{\rho-\phi_{6}}&&\\ {\rho}&{\rho}&{\rho}&{\rho}\end{array}\right]\left[\begin{array}{l l l l}{\rho(\phi_{1}-\phi_{2}-\phi_{3})}&&&{,}\\ &{\phi_{2}-\phi_{3}}&&{,}\\ &{\rho(\rho-\phi_{3}-\phi_{4})}&&{,}\\ {\rho(\rho-\phi_{4}-\phi_{5})}&&&{\rho-\phi_{6}}\end{array}\right]}\\ {\left[\begin{array}{l l l l}{\rho(\rho-\phi_{1}-\phi_{3})}&&&{,}\\ &{\rho(\rho-\phi_{2}-\phi_{3})}&&\\ &{\rho(\rho-\phi_{4}-\phi_{5})}&&\\ &{\rho(\rho-\phi_{5}-\phi_{6})}&&\\ &{\rho(\rho-\phi_{6}-\phi_{7})}&&\end{array}\right]\left[\begin{array}{l l l l}{\rho(\rho-\phi_{1}-\phi_{3})}&&&{,}\\ &{\rho(\rho-\phi_{1}-\phi_{4})}&&{,}\\ {\rho(\rho-\phi_{3}-\phi_{4})}&&&{\rho(\rho-\phi_{4}-\phi_{5})}\end{array}\right]\left[\begin{array}{l l l l}{\rho(\rho-\phi_{1}-\phi_{2}-\phi_{3})}&&&{,}\\ &{\rho(\rho-\phi_{3}-\phi_{4})}&&{,}\\ {\rho(\rho-\phi_{4}-\phi_{5}-\phi_{4})}&&&{\rho(\rho-\phi_{5}-\phi_{6})}\end{array}\right]}\\ {\left[\begin{array}{l l l l}{\rho(\rho-\phi_{1}-\phi_{3})}&&&{,}\\ & 
+$$  
+
+$$
+\begin{array}{r l}&{C^{\scriptscriptstyle{D M C}}=\left[\begin{array}{c c c c}{C^{\scriptscriptstyle{D M C}(i)}}&{\theta}&{\cdots}&{\theta}\\ {\theta}&{C^{\scriptscriptstyle{D M C}(2)}}&{\theta}&{\theta}\\ {\vdots}&{\ddots}&{\ddots}&{\ddots}\\ {\theta}&{\theta}&{C^{\scriptscriptstyle{D M C}(N)}}\end{array}\right]-\left[\begin{array}{c c c c}{D^{(i)}}&{\theta}&{\cdots}&{\theta}\\ {\theta}&{D^{(2)}}&{\theta}&{\theta}\\ {\vdots}&&{\ddots}&\\ {\theta}&{\theta}&&{D^{(N)}}\end{array}\right]\left[G\right]_{\scriptscriptstyle{\varphi}}\right]^{-i}\frac{\partial U}{\partial y}\Bigg|_{\varphi}\left[\begin{array}{c c c c}{C^{\scriptscriptstyle{D M C}(i)}}&{\theta}&{\cdots}&{\theta}\\ {\theta}&{C^{\scriptscriptstyle{D M C}(2)}}&&{\theta}\\ {\vdots}&&{\ddots}&\\ {\theta}&{\theta}&&{C^{\scriptscriptstyle{D M C}(N)}}\end{array}\right],\;a^{\scriptscriptstyle{D M}}}\\ &{D=\left[\begin{array}{c c c c}{D^{(i)}}&{\theta}&{\cdots}&{\theta}\\ {\theta}&{D^{(2)}}&&{\theta}\\ {\vdots}&&{\ddots}&\\ {\theta}&{\theta}&{D^{(N)}}\end{array}\right]\left[G\right]_{\scriptscriptstyle{\varphi}}\right]^{-i}\frac{\partial U}{\partial\tilde{u}}\Bigg|_{\varphi}\;.}\end{array}
+$$  
+
+Equation (18) is the linearized form of the nonlinear coupled system equations of Eq. (7) with additional input perturbations, $\boldsymbol{\varDelta u}^{+}$ .  Similar to how the constraint-state equations have been eliminated from the linearized system of Eq. (12), the inputs that act as additional constraint states in Eq. (7c) have been eliminated from Eq. (18).  This was possible because, once linearized, the additional constraint-state equations (i.e., the linearized input-output transformation equations) can be easily solved, essentially eliminating $\varDelta\tilde{u}$ as separate variables.  Eliminating $\varDelta\tilde{u}$ causes each of the matrices of Eq. (19) to depend on the matrix $G\Big\vert_{o p}$ , which is the matrix $G$ from Eq. (9) evaluated at the OP,  
+
+$$
+G\Bigr|_{o p}=\frac{\partial U}{\partial\tilde{u}}\Bigr|_{o p}+\frac{\partial U}{\partial y}\Biggr|_{o p}\left[\begin{array}{c c c c}{D^{(I)}}&{0}&{\cdots}&{0}\\ {0}&{D^{(2)}}&&{0}\\ {\vdots}&&{\ddots}&\\ {0}&{0}&&{D^{(N)}}\end{array}\right]\mathrm{with}\,\left|G\right|_{o p}\Bigr|\neq0\,.
+$$  
+
+Although they appear nearly identical, Eq. (18) should not be confused with Eq. (12).  Equation (18) applies to the coupled linear solution of all modules whereas Eq. (12) applies only to an individual linearized module. Effectively, the coupling illustrated in Figure 13 has been reduced to the system illustrated  
+
+in Figure 14.  Please note that if the discrete time step (interval), Δt, differs between individual modules, separate discrete-state-perturbation and linear discrete-state-equation groupings are needed, even though they have been shown grouped together in Eq. (18).  
+
+Equation (19) shows that all matrices of Eq. (18) are derived from the linearized inputoutput transformation functions evaluated at the OP from Eq. (15) and from the linearized matrices of each individual module from Eq. (13).  So, once all individual modules have been linearized, the linearized system model of the complete coupled system can be assembled.  All matrices of Eq. (18) are time invariant.  
+
+![](images/49b1113b4c971855509fd473c771195d706988aaccf0e14a9fcc5110cd0e7031.jpg)  
+Figure 14.  Effective coupled linear system.  
+
+Like Eq. (12), Eq. (18) is readily identifiable as a general state-space representation of an LTI system with a combination of continuous and discrete states (i.e., a “sampled” or “hybrid” LTI system).  A coupled linearized system with continuous states but without discrete states reduces to a standard continuous LTI state-space model characterized by matrices $A,\,B,\,C,$ , and $D$ .  A coupled linearized system with discrete states but without continuous states reduces to a discrete LTI state-space model, but with continuous inputs and outputs, characterized by matrices $A^{d},B^{d D C},C^{D\bar{A}C}$ , and $D$ .  
+
+The matrices of the linearized coupled system model are very illustrative of the coupling properties inherent in coupled system modeling.  While many of the matrices in Eq. (19) are block-diagonal, matrices $\left[G\right|_{o p}]^{-I},\left.\frac{\partial U}{\partial\tilde{u}}\right|_{o p},$  
+
+and $\frac{\partial U}{\partial y}\rvert_{c}$ may be full in general, so, the matrices of the linearized coupled system model may also be full in o   
+general.  The module-to-module coupling is apparent.††† −D( 2)   
+$\dagger\dagger\dagger\dagger_{\mathrm{In}}$ the two module example of footnotes †† and ‡‡, G with D( 1)D( 2) ≠0. D 1 I   
+The matrices of the linearized coupled system model for this two module example are given below and clearly show   
+off-diagonal terms and coupling between the two modules:  
+
+$$
+\begin{array}{r l}{{\cal A}=\left[\begin{array}{l}{\dot{\mathcal{A}}^{(\mu+\frac{1}{2})\bar{\mathcal{A}}^{(\mu)}}[\mathcal{I}-\mathcal{B}^{(\mu)}]^{2}\left[\mathcal{A}^{(\mu)}\left[\mathcal{A}^{(\nu)}\left[\mathcal{F}^{(\mu+\frac{1}{2})}\mathcal{F}^{(\nu)}\right]^{2}\right]^{-\ell(\nu)}}\\ {~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}\mathcal{A^{A}(\nu-\mu^{-\nu})^{\alpha}\left[\mathcal{F}^{(\mu)}\left[-{D}^{(\nu)}\right]^{2}\right]^{\frac{1}{2}}\mathcal{P}^{(\nu)\alpha}}\\ {~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}\end{array}\right]\left[\begin{array}{l}{\dot{A}^{(\mu)}[\mathcal{F}^{(\nu)}]^{\alpha^{\alpha}}[-\mathcal{B}^{(\mu)}]^{2}[-\mathcal{F}^{(\nu)}]^{2}\right]^{-\ell(\nu)^{\alpha}}}\\ {~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}\end{array}\right]\left[\begin{array}{l}{\dot{A}^{(\mu)}[\mathcal{F}^{(\nu)}]^{2}[-\mathcal{B}^{(\mu)}]^{2}[-\mathcal{F}^{(\nu)}]^{2}\right]^{\frac{1}{2}}\mathcal{P}^{(\nu)\alpha^{\alpha}}}\\ {~~~~~~~~~~~~~~~~~ 
+$$$$
+\begin{array}{r l}&{a=\left[{\begin{array}{l}{\rho^{10}+\Delta^{11}\rho^{10}\left({I-2.979}\right)^{-1}\rho^{11}}\end{array}}\right]^{-1}\rho^{11}\rho^{11}\left({I-2.979}\right)^{-1}\Bigg[\left({I-2.979}\right)^{-1}\rho^{11}\Bigg]^{-1},}\\ &{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}\rho^{12}\left[\left({I-2.979}\right)^{-1}\rho^{12}\right]^{-1},}\\ &{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}\rho^{-1}\left({I-2.97}\right)^{-1}\rho^{-1}\left({I-2.979}\right)^{-1},}\\ &{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}\rho^{-1}\left({I-3.91}\right)^{-1}\rho^{-1}\nabla^{12}\rho^{-1}\left({I-3.931}\right)^{-1}\rho^{-1}}\\ &{~~~~~~~~~~~~~~~~~ 
+$$  
+
+In the example of footnote §§ where the output of one of the modules is the conjugate quantity of the output of the other module, $\boldsymbol{D}^{(l)}=\boldsymbol{0}$ , so all terms in the equations above involving $\left[I-D^{(l)}D^{(2)}\right]^{-I}$ reduce to $I.$ and the  
+
+The continuous-state matrix for the linearized continuous-state equations of the coupled system, $A$ , depends on the continuous-state and input matrices for the linearized continuous-state equations from each individual module, the continuous-state and input-transmission matrices for the linearized output equations from each individual module, and the linearized input-output transformation functions evaluated at the OP, as shown in Eq. (19a) and Eq. (20).  These dependencies are similar for the discrete-state matrix for the linearized continuous-state equations of the coupled system, $A^{D A C}$ ; the continuous-state and discrete-state matrices for the linearized discrete-state equations of the coupled system, $A^{A D C}$ and $A^{d}.$ , respectively; and the continuous-state and discrete-state matrices for the linearized output equations of the coupled system, $C$ and $C^{D A C}$ , respectively.  The input matrix for the linearized continuousstate equations of the coupled system, $B$ , depends on the input matrix for the linearized continuous-state equations from each module, the input-transmission matrices for the linearized output equations from each individual module, and the linearized input-output transformation functions evaluated at the OP, as shown in Eq. (19c) and Eq. (20). These dependencies are similar for the input matrix for the linearized discrete-state equations of the coupled system, $B^{A D C}$ , and input-transmission matrix for the linearized output equations of the coupled system, $D$ .  It is interesting to note that the input-transmission matrices of each individual module, $D$ , impact all matrices of the linearized coupled system as illustrated in footnote †††, further highlighting the importance of the role played by direct feedthrough of input to output in the coupled system response.  
+
+# IV. Conclusion  
+
+NREL recently has put considerable effort into improving the overall modularity of its FAST wind turbine aerohydro-servo-elastic tool to (1) improve the ability to read, implement, and maintain source code; (2) increase module sharing and shared code development across the wind community; (3) improve numerical performance and robustness; and (4) greatly enhance flexibility and expandability to enable further developments of functionality without the need to recode established modules.  The new FAST modularization framework supports moduleindependent inputs, outputs, states, and parameters; states in continuous-time, discrete-time, and constraint form; loose and tight coupling; independent time and spatial discretizations; time marching, operating-point determination, and linearization; data encapsulation; dynamic allocation; and save/retrieve capability.  This paper explains the features of the new framework, as well as the concepts and mathematical background needed to understand and apply them correctly.  Table 3 summarizes the features of the new FAST modularization framework dependent on whether the modules are loosely or tightly coupled.  
+
+It is envisioned that the new modularization framework will transform FAST into a powerful, flexible, and robust 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.  
+
+At its core, the new FAST modularization framework is a means by which various mathematical systems are implemented in distinct modules and interconnected to solve for the global, coupled, dynamic response of a system. While the intent was to establish a framework to better model wind turbines, in reality, the new FAST modularization framework is quite general and could be applied to systems other than wind turbines, as long as the mathematical models of those systems are implemented in modules developed under the same framework.  
+
+# Future Work  
+
+Work is ongoing to convert the existing modules of FAST, including AeroDyn and HydroDyn, to the format required of the new modularization framework, including the data structures and interface procedures.  The loose coupling functionality and mesh library will be completed first, followed by the tight coupling functionality for time  
+
+<html><body><table><tr><td>matrices A()</td><td>of the</td><td></td><td>linearized</td><td></td><td>coupled</td><td>system</td><td>model</td><td>simplify</td><td></td><td>greatly</td><td>as</td><td>follows:</td></tr><tr><td colspan="2" rowspan="4">A=</td><td rowspan="4">+ B()D(2)c(l) B(2)c(l)</td><td rowspan="4">B()C(2) A(2)</td><td rowspan="4"></td><td rowspan="4">ADAC</td><td colspan="2" rowspan="4">ADAC(1) + B()D(2)CDAC(I)</td><td rowspan="4">B()CDAC(2)</td><td rowspan="4">ADAC(2)</td><td rowspan="4">B(l) B= 0</td><td rowspan="4"></td><td rowspan="4">B()D(2)] B(2)</td></tr><tr><td>B(2)CDAC(1) Ad(1)</td><td>BADC(1)CDAC(2)</td></tr><tr><td rowspan="2">+ BADC(I) D(2)C(l)</td><td rowspan="2"></td><td rowspan="2">+ BADC(1) D(2)CDAC(1)</td></tr><tr><td rowspan="2">ADC(1) A AADC BADC</td><td rowspan="2">BADC(1)C(2) AADC(2)</td></tr><tr><td>BADC(2)C(l) BADC(1) BADC(1) D(2) 0 BADC(2)</td><td>C =</td><td>C() D(²)c(l)</td><td>0 C(2)</td><td>Ad D(2)CDAC(1)</td><td>BADC(2)CDAC(1) CDAC(1) 0 CDAC(2)</td><td>and D</td><td>0</td><td>Ad(2) 0 0 D(2)</td></tr></table></body></html>  
+
+Table 3.  Features of loose and tight coupling.   
+
+
+<html><body><table><tr><td>Features</td><td>Loose</td><td>Tight</td></tr><tr><td colspan="3">Module-lndependentVariables</td></tr><tr><td>·Inputs</td><td></td><td></td></tr><tr><td>·Outputs</td><td></td><td>√</td></tr><tr><td>·Parameters</td><td>√</td><td>√</td></tr><tr><td>· Continuous states</td><td></td><td>√</td></tr><tr><td>· Discrete states</td><td></td><td>√</td></tr><tr><td>· Constraint states</td><td>√</td><td>√</td></tr><tr><td>System Formulation</td><td>√</td><td></td></tr><tr><td colspan="3">Explicit continuous-time ODEs</td></tr><tr><td></td><td>√</td><td>√</td></tr><tr><td>· Explicit discrete-time updates</td><td></td><td>√</td></tr><tr><td>·Constraint equations of index 1</td><td></td><td>√</td></tr><tr><td>Output equations with direct feedthrough</td><td></td><td>√</td></tr><tr><td>Semi-explicit DAEs of index 1</td><td>√</td><td></td></tr><tr><td>·Systems of any form</td><td></td><td></td></tr><tr><td colspan="3">Independent Spatial Discretizations</td></tr><tr><td>·Available</td><td></td><td></td></tr><tr><td colspan="3">Operating-Point Determination</td></tr><tr><td>· Static equilibrium</td><td></td><td>√</td></tr><tr><td>· Steady state</td><td></td><td>√</td></tr><tr><td>· Periodic steady state</td><td></td><td>√</td></tr><tr><td>· With trim of inputs</td><td></td><td>√</td></tr><tr><td colspan="3">Linearization</td></tr><tr><td>· About given initial conditions</td><td></td><td></td></tr><tr><td>· About given time</td><td></td><td>√</td></tr><tr><td>· About operating point</td><td></td><td>√</td></tr><tr><td colspan="3">Time Marching</td></tr><tr><td>·From given initial conditions</td><td></td><td></td></tr><tr><td>·From operating point</td><td></td><td></td></tr><tr><td>· Independent time steps for continuous states between modules</td><td></td><td></td></tr><tr><td>·Independent time steps for discrete states between modules</td><td></td><td></td></tr><tr><td colspan="3">Solution</td></tr><tr><td>·Solver implementation is up to the module developer</td><td></td><td></td></tr><tr><td>·Solveris selectable from those available in the glue</td><td></td><td></td></tr><tr><td>· Overall solvability, numerical stability, and convergence verifiable</td><td></td><td>√</td></tr><tr><td colspan="3">Data Encapsulation and No Global Data</td></tr><tr><td>· Required</td><td></td><td></td></tr><tr><td colspan="3">DynamicAllocation of InstancesofModules</td></tr><tr><td>·Available</td><td></td><td></td></tr><tr><td colspan="3">Save/Retrieve Capability</td></tr><tr><td>·Available</td><td></td><td>√</td></tr></table></body></html>  
+
+marching,  OP  determination,  and linearization,  respectively. The recently initiated assessment of the numerical stability, numerical accuracy, and computational performance  of various  coupling schemes12 will also continue.  
+
+The development of  new framework-compatible  modules  of higher fidelity is ongoing and will continue in the future.  Near- and long-term  developments  include implementing higher-fidelity models of the wind inflow (e.g., based on site-specific measurements), aerodynamics (e.g., vortex-wake and dynamic meandering wake (DWM) models), hydrodynamics (e.g., multimember support-structure hydrodynamics, high-order wave and loading theories, and ice loading), control and  electrical system dynamics (e.g., Type 1-4 generator topologies, deformable trailing edges, and wind farm control) and structural dynamics (e.g., multi-member support structures; mooring dynamics;  blades  with  composite cross  sections, precurve and presweep,  large deflection, and torsion; and drivetrain dynamics).  
+
+# Acknowledgments  
+
+The  contributions  of  Bonnie Jonkman,  John  Michalakes,  Mike Sprague, and Amir Gasmi of NREL to the overall development of the new FAST modularization framework are gratefully  
+
+acknowledged.  Specifically, Bonnie Jonkman led the development of the programmer’s handbook and source-code template that support the framework.  The source-code template is supported by a registry for automatic source-code generation developed by John Michalakes.  John Michalakes and Mike Sprague led the development of a library supporting independent spatial discretizations.  Mike Sprague and Amir Gasmi led the assessment of loose and coupling schemes.  Thanks also to Kathryn Ruckman and Bonnie Jonkman of NREL for editing this paper to make it much more readable.  
+
+This work was performed at NREL in support of the U.S. Department of Energy under contract number DEAC36-08GO28308 and under a project funded through topic area 1.1 of the Funding Opportunity Announcement (FOA) number DE-FOA-0000415.  
+
+# References  
+
+1Jonkman, J. M. and Buhl Jr., M. L. FAST User’s Guide. NREL/EL-500-38230. Golden, CO: National Renewable Energy Laboratory, August 2005.  
+
+2Jonkman, J. FAST Theory Manual. NREL/TP-500-32449. Golden, CO: National Renewable Energy Laboratory (to be published).  
+
+3Laino, D. J. and Hansen, A. C. User’s Guide to the Wind Turbine Dynamics Aerodynamics Computer Software AeroDyn. Salt Lake City, UT: Windward Engineering LLC. Prepared for the National Renewable Energy Laboratory under Subcontract No. TCX-9-29209-01, December 2002. 4Moriarty, P. J. and Hansen, A. C. AeroDyn Theory Manual. NREL/EL-500-36881. Golden, CO: National Renewable Energy Laboratory, December 2005. 5Jonkman, 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; NREL/TP-500-41958. Golden, CO: National Renewable Energy Laboratory. 6Jonkman, J. M. “Dynamics of Offshore Floating Wind Turbines—Model Development and Verification.” Wind Energy. Vol. 12, No. 5, July 2009, pp. 459-492; NREL/JA-500-45311. Golden, CO: National Renewable Energy Laboratory; DOI: 10.1002/we.347. 7Web page: http://wind.nrel.gov/designcodes/ (accessed June 1, 2012). 8Laino, D. J. and Hansen, A. C. User’s Guide to the Computer Software Routines AeroDyn Interface for ADAMS®. Salt Lake City, UT: Windward Engineering LLC. Prepared for the National Renewable Energy Laboratory under Subcontract No. TCX-9- 29209-01, September 2001. 9van Garrel, A. Development of a Wind Turbine Aerodynamics Simulation Module. ECN-C--03-079. Petten, The Netherlands: Energy research Centre of the Netherlands, August 2003. 10Jonkman, B. J.; Michalakes, J.; Jonkman, J. M.; Buhl Jr., M. L.; Platt, A.; and Sprague, M. A. NWTC Programmer’s Handbook: A Guide for Software Development within the FAST Computer-Aided Engineering Tool. Golden, CO: National Renewable Energy Laboratory (to be published). 11Cordle, A. and Jonkman, J. “State of the Art in Floating Wind Turbine Design Tools.” The Proceedings of the Twenty-First (2011) International Offshore and Polar Engineering Conference, 19–24 June 2011, Maui, HI [CD-ROM]. Vol. 1, pp. 367-374. 2011-TPC-1059. Cupertino, CA: International Society of Offshore and Polar Engineers (ISOPE), June 2011; NREL/CP-500- 50543. Golden, CO: National Renewable Energy Laboratory. 12Gasmi, A.; Sprague, M. A.; Jonkman, J. M.; and Jones, W. B. “Numerical Stability and Accuracy of Temporally Coupled Multi-Physics Modules in Wind-Turbine CAE Tools.” $5I^{s t}$ AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, 7–10 January 2013, Grapevine (Dallas/Ft. Worth Region), TX. Reston, VA: American Institute of Aeronautics and Astronautics, January 2013 (to be published); NREL/CP-500-57298. Golden, CO: National Renewable Energy Laboratory. 13Bendtsen, C. and Thomsen, P. G., Numerical Solution of Differential Algebraic Equations.  IMM-REP-1999-8. Lyngby, Denmark: Technical University of Denmark, May 1999. 14Shampine, L. F.; Reichelt, M. W.; and Kierzenka, J. A., “Solving Index-1 DAEs in MATLAB and Simulink.” SIAM Review, Vol. 41, No. 3, 1997, pp. 538-552. 15Felippa, C. A. Introduction to Finite Element Methods, Boulder, CO: University of Colorado, Fall 2004. 16Lee, S.; Churchfield, M.; Moriarty, P.; Jonkman, J.; and Michalakes, J. “Atmospheric and Wake Turbulence Impacts on Wind Turbine Fatigue Loadings.” $50^{t h}$ AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition,  9–12 January  2012,  Nashville, TN  [online  proceedings]. URL: http://pdf.aiaa.org/getfile.cfm?urlX=527I%276D%26X%5BRO%2FSPWIP4S%5EQ%3AG%224%3A%5C%25%0A.  AIAA2012-0540. Reston, VA: American Institute of Aeronautics and Astronautics, January 2012; NREL/CP-500-53567. Golden, CO: National Renewable Energy Laboratory.