Cooling of a channeled lava flow with non-Newtonian rheology : crust formation and surface radiance

We present here the results from dynamical and thermal models that describe a channeled lava flow as it cools by radiation. In particular, the effects of power-law rheology and of the presence of bends in the flow are considered, as well as the formation of surface crust and lava tubes. On the basis of the thermal models, we analyze the assumptions implicit in the currently used formulae for evaluation of lava flow rates from satellite thermal imagery. Assuming a steady flow down an inclined rectangular channel, we solve numerically the equation of motion by the finite-volume method and a classical iterative solution. Our results show that the use of power-law rheology results in relevant differences in the average velocity and volume flow rate with respect to Newtonian rheology. Crust formation is strongly influenced by power-law rheology; in particular, the growth rate and the velocity profile inside the channel are strongly modified. In addition, channel curvature affects the flow dynamics and surface morphology. The size and shape of surface solid plates are controlled by competition between the shear stress and the crust yield strength: the degree of crust cover of the channel is studied as a function of the curvature. Simple formulae are currently used to relate the lava flow rate to the energy radiated by the lava flow as inferred from satellite thermal imagery. Such formulae are based on a specific model, and consequently, their validity is subject to the model assumptions. An analysis of these assumptions reveals that the current use of such formulae is not consistent


Introduction
Laboratory studies on basaltic melts show that lava rheology can show nonNewtonian behavior under certain conditions, which include vesicularity [Spera et al. 1988, Stein and Spera 1992, Badgassarov and Pinkerton 2004], crystal concentration [Pinkerton and Stevenson 1992, Smith 2000, Sonder et al. 2006, Champallier et al. 2008], and temperature and shear rates [Shaw et al. 1968].The literature indicates that magmas have nonNewtonian, pseudoplastic behavior, with the exception of Smith [2000], who attributes a dilatant rheology to lava at high crystal concentrations.Pseudoplasticity and dilatancy can be described by power-law rheology.
When solving the problem of gravity-driven lava flow, power-law rheology introduces nonlinearity into the diffusion term of the momentum equation, and an analytical solution of the governing differential equations is not possible.This has given rise to various approximate solutions to the problem.The fully developed laminar flow of powerlaw fluids has been studied numerically using the finite element method [Syrjala 1995] and the finite volume method [Capobianchi 2008] for a pressure driven flow in a horizontal rectangular duct.
The equation of motion for gravity-driven lava flow with a power-law rheology was solved by Filippucci et al. [2010] using the finite volume method for the cases in which the rheology is temperature independent.Several models describing lava tube formation have been proposed [Greeley 1971, Peterson and Swanson 1974, Rowland and Walker 1990, Hon et al. 1994, Dragoni et al. 1995, Valerio et al. 2008, Valerio et al. 2010].In the present study, the cooling of a channel of lava with power-law rheology is considered.A two-dimensional thermal model with heat flux assigned at the upper surface is introduced to describe lava cooling due to radiation and convection into the atmosphere.The lava crust is considered as a plastic body, and its rheology is described through the introduction of a yield strength as a function of temperature.
Bends in lava flow are often observed, and these can strongly affect the flow dynamics and formation of lava tubes [Greely 1971, Peterson et al. 1994].In the present study, we consider the effects of the curvature of a channel on the lava velocity and shear stress, and on the formation of the solid crust at the flow surface.To allow an analytical solution of the Navier-Stokes equation in the presence of bends, the lava is regarded as a Newtonian, homogeneous, isotropic and incompressible fluid [Valerio et al. 2011].Heat radiation and convection into the atmosphere are considered as the main cooling mechanisms.The model analyses the effects of curvature on the development of surface solid plates.
The availability of high-resolution thermal imagery of active lava flows has stimulated the use of radiance maps for evaluation of lava effusion rates [Mouginis-Mark et al. 1994, Harris et al. 1995, Calvari et al. 2010].This has been made possible by the use of simple formulae that relate the lava flow rate to the energy radiated per unit time from the planimetric surface of the flow [Pieri and Baloga 1986, Harris et al. 1997, Harris et al. 2005a].Such formulae are based on a specific flow model, and consequently, the results depend on the flow-model assumptions.An analysis of these assumptions was performed by Dragoni and Tallarico [2009], which revealed that the current use of these formulae is not consistent with the model.

Crust formation in a lava channel with power-law rheology
The constitutive equation for a power-law fluid is: (1) where v ij is the stress tensor, ė ij is the strain rate tensor, k is the fluid consistency, n is the power-law exponent, which is a measure of the nonlinearity, and: (2 where I 2 is the second invariant of the strain rate tensor.The apparent viscosity of the fluid is given by: (3) If n < 1, the fluid is pseudoplastic and it thins (i.e. it deforms more rapidly) with an increase in stress.If n = 1, the fluid is Newtonian.If n > 1, the fluid is dilatant and it thickens with an increase in stress.
An illustration of the model with the coordinate system is shown in Figure 1.The equation of motion for a gravitydriven flow down an inclined rectangular channel is: (4) where v x is the x component of velocity, g is the gravity acceleration, and a is the slope angle.The apparent viscosity is: (5)  2ke e The boundary conditions are the nonslip at the walls, and free surface at the top of the flow: The dynamic problem has been solved using the finite volume method with an iterative solver [Patankar 1980].The numerical tests carried out by Filippucci et al. [2010] indicated that the numerical solution, with respect to the analytical solution, that is available for the Newtonian rheology produces an error of about 0.03% with a grid of 51 × 51 control volumes (CV).The error due to the mesh refinement tested with the power-law rheology, for which the analytical solution is not available, is about 0.03% with a grid of 51 × 51 CVs, so the grid of 51 × 51 CVs was adopted for all of the following computations, as a good compromise between precision and time for calculation.
The effects of the slope a on the average velocity and the flow rate can be remarkably different for Newtonian and nonNewtonian fluids.Figure 2a shows the effects of a on the average velocity, with varying n, and fixed k, channel width a, and thickness h.Differences between the linear and the nonlinear rheology are enhanced by increases in a.The effects of the slope have been evaluated also for a constant volume flow rate q (Figure 2b).The effects of nonlinearity on v x vary with a and k.Finally, we evaluated the dependence of the volume flow rate q on h for fixed values of a and a (Figure 2c).
We simulate two possible scenarios of nonlinearity (pseudoplastic, with n = 0.5, and dilatant, with n = 1.5) and compare these with the Newtonian case (n = 1) for the Laghetto lava channel (Figure 3a, b) and for another active lava channel, both of which were formed during the 2001 eruption of Mount Etna (Figure 4).Our results show that the COOLING OF CHANNELED LAVA FLOW  ( , velocity profiles and average velocities can vary significantly, depending on the assumption of linear or nonlinear rheology. We have also considered the May 1997 flow from Pu'u 'O'O (Kilauea, Hawai'i).The apparent viscosity varies from 9 to 879 Pa s, and decreases with an increase in the shear rate, so the system is pseudoplastic [Weed et al. 1986].The channel geometry was described by Harris and Rowland [2001].We simulated the flow with pseudoplastic rheology, setting n = 0.9, 0.7 and 0.5.Given the large variations in channel width and the possible variations in channel thickness, we considered two possible geometries for the channel section.In this way, we can observe how the average velocities are modified by a change in the aspect ratio b = h/a of the channel.As a result, with an increase in b the average velocity increases, and this effect is greater as n decreases (Figure 5a).Moreover, the volume flow rate q, which we hypothesized as constant with respect to b, increases very rapidly as n decreases (Figure 5b).
Finally, we considered the 1984 Mauna Loa lava channels.We used power law rheology to fit the flow velocity data measured from videotapes by Sakimoto and Gregg [2001].As a result, we found n in the dilatant range (Figure 6a, b).
Cooling is modeled for a lava flow exiting a vent with an effusion temperature T 0 = 1,273 K.We assume that the flow is laterally isothermal and that the cooling process occurs from the upper surface with a fixed value of the heat flux q 0 .This describes the average heat loss from the surface that is covered by a solid crust and from the crust-free shear regions [Valerio et al. 2008].The energy equation governing the heat transfer is: where T is the temperature and | is the thermal diffusivity: (9) TALLARICO ET AL.
where c p is the specific heat, and K is the thermal conductivity.
The boundary conditions are the adiabatic conditions at the walls and the constant heat flux at the upper surface.Initially, the fluid has a uniform temperature T 0 .At the time t > 0, a constant heat flux q 0 is set and the thermal boundary conditions become: (10) (11) In addition to the spatial integration, the discretization of the energy equation involves integration of Equation ( 8) over the time interval from t to t + Dt.This integration is operated by using a fully implicit scheme, which ensures simplicity and physical consistency [Patankar 1980].Details of the numerical test of this solution are given below.
We substituted the time dependence of the solution in T with the space dependence by the relation: (12) where: (13) We first considered the conductive cooling of a steady lava flow, and the consequent surface crust formation.Indeed, many basaltic lava channels have a surface crust that appears dark due to its lower radiative temperature, and that is carried along by the flow.We evaluated the effects of the nonlinearity of the lava rheology on crust formation, by modeling the plastic behavior of the crust using a temperature-dependent yield strength x, given by: where T s is the solidus temperature, x 0 is the maximum yield strength, and b is a constant.
We call d the width of the flow surface covered by crust.For temperatures below the solidus, crust formation and growth are possible in the region |y|< d/2, where the horizontal shear stress v xy < x.Shearing movements occur and crust formation is inhibited where v xy > x.As the surface temperature reaches the solidus, x starts to increase up to the threshold.The distance from the vent where this happens strongly depends on the assumption on the rheology (Figure 7a).Simultaneously, the crust width d/a also grows, and once it attains its maximum width, it stops widening.In particular, for the dilatant fluid, the crust stops widening earlier in space with respect to the Newtonian and the pseudoplastic fluid, which will cover a much longer distance before the crust reaches its maximum width (Figure 7b).
Models of formation of lava tubes as a consequence of the crust welding to the channel levees have been proposed that consider lava as a Bingham liquid [Dragoni et al. 1995] or a Newtonian liquid [Cashman et al. 2006, Valerio et al. 2008].We also explored how variations in the channel width, ground slope and volume flow rate can affect the formation of a lava tube, assuming a power-law rheology.As a result, we find that for topographic and morphological variations of the lava channel, the nonlinearity of the rheology has only a slight influence on the formation of a lava tube, and the results for the Newtonian fluid appear to also be valid for power-law fluids.

Numerical test on the solution of the heat equation
The solution of the heat equation was tested to study the accuracy, by comparing it with the analytical solution with constant surface flux q 0 [Turcotte and Shubert 1982]: ( , We evaluated the difference DT u between the upper temperatures of the two solutions and the difference DT between the average temperatures [Syrjala 1996, Capobianchi andWagner 2010], as given by: (16) As we expected, both and DT u always decrease as the number N CV of the control volumes increases, i.e. by thickening of the computational mesh, which means that the numerical temperature gradually approaches the analytical value (Figures 8, 9, respectively).In this test, is always about two orders of magnitude less than DT u , and this difference is a consequence of the cooling of the lava channel, which at the time considered, involves only the upper layers of the flow; the lower part remains at a constant temperature.It can be seen that is very low also for a coarse grid, as it reaches 0.01% for a grid 61 × 61 CVs (Figure 8), while DT u needs a grid of 201 × 201 CVs to drop under 0.05% (Figure 9).In Figures 8 and 9, the errors on and T u due to the mesh refinement are also plotted.It can be seen that due to mesh refinement, is a decreasing function of N CV , starting from the grid of 51 × 51 CVs.For Figure 10, we evaluated the trend of with N CV at different times.This test indicates that decreases as the time increases, which indicates that the discretization error does not propagate with time and the solution is convergent.

Crust formation in the presence of a bend in the channel
Bends in lava flows arise from local variations in the topography, and they are commonly observed in volcanic fields.The curvature of a channel affects the flow dynamics and morphology [Greeley 1971].
An analytical model was proposed by Valerio et al. [2011] to explain the effects of curvature on the velocity and viscous stress.In a bending channel, a cylindrical coordinate system is assumed (Figure 11).A flow segment with a constant curvature is limited by arcs of concentric circumferences, with their centers in the origin of the coordinate system.We studied the equation of motion for a viscous, homogenous, isotropic, Newtonian fluid in the gravity field.In the near-vent portion of the flow, neglecting nonNewtonian effects is a reasonable assumption, as a consequence of the high temperature.The constant density t and viscosity h are used.Lava flows in a rectangular channel with a constant slope.The steady-state motion is assumed to be unidirectional in the azimuthal direction: the flow is parallel to the levees.This implies that no radial component of the gravity force acts on it.This condition often applies when a channel changes direction as a consequence of local topography.Further assumptions are introduced: the velocity only depends on the radial    where: (20) The flow is not symmetric with respect to r c , the center of the channel.The gradient of the velocity is greater close to the levee with the higher curvature.This behavior is particularly noticeable in comparison with the velocity in a rectilinear channel (Figure 12), which is calculated considering analogous hypotheses: the unidirectional flow of a Newtonian fluid in a gravity field that depends on the horizontal coordinate, while dependence on depth and direction of flow is neglected.In this case, the velocity is symmetric with respect to the center of the channel, where it reaches its maximum.The maximum velocity in the curved channel is shifted toward the internal levee.The extent of the displacement from r c is a decreasing function of the radius of curvature, which shows that the asymmetry of the velocity is higher for sharp, rather than smooth, bends, and it increases with channel width (Figure 13).
We calculate the viscous stress as a function of the radius from the constitutive equation of a Newtonian fluid, obtaining: (21) In analogy with the velocity gradient, the shear stress reaches greater values in the region close to the internal levee.In particular, it reaches its maximum in the internal levee.We compared the viscous stress calculated in a curvilinear channel with the stress in a rectilinear channel (Figure 14).Here, the curvature produces a reduction in stress close to the external levee, and an increase close to the internal levee.The maximum shear stress is calculated as a function of r c : it is a decreasing function of the radius of COOLING OF CHANNELED LAVA FLOW    We consider radiation and convection into the atmosphere to be the dominant mechanisms in the cooling process.We calculate iteratively the surface heat-flow density as a function of temperature along the down-flow direction.( 22) We use the model of cooling of a conductive half-space, described by the time-dependent heat equation for the heat flux [Turcotte and Schubert 1982].
Substituting the Fourier law in the solution, and on the assumption that the temperature is initially uniform in the medium, we obtain the temperature.Substituting the time t according to Equation (12), we obtain the flow temperature as a function of the coordinates z and L, which represent the vertical coordinate and the distance covered, respectively: where T 0 is the initial uniform temperature, and ū is the average velocity.
As long as the lava cools, crust platforms form at the flow surface.These platforms of solidified lava are driven by the underlying fluid, and laterally they are limited by the action of the shear stress.Greeley [1971] described some effects of the curvature of a channel on the lava surface morphology in Hawaiian volcanic fields.The crust formed at a distance of a few meters from the eruptive vent.At the first bend in the channel, this crust breaks into plates, which float on the flow surface.The increased velocity gradient provokes fragmentation of the plates and their inclusion in the flow through remelting.Peterson et al. [1994] also noted that the crust plates break in the lateral flow region when they run across a curve.
A common phenomenon in bends is the presence of a different level of deposit on the two levees of the channel, which can be observed when the eruption ends and the lava drains.The level of deposits is considerably higher on the external levee rather than on the internal levee, as observed on Kilauea volcano, Hawaii [Heslop et al. 1989] and on Pico Partido volcano, Lanzarote [Woodcock and Harris 2006].The centripetal force provokes only a slight difference in height between the two lateral walls.Its effects are evaluated from the radial component of the Navier-Stokes equation, which considers the variation in pressure due to the variation in height in the flow surface.The calculated superelevation explains an increase in the height lower than 10% of the superelevation inferred from the data.The present model applies to cases where the flow is parallel to the levees, whereas the measured different amounts of deposits occur in flows where the radial component of the gravity force acts in bends, due to the local topography.

Surface radiance and implications for satellite thermal imagery
The effusion rate of the lava from an eruption vent is the primary quantity that controls the evolution of the lava flow that ensues.For this reason, a lot of effort has been devoted to the evaluation of this quantity, which involves measurement of the flow velocity and the cross-sectional area [e.g.Calvari et al. 2002].The knowledge of effusion rates also has a major role in real-time simulations of lava flow paths that are carried out when lava flows threatens inhabited areas [e.g.Ishihara et al. 1990, Miyamoto and Sasaki 1997, Vicari et al. 2007, Rongo et al. 2008, Hérault et al. 2009].
Direct measurements of effusion rates in the field are difficult, and calculations from other flow parameters would require the measurement of such parameters [Tallarico et al. 2006], which can be equally problematic.In recent years, the availability of thermal images of volcanoes during eruptions has stimulated the evaluation of lava effusion rates from the measurement of their heat radiation.The thermal data obtained from high-resolution radiometers mounted in Earth satellites were first used to monitor active lavas by Mouginis-Mark et al. [1994] and Harris et al. [1995].Thermal images from hand-held radiometers on the ground have also been used [Harris et al. 2005b].
The use of radiance maps for the evaluation of lava effusion rates is made possible by simple formulae that relate the lava flow rate to the energy radiated per unit time from the planimetric surfaces of the flow.Such formulae are based on a specific flow model, and consequently, their validity is subject to the model assumptions.The evaluation of lava effusion rates is based on a formula originally proposed by   Pieri and Baloga [1986] that relates the flow rate to the planimetric area of a flow: (24) where q is the volume flow rate, W is the heat radiated per unit time from the flow surface, and DT is the difference between the lava temperatures at the vent and at the front.
If the flow rate is assumed to coincide with the effusion rate, this formula states that the effusion rate is proportional to the heat radiated per unit time by the surface of the flow.Harris et al. [1997] modified the formula, by including the thermal contribution W 1 of convection in the air, and the degree { 0 of crystallization of the lava.The contribution W 2 of heat conduction to the ground was included later, obtaining [Harris et al. 2005a]: (25 where c L is the latent heat of crystallization. The radiance of a hot body is related to the surface temperature at a given instant of time.If the body is a flowing liquid, the connection between the radiated heat and the flow rate is not straightforward.The singularity of Equation ( 24) was noted by Wright et al. [2001].They suggested that the apparent success of the formula in giving reasonable values of effusion rates is due to numerical coincidence between a constant by which Harris et al. [1997] multiply space-based estimates of lava flow area, and the empirical ratio between effusion rates and flow areas reported by Pieri and Baloga [1986] for a suite of Hawaiian flows.Harris et al. [2007] asserted that the technique gives time-averaged discharge rates.Dragoni and Tallarico [2009] proposed a different explanation.Indeed, the formula of Pieri and Baloga [1986] can only be obtained in the framework of a lava flow model based on several simplifying assumptions.They explicitly addressed such assumptions, to ascertain which of them are acceptable approximations of real lava flows, and which instead impose strong limits to the applicability of the model.They also considered whether current use of the formula is consistent with the model itself.The analysis of the thermal model from which the formulae by Pieri and Baloga [1986] and Harris et al. [1997] are obtained shows that the measurement of only the total instantaneous heat flow from the flow surface is not sufficient to calculate the lava flow rate.It is thus necessary to provide further information, i.e., the difference between the lava temperatures at the vent and at the front.If this difference is assumed to be the same for all lava flows, as is currently assumed, then the formula states that the effusion rate is an increasing function of flow length and a decreasing function of flow thickness, which is unrealistic.Apparently reasonable results are obtained as the assumption of a uniform crust temperature is compensated for by the inconsistent use of the formulae.In the real world, the average surface temperature is lower for longer flows, and is higher for thicker flows.The measured heat flow incorporates these effects, which happen to counterbalance the use of a constant temperature difference between the vent and the front of the flow.

Conclusions
The models presented in the present study indicate the relevance of thermal and rheological processes on the dynamics of channeled lava flows.In particular, the effects of nonlinearity on the velocity and the flow rate of a gravitydriven lava flow are not negligible, especially for steep slopes.Moreover, the effects of nonlinearity on the average velocity are enhanced by an increase in the fluid consistency parameter.Generally, the average velocity tends to increase with the aspect ratio of the channel, and this effect grows by decreasing the exponent n of the power-law rheology.
The rapidity of crust formation due to the constant heat flow at the channel surface strongly depends on the degree of non-linearity of the rheology.The crust grows, attains its maximum width, and stops widening, depending on the rheology.While for the dilatant fluid, the crust stops widening very close to the point where the temperature reaches solidus, a pseudoplastic fluid will cover a much longer distance before the crust reaches its maximum width.This result is a consequence because with our assumptions, the pseudoplastic fluid flows with a greater velocity than the Newtonian fluid, which in turn flows with a greater velocity than the dilatant fluid.
Numerical tests were carried out with the purpose of testing the stability and accuracy of the solution of the thermal equation.We find that the numerical error decreases as the number of cells of the mesh increases, and it does not propagate with time.
The presence of bends in the lava path significantly affects the flow dynamics.A simplified form of the equation of motion for a Newtonian fluid has been solved, which analytically obtained the velocity and stress fields in a cylindrical coordinate system as a function of the radial coordinate.Fluid velocity and viscous stress show asymmetric behavior with respect to the center of the channel.The maximum of the velocity is shifted towards the internal levee, in proximity of which the shear stress also reaches higher values.Although we use strong assumptions to simplify the dynamic, rheological and thermal aspects of lava flows, the model provides a possible explanation of some field observations.An analysis was carried out of the formulae currently used to evaluate the effusion rate of lava flows on the basis of the measurement of the instantaneous heat flow from the flow surface.The analysis of the model from which the formulae are derived shows that further information should be supplied, i.e. the difference between the lava temperatures at the vent and at the front.If this difference is assumed to be the same for all lava flows, as is currently assumed, the formula states that the effusion rate is an increasing function of the flow length and a decreasing function of flow thickness, which is not realistic.
Consistent application of the models presented here to actual eruptions needs reliable data from both field and laboratory measurements of viscosity, which deserve much more effort in the future.

Figure 3 .
Figure 3. Velocity profiles of the Laghetto channel (Mount Etna, 2001).a) Horizontal profiles at z = 0. b) Vertical profiles, at y = 0.The velocity profiles are plotted for three different values of the power-law exponent n (0.5, 1, 1.5).The parameters of the Etna channel are: a = 12 m, h = 6 m, a = 5˚(data from Ferlito and Sievert [2006]); the parameters of the lava are: t = 2030 kg m −3 and k = 8625 Pa s n (data fromPinkerton and Sparks [1978]).

Figure 4 .
Figure 4. Velocity profiles of the channel flow of the 2001 Mount Etna eruption.a) Horizontal profiles at z = 0. b) Vertical profiles at y = 0.The velocity profiles are plotted for three different values of the power-law index n (0.5, 1, 1.5).The parameters of the Etna channel are: a = 3 m, h = 1.2-1.5 m, a = 7.5˚; the parameters of the lava are: t = 2030 kg m −3 and k = 1.1 × 10 4 Pa s n (data from Harris et al. [2005a]).

Figure 5 .
Figure 5. Pu'u 'O'O (Kilauea, Hawai'i) case study.a) Average velocity v a vs n for two values of b. b) Volume flow rate q vs n, with q equal for the two values of b.The parameters of the channel are: a = 5 m, h = 3 m, a = 10˚; the parameters of the lava are: t = 2600 kg m −3 and k = 895 Pa s n .

Figure 6 .
Figure 6.Matches between the measurements of the flow velocity (crosses) and the velocity profiles (lines) computed with a dilatant model (n > 1) for the channel flow of April 2, 1984, Mauna Loa eruption.a) Horizontal profiles at z = 0. b) Vertical profiles at y = 0.The parameters of the channel are: a = 14.6 m, h = 4.1 m, a = 8˚; the parameters of the lava are: t = 1200 kg m −3 and k = 850-1400 Pa s n (data from Sakimoto and Gregg [2001]; viscosity from Harris and Allen III [2008]).

Figure 7
Figure 7. a) Yield strength x as a function of x´for different values of n. b) Crust coverage on the lava surface d/a as function of x´for different values of n (a = 10 m, h = 3 m, a = 0.2 rad, t = 2800 kg m −3 , k = 10 4 Pa s n , c p = 837 J kg −1 K −1 , K = 3 W K −1 m −1 , q 0 = 10 4 W m −2 , x 0 = 8100 Pa, b =170).x´is the point on the flux direction where the surface temperature reaches the solidus, and it varies with varying n.

Figure 8 .
Figure 8. Numerical tests on the heat conduction equation.Percentage errors as function of the number N CV of control volumes of the computational mesh at 301 s of cooling.Solid line, between the numerical and analytical solutions.Dashed line, between the actual grid and the coarser grid computed with a grid step of 10 × 10 CVs.

Figure 9 .
Figure 9. Numerical tests on the heat conduction equation.Percentage errors DT u as function of the number N CV of control volumes of the computational mesh at 301 s of cooling.Solid line, DT u between the numerical and analytical solutions.Dashed line, DT u between the actual grid and the coarser grid computed with a grid step of 10 × 10 CVs.

Figure 10 .
Figure 10.Numerical tests on the heat conduction equation.Percentage errors as a function of the number N CV of control volumes of the computational mesh at different times t = 101, 201, 301 s of cooling.
the pressure gradient is negligible.This obtains equation of motion as: (17) where u is the velocity, and: (18) By solving the equation with no slip boundary conditions at the lateral walls of the channel, we obtain the velocity as a function of the radius at the surface of the flow: (19)

Figure 11 .
Figure 11.Sketch of the model for bends and the coordinate system.

Figure 12 .
Figure 12.Comparison between the velocity as a function of the horizontal coordinate in a rectilinear channel (dashed line), and in a curvilinear channel (solid line).The vertical line indicates the center of the channel r c .

Figure 13 .
Figure 13.Distance between the point where the velocity is maximum and the center of the channel r c , as a function of the radius of curvature of the channel.Each curve is calculated for a different value of channel width a.
reaches higher values with sharp bends.

Figure 14 .
Figure 14.Comparison between the viscous stress as a function of the horizontal coordinate in a rectilinear channel (dashed line), and in a curvilinear channel (solid line).The vertical line indicates the center of the channel r c .