Modelling, simulation and optimisation of parabolic trough power plants

Abstract

We present a mathematical model built to describe the fluid dynamics for the heat transfer fluid in a parabolic trough power plant. Such a power plant consists of a network of tubes for the heat transport fluid. In view of optimisation tasks in the planning and in the operational phase, it is crucial to find a compromise between a very detailed description of many possible physical phenomena and a necessary simplicity needed for a fast and robust computational approach. We present the model, a numerical approach, simulation for single tubes and also for realistic network settings. In addition, we optimise the power output with respect to the operational parameters.

1 Introduction

It is well known that the worldwide energy consumption is increasing rapidly both in industrial and in developing countries. The classical way to produce energy is co-responsible for the climate change and for environmental pollution. Therefore, there is a strong need to reduce the use of traditional non-renewable resources and to strongly increase the use of renewable resources. This refers to all types of energy, the electrical current as well as to energy in form of heat.

A important renewable source is given by the solar radiation. There are different possibilities to use the solar radiation. The radiation is typically used to heat a medium. In some cases, the heat is used directly, e.g. for heating or hot water production. However, in most of the cases on a industrial scale of bigger power plants, the heat is used to produce electrical energy, which can be transported easily over long distances.

One of the most promising types of solar power plants are the so called parabolic trough power plants. In parabolic trough’s, a special fluid in tubes positioned in the focus of big parabolic mirrors is heated by the solar radiation (see Figure 2–3). This fluid circulates and heats up a big heat storage. From there on, demand heat is extracted and used to power a steam turbine to finally produce electric current in an electric generator. This technique is realised in many places all over the world (see [26] as a database on existing and planned powerplants). Today we have already a few thousands of GW power available by parabolic trough power stations worldwide. A typical modern parabolic trough power plant has a pike power in the order of magnitude $10^2 -10^3\,$ MW.

As an example in this paper, we will refer often to the power plant NOOR I in Morocco [6, 9], see Figure 1. The Moroccan solar plan implemented by Moroccan Agency For Solar Energy (MASEN) consists of creating various solar power plants to be located in different regions [23]. NOOR I is one of the first power complex constructed by MASEN in Ouarzazate-Morocco in 2016. The plant generates solar power using parabolic trough collectors that uses thermal oil as a heat transfer fluid (HTF) [2].

For a deeper and profound understanding of the possibilities and limitations of parabolic trough systems, accurate models have to be available. We describe the complete power plant by including the main components, e.g. the thermo-fluid dynamics in the tube network and the pumps as driving forces for the fluid. This is done by a set of PDE’s describing the thermo-fluid dynamics in a network of tubes. However, to make such an approach realisable a series of simplifications are necessary. First we reduce the equations to a one (spacial and longitudinal) dimensional system on the tube network by appropriately describing the relevant effects from the (omitted) cross-sectional dimensions. Then we apply an asymptotic analysis with respect to certain small parameters to further reduce the system. Finally we obtain a set of 3 nonlinear coupled first-order PDE’s on every tube complemented with a set of physical meaningful conditions in the nodes of the network. This model allows for fast and robust simulations and thus for optimisation approaches. On one side, we optimise with respect to the various input data such as applied pressures or inflow temperature. On the other side, we can optimise with respect to various system parameters, e.g. network structure, dimensions and properties of the tubes and the thermo fluid in use. Therefore, this model can be applied both in the planning phase and in the operational phase of a parabolic trough power station. Also it is a good starting point for modelling more complex hybrid power plants which will become more and more important.

Figure 1. Parabolic trough solar power plant (Noor I) in Ouarzazate, Morocco (from Google Earth (left) and from https://commons.wikimedia.org/wiki/File:Mirrors.JPG (right)).

Figure 2. Schematic of a parabolic trough plant [5, dredit to NREL/SAM].

Figure 3. Schematic diagram of the different aspects of the operation of a parabolic solar collector.

There is existing literature on various aspects of parabolic trough power plants. There is overview articles on collectors [12], there is literature comparing the parabolic trough and the heliostat technology [31] using space independent models for a single pipe. There is literature studying and comparing the parabolic trough plants with thermo fluids to direct steam generating power plants [35, 8]. Another direction is literature considering in particular alternative thermo fluids in [21, 1]. There is models which focus on the detailed cross-sectional physics and dynamics in a tube, i.e. in [13, 28, 10]. There is literature doing a space independent thermal analysis for a single pipe in [25] with related optimisation tasks in [24]. There is other Ansatzes based on simple algebraic models or space independent models aiming to optimise economic [3], structural [18] or network structural [22] aspects. Control aspects can be found in [34]. Network simulations based on algebraic models can also be found in [32]. Optimisation based on complex 3D models can be found in [38]. Energy balancing models for the power plant Noor I in Marocco can be found in [1]. Then, there is literature doing thermo fluid dynamic simulations with optimisation tasks on tailored software tools like ColSimCSP in [29]. Unfortunately, we could not find out which model in detail is coded in the software package ColSimCSP. However, not considering fluid dynamics aspects (apart from the temperature) is a very rough approximation since in a HTF fluid in the relevant temperature range from 300 to 700 K the density is varying up to 50% (see Figure 5).

To summarise, our approach represents a complete thermo-fluid model on a network with spatial and time dependency. The model is a result of various asymptotic processes simplifying the model significantly. It includes appropriate node conditions for the network. The model is robust and fast in simulations. Therefore, it allows us to optimise with reasonable effort, i.e. the net power output.

To our knowledge, the model presented here is unique and one of the most complete models available to describe a parabolic trough power station. In particular, it is based on first principal thermo–fluid equations and does not only consider energy balances. It promises a wide range of applicability.

In Section 2, we set up and analyse a one-dimensional mathematical model on the network, including the scaling and the asymptotic analysis. In Section 3, we deal with the numerical simulation of the model and present some examples. Finally, in Section 4 we select real world power plants models and present optimisation results for both the planning and the operational phase.

2 Mathematical model

Here, we introduce a system of balance laws for the relevant quantities under consideration that we will use to model, to simulate and to optimise a parabolic trough power plant described in the previous section. Since the power plant consists of a network of tubes conducting the heated flow, we start with the description of a single pipe in such a network.

Models of the type we are going to use here are widely used in similar applications where the heat transport is crucial. We refer to applications in gas pipelines [7], exhaust pipes [16, 17], tunnel ventilation [15], solar updraft towers [14, 37], energy towers [4] etc. In particular, the application for tubes in solar parabolic trough’s were studied in [27, 30].

In the following, a mathematical model for the heat-transfer fluid flow in a collector is introduced. Different fluid quantities including fluid density, velocity, temperature and pressure are described as functions of time and space.

2.1 The unidirectional flow system

Generally, the cross-sectional area of the collector flow tube is very small compared to its length. Based on that, a one-dimensional flow model is derived by averaging variables over the cross-sectional variables $\tilde{y}$ and $\tilde{z}$ directions. Namely, for a given quantity $\tilde{F}$ , we have

(2.1) \begin{equation}\tilde{A}(\tilde{x}) \tilde{f}(\tilde{x},\tilde{t}) = \iint \tilde{F}(\tilde{x},\tilde{y},\tilde{z},\tilde{t}) d\tilde{y} d \tilde{z},\end{equation}

where the $\;\tilde{}\;$ indicates unscaled variables (with physical dimensions). $\tilde{f}$ denotes the cross-sectional averaged quantity. Here, $\tilde{x}$ , $\tilde{t}$ , $\tilde{A}(\tilde{x})$ are the longitudinal space variable (along the tube), the time and the cross-sectional area of the tube, respectively.

The governing equations are derived from the conservation and balance laws for the mass, momentum and energy, where the unknown variables are the density $\tilde{\rho} = \tilde{\rho}(\tilde{x},\tilde{t})$ , the velocity $\tilde{u} = \tilde{u}(\tilde{x},\tilde{t})$ , the temperature $\tilde{T} = \tilde{T}(\tilde{x},\tilde{t})$ and the pressure $\tilde{p} = \tilde{p}(\tilde{x},\tilde{t})$ , all intended as cross-sectional mean values. Thus, the governing equations are given by

1. (i) Conservation of Mass

   (2.2) \begin{equation}(\tilde{A}\tilde{\rho})_{\tilde{t}} + (\tilde{A}\tilde{\rho}\tilde{u})_{\tilde{x}} = 0\end{equation}

2. (ii) Momentum balance

   (2.3) \begin{equation}\left(\tilde{A}\tilde{\rho}\tilde{u}\right)_{\tilde{t}} + \left(\tilde{A}(\tilde{\rho}\tilde{u}^2+\tilde{p})\right)_{\tilde{x}} = -\tilde{U}_i\tilde{\tau}_w + \tilde{A}\left(\tilde{\tau}_{11}\right)_{\tilde{x}}\end{equation}

3. (iii) Energy balance

   (2.4) \begin{align}& \left(\tilde{A}(\tilde{c}_v\tilde{\rho}\tilde{T}+\frac{1}{2}\tilde{\rho}\tilde{u}^2)\right)_{\tilde{t}} + \left\{\tilde{A}\tilde{u}\left(\tilde{c}_v\tilde{\rho}\tilde{T}+\frac{1}{2}\tilde{\rho}\tilde{u}^2+\tilde{p}\right)\right\}_{\tilde{x}} \nonumber \\[5pt]&\quad =\tilde{U}_i\tilde{u}\tilde{\tau}_w + \tilde{A}\left(\tilde{u}\tilde{\tau}_{11}\right)_{\tilde{x}} + \tilde{A}\left(\tilde{k}\tilde{T}_{\tilde{x}}\right)_{\tilde{x}}+ \tilde{\dot{q}}_{s} - \tilde{U}_{o}\tilde{\dot{q}}_{rad} - \tilde{U}_{o}\tilde{\dot{q}}_{conv}, \end{align}

where $\tilde{\tau}_{11} = \tilde{\tau}_{11}(\tilde{x},\tilde{t})$ is the internal shear stress, $\tilde{\tau}_w = \tilde{\tau}_w(\tilde{x},\tilde{t})$ the wall friction, $\tilde{\dot{q}}_{s} = \tilde{\dot{q}}_{s}(t)$ the beam solar radiation, $\tilde{\dot{q}}_{conv} = \tilde{\dot{q}}_{conv}(\tilde{x},\tilde{t})$ the convective heat loss and $\tilde{\dot{q}}_{rad} = \tilde{\dot{q}}_{rad}(\tilde{x},\tilde{t})$ the radiation heat exchange between the absorber and the atmosphere, $\tilde{U}_i = \tilde{U}_i(\tilde{x}) = \pi \tilde{D}_i(\tilde{x})$ is the inner circumference (diameter) of the tube and $\tilde{U}_o = \tilde{U}_a(\tilde{x}) = \pi \tilde{D}_a(\tilde{x})$ is the absorber circumference (diameter). The physical parameters $\tilde{c}_v$ and $\tilde{k}$ denote the specific heat capacity and the heat conductivity of the fluid, respectively.

The shear stress of the turbulent fluid on the wall is described classically as follows

(2.5) \begin{equation}\tilde{\tau}_w = \frac{\xi}{8}\tilde{\rho} \tilde{u}|\tilde{u}|,\end{equation}

where $\xi$ denotes the coefficient of friction calculated according to the Colebrook equation [20, 1].

For the unidirectional flow of a Newtonian fluid the shear stress $\tilde{\tau}_{11}$ is given as follows [27]:

(2.6) \begin{equation}\tilde{\tau}_{11} = \frac{4}{3} \tilde{\mu} \tilde{u}_{\tilde{x}},\end{equation}

where $\tilde{\mu}$ is the viscosity of the fluid.

The solar radiation received is considered as a heat flux/rate and can be expressed as [1]:

(2.7) \begin{equation}\tilde{\dot{q}}_{s}(\tilde{t}) = \gamma \alpha_{g} r_m \tilde{W}_a k_\theta \tilde{G}_{bt}(\tilde{t}),\end{equation}

where $\tilde{G}_{bt}$ is the beam solar radiation, $W_a$ is the collector width, $\gamma$ is the intercept factor, $\alpha_g$ is the absorbance of glass cover, $r_m$ is the specular reflectance of the mirror and $k_\theta$ is the incident angle modifier.

The heat loss by convection is given by [1]:

(2.8) \begin{equation}\tilde{\dot{q}}_{conv} =\tilde{h}_w \left(\tilde{T} - \tilde{T}_a\right),\end{equation}

where $\tilde{T}_a$ is the ambient temperature and $\tilde{h}_w$ is the convective heat transfer coefficient.

The radiation heat loss is described by [1]:

(2.9) \begin{equation}\dot{q}_{rad} = \epsilon_{rad} \tilde{\sigma}\!\left(\tilde{T}^4 - \tilde{T}_{sky}^4\right),\end{equation}

where $\tilde{T}_{sky}$ is the sky temperature assumed to be [1]:

(2.10) \begin{equation}\tilde{T}_{sky} = 0.0552 \, \tilde{T}_a\end{equation}

with $\tilde{\sigma}$ as Stefan–Boltzmann constant and $ \epsilon_{rad}$ as the emittance.

The equations (2.2)–(2.4) form a set of 3 equations for the 4 unknowns $\tilde{\rho}, \tilde{u}, \tilde{T}, \tilde{p}$ . We still need to define a constitutive relation for the medium under consideration. The fluid under consideration is a special oil called Terminol VP-1 [36]. The constitutive relation is only known by measured data (given in [36]). The temperature dependence of the data (for a given pressure) can be approximated by polynomials (as in [1]). In Figure 5(a), we can see that the dependence of the density on the temperature is almost linear. There is also a theoretical pressure dependence of the density, which was studied in detail in [27], combining the data at different pressure values from [36] with additional data for a very similar fluid, called Dow-Corning 210 H [11]. The resulting approximative relation derived in [27] reads as

(2.11) \begin{equation}\tilde{\rho}(\tilde{p},\tilde{T}) = \tilde{G} + \tilde{B} \tilde{p} + \tilde{C} \tilde{T} + \tilde{D} \tilde{p} \tilde{T}.\end{equation}

However, as we will see after scaling in Section 2.3, the pressure dependence is neglectable. Thus finally we work with a linear approximation of the density (as function of the temperature). See also Table A.2 for the parameters, Table A.3 and Figure 5 for the properties of the fluid.

2.2 Scaling

The governing equations are scaled in order to obtain a dimensionless mathematical model, to reduce the number of (dimensionless composed) parameters and to identify order of magnitudes in the various terms. The reference values are denoted by an index r, the dimensionless quantities have no index and no tildes, i.e. for a general quantity f

\begin{equation*}f(x,t) = \frac{\tilde{f}(\tilde{x},\tilde{t})}{f_r} = \frac{\tilde{f}(x x_r,t t_r)}{f_r}.\end{equation*}

Moreover, the scaling of the Therminol VP-1 constitutive law leads to the following nondimenstional parameters:

(2.12) \begin{equation}G = \frac{\tilde{G}}{\rho_r}, \quad B = \tilde{B} u_r^2, \quad C = \tilde{C}\frac{T_r}{\rho_r}, \quad D = \tilde{D} T_r u_r^2.\end{equation}

The reference values are associated with the real-world data of the considered power plant, in our case the Ouarzazate Noor I [1]. Table A.1 shows the typical reference values that we use for scaling and Table A.2 contains the Noor I model parameters.

From now on, we assume the cross-sectional area A to be constant in a single tube. After scaling the governing equations (2.2)–(2.4) together with (2.11) can be written in dimensionless form as follows

(2.13) \begin{align} \rho_t + \left(\rho u\right)_x =0\end{align}

(2.14) \begin{align} \left(\rho u\right)_t + (\rho u^2+\frac{1}{\epsilon}p)_x =-\frac{\eta}{D_i}\rho u|u| +\left(\mu u_x\right)_{x} \end{align}

(2.15) \begin{align} (\rho T+ \epsilon \delta \frac{1}{2}\rho u^2)_t + (\rho u T+ \epsilon \delta \frac{1}{2} \rho u u^2 + \delta p u)_x&=\epsilon\delta\frac{\eta}{D_i}\rho u^2|u| + \left(\mu u u_x \right)_x+ \left(k T_x \right)_{x} + \frac{\kappa}{A}\dot{q}_{s} \nonumber \\[5pt]&\quad- \kappa\!\varsigma\frac{U_o}{A}\left(T^4 - T_{sky}^4\right) - \kappa \varrho \frac{U_o}{A}\left(T - T_{a}\right)\end{align}

(2.16) \begin{align} \rho(p,T) = G + \epsilon B p + C T + \epsilon D p T,\end{align}

where the dimensionless coefficients and its order of magnitude are expressed in Table 1. We see that we have series of small and very small parameters. This will be used to simplify the model.

Table 1. Dimensionless composed parameters with their order of magnitude

2.3 Asymptotic analysis

In this section, we look at the order of magnitudes (see also [27]). We see that the parameters $\mu$ , k, B and D are of order $\mathcal{O} (10^{-7})$ and below. Therefore, we skip the corresponding terms. Not surprisingly, this refers to viscous and to heat conducting terms (in this 1d approach). And it simplifies the constitutive law to a linear relation. With this we are left to

(2.17) \begin{align} \rho_t + \left(\rho u\right)_x =0\end{align}

(2.18) \begin{align} \left(\rho u\right)_t + (\rho u^2+\frac{1}{\epsilon}p)_x = -\frac{\eta}{D_i}\rho u|u|\end{align}

(2.19) \begin{align}(\rho T+ \epsilon \delta \frac{1}{2}\rho u^2)_t + (\rho u T+ \epsilon \delta u \frac{1}{2} \rho u^2 + \delta p u)_x&= \epsilon \delta \frac{\eta}{D_i} \rho u^2|u| + \kappa\frac{1}{A}\dot{q}_{s} - \kappa\varsigma\frac{U_o}{A}\left(T^4 - T_{sky}^4\right) \nonumber \\[5pt]&\quad - \kappa\varrho\frac{U_o}{A}\left(T - T_{a}\right)\end{align}

(2.20) \begin{align} \rho(p,T) =G + C T.\end{align}

Now we take advantage of the fact that $\epsilon$ and $\delta$ are always small in this application, everywhere and at any time. This is used to do an asymptotic approximation in terms of the small parameters $\epsilon$ and $\delta$ . We develop the quantities asymptotically with respect to $\epsilon$ and $\delta$ , i.e. for $f=(\rho, u, T, p)$ we write

(2.21) \begin{equation}f(x,t) = f_0 (x,t) + \epsilon f_{10} (x,t) + \delta f_{11} (x,t) + \mathcal{O}(\epsilon^2, \epsilon \delta, \delta^2).\end{equation}

In equation (2.14), the leading order term on the left $(p_0)_x$ (order $\epsilon^{-1}$ ) has as only relevant counterpart on the right the term $-\epsilon \frac{\eta}{D_i} \rho_0 u_0|u_0|$ . For big longitudinal extensions $x_r$ , small diameters $D_i$ and/or big friction coefficients $\xi$ we have

(2.22) \begin{equation}\epsilon \frac{\eta}{D_i} \approx \mathcal{O}(1),\qquad \mbox{ or } \qquad\frac{\eta}{D_i} \approx \mathcal{O}(\epsilon^{-1}).\end{equation}

This reflects the fact that friction effects cannot be neglected in this application. Therefore in leading order (in $\epsilon$ and $\delta$ ), we have the equations

(2.23) \begin{align} (\rho_0)_t + \left(\rho_0 u_0\right)_x = 0 \end{align}

(2.24) \begin{align} (p_0)_x = -\varepsilon \frac{\eta}{D_i} \rho_0 u_0|u_0| \end{align}

(2.25) \begin{align} (\rho_0 T_0)_t + (\rho_0 u_0 T_0)_x =\kappa\frac{1}{A}\dot{q}_{s}- \kappa\!\varsigma\frac{U_o}{A}\left(T_0^4 - T_{sky}^4\right)- \kappa\varrho\frac{U_o}{A}\left(T_0 - T_{a}\right)\end{align}

(2.26) \begin{align} \rho_0 = G + C T_0.\end{align}

This model includes as physical most relevant phenomena friction losses in the momentum balance, the solar input and the radiative and convective heat losses in the energy balance. We can rewrite and simplify the model using (2.26) in (2.25). This gives (dropping the index $()_0$ )

(2.27) \begin{eqnarray}\rho_t+(\rho u)_x =0\end{eqnarray}

(2.28) \begin{eqnarray}p_x =-\frac{\epsilon \eta}{D_i} \rho u |u|\end{eqnarray}

(2.29) \begin{eqnarray}u_x =-\frac{C \kappa}{\rho^2}\left(\frac{1}{A}\dot{q}_{s}- \varsigma\frac{U_o}{A}\left(T^4 - T_{sky}^4\right)- \varrho\frac{U_o}{A}\left(T - T_{a}\right)\right) \end{eqnarray}

(2.30) \begin{eqnarray}T = \frac{\rho-G}{C}. \end{eqnarray}

This is a coupled system of 3 first order nonlinear PDE’s and a linear algebraic relation. Therefore, we need initial data for the quantities

(2.31) \begin{equation}p(x,0)=p_I, \qquad T(x,0)=T_I,\qquad u(x,0)=u_I, \end{equation}

where the initial condition for the density $\rho$ is then given by (2.37). And we need 3 physically meaningful boundary conditions. Typically a certain pressure is applied to the tube and as a consequence a flow is induced. This is realised by prescribing the pressures at the inlet and the outlet and consequently the density at the inlet (inflow condition)

(2.32) \begin{eqnarray}p(0,t)=p_l,\quad p(1,t)=p_r\end{eqnarray}

(2.33) \begin{eqnarray}\rho(0,t)=\rho_l \mbox{ if }u(0,t) >0\quad \mbox{ and/or }\qquad \rho(1,t)=\rho_r\mbox{ if } u(l,t) <0.\end{eqnarray}

2.4 Network

Since a typical parabolic trough power plant is composed as a network of pipes, we now consider the network setting (see also [30]). We have – including an inlet and an outlet tube – $n_p$ tubes connected in $n_v$ nodes. For each tube $i=1,...,n_p$ we have density $\rho^i$ , velocity $u^i$ , temperature $T^i$ and pressure $p^i$ . In each of the tubes various parameters may differ, e.g. the length $L^{i}$ , the cross-sectional area $A^i$ and the diameter $D^i$ .

The fluid in each tube is described by (2.34)–(2.37)

(2.34) \begin{eqnarray}\rho^i_t+u^i\rho^i_x = {\frac{C \kappa}{(\rho^i)}}\left(\frac{1}{A}\dot{q}_{s}- \varsigma\frac{U_o}{A}\left((T^i)^4 - T_{sky}^4\right)- \varrho\frac{U_o}{A}\left(T^i - T_{a}\right)\right) \end{eqnarray}

(2.35) \begin{eqnarray}\!\!\!\!\!\!\!\!\!\!\!\!p^i_x= -\frac{\epsilon \eta}{D^i_i} \rho^i u^i |u^i| \qquad\qquad\qquad\qquad\qquad\quad \end{eqnarray}

(2.36) \begin{eqnarray} u^i_x =- {\frac{C \kappa}{(\rho^i)^2}}\left(\frac{1}{A}\dot{q}_{s}- \varsigma\frac{U_o}{A}\left((T^i)^4 - T_{sky}^4\right)- \varrho\frac{U_o}{A}\left(T^i - T_{a}\right)\right) \end{eqnarray}

(2.37) \begin{eqnarray}T^i = \frac{\rho^i-G}{C}.\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad \end{eqnarray}

The parameters are given in Table A.4. We have boundary conditions for the network at the inlet and the outlet tube, i.e. we add boundary conditions for the pressure $p_l$ and the density $\rho_l$ at the inlet tube and for the pressure $p_r$ at the outlet tube

(2.38) \begin{eqnarray}& p_{l} = p_{in}(0,t) & \qquad p_r = p_{out}(L_{out},t)\end{eqnarray}

(2.39) \begin{eqnarray} \;\;\rho_l = \rho_{in}(0,t). \qquad\qquad\qquad\qquad \end{eqnarray}

The inflow temperature $T_{in}$ depends via (2.37) from the inflow density $\rho_{in}$ . The subscripts in and out refer to the tubes connected to the outside world with in or outgoing flow.

We complete the system with intial condition for temperature, pressure and velocity

(2.40) \begin{eqnarray}p^i(x,0)=p^i_I(x), \ \ \ T^i(x,0)=T^i_I(x),\ \ \ u^i(x,0)=u^i_I(x),\end{eqnarray}

where the initial density is again given by (2.37). In the network, we have to define node conditions which act as boundary conditions for the internal (or non-external) ends of the tubes, i.e. those parts ending in a network node. We need 3 boundary conditions per tube – e.g. 2 pressure and 1 density inflow condition – and thus in total $3 \cdot n_p$ boundary conditions in the whole network.

For the node dynamics, we make the assumptions that the fluid is homogeneoulsy mixed such that

1. (i) we have a single pressure value $P_j = P_j(t)$ at each node j,

2. (ii) we have a single density value $\rho_{inflow,j}=\rho_{inflow,j}(t)$ for each tube with (node-)outgoing flow at node j.

This reduces the unknowns to a single pressure values and a single density values at each node.

We formulate physical reasonable and necessary node conditions (like in [15] or [7]):

1. (iii) mass is conserved in the nodes, i.e. the in and outgoing mass fluxes are balanced

   (2.41) \begin{equation}\sum_{sign(i,j)\neq 0}{sign(i,j)\rho_i u_i A_i}=0\end{equation}

2. (iv) the inner energy is conserved in the nodes, i.e. the in and outgoing fluxes for the inner energy are balanced

   (2.42) \begin{equation}\sum_{sign(i,j)\neq 0}sign(i,j)\rho_i u_i T_i A_i = 0.\end{equation}
   Remember that the kinetic energy parts were of orders of magnitude smaller and not relevant in this leading order approximation.

where the function $sign(\text{tube } i, \text{node } j)$ indicates at which end of the tube i the node j is located

\begin{eqnarray*}sign(\text{tube } i \text{, node }j)=\begin{cases} 1 & \text{ left end of the tube} \; {i} \; \text{ends in node} \; {j} \\[5pt]-1&\text{ right end of the tube} \; {i} \; \text{ends in node} \; {j}\\[5pt]0&\text{ no connection of tube} \; {i} \; \text{with node} \; {j}.\end{cases}\end{eqnarray*}

Using (2.37) in the inner energy balance at the nodes (2.42), then using (2.41) and the fact all outflowing tubes have the same (well mixed) denisty, we obtain at every node $j=1,...,n_v$

(2.43) \begin{eqnarray} \sum_{sign(i,j) u_i <0} sign(i,j)(\rho_i)^2 u_i A_i \end{eqnarray}

(2.44) \begin{eqnarray}= - \sum_{sign(i,j) u_i >0} sign(i,j)(\rho_i)^2 u_i A_i \end{eqnarray}

(2.45) \begin{eqnarray}= - \rho_{inflow,j}\sum_{sign(i,j) u_i >0} sign(i,j) \rho_i u_i A_i.\end{eqnarray}

Using again (2.41), we obtain coupling conditions for the outflow density in node j

(2.46) \begin{eqnarray}\rho_{inflow,j}= \frac{\sum\limits_{sign(i,j)u_i<0}\left.\left(\rho_i^2 u_i A_i\right)\right|_{x=L_i}}{\sum\limits_{sign(i,j)u_i<0} \left.\left(\rho_{i} u_i A_i\right)\right|_{x=L_i}}, \ \ j=1,...,n_v\end{eqnarray}

which only depends on inflow data at the related node. We remark that (2.46) by construction guarantees the conservation of mass in every node. With this the solution of the network problem is prepared. In Section 3, the details of the numerical realisation of the solution of the network problem are presented.

We start with given data at a certain time (e.g. initial data) on every tube of the network and a set of pressures $P_j$ at the nodes. We perform a step forward (in time) in the continuity equation (2.34). Then we solve the combined network problem

(2.47) \begin{gather}\Delta p^i =-\frac{\epsilon \eta}{D^i_i} \int_0^1{\rho^i u^i|u^i| dx}\ \ i=1,...,n_p\end{gather}

(2.48) \begin{gather}\sum\limits_{sign(j,i)u^i>0}\left(\rho^i u^i A^i\right)_{x=0}=\sum\limits_{sign(j,i)u^i<0}\left(\rho^i u^i A^i\right)_{x=L_i} \ \ j=1,...,n_v,\end{gather}

where $\Delta p_i$ denotes the pressure differences of right an left node pressure acting on tube i. This problem promises to behave well since we have – at least formally – a ‘monotonicity’ behavior of the node pressures at every node. A strategy o solve this problem is in a first step to guess the pressures $P_j$ at all the nodes. The higher we chose the pressure $P_j$ in a node j the higher is the outflow in the outflowing tubes and the lower is inflow in the inflowing tubes (at node j). As long as there is no change in the flow direction with respect to the tubes linked to node j, there is exactly on pressure value $P_j$ such that outflow balances inflow at that node. This monotinicity remains even true if by lowering (highering) the node pressure $P_j$ an outflowing (inflowing) tube switches and becomes an inflowing (outflowing) tube. Insofar we expect a unique solvability with respect to the pressures. This algorithm will be used in Section 3 to solve the network problem.

2.5 Power optimisation

In the previous section, we studied and performed direct numerical simulations on a network for a given set of parameters and initial and boundary data. As already mentioned, we would like to go further and – as required for the application – optimise the output. A natural quantity to be optimised is the power output. In an operational phase, an important control variable is the pressure drop between the outlet and inlet to the solar field, $\Delta p$ .

The power P is the energy converted per unit of time, $P=\frac{\Delta E}{\Delta t}$ , where E is the energy. The thermal (gross) output of a power plant is the amount of heat generated per unit of time. On the other hand, we have to employ pumping power to push the fluid trough the network. At the end the most interesting quantitiy is the net power, i.e. the gross power minus the necessary power for running the system.

We start with the thermal output. The amount of heat $\tilde{P}_{thermal}$ , which is generated by a temperature difference is calculated by

(2.49) \begin{equation}\tilde{P}_{thermal} = \tilde{A} \tilde{\rho} \tilde{u} \cdot (\tilde{c}_v(\tilde{T}_{out}) \tilde{T}_{out} - \tilde{c}_v(\tilde{T}_{in}) \tilde{T}_{in}),\end{equation}

where $\tilde{T}_{in},\tilde{T}_{out}$ are the in and outlet temperatures, $\tilde{c}_v$ the (temperature-dependent) specific heat and $\tilde{A} \tilde{\rho} \tilde{u}$ the fluid mass push through the tubes per time. Here, we assume incompressibility and thus the inflowing mass is equal to the outflowing mass. Applied to our network we get for the relevant outgoing tube or tubes (if more than one outgoing tube)

(2.50) \begin{equation}\tilde{P}_{thermal}=\left.\left(\tilde{A}_{out}\tilde{\rho}_{out}\tilde{u}_{out}\cdot\left(\tilde{c}_v(\tilde{T}_{out})\tilde{T}_{out}-\tilde{c}_v(\tilde{T}_{in}) \tilde{T}_{in}\right)\right)\right|_{\tilde{x}=\tilde{L}_{out}}.\end{equation}

Now we integrate $\tilde{P}_{thermal}$ over time to account for the performance over the entire simulation period. In addition, we take into account the fact that there can be several output pipes and form the sum over the corresponding pipes

(2.51) \begin{equation}\tilde{J}_{thermal}=\int_0^{\tilde{t}_{end}}\sum_{i=out}\left.\left(\tilde{A}^i_{out}\tilde{\rho}^i_{out}\tilde{u}^i_{out}\cdot\left(\tilde{c}_v(\tilde{T}^i_{out})\tilde{T}^i_{out}-\tilde{c}_v(\tilde{T}_{in}) \tilde{T}_{in}\right)\right)\right|_{\tilde{x}=\tilde{L}_{i}} dt\end{equation}

representing the gross power generated in the solar field by solar thermal heating. The power functional is still subject to dimensions.

Now we take into account in our functional the needed pumping power to generate the pressure. This is calculated according to Sterner et al. [33], as follows

(2.52) \begin{equation}P_{pump} = \left.\tilde{A}_{in}\tilde{u}_{in}\right|_{\tilde{x}=\tilde{0}}\frac{\Delta \tilde{p}}{\eta_{pump}},\end{equation}

where $0 < \eta_{pump}< 1$ designates the efficiency of the pump and $\Delta p = p_{in} - p_{out}$ the produced overpressure between inlet and outlet pipes. Since pumps act always at the inlet(s), we have (again integrated over the simulation period)

(2.53) \begin{equation}\kern-1pt\tilde{J}_{pump}=\int_0^{\tilde{t}_{end}}\sum_{i=in}\left.\tilde{A}^i_{in}\tilde{u}^i_{in}\right|_{\tilde{x}=\tilde{0}}\frac{\tilde{p}^i_{in}-\tilde{p}_{out}}{\eta_{pump}}dt.\end{equation}

We assume in our model that the pump is installed at the entrance of the first pipe (only one entrance pipe). Overall, we can now set up the following net power functional $\tilde{J}$ :

(2.54) \begin{eqnarray}\kern-1pt\tilde{J} & = &\int_0^{\tilde{t}_{end}}\left(\sum_{i=out}\left.\left(\tilde{A}^i_{out}\tilde{\rho}^i_{out}\tilde{u}^i_{out}\cdot\left(\tilde{c}_v(\tilde{T}^i_{out})\tilde{T}^i_{out} -\tilde{c}_v(\tilde{T}_{in}) \tilde{T}_{in}\right)\right)\right|_{\tilde{x}=\tilde{L}_{i}}\right.\nonumber \\[5pt]&& -\left.\left.\tilde{A}_{in}\tilde{u}_{in}\right|_{\tilde{x}=\tilde{0}}\frac{\tilde{p}_{in}-\tilde{p}_{out}}{\eta_{pump}}\right)dt.\end{eqnarray}

We mention that the (in time integrated) net power functional represents the amount of net energy produced by the power station over the time intervall $[0, \tilde{t}_{end}]$ .

It remains to scale the functional using the reference values from the Subsection 2.2 (the reference value for $\tilde{J}$ is given by $T_{r} \rho_{r} u_{r} A_r$ ). We obtain ( $c_v$ denotes a scaled version of the heat capacity)

(2.55) \begin{eqnarray}J & = &\int_0^{t_{end}}\left(\sum_{i=out}\left.\left(A^i_{out} \rho^i_{out} u^i_{out}\cdot\left(c_v(T^i_{out})T^i_{out} -c_p(T_{in}) T_{in}\right)\right)\right|_{x=L_{i}}\right.\nonumber \\[5pt]&& -\left.\left.A_{in} u_{in}\right|_{x=0}\frac{p_{in}- p_{out}}{\eta_{pump}}\right)dt.\end{eqnarray}

The optimisation (maximisation) of the presented net (integrated) power functional is done under the additional constraint $\max T^{out} < T_{c}$ since the HTF (in our case Therminol VP-1) should not go above the critical temperature $T_c = \frac{663.15 \ K}{T_r}$ ( $390^\circ\mathrm{C}$ ). In the next section, we show how the numerical realisations of this optimisation is done (in MATLAB based on optimisation toolbox).

3 Numerical algorithm

In this section, the algorithm for a single pipe is first introduced accompanied with the initial and boundary values. The coupling conditions are then used for the calculation of missing boundary values of all pipes.

In order to solve the system of equations (2.34)–(2.37), a temporal and spatial discretisation is required. The discrete grid points for space and time are denoted, respectively, by $x_1, \ \ldots, \ x_{n_x}$ and $t_1, \ ..., \ t_{n_t}$ . The discrete variables are accordingly provided with the indexing (n, j), where the first index is relating to the time coordinate, the second to the spatial coordinate.

In the following, the index for the pipe number is dropped as the algorithm is introduced for a single pipe. The explicit upwind method is used to solve the mass conservation, equation (2.34), where CFL condition is preserved [30]. Calculation of the discrete temperature values result from the linear relationship (2.37) between T and $\rho$ . On the other hand, the velocity solution can be obtained by integrating (2.36)

(3.1) \begin{equation}u(t,x) = v(t) + \int_0^x -\frac{C \kappa}{\rho^2}\left(\frac{1}{A}\dot{q}_{s}- \varsigma\frac{U_o}{A}\left(T^4 - T_{sky}^4\right)- \varrho\frac{U_o}{A}\left(T - T_{a}\right)\right) ds,\end{equation}

where $v=v(t)$ (space independent) represents the value of the velocity at the pipe inlet. To calculate $v=v(t)$ , we use the pressure equation (2.35). If we integrate equation (2.35) over the whole pipe (from $x=0$ to $x=L$ ), insert (3.1) and the two pressure boundary values, we get an equation for $v=v(t)$

(3.2) \begin{align} & p(1,t) - p(0,t) \nonumber\\[5pt] & =-\frac{\epsilon \eta}{D_i}\int_0^1 \rho \left(v(t) + \int_0^x -\frac{C \kappa}{\rho^2}\left(\frac{1}{A}\dot{q}_{s}- \varsigma\frac{U_o}{A}\left(T^4 - T_{sky}^4\right)- \varrho\frac{U_o}{A}\left(T - T_{a}\right)\right) ds\right) \nonumber \\[5pt]& \left|v(t) + \int_0^x -\frac{C \kappa}{\rho^2}\left(\frac{1}{A}\dot{q}_{s}- \varsigma\frac{U_o}{A}\left(T^4 - T_{sky}^4\right)- \varrho\frac{U_o}{A}\left(T - T_{a}\right)\right) ds\right| dx. \end{align}

After the density, temperature and velocity are solved for the respective time step, the pressure can also be calculated. In order to take into account the two given pressure boundary values, the pressure equation (2.35) is differentiated with respect to x

(3.3) \begin{equation}p_{xx} = -\frac{\epsilon \eta}{D_i} \left(\rho u |u|\right)_x.\end{equation}

Finally, the pressure values are obtained by solving the above equation using second-order central finite different scheme.

This is the algorithm for a single pipe. The network approach still requires the numerical implementation of the coupling conditions in order to derive the required boundary values for the individual pipes.

As described in Subsection 2.4, a boundary value for the density is required for the explicit upwind method. For internal pipes, the boundary value for the density is determined via the equation (2.46) at every node $j=1,...,n_v$ . With the input density determined in this way for each pipe, the continuity equation can be solved using the explicit upwind methods.

Furthermore, the $n_v$ pressure values at the nodes are required to calculate the speed and the pressure in the individual pipes. For this, the following $n_v$ equations at the nodes and the $n_p$ equations for the pipes form a non-linear system of equations for $i=1,\ldots,n_p$ , and $j=1,\ldots,n_v$ .

(3.4) \begin{align}& \Delta p^i = -\frac{\epsilon \eta}{D_i} \int_0^{L_i} \rho^i \left(v^i(t) + \int_0^x - {\frac{C \kappa}{(\rho^i)^2}}\left(\frac{1}{A}\dot{q}_{s}- \varsigma\frac{U_o}{A}\left((T^i)^4 - T_{sky}^4\right)- \varrho\frac{U_o}{A}\left(T^i - T_{a}\right)\right)ds\right) \nonumber\\[5pt]&\quad \times\left|\left(v^i(t) + \int_0^x - {\frac{C \kappa}{(\rho^i)^2}}\left(\frac{1}{A}\dot{q}_{s}- \varsigma\frac{U_o}{A}\left((T^i)^4 - T_{sky}^4\right)- \varrho\frac{U_o}{A}\left(T^i - T_{a}\right)\right) dx\right)\right| dx\end{align}

(3.5) \begin{align} & \sum_{sign(i,j) u^i > 0} \rho^i v^i A^i \nonumber \\[5pt]& = \sum_{sign(i,j) u^{i} < 0} \rho^i \left(v^i + \int_0^{L_i} - {\frac{C \kappa}{(\rho^i)^2}}\left(\frac{1}{A}\dot{q}_{s}- \varsigma\frac{U_o}{A}\left((T^i)^4 - T_{sky}^4\right)- \varrho\frac{U_o}{A}\left(T^i - T_{a}\right)\right)ds\right) A^i,\end{align}

where $n_v+n_p$ equations are solved to obtain the $n_v$ pressure values at the nodes $P_j$ and the $n_p$ input velocities $v^i$ of the pipes.

The required boundary values for each individual pipe within the entire network are then available and the quantities velocity, pressure and temperature can be determined according to the procedure explained above.

4 Numerical simulation

In the following, the presented model system is simulated for a set of different examples. A number of MATLAB codes are implemented in order to study the different numerical cases.

At first, a simple case is considered where the presented model is numerically solved for one single pipe and validated with existing results in literature. In the second case, a solar field collector row composed of 8 collector pipes is considered.

The solar radiations and ambiant temperature for a typical sunny day in the region of Ouarzazate (Morocco) are considered [1]. The hourly variation of the considered data are plotted in Figure 4 between the sunrise and the sunset, where the observed peak solar radiation is around $1000\,\mathrm{W/m}^2$ and the peak ambient temperature is around 308 K. Some of the coefficients involved in the governing equations are determined using the temperature-dependent properties of the Therminol VP-1 in order to obtain accurate numerical results at each time step. The thermal properties of the HTF for a range of temperatures from 285 to 693 K are depicted in Figure 5 alongside the linear interpolation expressed in Table A.3.

Figure 4. Hourly variation of solar radiation (left axis) and ambient temperature (right axis) in the region of Ouarzazate [1].

Figure 5. Linear fitting of Therminol VP-1 thermal properties from [36].

4.1 Validation of the single pipe model

First, both models should be examined with regard to their plausibility. In the first simulation, we want to represent a situation that is as simple as possible. A simple fluid transport through the pipe from left to right, the results obtained by the one pipe model are compared with existing results from Allouhi et al. [1]. We will always select the same data in the subsequent simulation series, i.e. model parameters, calculation parameters, initial and boundary conditions for the two systems in order to be able to compare them. The associated data is listed in the previous sections.

For the case of a single collector pipe, the obtained results are presented to validate the derived model by comparing the predicted temperature variations with existing values in the literature. The inlet temperature is set at 320 K ( $46,85^\circ\mathrm{C}$ ). Figure 6, shows the variation of the obtained results for the fluid properties, namely, the density, velocity, pressure and temperature as function of space (along the pipe) and time (hours of the day). It is noticed from the different plotted results, that at the peak value of the solar radiation during the day, the minimum value of the density and the maximum values of the velocity and temperature are observed. Moreover, the collector generates a maximum temperature value (375 K) at the outlet of the pipe when the solar radiation is at the peak value during the day, similar values are obtained in [1] for the same NOOR I model parameters. Compared to other models in the literature and the work of Allouhi et al. [1], the presented model gives more and new insight to the spacial and time variation of other fluid properties besides temperatures, such as velocity and density. This additional information can only be generated by a thermo fluid dynamics model as the one we use here in this paper.

Figure 6. Evolution of the different variables in the single pipe model.

4.2 Validation of the network model

In this case, a simple network model constituted of 8 collectors connected in series is assumed, where the same model parameters and properties are considered. In order to check the coupling conditions of the system (3.4)–(3.5), the numerical results for this model can be compared with the results of a one pipe model with the length of 8 collectors ( $8\times12.27$ m). It is observed from the results, that both systems achieve the same results, see Figure 7. The network model clearly keeps the continuity of the temperature variation from a pipe to an other. Similar to the single tube model, the maximum temperature value is at the outlet of the last pipe when the solar radiation reaches the peak value. This can be seen as a basic test for the network approach.

Figure 7. Temperature variation of collectors row using one pipe model (left) and network model (right).

5 Design and optimization: NOOR I power plant

5.1 In series design of the collector

The NOOR I power plant (Ouarzazate, Morocco) is considered as an application for the presented network model. The network case of eight pipes connected in series is a realistic case of the NOOR I plant, the model can be then used for further optimization of the power functional with respect to different model parameters. In this paper, the (optimal) pressure drop is the main output of the optimisation approach. For the usual design of each collector row, the pressure drop is optimised to maximise the net power output. Table 2 shows the optimisation iterations and the associated values of the pressure drop ( $\Delta p$ ) and power output (J). It is observed that the maximum output temperature exceeds the critical value of $663.15/\mathrm{K}$ ( $390/^\circ\mathrm{C}$ ) in the first 11 iterations, but then finally at iteration 12 it stays below the critical value and subsequently maximises the generated power under the additional constraint. Using the optimal value of the pressure drop $\Delta P = 8.53$ , the different fluid variables are plotted again with respect to time and space, see Figure 8. The results show a decrease in the minimum value of density and an increase in both maximum values of the velocity and temperature, especially at the output.

Table 2. Iterations to reach the optimal of the net power output satisfying also the upper limit for the temperature; iterations 1–11 do not satisfy the upper limit. The numbers in red shows the outlet temperatures that exceed the critical value of 663.15 K

Figure 8. Evolution of the different fluid variables for the obtained optimal pressure drop.

5.2 Complex design of collector rows

Usually collector rows are designed in a series form (Design 1). More complicated designs can be studied with respect to possible better net power output. Thus, three different additional designs are considered as shown in Figure 9. The first design (Design 2) consists of two parallel collector pipes at the beginning of the row and the rest are connected in series, while the second design (Design 3) consists of two parallel collector pipes at the beginning of the row and another two parallel collector pipes in the middle of the row, while the other pipes are connected in series. The third design (Design 4) consists of six pipes in parallel.

Figure 9. Skitches of the four studied network designs.

For the four designs, the same optimisation approach is used in order to obtain the optimal pressure drop with a maximum net power output. As shown in Table 3, the value of the pressure drop varies from one design to another with the minimum value observed for the D4 and the maximum value observed for the classical one d1.

Table 3. Optimal pressure drops, thermal, pumping and net power outputs for the different designs

The concluding observation from the obtained results is that a more complex designs of the power plant parabolic trough collectors can result in better output values under better conditions (e.g. lower temperatures). Therefore, further research and studies can be conducted in terms of network design and optimisation using simplified models similar to the one introduced in this paper.

In Figure 10, the fluid temperature variations for the different designs are presented. It is observed that for the parallel pipes, the temperature increases more than for a single pipe, while the output temperature is similar to design 1 with lower pressure drop.

Figure 10. Temperature variation through the collectors row for the different designs.

6 Conclusions

We consider the presented model as a significant step towards accurate modeling of parabolic trough power plants. The main physical phenomena are included. Nevertheless, the model is simple enough to allow for fast and robust simulations and optimisation tasks. The model was validated with data from existing parabolic trough power plants. The performed simulations show very good agreement with the available measurements.

However, we know that this is only a starting point. It is well known that modelling and simulation has become a significant tool for the design of complex technical systems. A parabolic trough power plant is such a complex technical system and its design and operation is all but trivial. In this context, mathematical models combined with fast and robust numerical and optimisaton approaches are needed for the design of the network (structure, substructure, pipe dimensions etc.), to evaluate non planar collector fields, to evaluate the losses and possible isolation measures, for introducing alternative fluids (e.g. nanofluids), to identify optimal settings in the operational phase etc.

Clearly the presented model is far from being overarching. It might be that phenomena which are not yet identified to be relevant are not yet described. This model cannot capture certain details due to its 1 dimensional character, i.e. details of the flow in a junction. It might be that at some point the model is not precise enough in the description, i.e. the heat transfer mechanism from the radiation into the fluid.

And there are mathematical and numerical challenges to be faced. So far there is no well-posedness theory for the model. This not surprising since we have a system of nonlinear coupled PDEs on a network. However, first steps have to be done. The numerical algorithm is kept as simple as possible, there is room of improvement. Adaptive methods could help to accelerate in particular in case of the optimisation task. Up to now the optimisation is done using standard tools. This works here for rather simple examples, it might be that i.e. for a large multiparameter optimisation examples we run into problems. Thus, more sophisticated approaches have to be considered, i.e. a Lagrangian based method.

Nevertheless, we hope that our approach can become a valuable and reliable tool in designing and operating parabolic trough power plants.