2.1 Node element

The concept of node originated from the lumped system analysis. The data structure of node element is defined as an abstract entity that represents flow cavities having a finite volume. The physical location of a node is at the barycentre of the represented cavity. A node is considered to represent the fluid properties in cavities in well-mixed condition [Reference Kulkarni and di Mare31, Reference Kulkarni and di Mare32], hence termed as the lumped-volume model. While this assumption is valid for small cavities, the gas states might vary significantly in large cavities. The node data structure is subdivided into several types. Figure 7 shows the taxonomy of node data structure in the present flow network modeller. It might be observed that this taxonomy structure allows instancing of specific types of nodes to represent specific types of flow cavities.

Figure 7. Taxonomy of node data structure.

As shown in Fig. 7, an axisymmetric node can be instanced to represent an axisymmetric cavity in flow network model. This type of node incorporates free vortex model to simulate the vortical flows in cavities. It can compute the losses generated by flow windage and due to the features protruding in cavity. The pre-swirl node is a sub-type of axisymmetric node and is specially developed to represent the turbine disc pre-swirl cavity [Reference Brillert, Reichert and Simon40]. The pre-swirl node has intrinsic methods to identify the links representing pre-swirl nozzles [Reference Benra, Dohmen and Schneider41] and turbine blade-root entry holes and to compute the pre-swirl flow losses.

A non-axisymmetric node is developed to represent the non-axisymmetric cavities in engine, wherein the tangential component of velocity is supressed. Junction node inherits the attributes of non-axisymmetric node data structure. It can compute both combining and dividing flow losses depending on the flow direction. Boundary nodes are generated at the boundaries of flow network model. All these nodes acquire their global number from the connectivity table.

The automated pre-processing methodology [Reference Kulkarni17, Reference Kulkarni and di Mare31, Reference Kulkarni and di Mare32] of node data structure calculates various internal parameters for axisymmetric node elements such as, the coordinates of cavity barycentre, number of solid walls, number of speed-frames and the description of protruding bolts etc. It examines the node cavity polygon to find the segments created from aerodynamic boundaries. These segments are interrogated to obtain the global number, headings, locations and pointers to the links connected to that node. The junction type node has specialised pre-processing methods to acquire their types [Reference Kulkarni and di Mare31, Reference Kulkarni and di Mare32] from the engine geometry model such as, a T-junction node, a Y-junction node, a valve node, a sudden or gradual expansion and contraction type node etc. By identifying the type of junction node, a representative flow loss model is explicitly allocated to that specific node instance. The boundary node data structure incorporates special database enquiry methods, which are designed to fetch and apply the flow boundary conditions to the network model.

2.2 Link element

Link is an abstract data structure that represents the flow paths connecting nodes of a flow network model. These flow-paths can include a variety of air system components and, therefore, the link data structure and its subtypes are primarily designed to contain the description of 1D flow loss models [Reference Kulkarni and di Mare31, Reference Kulkarni and di Mare32]. Flow paths are considered to have negligible volume and hence the time rate change of flow quantities across a link is very small.

Figure 8 depicts the taxonomy of link data structure. The link data structure is categorised into flow links and thermal links, although this paper describes only the flow links. The flow link data structure is further divided into two major subtypes, termed as the loss links and the mass links. Their data structures inherit all properties, data members and methods of the flow link data structure. The mass-type links calculate the mass flow rate passing through a flow link; whereas loss-type links estimate the total pressure drop across a link. For example, a mass-type orifice link computes the mass flow rate based on an implicit measure of loss such as the coefficient of discharge, Cd, whereas the loss-type links such as labyrinth seals, pipes, junctions etc. explicitly calculate the total pressure loss across their loss stations and thereby obtain the mass flow rate passing through those links. Together these link elements can represent essentially all kinds of SAS components in a gas turbine engine.

Figure 8. Taxonomy of the library of link data structure.

The link elements are instanced by link-set data structure during the automated flow network model generation process. Its automated pre-processing method supplies various input variables to the link data structures [Reference Kulkarni17, Reference Kulkarni and di Mare31, Reference Kulkarni and di Mare32]. The individual pre-processing methods of flow links then fetch the geometric and non-geometric metadata from the geometry modeller and initialise various internal flow variables used by the flow loss models. During the flow computations, various link-loss elements exchange the information with the flow network solver through a common interface mechanism. The design of this interface mechanism ensures the consistency and the ease of adding more loss elements into the loss model library of Virtual Gas Turbines design-simulation environment.

2.3 Link-set element

Link-set is a container-type data structure that is designed to represent both short and long flow paths [Reference Kulkarni and di Mare31, Reference Kulkarni and di Mare32]. A short flow path contains only one flow link representing an air system component such as an orifice or a seal etc., whereas a long flow-path such as a long pipe usually includes multiple SAS components such as, bends, valves, manifolds, sudden or gradual contractions and expanders and internal constrictions etc. The link-set data structure assembles all links in a flow-path that are connected in series, registers the nature of each link, maintains the sequence of those links, associates a loss model to each link and supplies the boundary conditions to each link element, thereby generating a unified flow-path description.

The link-set data structure holds various important variables. It gets the global number of connected nodes, the number of links contained in it and the number of link-set repetitions along circumferential direction. The detailed process of link-set assembly is explained in Refs. [Reference Kulkarni and di Mare31, Reference Kulkarni and di Mare32]. Such a loss-link assembly generates a system of tightly coupled simultaneous linear equations for each flow path. The link-set instances the appropriate type and number of link data structures during the pre-processing phase and maintains the pointers to flow links and adjacent nodes. Link-set then communicates with each contained loss link and starts loading the relevant data in them.

The link-set data structure establishes connectivity in the flow network model by virtue of its association properties. It provides a common interface to all types of link elements; thereby facilitating the information exchange between the flow network solver and the link elements. The pre- and post-processing methods of link-set data structure include data enquiry methods, and thus, the flow network modeller can fetch and return the geometry parameters and the boundary conditions parameters to other application spaces of Virtual Gas Turbines design-simulation environment.

The link-set data structure also has certain limitations. For example, it cannot include both mass-type links and loss-type links in the same link-set. The pre-processing method of link-set data structure ensures that these flow link types cannot be assembled as they employ different entropy generation mechanisms, which prevents the generation of the tightly coupled system of equations in flow network solver.

2.4 Solver connectivity

The present flow network solver has been designed to operate with any kind of node connectivity pattern. The connectivity in the flow network modeller is maintained by both node and link-set data structures. The node data structure establishes the connectivity with link-sets by recording the number of link-set connections and the identity of each link-set. Similarly, the link-set data structure records the global numbers of adjacent nodes. Both node and link-set data structures retain the pointers to the connected entities; thus, allowing the communication of their internal variables through shared methods.

Figure 9 shows the hierarchical twin-layered arrangement of the present flow network solver. The flow network solver iterates on the node elements to compute the state of the gas in the SAS cavities, and it communicates with the link-set solver to obtain the state of fluxes across the link-sets. The link-set solver collects loss variables from different types of loss link elements and solves gas dynamic equations for each flow path. It then returns the state of fluxes at each node–link-set interface to the flow network solver.

Figure 9. Hierarchical structure of network (node) and link-set (loss links) solvers.

2.5 Application of boundary conditions

The application of flow boundary conditions is one of the most tedious tasks in SAS flow network modelling. In the present flow network modeller, the boundary node data structure is furnished with an automated pre-processing method. The boundary nodes send enquiries to the engine performance model along with their headings and location coordinates in the SAS flow network model. The associative object-oriented framework of Virtual Gas Turbines design-simulation environment provides monitored channels for such information exchange between its application spaces or devices. After confirming the matching of both heading and the location coordinates, the engine performance model returns the performance data to the flow network model for each operation point of interest.

Similarly, the application of rotational speed boundary conditions to SAS flow network model is another labourious, error-prone and non-value-added activity. In the present flow network modeller, the cavity node data structure holds an automated pre-processing method that enquires, fetches and applies the rotational speed boundary conditions to the cavity walls. Each node sends a data request to the feature dynasty-based geometry model [Reference Kulkarni17, Reference Kulkarni, Lu, Wang and di Mare28, Reference Kulkarni, Lu, Wang and di Mare29] of the gas turbine engine. The geometry modeller acquires the rotational speed from the speed-dynasty or the ancestor spool feature and returns it to the node pre-processor. The node pre-processor then calculates the average tangential velocity of each cavity wall.

Figure 10 displays the flow network model of the secondary air system of a three-spool gas turbine engine. The automated process of generating such a whole-engine flow network model requires minimal computing resources compared to that required for generating a full-scale CFD model. If the engine geometry model already exists, it takes only a few minutes to transform the features, to fetch and apply the flow and speed boundary conditions and to initiate the SAS analysis at multiple operating points. This comprehensive methodology significantly reduces the overall SAS design-analysis time-cost and the dependency on human expertise.

Figure 10. Generation of flow network model consisting of nodes (red) and flow links (green) [Reference Kulkarni17].

3.1 The node solver

This section describes the mathematical formulation of the present 1D flow network solver, its flow variables, and its communication protocols. This flow network solver solves the conservation equations of mass flow rate ( $\dot w$ ), meridional momentum $\left( {{M_m}} \right)$ , tangential momentum $\left( {{M_\theta }} \right)$ and internal energy (e) to find the steady state solution for four primary flow variables namely, the meridional velocity (u _m ), tangential velocity (u _θ ), static temperature (T) and static pressure (P) at each node in the network. In its most basic form, the conservation equation for any node k can be written as shown in Equation (5),

(5) \begin{align}\frac{{\partial {\mathbb{q}^{t + 1}}}}{{\partial t}} = \; - \mathfrak{d}\dot \Re \end{align}

where, $\mathbb{q}$ is the vector of conservation quantities and $\mathfrak{d}\dot{\mathfrak{R}}$ is the vector of time rate change of derivatives of residuals of conservation quantities. The conservation quantities are related to the primary flow variables through the time rate change of residuals at interfaces of links, losses in cavities and volumetric forces within cavities. The fundamental nature of this mathematical scheme is to drive the residuals of conservation equations to zero [Reference Kulkarni17].

Apart from the primary flow variables, the solver also computes several secondary and auxiliary flow variables. A total of six secondary variables are evaluated, namely, the total temperature (T ⁰), total pressure (P ⁰), meridional sonic velocity ( $u_m^{\ast}$ ), tangential sonic velocity ( $u_\theta ^{\ast}$ ), static temperature at sonic condition ( ${T^{\ast}}$ ) and static pressure at sonic condition ( ${P^{\ast}}$ ). The auxiliary variables include fluid density ( $\rho$ ), entropy at static condition (s) and at total condition (s ⁰), enthalpy at static condition (h) and at total condition ( ${h^0}$ ) and the sonic velocity (a). These variables are arranged as shown in Equations (6)–(9),

(6) \begin{align}q{n^{prim}} = \;\left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{u_m}} & {{u_\theta }} & T & P\end{array}} \right\}\quad aux{n^{prim}} = \;\left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\rho & h & s & a\end{array}} \right\}\end{align}

(7) \begin{align}q{n^{sonic}} = \;\left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{u_m^*} & {u_\theta ^*} & {{T^*}} & {{P^*}}\end{array}} \right\}\quad aux{n^{sonic}} = \;\left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{\rho ^*}}& {{h^*}}& {{s^*}}& a\end{array}} \right\}\end{align}

(8) \begin{align}qn{0^{sec}} = \;\left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 0& 0& {{T^0}}& {{P^0}}\end{array}} \right\}\quad auxn{0^{sec}} = \left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{\rho ^0}}& {{h^0}}& {{s^0}}& {{a^0}}\end{array}} \right\}\end{align}

(9) \begin{align}qn = \;\left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}w & {{M_m}}& {{M_\theta }}& e\end{array}} \right\}\quad \mathfrak{d}\dot \Re = \;\left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{\dot w} & F& {{T_q}}& {\dot e}\end{array}} \right\}\end{align}

The variable tags ‘q’ and ‘aux’ indicate primary and auxiliary quantities, respectively, whereas ‘n’ indicates node. This variable structuring methodology is useful for linking the flow network solver with the Virtual Gas Turbines gas dynamic library [Reference di Mare42, Reference di Mare43]. This library consists of several hodograph transformations, and these are extensively used in the current formulation [Reference Kulkarni17, Reference Shapiro44, Reference Ames45]. Moreover, this solver incorporates several functions and methods for receiving the constructed datasets from plugins [Reference Kulkarni and di Mare31, Reference Kulkarni and di Mare32], for instancing and initialising the derived datasets, for visualising the network model [Reference Kulkarni and di Mare31, Reference Kulkarni and di Mare32], for evaluating the gas dynamic equations and for the post-processing of results etc.

At node k, $\mathfrak{d}\dot \Re $ can be expanded as shown in Equation (10). It is a vector of summation of all mass flow rates, the fluid forces at link-set openings, the torques applied by or on flow and the rate of change of flow energy at node k.

(10) \begin{align} \mathfrak{d}\dot{\mathfrak{R}}_{k}={\ \mathfrak{d}\left[ \begin{array}{c}\dot{w} \\F \\T_q \\\dot{e} \end{array}\right]}_k\ =\ {\mathfrak{d}\left[ \begin{array}{c}\rho Au_m \\\rho Au^2_m+\ PA \\\rho Au_{\theta }r \\\rho Au_mh \end{array}\right]}_k \end{align}

The time rate change of all residuals at k ^th node are obtained using existing instantaneous values of primary flow variables at adjacent nodes and links. $\mathfrak{d}\dot{\mathfrak{R}}_{k}$ could be expressed as shown in Equation (11), where the residuals of conservation quantities are partially differentiated with respect to time t.

(11) \begin{align}\frac{\partial }{{\partial t}}{\left[ {\begin{array}{c}{\rho V}\\{\rho V{u_m}}\\{\rho V{u_\theta }r}\\{\rho Ve}\end{array}} \right]_k} = \;\frac{\partial }{{\partial t}}\left( {\mathfrak{d}\Re _k^{cavity} + \;\mathop \sum \limits_{n = 1}^{n = nlinks} \mathfrak{d}\Re _k^{links} + \;\mathfrak{d}\Re _k^{loss}} \right)\end{align}

The term at RHS, $\sum \mathfrak{d}\Re _k^{links}$ , represents the summation of residuals of conservation quantities at the interfaces of node k and flow links that connect to it. The term $\mathfrak{d}\Re _k^{cavity}$ represents the residuals of volumetric forces applied by walls of cavity and $\mathfrak{d}\Re _k^{loss}$ represents the residuals of conservation quantities generated due to flow losses in the cavity node k. These vector notations may be further expanded as shown in Equation (12) [Reference LeVeque46],

(12) \begin{align} \frac{\partial }{\partial t}{\left[ \begin{array}{c}\rho V \\\rho Vu_m \\\rho Vu_{\theta }r \\\rho Ve \end{array}\right]}_k=\ \frac{\partial }{\partial t}\left\{\ \ \ \ \ {\left[ \begin{array}{c}\dot{w} \\F \\T_q \\\dot{e} \end{array}\right]}^{cavity}_k+\ \sum^{n=nlinks}_{n=1}{{\left[ \begin{array}{c}\dot{w} \\F \\T_q \\\dot{e} \end{array}\right]}^{links}_k}+\ {\left[ \begin{array}{c}\dot{w} \\F \\T_q \\\dot{e} \end{array}\right]}^{loss}_k\right\} \end{align}

Flow links have negligible volume compared to that of nodes and thus flow perturbations can transmit through them instantaneously. Consequently, the residuals computed at the interfaces of links are always instantaneous in nature. If each of these quantities is identified separately owing to its origin, Equation (12) might be stated in a more compact and generic form as shown in Equation (13),

(13) \begin{align}\frac{{\partial {\mathbb{q}^{t + 1}}}}{{\partial t}} = {\left( {A + B + L} \right)^t}\end{align}

where, A represents the vector of time rate change of conservation quantities at the interfaces of connecting links, B represents the vector of time rate change of volumetric quantities in the node and L represents the vector of time rate change of conservation quantities due to the losses in cavity. In finite difference form, the conservation equation at a node can be expressed as shown in Equation (14),

(14) \begin{align}\frac{{{\mathbb{q}^{t + 1}} - {\mathbb{q}^t}}}{{\Delta \tau }} = {\left( {A + B + L} \right)^{t + 1}}\end{align}

For obtaining the steady state solution, conservation quantities should not vary with time and hence RHS must diminish to zero. This means that the sum of all residuals ${\left( {A + B + L} \right)^{t + 1}}\; \cong 0$ . The aforementioned formulation thus represents a well-configured system of simultaneous equations. As the volumes of all cavities remain constant, Equation (14) could be written as,

(15) \begin{align}\frac{{V\!\left( {{\mathbb{q}^{t + 1}} - \;{\mathbb{q}^t}} \right)}}{{\Delta \tau }} = -\mathfrak{d}\Re \end{align}

Equation (15) can be conveniently represented as shown in Equation (16), where ${\mathbb{q}^{t + 1}}$ represents the solution of Equation (15) based on the geometry parameters of node k and the residuals of conservation quantities.

(16) \begin{align}{\mathbb{q}^{t + 1}} = \; - \frac{{\sigma \Delta \tau }}{V}{\textrm{ }}\mathfrak{d}\Re \end{align}

The term $\frac{{\sigma \;\Delta \tau }}{V}$ is used for scaling the residuals of conservation quantities at node k having volume V. This factor also converts time rate change of conserved quantities into actual conservation quantities such as fluid mass, meridional momentum, tangential momentum and internal energy of fluid. $\sigma $ is Courant number for flow network model, which scales the speed of information propagation within the cavity of node k. $\Delta$ τ is the time required to form homogeneous flow mixture in cavity. This time period is important to fulfil the assumption that conservation quantities computed at the barycentre of cavity represent the fully mixed flow conditions. This residence time is not known a priory, and it needs to be computed considering the volume of each cavity and the total time steps required to achieve steady state in complete network. The residence time would be shorter for the nodes having smaller volumes, whereas it is longer for the nodes having larger volumes. The solver always uses the shortest possible time step for time marching to ensure that the link perturbations are not transmitted through a node before the end of time step.

The nonlinear unsteady flow equation, shown in Equation (16), states the conservation of flow quantities in node k. In order to obtain the change in primary flow variables, the Equation (16) is converted into locally linear simultaneous algebraic equation as shown in Equation (17). The transformation Jacobian [Reference Laney47] has been defined as,

(17) \begin{align}\left[ {\partial {\mathbb{q}^{t + 1}}} \right]\left[ {\mathfrak{d}\mathfrak{q}} \right] = -\mathfrak{d}{\mathbb{q}^t} = - \frac{{\sigma \;\Delta \tau }}{V}\mathfrak{d}\Re \end{align}

The terms at LHS are Jacobian matrix of partial derivatives of conservation quantities $\left[ {\partial {\mathbb{q}^{t + 1}}} \right]$ with respect to primary flow variables and the vector of incremental change in primary flow variables $\left[ {\mathfrak{d}\mathfrak{q}} \right]$ , respectively. RHS represents residuals of conservation quantities, expressed by symbol $ -\mathfrak{d}{\mathbb{q}^t}$ and computed by term $ - \dfrac{{\sigma \;q\tau }}{V}\mathfrak{d}\Re $ . This equation represents the system of simultaneous equations as shown in Equation (18),

(18) \begin{align}\mathfrak{A}{\textrm{ }}\mathfrak{d}\mathfrak{q} = - \frac{{\sigma \;\Delta \tau }}{V}\mathfrak{d}\Re \end{align}

where, $\mathfrak{A}$ is the Jacobian matrix comprising local slopes (partial derivatives) of conservation quantities with respect to the primary flow variables. $\mathfrak{d}\mathfrak{q}$ is the vector comprising differential changes in the primary flow variables. The flow mixing time scale, $\Delta \tau $ , in Equation (18) is calculated based on the maximum speed of information propagation through nodes. The speed of information propagation is the highest possible magnitude of velocity in the direction of flow in links. For simplicity, the highest eigenvalue of the system of simultaneous equations for connecting link-sets is chosen, which is also known as spectral radius, ${\mathfrak{S}^R}$ . In the present flow network solver, the spectral radius is calculated as $A\;\left| {{u_m} + \;a} \right|$ , where A is area of opening of a link-set into cavity, ${u_m}$ is the highest flow velocity in connecting link-set and $a$ is speed of sound at existing static temperature in that link-set. This system of equations is represented by the conventional notation shown in Equation (19),

(19) \begin{align}\mathfrak{A}{\textrm{ }}\mathfrak{d}\mathfrak{q} = -\mathfrak{d}\Re \end{align}

And thus,

(20) \begin{align}\mathfrak{d}\mathfrak{q} = -\!\mathfrak{A}{{\textrm{ }}^{ - 1}}\mathfrak{d}\Re \end{align}

For obtaining the solution of system of simultaneous equations, it is necessary to invert Jacobian matrix $\mathfrak{A}$ assuming that it is not singular and its inverse does exist. Equation (21) can be obtained after substituting the terms in Jacobian matrix $\mathfrak{A}$ and the vector $\mathfrak{d}\mathfrak{q}$ ,

(21) \begin{align}\left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{\dfrac{{\partial w}}{{\partial {u_m}}}} & {\dfrac{{\partial w}}{{\partial {u_\theta }}}} & {\dfrac{{\partial w}}{{\partial T}}} & {\dfrac{{\partial w}}{{\partial P}}}\\[9pt] {\dfrac{{\partial {M_m}}}{{\partial {u_m}}}} & {\dfrac{{\partial {M_m}}}{{\partial {u_\theta }}}} & {\dfrac{{\partial {M_m}}}{{\partial T}}} & {\dfrac{{\partial {M_m}}}{{\partial P}}}\\[9pt] {\dfrac{{\partial {M_\theta }}}{{\partial {u_m}}}} & {\dfrac{{\partial {M_\theta }}}{{\partial {u_\theta }}}} & {\dfrac{{\partial {M_\theta }}}{{\partial T}}} & {\dfrac{{\partial {M_\theta }}}{{\partial P}}}\\[9pt] {\dfrac{{\partial e}}{{\partial {u_m}}}} & {\dfrac{{\partial e}}{{\partial {u_\theta }}}} & {\dfrac{{\partial e}}{{\partial T}}}& {\dfrac{{\partial e}}{{\partial P}}}\end{array}} \right]\;\;\;\;{\left[ {\begin{array}{c}{\begin{array}{c}{\mathfrak{d}{u_m}}\\[5pt] {\mathfrak{d}{u_\theta }}\\[5pt] {\mathfrak{d}T}\\[5pt] \end{array}}\\[5pt] {\mathfrak{d}P}\end{array}} \right]_k} = \; - \left[ {\begin{array}{c}{\mathfrak{d}w\;}\\[5pt] {\mathfrak{d}{M_m}}\\[5pt] {\mathfrak{d}{M_\theta }}\\[5pt] {\mathfrak{d}e}\end{array}} \right]_k^\Re \end{align}

The $\mathfrak{A}$ matrix can be represented as shown in Equation (22) by substituting the flux variables for conservation variables [Reference LeVeque46]. It is possible to write the exact analytical forms of partial derivatives of flux variables, which simplify the computations of partial derivatives and hence the actual matrix inversion is not necessary. This reduces the number of computations substantially.

(22) \begin{align}\frac{V}{\tau }\;\left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{\dfrac{{\partial \rho }}{{\partial {u_m}}}}& {\dfrac{{\partial \rho }}{{\partial {u_\theta }}}}& {\dfrac{{\partial \rho }}{{\partial T}}}& {\dfrac{{\partial \rho }}{{\partial P}}}\\[9pt] {\dfrac{{\partial\!\left( {\rho {u_m}} \right)}}{{\partial {u_m}}}} & {\dfrac{{\partial\!\left( {\rho {u_m}} \right)}}{{\partial {u_\theta }}}} & {\dfrac{{\partial\!\left( {\rho {u_m}} \right)}}{{\partial T}}} & {\dfrac{{\partial\!\left( {\rho {u_m}} \right)}}{{\partial P}}}\\[9pt] {\dfrac{{\partial\!\left( {\rho {u_\theta }r} \right)}}{{\partial {u_m}}}} & {\dfrac{{\partial\!\left( {\rho {u_\theta }r} \right)}}{{\partial {u_\theta }}}} & {\dfrac{{\partial\!\left( {\rho {u_\theta }r} \right)}}{{\partial T}}} & {\dfrac{{\partial\!\left( {\rho {u_\theta }r} \right)}}{{\partial P}}}\\[9pt] {\dfrac{{\partial\!\left( {\rho e} \right)}}{{\partial {u_m}}}} & {\dfrac{{\partial\!\left( {\rho e} \right)}}{{\partial {u_\theta }}}} & {\dfrac{{\partial\!\left( {\rho e} \right)}}{{\partial T}}} & {\dfrac{{\partial\!\left( {\rho e} \right)}}{{\partial P}}}\end{array}} \right]{\left[ {\begin{array}{c}{\begin{array}{c}{\mathfrak{d}{u_m}}\\[5pt] {\mathfrak{d}{u_\theta }}\\[5pt] {\mathfrak{d}T} \end{array}}\\[5pt] {\mathfrak{d}P}\end{array}} \right]_k} = - \dfrac{V}{\tau }\left[ {\begin{array}{c}{\mathfrak{d}\!\left( \rho \right)\;}\\[5pt] {\mathfrak{d}\!\left( {\rho {u_m}} \right)}\\[5pt] {\mathfrak{d}\!\left( {\rho {u_\theta }r} \right)}\\[5pt] {\mathfrak{d}\!\left( {\rho e} \right)}\end{array}} \right]_k^\Re \end{align}

After inverting Jacobian matrix and multiplying it by RHS, a vector containing differential change in primary flow variables is obtained. This vector is used to update the values of primary and auxiliary flow variables at node k as shown in Equations (23)–(26),

(23) \begin{align}u_m^{t + 1} = u_m^t + \mathfrak{d}{u_m}\qquad {\rho ^{t + 1}} = {\rho ^t} + \mathfrak{d}\rho \end{align}

(24) \begin{align}u_\theta ^{t + 1} = u_\theta ^t + \mathfrak{d}{u_\theta }\qquad {h^{t + 1}} = {h^t} + \mathfrak{d}h\end{align}

(25) \begin{align}{T^{t + 1}} = {T^t} + \mathfrak{d}T\qquad {s^{t + 1}} = {s^t} + \mathfrak{d}s\end{align}

(26) \begin{align}{P^{t + 1}} = {P^t} + \mathfrak{d}P\qquad {a^{t + 1}} = {a^t} + \mathfrak{d}a\end{align}

These new values of flow variables at each node are used as instantaneous guess values in the next iteration of solver. The network solver continues iterations to compute the conservation quantities at the nodes until a convergence criterion is satisfied.

3.2 Solver interfaces

The RHS terms in Equation (12) represent the summation of residuals of conservation quantities within the node k and at the interfaces of node k and flow links that connect to it. The flow network solver obtains the residual vectors $\mathfrak{d}\Re _k^{cavity}$ and $\mathfrak{d}\Re _k^{loss}$ from the representative loss model of cavity node k. Depending on the type of the node such as, axisymmetric node, non-axisymmetric node, junction node etc., the specific loss model for node k generates these residual vectors and returns those to the flow network solver.

Further, the flow network solver obtains the RHS term $\sum \mathfrak{d}\Re _k^{links}$ from each link-set connected to node k. It supplies the global node numbers and flow boundary variables to a link-set j. The link-set j acquires the partial derivatives of flux and loss variables from the loss models of each contained flow links. It also provides the spectral radius, ${\mathfrak{S}^R}$ , of its Jacobian matrix to calculate the flow residence time in node k. These variables are substituted in Equation (12). More explanation is provided in Section 3.3.

3.3 The link-set solver

The fundamental mathematical scheme of link-set solver is similar to that of flow network solver. It uses linearised gas dynamic equations, which are coupled to the linearised loss correlations in flow loss models. These loss models represent the flow physics in air system components.

Figure 11. The conservation variables computed by link-set solver.

Figure 11 shows the schematic representation of link-set j connecting nodes k and m. The flow variables at the interfaces of connected nodes, k and m, are used as the boundary conditions. The link-set j can have multiple internal stations that are physically spaced in meridional plane. The number and the locations of internal stations are determined during the pre-processing stage of the link-set [Reference Kulkarni and di Mare31, Reference Kulkarni and di Mare32]. The link-set solver computes four primary variables namely, the meridional velocity $\left( {{u_m}} \right)$ , tangential velocity $\left( {{u_\theta }} \right)$ , static temperature (T) and static pressure (P). It also consists of two secondary variables and four sonic state variables, which are total temperature $\left( {{T^0}} \right)$ , total pressure $\left( {{P^0}} \right)$ , meridional sonic velocity $\left( {u_m^*} \right)$ , tangential sonic velocity $\left( {u_\theta ^*} \right)$ , enthalpy at sonic condition $\left( {{h^*}} \right)$ and entropy at sonic condition $\left( {{s^*}} \right)$ , respectively. Apart from these, auxiliary variables are also included in link-set solver namely, fluid density $\left( \rho \right)$ , enthalpy (h), entropy (s) and speed of sound (a). The auxiliary variables are also computed at total and sonic state conditions. All these flow variables are coupled through gas dynamic equations and geometry variables of link-set such as cross-section area (A). The link-set solver computes the mass flow rate $\left( {\dot w} \right)$ through all links in link-set j based on total pressure loss $\left( {{P^{0,\;\;loss}}} \right)$ between consecutive loss stations. These variables are arranged as shown in Equations (27)–(29),

(27) \begin{align}q{l^{\,primary}} = \;\left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{u_m}} & {{u_\theta }}& T & P\end{array}} \right\}\qquad aux{l^{primary}} = \;\left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\rho & h & s & a\end{array}} \right\}\end{align}

(28) \begin{align}q{l^0} = \;\left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 & 0& {{T^0}}& {{P^0}}\end{array}} \right\}\qquad aux{l^0} = \left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{\rho ^0}}& {{h^0}}& {{s^0}}& {{a^0}}\end{array}} \right\}\end{align}

(29) \begin{align}q{l^{sonic}} = \;\left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{u_m^*}& {u_\theta ^*}& {{T^*}}& {{P^*}}\end{array}} \right\}\qquad aux{l^{\,sonic}} = \left\{ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{\rho ^*}}& {{h^*}}& {{s^*}}& a\end{array}} \right\}\end{align}

The names ‘q’ and ‘aux’ indicate arrays of primary and auxiliary variables, respectively, whereas ‘l’ indicates link or link-set. Such a variable arrangement has been useful to establish a common interface between the link-set solver and the loss model library of Virtual Gas Turbines design-simulation environment [Reference Kulkarni and di Mare48, Reference Kulkarni and di Mare49].

Equation (30) shows the expanded view of the system of simultaneous linear algebraic equations for station k of the link-set j. The LHS of Equation (30) includes a Jacobian matrix containing the partial derivatives of conservation quantities with respect to primary, total and sonic flow variables and a vector representing incremental changes $\left\{ {\Delta {\textrm{X}}} \right\}$ in primary flow variables at station k. It might be observed that the link-set Jacobian matrix is relatively sparse, and it is not symmetric about its diagonal. Some of the loss models generate partial derivatives of conservation quantities with respect to geometry variables. These partial derivatives occupy several diverse locations in the link-set Jacobian matrix. During the matrix assembly, each conservation equation is partially differentiated with respect to a flow variable. This process is repeated for all four primary, two total and four sonic state flow variables.

(30)

The RHS of Equation (30) consists of a vector representing the residuals of primary, total and sonic state flow variables. Equation (31) shows the generalised form of residual vector at station k.

(31) \begin{align}{\left[ {\begin{array}{c}{\mathfrak{d}\dot w}\\[5pt] {\mathfrak{d}{u_\theta }}\\[5pt] {\mathfrak{d}h}\\[5pt] {\mathfrak{d}s}\\[5pt] {\mathfrak{d}{T^0}}\\[5pt] {\mathfrak{d}{P^0}}\\[5pt] {\mathfrak{d}u_m^*}\\[5pt] {\mathfrak{d}u_\theta ^*}\\[5pt] {\mathfrak{d}{h^*}}\\[5pt] {\mathfrak{d}{s^*}}\end{array}} \right]_k} = {\left[ {\begin{array}{c}{{{\left( {\rho \;{u_m}A} \right)}_k}\; - \;{{\left( {\rho \;{u_m}A} \right)}_m} + \;\Delta \dot w}\\[5pt] {{{\left( {{u_\theta }} \right)}_k} - \;{{\left( {{u_\theta }} \right)}_m} + \;\Delta {u_\theta }}\\[5pt] {{{\left( h \right)}_k} - \;{{\left( h \right)}_m} + \;\Delta h}\\[5pt] {{{\left( s \right)}_k} - \;{{\left( s \right)}_m} + \;\Delta s}\\[5pt] {{{\left( {{T^0}} \right)}_k} - \;{{\left( {{T^0}} \right)}_m} + \;\Delta {T^0}}\\[5pt] {{{\left( {{P^0}} \right)}_k} - \;{{\left( {{P^0}} \right)}_m} + \;\Delta {P^0}}\\[5pt] {{{(u_m^*)}_k} - {{(u_m^*)}_m} + \Delta u_m^*}\\[5pt] {{{(u_\theta ^*)}_k} - {{(u_\theta ^*)}_m} + \Delta u_\theta ^*}\\[5pt] {{{\left( {{h^*}} \right)}_k} - \;{{\left( {{h^*}} \right)}_m} + \;\Delta {h^*}}\\[5pt] {{{\left( {{s^*}} \right)}_k} - \;{{\left( {{s^*}} \right)}_m} + \;\Delta {s^*}}\end{array}} \right]_k}\end{align}

Figure 12 shows the generalised assembly of the Jacobian matrix and the vector of residuals for link-set j having two stations, k and m. The loss model of the SAS component in link-set j provides the partial derivatives of its linearised loss correlations to the Jacobian matrix and the incremental variables to the vector of residuals. Both the Jacobian matrix and the vector of residuals are automatically assembled in block column-wise manner. The partial derivatives of loss correlations occupy various lower diagonal positions in the Jacobian matrix; thereby coupling the flow conditions at stations, k and m.

Figure 12. Assembly of the Jacobian matrix and the residual vector for link-set j.

Further the static pressure boundary conditions are applied by adding one equation and one unknown to the system of equations. The extra equation represents a mass continuity constraint, whereas the extra variable represents the exit pressure. The additional equations and unknowns appear in the Jacobian matrix as one extra row and column each, which are highlighted in Fig. 12.

To obtain the steady state solution for link-set j, the vector of residuals at the RHS (see Fig. 12) should diminish to zero. The link-set solver inverts the Jacobian matrix [Reference Golub and van Loan50] and multiplies it to the RHS using the standard routines of LAPACK library [51]. It solves the system of tightly coupled linearised simultaneous equations by balancing the mass flow rate, tangential velocity, total pressure and total temperature across the link-set j and then computes the residuals and the spectral radius of Jacobian matrix.

3.4 Interface with loss models

The link-set solver is coupled to the link data structures and they exchange data through a common interface irrespective of the type of the link. The solver provides the banded vectors of primary, total and auxiliary flow variables (see Equations (27)–(29)) to the link loss models. The link-set flow variables are associated through gas dynamic relations and hodograph transformations [Reference Kulkarni17, Reference di Mare42, Reference Ames45] and hence cannot be modified by the loss models. The loss models generate loss estimates using these flow variables and return the loss variables and their partial derivatives to the link-set solver through the loss vectors, $q{L^{vect}}$ and $dq{L^{vect}}$ . Figure 13 illustrates the variable exchange between link-set solver and loss models.

Figure 13. Interface between link-set solver and the SAS link loss models.

3.5 The SAS loss models

The present loss model library [Reference Kulkarni17] incorporates two forms of SAS loss models namely,

1. a) the loss-type links, which compute the mass flow rate based on inlet total pressure, exit static pressure and drop in total pressure (or rise in entropy)

2. b) the mass-type links that compute the rise in entropy (or drop in total pressure) based on inlet total pressure, exit static pressure and the mass flow rate through the link

The linearisation methods for both loss-type and mass-type loss models evaluate the partial derivatives of loss correlations with respect to flow and geometry parameters [Reference Kulkarni17, Reference Kulkarni and di Mare48, Reference Kulkarni and di Mare49] and generate four incremental loss variables, namely mass flow rate difference ( $\Delta \dot w$ ), total pressure drop ( $\Delta {P^0}$ ), total temperature rise ( $\Delta {T^0}$ ) and the difference in tangential velocity ( $\Delta {u_\theta }$ ). The vector of these loss quantities and the vector containing their partial derivatives are shown in Equation (32).

(32) \begin{align}q{L^{vect}} = \;\left[ {\begin{array}{c}{\Delta \dot w}\\{\Delta {P^0}}\\{\Delta {T^0}}\\{\Delta {u_\theta }}\end{array}} \right]\quad dq{L^{vect}} = \;\left[ {\begin{array}{c}{\mathfrak{d}\dot w}\\{\mathfrak{d}{P^0}}\\{\mathfrak{d}{T^0}}\\{\mathfrak{d}{u_\theta }}\end{array}} \right].\end{align}

Depending on the type of the loss model, either the total pressure drop or the mass flow rate can be predicted. Therefore, for a loss model that predicts total pressure drop, the mass flow rate through that link is considered as a free variable and vice-a-versa. This difference in predicting the loss variables based on loss-type approach and mass-type approach necessitates the loss vectors to be modified as shown in Equations (33) and (34).

For loss-type link-set:

(33) \begin{align}q{L^{vect}} = \;\left[ {\begin{array}{c}{freevar}\\{\Delta {P^0}}\\{\Delta {T^0}}\\{\Delta {u_\theta }}\end{array}} \right]\quad dq{L^{vect}} = \;\left[ {\begin{array}{c}{\mathfrak{d}freevar}\\{\mathfrak{d}{P^0}}\\{\mathfrak{d}{T^0}}\\{\mathfrak{d}{u_\theta }}\end{array}} \right].\end{align}

For mass-type link-set:

(34) \begin{align}q{L^{vect}} = \;\left[ {\begin{array}{c}{\Delta \dot w}\\{freevar}\\{\Delta {T^0}}\\{\Delta {u_\theta }}\end{array}} \right]\quad dq{L^{vect}} = \;\left[ {\begin{array}{c}\mathfrak{d}{\dot w}\\{\mathfrak{d}freevar}\\{\mathfrak{d}{T^0}}\\{\mathfrak{d}{u_\theta }}\end{array}} \right].\end{align}

The loss variables are returned to the RHS vector of Equation (31) through $q{L^{vect}}$ , whereas the partial derivatives of loss correlations are returned to the link-set Jacobian matrix through $dq{L^{vect}}$ .

3.6 Handling of flow reversal and flow choking

Flow reversal and flow choking are commonly observed in SAS flow network models. Flow reversal might occur either due to poor SAS design or due to numerical instability, whereas flow choking is intentionally used to meter and control the mass flow rates in specific regions of SAS. Various flow paths in the SAS are designed to choke during certain operating conditions. Therefore, the present network solver incorporates both the reversed flow and choked flow handling capabilities.

Due to the tightly coupled nature of equations in the present solver, it can handle simultaneous flow reversal conditions in multiple link-sets. While iterating, the link-set solver monitors pressure, temperature and meridional velocity at both source and sink nodes to estimate the possibility of flow reversal. If flow seems to be reversing, the link-set solver swaps the pointers to the source and sink nodes; thereby updating the link-set boundary conditions.

Flow choking can be detected, if meridional flow velocity at any station in a link-set equals the local sonic velocity. In the present link-set solver, local sonic conditions are computed at each internal station of link-set and thus flow choking can be easily detected. At present, the loss models in the library cannot predict flow losses at sonic or supersonic flow conditions; thus keeping the mass flow rate constant through a choked link-set.

Figure 14. The overall workflow chart for SAS flow network analyses.