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 [Reference Kalogirou20, Reference Allouhi, Benzakour Amine, Saidur, Kousksou and Jamil1].

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

(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 [Reference Allouhi, Benzakour Amine, Saidur, Kousksou and Jamil1]:

(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 [Reference Allouhi, Benzakour Amine, Saidur, Kousksou and Jamil1]:

(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 [Reference Allouhi, Benzakour Amine, Saidur, Kousksou and Jamil1]:

(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 [Reference Allouhi, Benzakour Amine, Saidur, Kousksou and Jamil1]:

(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 [Reference Allouhi, Benzakour Amine, Saidur, Kousksou and Jamil1]). 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 [Reference Parfenov27], combining the data at different pressure values from [36] with additional data for a very similar fluid, called Dow-Corning 210 H [Reference Fakhreddine and Zoller11]. The resulting approximative relation derived in [Reference Parfenov27] 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 [Reference Allouhi, Benzakour Amine, Saidur, Kousksou and Jamil1]. 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 [Reference Parfenov27]). 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 [Reference Schuster30]). 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 [Reference Gasser and Kraft15] or [Reference Brouwer, Gasser and Herty7]):

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. [Reference Sterner33], 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).

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. [Reference Allouhi, Benzakour Amine, Saidur, Kousksou and Jamil1]. 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 [Reference Allouhi, Benzakour Amine, Saidur, Kousksou and Jamil1] for the same NOOR I model parameters. Compared to other models in the literature and the work of Allouhi et al. [Reference Allouhi, Benzakour Amine, Saidur, Kousksou and Jamil1], 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.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.