Forced Laminar Flow in Pipes Subjected to Asymmetric External Conditions: The HEATT© Platform for Online Simulations

Abstract

This chapter studies the fluid flow within pipes subjected to thermal asymmetrical boundary conditions. The phenomenon at hand takes place in many real-world industrial situations, such as solar thermal devices, aerial pipelines. A steady-state analysis of laminar forced-convection heat transfer for an incompressible Newtonian fluid is studied. The fluid is considered to flow through a straight round pipe provided with straight fins. For the case studied, axial heat conduction in the fluid has been considered and the effects of the forced convection have been considered to be dominant. A known uniform temperature field is applied at the upper external surface of the assembly. The 3D assembly has been created combining cylindrical and Cartesian coordinates. The governing differential equation system is solved numerically through suitable discretization in a set of different finite volume elements. The results are shown through the thermal profiles in respect of longitudinal and radial-azimuthal coordinates and the problem characteristic length. To facilitate the resolution of this phenomenon, an open computing platform called HEATT©, based on this model, has been developed, and it is also shown here. The platform is now being built and is expected to be freely available at the end of year 2022.

1. Introduction

Laminar forced convection mechanisms take place in many industrial installations where the processes involve the use of fluids in pipes. Numerous applications exist of such flows, including flat solar thermal collectors [1, 2, 3, 4], solar trough devices [5], nanofluids [4, 6], mini and micro channels [7] and a wide variety of heat exchangers [8]. Currently, the interest on the studies of thermal behaviour in oil or gas pipelines is growing due to the international context, as well as the study of hydrogen through liquefied petroleum gases (LPG) pipelines that are already installed [9, 10, 11, 12].

From the point of view of the analysis of fluids’ heat exchange on pipelines, it is technically relevant to consider that flow and heat processes occur simultaneously, i.e. they are coupled, which increases the complexity of the process. Understanding flow behaviour under these conditions is key to pipeline design and device efficiency.

The fundamentals of the thermal mechanisms involved on fluid flow have been extensively studied [13, 14, 15]. The complexity in their analysis and the geometry of practical applications has made common its resolution using experimental and numerical studies, examples of them are [16, 17].

One of the more significant models of fluid behaviour within a pipe is the one known as Graetz Problem (GP), stated by Graetz in 1882 [18], where in a given point of the pipe a fluid flowing in laminar forced flow is subject to a sudden change in its external boundary conditions, either temperature or heat flow. Graetz proposed a bi-dimensional approach, neither considering pipe nor axial fluid conduction, which was analytically solved. More than a 100 years after his work, Graetz’s problem continues to receive the attention of researchers. In the present century, some researchers have extended this known problem to take into account both the physical presence of the pipe and the axial fluid conduction (conjugate extended Graetz problem), finding their results by analytical [19] or numerical procedures [20, 21], also some have included transient process [22]. Other researchers have studied the flow with periodically varying inlet temperature in pipes of different shapes [23], pipes subjected to a sudden [24] or periodical change [25] in external heat or ambient temperature, etc. Also significant is the formulation of the concept of characteristic length of the process, carried out by discriminated dimensional analysis by Seco-Nicolás et al. [26]. All these studies, based on radial symmetry, assume the 2D hypothesis. However, in certain cases such relatively simplified models do not provide results as accurate as those obtained through tri-dimensional numerical models.

This chapter faces the problem of fluid flow within pipes subjected to thermal asymmetrical boundary conditions which take place in many real industrial situations such as those related to solar thermal devices, aerial pipelines subjected to external temperatures, etc. Other examples and attempts to solve this problem can be seen on: [27, 28, 29, 30, 31]. Despite the asymmetry of the problem, much simpler bi-dimensional models are currently used for pipe design purposes, ignoring the important consequences of the asymmetry that it exists.

This work presents a steady-state analysis of the laminar forced-convection heat transfer process for a liquid flowing through a straight round pipe when radially asymmetrical external conditions are applied to the tube’s external surface (a known uniform temperature to the upper surface and adiabatic condition to the lower) and taking into consideration axial heat conduction in the fluid.

A governing differential equation system is coupled to the Laplace equation for the solid and is solved numerically through suitable discretisation in a set of different finite volume elements, considering the axial heat conduction in the fluid, but neglecting the heat generation by viscous dissipation, the buoyancy effects or the variation of the thermal properties of the materials.

Many techniques have been developed to simulate convective flows using finite element techniques [32], finite difference solvers [33], method of lines [34] and many others. In the present case to evaluate the proposed model Network Simulation Method (NSM) [35, 36, 37], a powerful numerical methodology has been chosen. This method, based on a finite-difference scheme, can virtually solve any ordinary and partial differential equation. The method starts from suitable discretisation of the problem [37, 38, 39]; after which an equivalent electric circuit of the process is built, including boundary conditions, the so-called network model. Finally, appropriate software is used to solve the circuit, from whose results the thermal response is obtained. NSM has been chosen in this case because it yields results that have been seen to be as accurate as those obtained with CFD software using significantly fewer computational resources. Open-source software NGSpice [40], originally created for electric circuit analysis, has been used in our case to solving this problem; also proprietary software such as Pspice© [41] can be used for this purpose. The model has been validated by experimental results obtained in an experimental solar thermal plant [42].

Related to this model, an open computing platform called HEATT©, based on this model and solved using the NTM, is now being built, which will allow online calculation of flow within pipes subject to complex thermal conditions. The platform is expected to be freely available to the public before the end of 2022.

2. Physical and mathematical model

The set-up studied here is composed by a flat plate welded to a round duct as shown in Figure 1. The flat plate acts as a fin. Its addition to the model is due to the fact that many of the pipes subjected to thermal stress incorporate fins to improve heat exchange.

Figure 1.

Simulation model outline.

Consider an isotropic fluid at uniform temperature, T₁, and a certain velocity profile inside a tube whose whole external surface is also at temperature T₁. At z = 0 the fluid enters a long duct, whose upper half surface is maintained from this point at a constant temperature T₂ (T₂ > T₁), while its lower surface is thermally isolated (Figure 1); these conditions are maintained along the whole duct length (z > 0). The plate is maintained at a constant temperature T₂ on its upper surface, and it is insulated on its lower surface and extremes in the same conditions of the pipe. Heat is transferred throughout the plate to the duct by conduction. Therefore, the plate acts as a fin.

Consequently, at z = 0, the fluid suffers a sudden change in its boundary conditions and heat is transferred through the tube by conduction, and then by convection to the fluid, in which conduction takes place in both radial and axial directions.

The fluid working conditions are similar to those of the bi-dimensional problem stated by Graetz [18], except for the radial asymmetry of the temperature boundary condition caused by the radially asymmetric thermal condition in the external surface of the pipe and the presence of the fin. These boundary conditions make it necessary to formulate a 3D model.

The geometry of the problem requires incorporating both cylindrical and Cartesian coordinate systems. Regarding cylindrical coordinates for the round duct, direction z is located parallel to the axis of the pipe, the r vector is normal to it, and the third dimension is described by angle φ. Regarding Cartesian coordinate system for the fin, direction z is also located parallel to the axis of the pipe, and directions x and y are orthogonal to it, as can be seen in Figure 2.

Figure 2.

3D Duct problem geometry. (a) Section planes outline. (b) 2D longitudinal section β including velocity profile. (c) 2D transversal section α.

As regards the material parts, the assembly consists of a round pipe, with k_s being the constant conductivity and e_p the constant thickness. The pipe is assumed straight and non-deformable. Heat generation by viscous dissipation or other sources is not considered. As regards the time domain, the model is considered to be at stationary state.

It is assumed that the fluid flows in laminar-forced convection in stationary regime, while its thermal properties (density, ρ_f, specific heat, c_f, and thermal conductivity, k_f) remain unchanged in their values at temperature T₁.

For the velocity profile, a polynomic profile of grade 10 has been chosen. This kind of profiles can be found in certain devices, such as inclined tubes [43]. The velocity will be considered invariable along the whole pipe (Figure 2b).

The viscous buoyancy-driven heat transfer is considered negligible due to dominant effects of the studied forced convection in a laminar and incompressible Newtonian fluid flow, which is the case of water and other fluids in certain conditions [44].

Under these conditions, the tri-dimensional equations that govern the coupled system are (on cylindrical coordinates) [45, 46]:

Equation of the solid (pipe) region. Cylindrical coordinates:

$$\frac{1}{r}\left[\frac{\partial }{\partial r}\left(r{k}_s\frac{\partial T_s}{\partial r}\right)\right]+\frac{1}{r^2}\frac{\partial }{\partial \phi }\left(k_s\frac{\partial T_s}{\partial \phi }\right)+\frac{\partial }{\partial z}\left(k_s\frac{\partial T_s}{\partial z}\right)={\left(\mathit{\rho c}\right)}_s\frac{\partial T_s}{\partial t}$$E1

Equation of the fluid (inside pipe) region. Cylindrical coordinates:

$$\frac{1}{r}\left[\frac{\partial }{\partial r}\left(r{k}_f\frac{\partial T_f}{\partial r}\right)\right]+\frac{1}{r^2}\frac{\partial }{\partial \phi }\left(k_f\frac{\partial T_f}{\partial \phi }\right)+\frac{\partial }{\partial z}\left(k_f\frac{\partial T_f}{\partial z}\right)={\left(\rho u_{z}c\right)}_f\frac{\partial T_f}{\partial z}+{\left(\mathit{\rho c}\right)}_f\frac{\partial T_f}{\partial t}$$E2

Equation of the fin. Cartesian coordinates:

$$\frac{\partial }{\partial x}\left(k_a\frac{\partial T_a}{\partial x}\right)+\frac{\partial }{\partial y}\left(k_a\frac{\partial T_a}{\partial y}\right)+\frac{\partial }{\partial z}\left(k_a\frac{\partial T_a}{\partial z}\right)={\left(\mathit{\rho c}\right)}_a\frac{\partial T_a}{\partial t}$$E3

This model can be considered as an extension of the classical conjugate-extended Graetz problem and has been widely used in the literature [20, 21, 22].

Time-dependent terms will be omitted in the current model due to fact that the study focuses on the stationary phenomenon, and conduction coefficients are considered equal in all directions as the fluid is considered isotropic and with invariable thermal properties. The rest of boundary conditions that define the problem are detailed in Table 1.

At z = 0, r < R+ ep	Tf = Ts = T1; ∂Ta/∂z = 0	(4)
At z = +∞	∂Tf/∂r = ∂Ts/∂r = 0; ∂Ta/∂z = 0	(5)
At r = R	Tf = Ts; (ks/kf)∂ Ts / ∂r = ∂Tf / ∂r	(6)
At r = R+ ep, φ ∈ (0, π]	Ts = Ta = T2	(7)
At r= R+ ep, φ ∈ 0	Ts = Ta	(7 bis)
At r = R+ ep, φ ∈ (π, 2π)	∂Ts/∂r = 0	(8)
At y = e/2	Ta = T2	(9)
At y = −e/2	∂Ta/∂y = 0	(10)
At x = a	∂Ta/∂x = 0	(11)

Table 1.

Analytical and dimensional boundary conditions.

3. Numerical model

Numerical solutions of coupled governing Eqs. (1)–(3) under boundary conditions (4)–(11) have been reached using Network Simulation Method (NSM) and the circuit solver NGSpice.

3.1 The network model

In NSM terminology, the equivalent, or analogous, electrical circuit of the considered process is called the network model. The construction of the network model begins with the discretisation of the mentioned system of differential equations in a finite volume elements mesh which transforms the governing partial differential equations of balance Eqs. (1)-(3) into the set of algebraic Eqs. (12)–(14).

Discretisation of the equation of the solid (duct) region. Cylindrical coordinates:

$$\begin{array}{l}\frac{1}{r}·\left[\frac{T_{j+∆r/2}-T_{j-∆r/2}}{∆r}\right]+\frac{1}{∆r}·\left[\frac{T_{j+∆r/2}-T_{j-∆r/2}}{∆r}\right]+\frac{1}{r^2\phi }\left[\left(T_{\phi +∆\phi /2}-T_{\phi -∆\phi /2}\right)/∆\phi \right]\\ +1/∆z·\left[\left(T_{i+∆z/2}-T_{i-∆z/2}\right)/∆z\right]\\ =0\end{array}$$E12

Discretisation of the equation of the fluid (duct) region. Cylindrical coordinates:

$$\begin{array}{l}1/r·\left[\left(T_{j+∆r/2}-T_{j-∆r/2}\right)/∆r\right]+1/∆r·\left[\left(T_{j+∆r/2}-T_{j-∆r/2}\right)/∆r\right]\\ +1/r^{2}1/\phi ·\left[\left(T_{\phi +∆\phi /2}-T_{\phi -∆\phi /2}\right)/∆\phi \right]+1/∆z·\left[\left(T_{i+∆z/2}-T_{i-∆z/2}\right)/∆z\right]\\ -{\left(\left(\rho ·c·u_z\right)/k\right)}_f\left[\left(T_{i+∆z/2}-T_{i-∆z/2}\right)/∆z\right]\\ =0\end{array}$$E13

Discretisation of the equation of the fin. Cartesian coordinates:

$$\begin{array}{l}1/∆x·\left[\left(T_{i+∆x/2}-T_{i-∆x/2}\right)/∆x\right]+1/∆y·\left[\left(T_{i+∆y/2}-T_{i-∆y/2}\right)/∆y\right]\\ +1/∆z·\left[\left(T_{i+∆z/2}-T_{i-∆z/2}\right)/∆z\right]\\ =0\end{array}$$E14

Based on the governing equations and boundary conditions, an electrical circuit for each equation has been created (Figure 3). Each term in these equations becomes an electric component. Consequently, each finite volume element (elementary cell) is composed of a set of electrical elements according to the thermo-electric analogy corresponding to the different terms of the previous equations. In this temperature is equivalent to voltage, and heat fluxes (∂T/∂x, ∂T/∂y, ∂T/∂z, ∂T/∂r and ∂T/∂φ) are equivalent to electric currents. More details of the fundamentals of the NSM can be found in Ref. [42].

Figure 3.

Electrical circuit scheme of the fluid basic cell.

From the point of view of the NSM, the terms contained in each of the above equations can be considered as currents according to the currents Kirchhoff’s law (because its summation over a node needs to be zero). As an example, the term $1/r·\left[\left(T_{j+∆r/2}-T_{j-∆r/2}\right)/∆r\right]$ introduces heat conduction in cylindrical co-ordinates both on the fluid and on the pipe equation. The step-by-step description of the elementary cell building (Figure 3) has been detailed in the Appendix of this chapter.

Figure 4.

Numerical model discretisation mesh.

Eq. (13) describes the three-dimensional flow behaviour which is dominated by the axial velocity u_z, assumed as a velocity function of order 10, profile equation of which is done by Eq. (15):

$$u_z=\begin{array}{l}1.06·{10}^{-2}+\left(\frac{r}{R}\right)·4.59·{10}^{-1}+3.75·{\left(\frac{r}{R}\right)}^2+19.4·{\left(\frac{r}{R}\right)}^3+15.6·{\left(\frac{r}{R}\right)}^4\\ -92.1·{\left(\frac{r}{R}\right)}^5+40.6·{\left(\frac{r}{R}\right)}^6+843·{\left(\frac{r}{R}\right)}^7+577·{\left(\frac{r}{R}\right)}^8-{\left(\frac{r}{R}\right)}^9·3.26·{10}^3\\ -{\left(\frac{r}{R}\right)}^{10}·4.93·{10}^3\end{array}$$E15

In total, the entire system has been discretised using a three-dimensional mesh of identical cells in a manner that the square-sectioned straight fin is divided into 200 cells in z direction, 10 cells in the x direction and 1 cell in the y direction (Figure 4). The pipe is also symmetrically divided into 200 cells in the z direction, 4 cells in φ direction and 7 cells in the r direction where 5 belongs to the fluid and 2 to the pipe thickness (Figure 4). This mesh has been considered to be sufficient for the purpose of this work.

Figure 5.

Longitudinal pipe section planes screened in simulations.

Figure 5 shows the five main planes of the pipeline (named similarly to the cardinal points) that will be referred to in Section 4 to identify the different points of the pipeline’s cross section. Although we have pointed out the asymmetry of the physical model, it should be noted that the pipe-fin assembly is radially asymmetric, but it presents symmetry with respect to the N-S plane.

Once the elementary circuit has been built, the boundary conditions must be implemented. In this case, some of the most relevant are voltage sources than fit the constant external temperature T₂ in the upper side of the assembly and infinite (very high) resistance in its lower side. Finally, elementary cells and devices corresponding to the boundary conditions are assembled building an electric circuit (network model), which can be solved using an appropriate software, such as NGSpice.

A complete description of the circuit construction as well as computing details of the simulation can be found in the Appendix and reference [42]. In order to know the numerical value of the different electrical components forming out network (resistors, voltage and current sources, etc.), it is necessary to build the circuit and to obtain realistic values of parameters and boundary conditions corresponding to the case study. In our case, a solar thermal collector. The parameters are listed on Table 2. As stated above, solar thermal devices are one of the practical applications of this model, with obvious similarities in their geometry.

Table 2.

Simulated model magnitudes.

In this case, when Reynolds number is high (41989) and Rayleigh number remains low (3.92·10³), the fluid flows can be considered in a laminar forced-convection regime, and in consequence, the buoyancy effects are negligible [44]. The minimum pipe length in which the thermal phenomena are developed came from the use of the Nusselt number approximations [47].

About the length of the pipe, it has been used the concept of characteristic length, l*, defined as the length needed for the fluid to fully develop the thermal process [26], in this case to achieve the external temperature T₂. Therefore, the simulated pipe should be long enough (L > l*), to ensure that the thermal phenomena due to the sudden change in temperature at z = 0 has ended. This hidden parameter must be previously determined. In this case, the preliminary simulation with a 2D model guided the determination of the characteristic length (afterwards checked). So, we found that the corresponding characteristic length in this case would be 2 m if the whole exterior pipe surface was kept isothermally surrounded. This is not exactly the case, because in the 3D model presented in this chapter, only upper half surface is at constant temperature and the rest is insulated; consequently, 4 m length duct has been simulated.

3.2 Model validation

The 3D numerical model built following these rules and using the parameters mentioned in the previous section was validated by comparing its results to literature and experimental data [42].

The case modelled in the present chapter is a complete innovation; therefore, no references for comparison can be found in the literature. Furthermore, no references of the 3D configuration of the conjugate extended Graetz problem could be found in the literature previously to the work of reference [42] of the same authors of this chapter. In comparison to the finned tube model of the present chapter, the above cited reference modelled a bare tube (without fins) with a parabolic velocity profile, both being models analogous in the rest of characteristics.

Figure 6 shows the comparison of the fluid temperature field at 0.5 m of the entrance in both models. Main differences are a consequence of the different fluid velocity and the presence of fins. Figure 6a depicts the isothermal lines of the finned tube, where a slight tendency towards horizontality can be appreciated in the surroundings of the tube-fin junction compared with that of the Figure 6b of the bare tube. Nevertheless, the bulk fluid is little affected by the presence of fins. This was to be expected as in the finned tube, the temperature of the fins is the same as that of the upper half-pipe in both models. Consequently, it can be said that the fluid temperature field of the finned tube model agrees with that of the bare tube, and it can be considered that the validation procedure carried out in Ref. [42] applies to the present model.

Figure 6.

Comparison of the cross-section temperature maps at 0.5 m in the finned tube (a) and in a bare tube (b).

In respect to the code, the results of the 3D model “acting as” a 2D model (radially symmetric boundary conditions) were compared with the 2D published results [21, 22], obtaining relative errors along the whole tube length below 2% and mean of 0.98%, the error standard deviation [42] being 0.56%. Analysis of the typified residuals of the 2D and 3D simulation results, error variance and a regression analysis were carried out. The coefficient of skewness (nearly 0) and kurtosis showed that the error data set was normally distributed. Finally, the relative error of the 3D model acting as a 2D model, and the 2D model was found to be 0.98 ± 7.74E-2% at a 95% confidence level.

On the other hand, the external temperatures (tube) of the cross section of the 3D simulation at different lengths were compared with the experimental measures yielded in a solar thermal experimental rig [42]. An analogous error study was conducted between numerical and experimental temperature data, concluding that the 3D numerical results were sufficiently close to those measured experimentally, relative errors being of 3.40 ± 0.601%, at 95% of confidence. More details of experimental rig, measured data and error study can be found in [42].

Consequently, those results confirmed the accuracy of the bare tube 3D model, and this conclusion can be extended to the finned tube model of this chapter, substantially equivalent to that of reference [42], especially as regards the fluid, and which can be considered as a continuation of the bare tube model.

4. Results and discussion

Typical fluid flow in laminar regime shows a parabolic profile. Nevertheless, there are situations in which velocity has a more complex profile [43]. Because of that a polynomial function of grade 10 (Figure 7) has been used in this work. Eq. (15) is the formula corresponding to this curve, and it is the value for u_z in Eqs. (2) and (13). Eq. (15) has been yielded from experimental results, and it corresponds to the case of the flow in tilted solar collectors [43].

Figure 7.

Axial velocity profile following [43].

4.1 Axial temperature profile and characteristic length of the process

Figure 8 shows the profile of the fluid temperature along the pipe in the different planes of the pipe (see Figure 4). In this case, the temperature rises faster than in the case of the parabolic profile of the velocity, reaching the external temperature T₂ even in the lower sections of the tube, i.e. flows with this speed profile mix much faster than with the laminar parabolic profile, yielding a higher convective coefficient.

Figure 8.

Temperature profiles for u_z given by Eq. (15).

If characteristic length, l*, is defined as the one where the temperature in the middle of the pipe (red line) reaches 90% of the total temperature jump, in this case 65.45 + 0.9⋅(100–65.45) = 96.55°C, from Figure 7 gives l* = 3.6 m. Beyond this point, the temperature of the fluid does not increase substantially, which indicates that the rest of the pipe is not efficient for heat transfer purposes, so the length could be reduced.

4.2 Fluid temperature maps in a cross section

Figure 9 shows the temperature field of the cross section of the fluid (perpendicular to the axis) and of some vertical and horizontal sections of it at 13 cm from the inlet. As expected, the temperature curves show symmetry about the vertical axis.

Figure 9.

Fluid temperature map of the cross section (a) and temperature profiles for different sections (b–e) at 13 cm from the entrance of the pipe.

Figure 9a represents the thermal map of the cross section. It can be observed that the fluid initially (z = 0) at T₁ (65.45°C) acquires at 13 cm from the entrance temperatures close to T₂ (100°C) in the layers near the top of the pipe (dark brown layers between the N (φ = π/2) and NE (φ = π/4) planes and nearby areas). At this point of the pipe, the upper half of the fluid shows almost parallel isotherms following the circular curve of the pipe, whose values decrease from plane N, (at the top) to the horizontal plane (plane E, φ = 0) and even up to almost plane SE (φ = –π/4). Sections A-A’ and D-D′ (Figure 8b and e, respectively) clearly show this situation, presenting a minimum at the N-S axis as a consequence of the shape of the isotherms. Compared with the case of velocity parabolic profile [42], it can be said that in the case of polynomial profile, the temperature increases much faster than in the former, reaching the external temperature T₂ in the lower sections of the tube relatively close to the entrance, i.e. flows with this speed profile mix much faster than with the laminar parabolic profile, yielding a higher convective coefficient. A consequence of this is the shorter characteristic length shown by the polynomial velocity process compared with that of parabolic velocity profile.

In contrast, near the south (S) and southeast (SE) planes, the isothermal curves take a U-shape, except near the N-S axis, where a loop is formed at 68°C approximately in the centre of the bottom half of the tube. This is due to the fact that the fluid is being heated from the top half of the tube (whose temperature has been imposed at a uniform value of T₂), while the lower half of the tube, which is externally insulated, is warmed by heat conduction from the solid pipe due to the fact that k_f > > k_s. As a result, in the S plane (φ = − π/2), the fluid near the tube is hotter than in the middle of the tube. This fact can also be detected in Section B-B′, which corresponds to the N-S plane (Figure 9c), where the minimum (66.98°C) is not found in the lower part of the section; you can see that the temperature profile is roughly flat along the bottom quarter of the section. This suggests that the thermal process is more dependent on pipe-wall conduction effects than on the velocity profile, even with such a complex and irregular profile, due to the low conductivity of the fluid and the laminar flux. It is relevant to note that the classic Graetz Problem does not consider the pipe through which the fluid flows, thus making it impossible to detect this temperature distortion.

Another issue is the influence of the fin on fluid temperatures. The fin-tube junction takes place in plane E, which corresponds to section A-A’ (Figure 9b). A small distortion of the isotherms can be observed at this point compared with the case of a bare duct [42]. In this case, this distortion is not very important because the external temperature T₂ is uniform, both in the pipe and in the fin, and a very thin fin with little heat conduction capacity has been modelled. Otherwise, it could be very relevant, as in the case of the heat boundary condition.

Although not studied in this work, consequences from the non-uniform temperature field, which affects most of the thermo-fluid properties (density, viscosity, etc.) among others, can be drawn. This behaviour could impact on the operation of the pipeline and the related equipment.

5. The HEATT© platform

The methodology explained above has been used to develop a 3D model to simulate a laminar flow under conditions of forced convection subject to asymmetric boundary conditions, such as those found in the tube grid of a low-temperature solar thermal collector, in which the upper half of each tube receives energy from the sun, while the bottom half of the tube remains embedded in a layer of insulation [42]. The simulated results could be compared with the previous bibliography and with the experimental results obtained from the Solar Laboratory of the University of Murcia under real operating conditions, obtaining great differences between the results obtained using 2D models versus 3D models. The model roughly coincides with the one presented in this chapter (except for the fact of considering tube without fins).

This model required the creation of a three-dimensional mesh with nearly 17,000 cells, in which an electrical equivalent circuit of approximately 140,000 elements was implemented. The model includes resistors, voltage and current sources and capacitors.

Once the problem has been solved for a sufficiently variety of real cases, the tool has been created for solving the problem using a web as a service, providing free service to any professional or researcher anywhere in the world without the need to acquire expensive software or to instal any application that becomes obsolete with the evolution of the Operating Systems.

The resulting web application has been named HEATT ©, acronym for Pipeline Thermal Analysis and Assessment Tool (in Spanish: Herramienta de Evaluación y Análisis Térmico de Tuberías).

This simple and friendly platform is currently being launched in its version 1.0 as a Proof of Concept, to be released “as a service” for general use and for researchers and professionals from all over the world to send their opinions and improvement proposals to gradually make it growing up.

6. Conclusions

A numerical physical–mathematical model is presented for a laminar forced-convection fluid flow within a finned round duct subjected to constant temperature on its upper side and insulated on its lower side. The governing heat transfer and fluid flow 3D partial differential equations, combining cylindrical and Cartesian coordinates, are solved for steady-state conditions. The numerical model is run using the accurate contrasted Network Simulation Method, a low-workload computational method.

In the present chapter, the Graetz Problem is extended to incorporate axial fluid conduction, 3D coordinates, wall thickness with attached longitudinal flat-fins, radially asymmetrical boundary conditions and highly non-linear velocity profiles. Realistic conditions corresponding to the flow in flat plate solar collectors have been used in the simulations. A high non-linearity velocity profile has been evaluated, confirming that solid conductivity and thickness effects are not negligible for the studied thermal phenomena.

Temperature evolution across the fluid is analysed in detail. Different temperature values were found for different angles within every plane, due to the radial asymmetry of the geometry. Temperature fluid maps were obtained at different distances from the entrance of the pipe. Isotherms show parallel-like shapes on the top half of the tube, which become distorted in the lower half, where some loops appear due to the conduction effects of the studied pipe wall thickness, illustrating the non-uniformity of the temperatures within the fluid.

In addition, 3D simulation reveals that, in cases of asymmetry, the thermal phenomena require much more length to completely develop the flow than the length yielded by the 2D radially symmetric model. This is relevant because the 3D simulation reveals that the pipe needs up to six times more length than that predicted by 2D model.

When the characteristic length of the problem was considered (a virtual dimension equivalent to the distance at which the thermal phenomena are fully developed), the 2D approach was found to be no longer valid. That occurs in most of the studied real cases, when the pipe is subjected to asymmetrical boundary conditions.

All those changes enhance the solution of the Graetz problem and bring it nearer to real pipe conditions. From this, a basis for future works, including heat boundary conditions, different union thermal resistances or a variety of complex velocity functions, among many other possibilities, may be mentioned in order to encourage future optimisation studies. Meanwhile, the findings of the present study have applications in solar energy collectors, thermal heat dissipators, oillines and heat exchangers among many other facilities.

In addition, an open computing platform called HEATT©, based on this model, is now being built. The platform is expected to be freely available to the public before the end of 2022.

Acknowledgments

This project has been possible thanks to the funds obtained in public concurrence in the 2021 call for Proofs of Concept from the Seneca Foundation – Agency for Science and Technology of the Region of Murcia (Spain). Ph.D, Ramallo-González would like to thanks the European Commission for their funding of project PHOENIX grant number 893079.

Conflict of interest

The authors declare no conflict of interest.

A.

The procedure that needs to be followed to reproduce the simulation results of this chapter is described below.

A.1 Discretisation of the spatial domain

In this case, we can account three different regions or sub-domains, as stated in the physical–mathematical model, every one described by its corresponding equation of thermal behaviour: solid (pipe) region, Eq. (1), fluid (inside pipe) region, Eq. (2) and the fin, Eq. (3). Each region is divided into an adequate number of bi-dimensional cells (the third dimension is the axial coordinate, which is common for all cells). In this problem, the number of discrete dimensions are:

1. Pipe: Two divisions in the radial coordinate and four divisions in the azimuthal coordinate (Figure 5). Taking into account the symmetry of the problem in respect to the N-S plane, only half section must be simulated.

2. Fluid: Five divisions in the radial coordinate and, in accordance to the pipe, four divisions in the azimuthal coordinate.

3. Fin: Ten divisions of the x-dimension and one (no discretisation) in the y-dimension (Figure 4).

The third dimension, i.e. z-dimension, has been divided into 200 parts (Figure 4) and is common for all the bidimensional cells previously accounted for.

In summary, an overall number of 38 divisions (8 for the pipe, 20 for the fluid and 10 for the fin) of the cross section have been done, which gives a total of 38x200 = 7600 cells for the whole domain. In this case, the results have shown that the discretisation is adequate.

A.2 Building the numerical model

This step is divided into different minor steps, each one has to be applied to the corresponding region. Due to the fact that the most complex region is that of the fluid, and in order not to be repetitive, the process of construction of its elementary cell. The starting point is the set of discretised Eqs. (12)–(14) of the problem.

A.3 Building the elementary cell

The Network Simulation Method is based in the well-known thermal-electrical analogy [35]. In accordance to this, temperature is equivalent to voltage, and ΔT/Δz, ΔT/Δr and ΔT/(rΔφ), which are related with heat fluxes, are equivalent to electric currents. Each of the Eqs. (12)-(14) represents the energy balance in the corresponding cell. These cells, at steady state must have a zero summation of currents. Consequently, each term of the differential equation must be converted into a current of the elementary circuit, which as a whole fulfils Kirchhoff’s Laws.

In order to arrange the equations in a more convenient way, Eqs. (12) and (13) are multiplied by Δz⋅Δr⋅Δφ; Eq. (14) is multiplied by Δx⋅Δy⋅Δz. Note that the elementary cell (Figure 3) is in fact divided in two parts, being the node at the centre of the cell.

Tables 3, 4 and 5 contain the formulae of the different devices that integrate the elementary cells of solid pipe, fluid and fin, respectively. As the most complex, the fluid region elementary circuit is explained. In this case, five terms are found:

Term	Discretisation of equations	Thermal-electrical analogy	Type of device
1	$\frac{∆\phi \bullet ∆z}{r}\left[T_{j+\frac{∆r}{2}}-T_{j-\frac{∆r}{2}}\right]$	$G_{r,s}=∆T\bullet \frac{∆z·∆\phi }{r}$	Voltage-dependent current source
2	$∆\phi \bullet ∆z·\left[\frac{T_{j+\frac{∆r}{2}}-T_j}{∆r_s/2}-\frac{T_j-T_{j-\frac{∆r}{2}}}{∆r_s/2}\right]$	$R_{r,s}=\frac{∆r_s}{2∆z·∆\phi }$	Resistor
3	$\frac{∆r\bullet ∆z}{r^2}·\left[\frac{T_{\phi +\frac{∆\phi }{2}}-T_{\phi }}{∆\phi /2}-\frac{T_{\phi }-T_{\phi -\frac{∆\phi }{2}}}{∆\phi /2}\right]$	$R_{\phi ,s}=\frac{r^2·∆\phi }{2∆r_s·∆z}$	Resistor
4	$∆r\bullet ∆\phi \left[\frac{T_{i+\frac{∆z}{2}}-T_i}{∆z/2}-\frac{T_i-T_{i-\frac{∆z}{2}}}{∆z/2}\right]$	$R_{z,s}=\frac{∆z}{2∆r_s·∆\phi }$	Resistor

Table 3.

Formulae and analogous electrical devices yielded from solid Eq. (12).

Term	Discretisation of equations	Thermal-electrical analogy	Type of device
1	$\frac{∆\phi \bullet ∆z}{r}\left[T_{j+\frac{∆r}{2}}-T_{j-\frac{∆r}{2}}\right]$	$G_{r,f}=∆T\bullet \frac{∆z·∆\phi }{r}$	Voltage-dependent current source
2	$∆\phi \bullet ∆z\left[\frac{T_{j+\frac{∆r}{2}}-T_j}{∆r_f/2}-\frac{T_j-T_{j-\frac{∆r}{2}}}{∆r_f/2}\right]$	$R_{r,f}=\frac{∆r_f}{2∆z·∆\phi }$	Resistor
3	$\frac{∆r\bullet ∆z}{r^2}·\left[\frac{T_{\phi +\frac{∆\phi }{2}}-T_{\phi }}{∆\phi /2}-\frac{T_{\phi }-T_{\phi -\frac{∆\phi }{2}}}{∆\phi /2}\right]$	$R_{\phi ,f}=\frac{r^2·∆\phi }{2∆r_f·∆z}$	Resistor
4	$∆r\bullet ∆\phi \left[\frac{T_{i+\frac{∆z}{2}}-T_i}{∆z/2}-\frac{T_i-T_{i-\frac{∆z}{2}}}{∆z/2}\right]$	$R_{z,f}=\frac{∆z}{2∆r_f·∆\phi }$	Resistor
5	${\left(\frac{-\rho ·c·u}{k}\right)}_f·\left[T_{i+\frac{∆z}{2}}-T_{i-\frac{∆z}{2}}\right]$	$G_{z,f}=-∆T·∆r_f·∆\phi \bullet {\left(\frac{\rho ·c·u}{k}\right)}_f$	Voltage-dependent current source

Table 4.

Formulae and analogous electrical devices yielded from fluid Eq. (13).

1. J₁ is a non-lineal term which comes from the increasing cell area with the radius. It must be implemented by a voltage-dependent current source, G_r,f (Figure 3), value of which is expressed by Eq. (15), as it can be seen in Table 4.

   $$J_1=\frac{∆z·∆\phi }{r}\left[\left(T_{j+\frac{∆r}{2}}-T_{j-\frac{∆r}{2}}\right)\right]$$E16

2. The second term of the Eq. (13) is a current which corresponds to the heat flux in r-direction, T_j + Δr/2 and T_j-Δr/2 being the temperatures at both sides of the elementary j-cell in this direction, which are made equivalent to the external voltages of the elementary circuit, Eq. (16).

   $$J_2=∆z·∆\phi \left[\frac{\left(T_{j+\frac{∆r}{2}}-T_{j-\frac{∆r}{2}}\right)}{\left(∆r/2\right)}\right]$$E17

   Consequently, the thermal resistance is represented by two electric resistances of value $R_{r,f}=\frac{∆r}{2∆z·∆\phi }$, as can be seen in Table 4. Similarly, the third and fourth terms become J₃ and J₄ currents corresponding to heat flux in the φ and z directions. The value of the thermal resistances R_φ,f and R_z,f yielded from these terms can be found in Table 4.

3. Finally, the fifth term represents the axial heat conduction due to the velocity field and has been implemented by a voltage-dependent current source, G_z,f, whose value is given by Eq. (17):

$$J_5=∆z·∆\phi {\left(\left(\rho ·c·u_z\right)/k\right)}_f·\left(T_{i+∆z/2}-T_{i-∆z/2}\right)$$E18

The location and connections of the different electric devices in the fluid elementary cell can be seen in Figure 3; plus or minus symbols in the different resistors distinguish both sides of the cell in respect to the central node.

This procedure must be followed for the different regions of the assembly, i.e. the fluid, the pipe and the fin. The pipe equation, Eq. (12) has only four terms, all of them coincide with those of the fluid except for the velocity term, and the fin, Eq. (14), has only three terms, one for each Cartesian direction. Tables 3–5 summarise the devices used in the elementary circuits.

Term	Discretisation of equations	Thermal-electrical analogy	Type of device
1	$∆y\bullet ∆z\left[\frac{T_{k+\frac{∆x}{2}}-T_k}{∆x/2}-\frac{T_k-T_{k-\frac{∆x}{2}}}{∆x/2}\right]$	$R_{x,a}=\frac{∆x}{2∆y·∆z}$	Resistor
2	$∆x\bullet ∆z\left[\frac{T_{l+\frac{∆y}{2}}-T_l}{∆y/2}-\frac{T_l-T_{l-\frac{∆y}{2}}}{∆y/2}\right]$	$R_{y,a}=\frac{∆y}{2∆x·∆z}$	Resistor
3	$∆x\bullet ∆y\frac{\left[T_{i+\frac{∆z}{2}}-T_i-T_i-T_{i-\frac{∆z}{2}}\right]}{\left(∆z/2\right)}$	$R_{z,a}=\frac{∆z}{2∆x·∆y}$	Resistor

Table 5.

Formulae and analogous electrical devices yielded from fin Eq. (14).

A.4 Building the network model or whole equivalent circuit

Once the elementary circuits of the different parts of the system have been built, the equivalent circuit of each region is assembled by adding as many cells as needed according to the discretisation carried out; then the circuits of the regions are joined according to the topology.

Finally, the problem boundary conditions, Eqs. (4)–(11), must be implemented. In this case, the most relevant are voltage sources for Eqs. (4), (7) and (9), very large resistance resistors are used to emulate the adiabatic conditions, Eqs. (4), (5), (8), (10) and (11) or just electrical continuity, Eqs. (6) and (7bis).

A.5 Writing the code

One of the advantages of NSM is the use of well-known and reliable software for the analysis of electronic and electrical circuits (NGSpice, PSpice or others), which are also very easy to program. Naturally, for the programming of the elementary (and complete) circuit, we direct the reader to the corresponding websites [40] and literature [41], where the documentation can be found.

B.

The whole electric circuit (network model) is solved using appropriate software, such as NGSpice, PSpice or others.

B.1 Obtaining results

The results obtained must be interpreted taking into account that in the thermoelectric analogy, voltages and currents are analogous to temperature and heat flow, respectively. Both temperature profile along the pipe-fin and temperature maps of the cross section are plotted using suitable software from the voltages at the appropriate points along the pipeline or the selected cross section using the software.

Formulae and analogous electrical devices yielded from Eqs. (12)–(14).

Nomenclature

a

absorber width (m)

c

specific heat (J·kg⁻¹·K⁻¹)

e

absorber thickness (m)

e_p

pipe thickness (m)

G

voltage-dependent current source (A)

k

thermal conductivity (W·m⁻¹·K⁻¹)

J

electric current (A)

L

pipe/absorber length (m)

l^*

characteristic length (m)

R

inner radius of the pipe (m)

R

Resistor (Ohm)

r

radial coordinate

T

Temperature (°C)

t

Time (s)

T₁

temperature in the entrance of the fluid (°C)

T₂

temperature in the exit of the fluid (°C)

u

Velocity (m·s⁻¹)

z

axial coordinate (m)

Superscripts

′

dimensionless variable

Subscripts

a

associated to absorber (fin)

ext

exterior

f

associated to fluid

i j, k, φ, l

discretisation counters

in

inner

m

medium/average

p

solid–liquid interface

r

associated to radius

s

associated to solid

z

associated to z axis

φ

assoc. to azimuthal coordinate

Greek characters

φ

azimuthal coordinate (rad)

μ

dynamic viscosity (kg·m⁻¹·s⁻¹)

ρ

density (kg·m⁻³)