MMIP Artemio Plana-Fattori
email: artemio.planafattori@agroparistech.fr

An Eulerian-Lagrangian Approach for coupling Fluid Flow, Heat Transfer and Liquid Food Product Transformation

Keywords:

Introduction

Computational fluid dynamics (CFD) is a powerful numerical tool that became widely used to simulate different processes in chemical engineering [START_REF] Kuipers | Computational fluid dynamics applied to chemical reaction engineering[END_REF][START_REF] Farmer | CFD analyses of complex flows[END_REF], including applications in the food industry [START_REF] Xia | Applications of computational fluid dynamics (CFD) in the food industry: A review[END_REF][START_REF] Norton | Computational fluid dynamics (CFD): An effective and efficient design and analysis tool for the food industry: A review[END_REF]. The design of reactors and development of products require the assessment of velocity, temperature and composition distributions on the entire equipment. Fluid flow, heat transfer and product transformation are coupled in a number of chemical engineering processes. Fluid flow and heat transfer determine the temperature field, which can influence the product transformation.

Transformation can in turn modify the thermal and rheological properties associated with the product. This two-way dependence has been taken into account in modeling of coupled processes.

Difficulties arise with the inclusion of complex transformation models.

Many problems of practical interest in chemical engineering imply particles and/or droplets and/or crystals whose size distribution evolves along the process equipment. The solution of such problems requires the inclusion of the population balance equation [START_REF] Ramkrishna | Population Balances Theory and Applications to Particulate Systems in Engineering[END_REF] into the CFD model. However, the computational burden placed by including the population balance equation as a dispersed phase modeling tool can be excessive even with state-of-the-art computers. This is the reason why engineering applications need to have a simple reduced model for discrete phases, without losing the detailed description of the phenomena inherently embedded in the population balance equation [START_REF] Attarakih | Integral formulation of the population balance equation: Application to particulate systems with particle growth[END_REF].

The numerical solvers used in CFD have become specialized, as a way for saving computer time and memory. In the present situation, only convection-diffusion-reaction equations can be included without depriving the solver efficiency. A first strategy for including the population balance equation into CFD models involves the discretization of the size distribution into a finite number of contiguous classes, as an attempt to represent the infinite range of sizes [START_REF] Kumar | On the solution of population balance equations by discretization -I. A fixed pivot technique[END_REF]. Methods of classes allow a discontinuous reconstruction of the size distribution over a truncated range [START_REF] Drumm | One-group reduced population balance model for CFD simulation of a pilot-plant extraction column[END_REF]; further, its inclusion into CFD models significantly increases the numerical cost because every class of size is represented using a convection-diffusionequation which needs to be solved with the fluid flow and heat transfer along the whole domain of interest. The use of numerous convection-diffusion equations in connection with those describing the fluid flow problem leads to large systems, which in turn are computationally heavy and time consuming to solve. A second strategy for including the population balance equation into CFD models involves the reduction of the size distribution to its moments [START_REF] Marchisio | Implementation of the quadrature method of moments in CFD codes for aggregation-breakage problems[END_REF][START_REF] Silva | Comparison of the accuracy and performance of quadrature-based methods for population balance problems with simultaneous breakage and aggregation[END_REF]. Methods of moments require relatively lower numerical cost but provide a small amount of information about the size distribution [START_REF] Drumm | One-group reduced population balance model for CFD simulation of a pilot-plant extraction column[END_REF].

There are product transformation processes of interest in chemical engineering whose coupling with CFD is even more difficult. For instance, agglomeration during spray drying is considered to be very difficult to control: it involves the atomization of the solution or suspension into droplets, then the mixing of spray and hot air, latter the drying of suspension droplets, and finally the collision of particles which might lead to coalescence or agglomeration. In this context, [START_REF] Verdurmen | Simulation of agglomeration in spray dryers: the EDECAD project[END_REF] have summarized the development of an Eulerian-Lagrangian strategy in order to combine elementary models for the main processes taking place in spray-dryers. Such a strategy combines the advantages of Eulerian and Lagrangian descriptions of fluid flow: on one hand, the velocity, pressure and temperature fields for the gas phase can be accurately predicted as continuous functions along the Eulerian frame; on the other hand, the microphysical properties associated with the particles under consideration can be evaluated along representative Lagrangian trajectories with the help of transformation models as complex as needed.

The present study combines the advantages of Eulerian and Lagrangian descriptions, focusing on the coupled representation of fluid flow, heat transfer and transformation of a liquid food product. Changes in physicochemical properties of liquid products are relevant in the dairy industry as well as in the preparation of almost any kind of liquid food product (soups, juices, sauces, etc). These products can be represented by a single phase fluid, in the sense that all the elements (water, solutes, macromolecules, small particles) move according to a same velocity field.

Fluid parcels containing all these elements are hereafter studied as the Lagrangian entities of the problem under consideration, rather than particles of a dispersed phase as usual.

The Eulerian-Lagrangian approach presented in Section 2 is conceived to be flexible enough for including any available methodology able to evaluate the transformation state of a liquid food product under heat treatment. The approach is illustrated by the evolution of an aqueous suspension of starch granules which progressively swell under continuous heat treatment; the transformation state of the suspension is characterized by the swelling degree of granules in water, whose variation with time can reasonably be represented through a kinetics equation (Section 3). Such a choice allows us to solve also the whole problem through a purely-Eulerian approach (Section 4). Despite the challenging appearance of the coupled problem under consideration, the algorithm implementing the Eulerian-Lagrangian approach is relatively simple (Section 5). Finally, the availability of a purely-Eulerian solution for the whole problem provides the opportunity to assess the consistency of the Eulerian-Lagrangian approach (Section 6).

Eulerian-Lagrangian Approach: Overview

The proposed Eulerian-Lagrangian approach is based on an iterative algorithm, as presented in Figure 1. As an example, the liquid food is assumed flowing in a tubular heat exchanger, which is represented by a rectangular domain: the dashed red and the bold black lines indicate the axis of symmetry and the heating wall, while the left and right sides represent the domain inlet and outlet respectively. Along this domain, fluid parcels running nearer to the wall move slower than those at the axis; dynamical and thermal histories differ accordingly.

The first iteration starts with the Eulerian calculation of velocity, pressure, and temperature fields using given thermal and rheological properties. For this first iteration, viscosity and thermal conductivity are assumed to exhibit simple behavior in the whole equipment, such as uniform value or simple thermal dependency.

The following step is the estimation of Lagrangian trajectories followed by representative fluid parcels of the liquid food product. This step involves also the retrieval of the shear rate and temperature T as functions of time along every trajectory.

Next, such shear rate and temperature histories are employed in evaluating the liquid food's transformation state along every trajectory. This is the key step regarding the flexibility of the Eulerian-Lagrangian approach proposed. Any available methodology can be employed in evaluating the degree of transformation of the liquid food product along the process unit: population balance models and even sub-scale models (for instance, involving transport phenomena inside particles of a liquid suspension), dynamic Bayesian networks, stochastic models, etc. The transformation state is evaluated independently from the fluid flow and heat transfer coupled problem; the model chosen for representing the transformation of the liquid food product has to be integrated over the time only, and not over time and space as in the case of the Eulerian approach. Once evaluated the transformation state, the thermal and rheological properties of the liquid food can be estimated along every trajectory with the help of suitable laws.

At this point, liquid food's thermal and rheological properties have been estimated along the domain of interest through a dense sampling with Lagrangian trajectories. The following step is devoted to the interpolation of such a large set of point estimates to the Eulerian frame, in which the variables of interest are continuous functions of position.

The final step is the comparison between updated and previous estimates of liquid food's thermal and rheological properties. If the former compare favorably to the latter, that is if their maximum relative difference reaches a value below a certain threshold, then the iteration series is stopped and the updated estimates are provided as answer of the problem. Otherwise, the updated estimates are employed in a new iteration.

The respective advantages of Eulerian and Lagrangian descriptions can be easily recognized.

On one hand, the Eulerian description of fluid flow and heat transfer can be applied for accurately obtaining velocity, pressure and temperature fields as continuous functions of position along the domain under consideration. On the other hand, the study of Lagrangian entities (the fluid parcels) offer the possibility of evaluating any product transformation of interest through any available methodology along a number of representative trajectories.

The main assumption regarding the application of the Eulerian-Lagrangian approach is that each fluid parcel represents a closed system. Such a hypothesis is most often valid under laminar liquid flow. Indeed, typically the mean residence time is often of the order of 100 s, and molecular diffusion (denoted D s ) is of the order of 10 -9 m 2 .s -1 ; hence the random walk of the molecules inside the fluid parcels between inlet and outlet of the equipment is typically less than 1 mm as estimated using the diffusion length l = t D 2 s . For micrometric particles, diffusivity is about 10 -12 m 2 .s -1 , so that the random walk is less than 20 µm.

Case Study: Starch Granule Swelling

Starch is the most widely employed thickener in the food industry. Most industrial processes for starch-thickened products involve temperatures that can be high enough for enabling starch phase transitions. The concept of "gelatinization" refers to the destruction of the crystalline structure in granules, which is an irreversible process that includes, in a broad sense and in time sequence, granular swelling, native crystalline melting and molecular solubilization [START_REF] Liu | Thermal processing of starch-based polymers[END_REF]. We focus our attention on the granule swelling, whose consequences on the viscosity of starch suspensions have been observed at least since the 1940's. Early contributions [START_REF] French | Physical properties of starch[END_REF] have identified a first reversible phase associated with limited granule swelling and weak viscosity increase, and a second phase within a small range of temperature at about 65 °C in which starch granules suddenly increase many times in size while the viscosity rapidly rises. More generally, the viscosity of the starch suspension relative to that of water was found to be dependent on the volume occupied by the particles, their shape and deformability; moreover, both swelling and deformability can depend on the heating rate employed to gelatinize the starch [START_REF] Ellis | A comparison of the viscous behaviour of wheat and maize starch pastes[END_REF]. Starch research is a very active subject, and much effort have been done in characterizing native and modified starch, including various genotypes and climates (e.g., [START_REF] Singh | Morphological, thermal and rheological properties of starches from different botanical sources[END_REF][START_REF] Xie | Rheological properties of starches with different amylose/amylopectin ratios[END_REF]. A number of starch swelling models have been proposed from experimental work. For instance, second-order kinetics have been successfully applied for describing the evolution of a) modified waxy maize starch until temperatures above 100 °C, b) native waxy maize starch until granules start to rupture, and c) native maize starch [START_REF] Lagarrigue | Swelling kinetics of wazy maize and maize starches at high temperatures and heating rates[END_REF].

The Eulerian-Lagrangian approach overviewed in Section 2 approach is hereafter illustrated in studying the evolution of modified waxy starch granules under hydrothermal treatment with the help of a thermo-rheological model based upon experimental work. Following [START_REF] Lagarrigue | Swelling kinetics of wazy maize and maize starches at high temperatures and heating rates[END_REF], the transformation state characterizing the evolution of starch granules is expressed in terms of the swelling degree defined by:

0 M 0 D D D D S (1)
where D is the mean volume-weighted diameter of starch granules at a given time along the heat treatment, where D 0 and D M are the mean diameter before any heat treatment and after complete heat treatment, respectively. The progressive increase of the swelling degree with time along the heat treatment is evaluated through a second-order kinetic equation,

2 S 1 V t d S) 1 ( d , ( 2
)
where the rate constant V is estimated through the Arrhenius law,

T R E exp V ) T ( V A 0 . ( 3
)
where E A is the activation energy. In the scope of the Eulerian-Lagrangian approach, equation ( 2) is integrated along Lagrangian trajectories of fluid parcels. After evaluating the swelling degree associated with the starch suspension, the mean volume-weighted diameter of its granules can be estimated from equation (1),

S ) D D ( D D 0 M 0 . ( 4
)
The granule volume fraction is the key parameter for assessing the relative viscosity increase, and it can be estimated as:

3 0 0 D D , ( 5
)
where 0 and D 0 values characterize the starch suspension before any thermal treatment.

The liquid food product under consideration is an aqueous suspension of starch granules.

The apparent viscosity associated with the starch suspension is obtained by multiplying the continuous phase (water) viscosity W and the relative viscosity R which is a function of the transformation state. The water viscosity can be suitably evaluated as a function of the temperature (e.g. IAPWS, 2008). Various approximations have been proposed for estimating the relative viscosity as a function of the volume fraction ( ) occupied by the particles under consideration (see for instance [START_REF] Cheng | Exponential formula for computing effective viscosity[END_REF]. We have retained the model proposed by [START_REF] Thomas | Transport characteristics of suspension VIII: A note on the viscosity of Newtonian suspensions of uniform spherical particles[END_REF],

16.6 exp 0.00273 10.05 2.5

1 2 R W , ( 6
)
because it provides an increase of relative viscosity with the particle volume fraction which is intermediate to those obtained with other approximations.

Selected starch dispersions have been studied by [START_REF] Lagarrigue | Swelling kinetics of wazy maize and maize starches at high temperatures and heating rates[END_REF]. Those authors followed the kinetics of swelling of starch granules in aqueous suspension when applying thermo-mechanical treatments in a laboratory simulator coupled to a rheometer. Data on granule swelling and suspension viscosity development were obtained under controlled temperature and shear conditions. In this study we focused the attention on waxy maize starch; its solubility under heat treatment is negligible because the low level of amylose molecules in its constitution. Further, we restrict the analysis to waxy maize starch which has been chemically modified, preventing the rupture of swollen granules. Hereafter, we consider the evolution of the 3.1 % w/w modified waxy maize starch dispersion studied by [START_REF] Lagarrigue | Swelling kinetics of wazy maize and maize starches at high temperatures and heating rates[END_REF]; it corresponds to D 0 15.5 and D M 39.6 m, 0 0.03, V 0 2.9 10 12 s , and E A 84.8 kJ.mol .

Purely-Eulerian Approach

Models and methods

Conservation equations for mass, momentum and energy can be written, under steady-state conditions, as:

0 v (7) T p v v I v v (8) T T C P v (9)
where , C P , and are the density, the heat capacity and the thermal conductivity, respectively.

Looking for a representation of the second-order kinetic equation ( 2) with the help of the Eulerian description of fluid flow, starch granule swelling can be evaluated through the following convection-diffusion equation:

S D S - 1 V S S 2 v , ( 10
)
where V is the rate constant expressed in terms of the Arrhenius law (3) and D S is the diffusion coefficient.

Coupled phenomena involving fluid flow, heat transfer and starch granule swelling are hereafter illustrated inside a tubular heat exchanger with length 1 m and radius 5 mm, which is represented through a two-dimension axial-symmetric domain. The liquid food product under consideration is an aqueous suspension of starch granules; it is assumed to be a Newtonian liquid, and its density and thermal properties are approximated by those associated with water. A uniform flux of 10 kW.m -2 is applied on the heating wall, and the temperature is assumed to be 20 °C at the domain inlet. A fully developed parabolic flow profile was assumed at the inlet. The flow rate was fixed at 10 liters per hour (mean velocity of 0.035 m.s ). According to the Reynolds number, the flow regime is laminar. These conditions enable a maximum temperature near the wall at the outlet of about 91 °C, which corresponds to a relevant degree of starch granule swelling (see Table 2 of [START_REF] Lagarrigue | Swelling kinetics of wazy maize and maize starches at high temperatures and heating rates[END_REF].

Equations (7-10) were solved under steady-state conditions with the finite-element method (simulation package COMSOL Multi-physics, version 3.5a). Second-order Lagrange finite elements were employed in all cases, excepting the pressure for which first-order finite elements were assumed. We choose the Parallel Sparse Direct Linear Solver (PARDISO) for solving the large systems resulting from the linearization of the coupled equations under consideration. PARDISO implements an efficient method for large systems of equations like those arising from the finite element method [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with PARDISO[END_REF]. Results discussed below were obtained by assuming a relative tolerance of 10 -6 . Looking for the solution of the whole problem, satisfying model convergence was reached after assuming the diffusion coefficient D S = 10 -12 m 2 .s -1 in equation ( 10).

Computations were carried out in a 64-bits calculator disposing of 24 Gb of RAM.

Mesh influence

Sensitivity tests were conducted with the purely-Eulerian modeling of the whole coupled problem in order to assess the influence of the mesh resolution. As [START_REF] Kuipers | Computational fluid dynamics applied to chemical reaction engineering[END_REF] have emphasized, only those computational results that possess invariance with respect to discretization could in a further step be confronted with experimental data. Following paragraphs discuss the influence of mesh resolution on mass-weighted results at the domain outlet (i.e. candidates for a future observation) and on local results along the outlet (i.e. the radial profiles).

Meshes under consideration are constituted of identical rectangular cells characterized by dimensions dz and dr along the domain's length (1 m) and radius (5 mm), respectively. In our case, the velocity, temperature, product transformation state, and viscosity gradients are mainly oriented along the domain's radius and therefore dr has to be small. On the other hand, dz can be relatively larger because the flow is mainly oriented along the domain's length. Selected aspect ratio values were assumed in order to compare results obtained with "squares" (dz/dr = 1), "short-" (= 5) and "long-rectangles" (= 10).

Figure 2 presents bulk (mass-weighted) estimates of temperature (left) and starch swelling degree (right) at the domain outlet as functions of the total number of cells under consideration (top displays). Bulk estimates exhibit a near asymptotic behavior as the total number of cells increases.

In the case of meshes with about 250000 cells, results obtained with distinct aspect ratios exhibit small differences: less than 0.002 °C in temperature and 10 -4 in swelling degree. Coarser resolutions (less than about 25000 cells) provided unreliable results somewhere in the domain, for instance swelling degree values higher than the unity. Finer resolutions (more than about 250000 cells) gave origin to LU factorization out of memory, indicating the limits in applying the PARDISO solver with our present calculation capabilities.

Bottom displays in Figure 2 show the same results as function of the number of cells in the domain's radial direction. The interest on such a presentation can be justified by the fact that the main variables of interest (velocity, temperature, transformation state, and suspension viscosity) are expected to exhibit more important variations from the axis of symmetry to the heating wall than from the inlet to the outlet. Both presentations in Figure 2 put in evidence that results obtained with "long rectangles" (dz/dr = 10) converge faster than for shorter ones with the increase of the number of cells. Among the meshes associated with dz/dr = 10, we retained as reference the one which is constituted of 100 rectangles along the domain's radius and 2000 rectangles along the length (hereafter, the 100x2000 mesh). This mesh is further referred as mesh 1A.

Table 1 compares, in its upper section, outlet bulk values of temperature and swelling degree after assuming selected meshes (bulk values are 20 °C and zero at the domain inlet). Relative differences are evaluated with respect to the range of bulk values between the outlet and the inlet;

for instance, the relative difference in temperature for test 1B (0.02 %) results from the difference between tests 1B and 1A (5.33 10 -3 °C) divided by the difference between the outlet (47.23 °C) and the inlet (20 °C) bulk temperatures from test 1A.

Results obtained with the 100x2000 mesh (test 1A) compare favorably to those obtained after assuming a 50 % smaller number of cells on both directions (50x1000, test 1B). Results become progressively poorer after assuming coarser meshes, like those constituted of 25 and 10 % of the number of cells previously assumed on both directions (tests 1C and 1D).

The lower section of Table 1 summarizes the range of values assumed along the outlet by temperature, swelling degree and suspension viscosity after assuming the 100x2000 mesh (test 1A).

It shows also the maximum absolute impact of replacing the mesh 1A by the meshes 1B to 1D. The adoption of progressively coarser meshes leads to larger impact, particularly in the case of the suspension viscosity.

In summary, a number of tests were performed regarding the mesh influence along the implementation of the purely-Eulerian approach for solving the whole problem under consideration.

Tests summarized in Figure 2 and Table 1 support the adoption of the 100x2000 mesh in looking for the solution of the whole problem. This mesh is also adopted in the Eulerian step (fluid flow and heat transfer only) of the Eulerian-Lagrangian approach.

Overview from purely-Eulerian results

Main results obtained with the numerical model simulations described above are presented in Figure 3. Its left section displays the temperature (top), swelling degree (center) and suspension viscosity (bottom) fields in the domain under consideration. As in Figure 1, the liquid food product is assumed flowing in a heat exchanger whose inlet and outlet are represented by the left and right sides of the rectangular domain; the liquid food product is heated from the wall (bottom of domain).

The suspension viscosity field was deformed according to the velocity, and inlet and outlet velocity profiles can be appreciated at domain extremes.

Fluid parcels move slower and reach higher temperatures when running near it. Starch granules swell after heating, and their volume fraction increases as well as the suspension viscosity.

In comparison with the conditions prevailing at the domain inlet, the suspension viscosity increases by a factor of about 7 near the corner of the domain outlet with the heating wall, slowing down the velocity in this region.

These results follow, in a partial extent, those obtained by [START_REF] Liao | Role of thermo-rheological behaviour in simulation of continuous sterilization of a starch dispersion[END_REF]. Those authors also observed in the heating section of their tubular exchanger an increase of viscosity near the heating wall due to gelatinization, leading to a velocity decrease near this wall and a velocity increase at the axis of symmetry. Nevertheless starch type and concentration were different (4 % of waxy rice starch), and heating conditions were stronger (145 °C); therefore a direct comparison of results seems difficult. Their maximum suspension viscosity was much higher (about 6 Pa.s), which caused velocity profiles almost flat were the food product was very viscous; this was not observed in our case. Other geometries than tubular exchangers have been adopted in numerical models coupling fluid flow, heat transfer and starch swelling, like axially rotating cans (e.g., [START_REF] Tattiyakul | Simulation of heat transfer to a canned corn starch dispersion subjected to axial rotation[END_REF]. Any tentative of comparison with such results seems far to be straightforward.

The evolution of liquid food product properties can also be studied with the help of hypothetical trajectories of fluid parcels released at the domain inlet. Each fluid parcel can be interpreted as a suspension droplet containing an ensemble of starch granules at similar thermodynamic and kinetic conditions. Figure 3 displays, in its right section, selected properties associated with fluid parcels released at the domain inlet, at 0.5 mm (black curves) and at 1 mm (gray curves) from the heating wall. Because the maximum velocity takes place at the axis of symmetry (exchanger center), fluid parcels released at 0.5 and 1 mm reach the domain outlet after about 86 and 36 s respectively. When reaching the domain outlet, fluid parcels released at 0.5 mm are associated with higher values of temperature, of mean diameter, and of suspension viscosity than those released at 1 mm from the heating wall. Other the granule starch swelling, the suspension viscosity depends on the continuous phase (water) viscosity. The latter slowly decreases with heating, and this fact explains why the suspension viscosity is slightly reduced during the first seconds after the fluid parcel release.

Positions A and B in Figure 3 were selected because they correspond to a same temperature value (about 66 °C). Position A represents the outlet for fluid parcels released at 1.0 mm from the heating wall; they reached this position after running about 36 s since the inlet. Position B is located before the domain outlet; fluid parcels released at 0.5 mm from the heating wall move slower and spent about 55 s to come there since the inlet. Despite the same temperature, position B corresponds to a mean diameter of starch granules that is about 10 % greater, and to a suspension viscosity about 45 % greater, than respective values at the position A. Such differences are explained by the fact that fluid parcels running at different distances from the heating wall move with different velocities and hence experience different temperature time series (or histories). We conclude that the influence of temperature on the suspension viscosity is not instantaneous. The liquid food product viscosity depends on the continuous phase viscosity as well as on the starch granule volume fraction. The former depends on the present local temperature whereas the latter depends on the whole thermal history of respective fluid parcels.

Eulerian-Lagrangian Approach

Models and methods

The Eulerian-Lagrangian approach has been overviewed in Section 2. A detailed algorithm is presented in Figure 4, taking into account the liquid food product transformation presented in Section 3. Rectangles summarize the most relevant steps, where variables indicated on the left (and right) sides constitute key inputs (and outputs).

Geometry and boundary conditions taken into account have been presented in the scope of the purely-Eulerian modeling of fluid flow and heat transfer. The 2000x100 mesh was retained (see Sub-section 4.2).

Before the first iteration (k = 1), a preliminary Eulerian modeling of fluid flow and heat transfer is needed. It provides a first estimate of the velocity v k=0 (x) and the temperature T k=0 (x)

fields along the Eulerian frame (positions x). Because this modeling step is performed before any liquid food product transformation, the viscosity is here assumed to be the one of the continuous phase (water).

The first step of a given k-th iteration is the estimation of Lagrangian trajectories associated with a number of representative fluid parcels. Each i-th trajectory is defined by a number of j-th positions, each position being described by axial (z) and radial (r) coordinates: x k = (r k i,j , z k i,j ), with i= max and j= max . Each of these positions is associated with a given time t i,j . The estimation of positions x k = (r k i,j , z k i,j ) from the Eulerian velocity field v k-1 (x) is performed with the help of a routine implementing a fourth-order Runge-Kutta method [START_REF] Abramowitz | Handbook of Mathematical Functions, with Formulas, Graphs, and Mathematical Tables[END_REF].

Fluid parcels are released at the domain inlet, where starting positions (r i,j=0 , z i,j=0 ) are uniformly distributed from the axis of symmetry towards the heating wall. No fluid parcels are released at the heating wall because the latter corresponds to vanishing velocities and hence to unphysical thermal histories. In order to represent the physical conditions prevailing at the heating wall, an additional trajectory is placed nearby; the radial component of the velocity is set to zero and the resulting trajectory is parallel to the heating wall. This additional trajectory is arbitrarily placed at a distance from the heating w width (5 mm / 1000 = 0.5 micrometer). In summary, the release positions of fluid parcels are the same for all the iterations that have to be considered, and they are defined as:

r i,j=0 = (i/i max ) R for i=1 max 1, r i,j=0 = R 0 -7 m for i=i max , and (11)
z i,j=0 = 0 for i=1 max ,
where R is the radius of the heat exchanger under consideration. A number of sensitivity tests were conducted regarding the number of trajectories as well as the number of positions along every trajectory (Sub-section 5.3).

Lagrangian trajectories are then employed in the reconstruction of thermal histories along the domain under consideration. Thermal histories are described by a series of positions (r k i,j , z k i,j ), times t k i,j , and temperature values T k i,j . These temperature values are estimated by interpolation from the last available Eulerian temperature field T k-1 (x).

The liquid food product transformation state is here assessed by the starch swelling degree, whose methodology of evaluation is based upon experimental work (Section 3). At k-th iteration, the thermal history associated with i-th trajectory (times t k i,j and temperature values T k i,j ) allows us to integrate equation (2) between two successive positions of interest:

t d T(t) R E exp V S 1 S d k j i, k 1 - j i, k j i, k 1 - j i, t t A 0 S S 2 . ( 12
)
After integrating the left-hand side and approximating the right-hand side with the help of the trapezoidal rule, the swelling degree at a given position can be evaluated as:

k 1 - j i, k 1 - j i, k 1 - j i, k j i, S 1 1 S 1 S S (13) t t T R E exp T R E exp V k 1 - j i, k j i, k 1 - j i, A k j i, A 0 .
These expressions apply for all positions beyond the domain inlet (j=0), where the starch swelling degree is assumed to be null. Suspension viscosity is obtained by multiplying the continuous phase (water) viscosity and the relative viscosity associated with the granules volume fraction. At k-th iteration, the water viscosity is straightforwardly evaluated from the temperature value T k i,j . The availability of the swelling degree S k i,j allows the evaluation of the mean diameter D (equation 4), then of the granule volume fraction (equation 5), and finally the relative viscosity (equation 6).

Suspension viscosity values as estimated above can appreciably vary between two successive iterations, turning longer the numerical procedure. The influence of large oscillations is here reduced by combining such new estimates

j , i k ~
with those resulting from the previous iteration:

j , i 1 - k j , i k j , i k 1 ~ , (14)
where the relaxation coefficient is set to 0.50. The application of such a strategy requires suspension viscosity values j , i 1 k

; they are estimated by retrieving, in the suspension viscosity field k-1 (x) corresponding to the Eulerian step performed at the previous iteration, those values which correspond to (r k i,j , z k i,j ) positions defining the present trajectories of interest. This step provides an ensemble of updated suspension viscosity values, k i,j , which improves the knowledge available regarding the transformation of the liquid food product.

Updated suspension viscosity values have to be informed to the Eulerian solver, as a previous task towards an updated solution for the fluid flow and heat transfer problem. Difficulties at this step come from the fact that Lagrangian trajectories provide updated viscosity values at a number of positions which are irregularly distributed along the Eulerian frame. The Delaunay triangulation was retained as strategy. The application of constraints relying successive positions along a same trajectory as well as those relying corresponding positions between successive trajectories prove useful. As a result, the suspension viscosity values k i,j available at Lagrangian positions are successfully interpolated into a suspension viscosity field k (x) needed at the nodes constituting the Eulerian frame.

Before the following iteration, an Eulerian modeling of fluid flow and heat transfer is performed. It provides updated fields for the velocity v k (x) and the temperature T k (x). Because this modeling step is performed by taking into account an improved suspension viscosity field, it is expected that these velocity and temperature fields represent a somewhat better picture of the fluid flow and heat transfer conditions under the occurrence of product transformation.

The final step is the comparison between updated k (x) and previous estimates k-1 (x) of the suspension viscosity along the Eulerian frame. If their maximum relative difference reaches a value below 0.1 %, then the iteration series is stopped and the updated estimates are provided as answer of the problem. Otherwise, the updated estimates become the previous ones for the following iteration of the numerical procedure.

In summary, the algorithm presented in Figure 4 for the Eulerian-Lagrangian approach is a particular case of that proposed in Figure 1. In studying the evolution of one starch dispersion under continuous heat treatment, we assumed negligible a) the influence of the shear rate on the swelling kinetics, and b) the influence of the solid fraction associated with progressively swollen starch granules on the suspension thermal conductivity.

The algorithm described above for the Eulerian-Lagrangian approach was implemented through a main program and a set of functions developed in the programming language MATLAB. Eulerian modeling of fluid flow and heat transfer was solved with the help of the COMSOL Multiphysics 3.5a simulation package (see Section 4). The implementation of both tasks, including the adoption of other mathematical solvers, seems feasible in other programming languages.

Following sub-sections shows how intermediate results of the Eulerian-Lagrangian approach evolve towards the solution reached for the whole problem including the product transformation (5.2), and then presents a number of sensitivity tests regarding the influence of the number of trajectories as well as that of the number of positions per trajectory on the solution reached for the problem (5.3). Section 6 compares selected results provided by the purely-Eulerian and Eulerian-Lagrangian approaches for solving the whole problem.

Evolution of selected results along the iterating procedure

A summary of the algorithm implementing the Eulerian-Lagrangian approach was presented in Figure 4. The solution reached for the coupled problem is obtained as the result of a number of iterations, each of them including the evaluation of the product transformation along Lagrangian trajectories as well as the Eulerian modeling of fluid flow and heat transfer. The algorithm assumes velocity and temperature fields associated with transformation-free conditions (k = 0) in obtaining the first evaluation of product transformation along Lagrangian trajectories (k = 1). The algorithm stops when the product transformation gives origin to a suspension viscosity field sufficiently close to the one obtained at the previous iteration.

The evolution of selected variables along the iterating procedure is discussed from results obtained with the Eulerian-Lagrangian approach by considering 400 trajectories with 10000 positions each. Under such configuration, 22 iterations are necessary for reaching model convergence (that is, successive suspension viscosity fields which are closer than 0.1 % over the whole Eulerian frame). Remaining iterations globally follow the tendencies identified above. Dark gray lines (k = 2) and bold black lines (k = 22) show that i) axial velocity values increase at the exchanger center and decrease in the vicinity of the heating wall, ii) temperatures increase more importantly at the heating wall, and iii) both the starch swelling degree and the suspension viscosity increase at distances progressively larger from the heating wall. The abovementioned inflexion point in the axial velocity profile roughly divide the computational domain in two regions; the first, towards the heating wall, where the product transformation and its effects on the suspension viscosity are progressively higher, while fluid parcels moving slower are progressively more exposed to heating and product transformation; and the second region, towards the exchanger center, where the fluid parcels move faster and hence are relatively less exposed to heating and product transformation.

Selected sensitivity tests

The Eulerian-Lagrangian approach for solving the whole coupled problem involves the Eulerian modeling of fluid flow and heat transfer as well as the evaluation of liquid food product transformation along Lagrangian trajectories. It is therefore expectable that the solution reached depends on both the Eulerian mesh resolution and the distribution of Lagrangian positions in the computational domain. In this study these two influences are discussed separately. The influence of the mesh resolution on results obtained with the purely-Eulerian modeling of the whole problem was discussed above (Subsection 4.2).

The distribution of Lagrangian positions depends on the velocity field taken into account as well as on the methodology adopted for computing them. The velocity field is obtained iteratively along the execution of the algorithm, hence it constitutes either a consequence of previous iterations or the starting condition for the next one. As explained above (Subsection 5.1), Lagrangian positions are computed from velocity components estimated with the help of a fourth-order Runge-Kutta method. The implementation of the Eulerian-Lagrangian approach includes an additional trajectory in the vicinity of the heating wall, whose positions are evaluated neglecting the radial component of the velocity field; such an additional trajectory significantly improves the representation of physical conditions prevailing near the heating wall. In summary, the distribution of Lagrangian positions depends on: i) the number of trajectories;

ii) the number of positions per trajectory;

iii) the distance from the heating wall at which is located the additional trajectory; and iv) the distribution of release starting positions at the domain inlet.

Regarding the latter, it would be convenient to distribute the starting positions in order to improve the sampling of regions associated with higher spatial variation; for instance, the number of trajectories might be greater in the vicinity of the heating wall, where temperature and hence suspension viscosity gradients are expected to be higher. In this study, looking for a strategy independent on the magnitude of physical conditions prevailing near the heating wall, the starting positions were uniformly distributed along the domain inlet (eq.11). Regarding the additional trajectory in the vicinity of the heating wall, it was located at 1/1000 of the computation domain's width; it was found to be a satisfying compromise. Larger distances from the heating wall progressively decrease our ability in properly representing the dynamical and thermal gradients prevailing there. Following paragraphs discuss the influence of the number of trajectories as well as the influence of the number of free positions per trajectory on the results obtained with the Eulerian-Lagrangian approach.

Tables 2 and 3 summarize the tests conducted, presenting results at the domain outlet only. A same reference set of results (test A) is adopted in both cases; it was obtained with the Eulerian-Lagrangian algorithm by considering 400 trajectories with 10000 positions each. All the tests in Table 2 assume 10000 positions per trajectory, while 400 trajectories are assumed in all tests in Table 3. In both tables, there are two tests which represent a worsening of the algorithm performance in solving the whole coupled problem: a decrease by a factor of 2 (test B) and later of 4 (test C) in the number of trajectories and in the number of positions per trajectory, respectively.

A few conclusions emerge from these tests. The role played by the number of trajectories seems to be more important than the number of positions per trajectory, as indicated by tests 2C and 3C: decreasing the number of trajectories by a factor of 4 has a larger impact on temperature values than decreasing the number of positions by a factor of 4. On the other hand, as indicated by tests 2B and 3B, there is narrow room for improvement. We can therefore conclude that reference results (test 2A) are robust regarding slight changes in both the number of trajectories and the number of positions per trajectories.

Comparison between Purely-Eulerian and Eulerian-Lagrangian Approaches

The consistency of the Eulerian-Lagrangian approach for solving the coupled problem involving fluid flow, heat transfer and liquid food product transformation is hereafter assessed in the case of the starch swelling, by comparing its results at the domain outlet with those obtained through the purely-Eulerian modeling of the whole problem. It must be stressed that, in comparing these approaches, we are looking for consistency rather than accuracy: although classical in solving coupled problems including kinetic-type transformations, the purely-Eulerian approach is not free of sources of uncertainty.

Previous sensitivity tests have indicated the degree of convergence of the purely-Eulerian approach with the mesh resolution (Subsection 4.2). Results obtained after assuming a mesh with 100x2000 rectangles (test 1A) can satisfactorily represent this approach. Similarly, the influence of the number of trajectories and that of the number of positions per trajectory on the results obtained with the Eulerian-Lagrangian approach have both been assessed (Subsection 5.3). Results obtained with this approach can be illustrated by those assuming an Eulerian mesh with 100x2000 rectangles and 400 Lagrangian trajectories with 10000 positions each (test 2A). Table 4 compares selected results at the domain outlet after tests 1A and 2A. As in previous tables, relative differences between two sets of results are evaluated with respect to the difference between bulk values at the outlet and at the inlet.

There is a good level of agreement between the purely-Eulerian and the Eulerian-Lagrangian approaches: differences between outlet bulk results are no larger than 0.4 % for temperature and swelling degree. Maximum differences between respective outlet profiles are smaller than 1 %.

Figure 6 presents in its left section the outlet profiles of temperature, swelling degree and suspension viscosity obtained after assuming the 100x2000 mesh in the purely-Eulerian approach, as functions of the distance from the exchanger center. In its right section, Figure 6 shows the impact of replacing the purely-Eulerian approach by the Eulerian-Lagrangian one for solving the whole problem. With respect to the results obtained with the former, the latter approach slightly underestimates the temperature along the whole outlet. Regarding the starch swelling degree and the suspension viscosity, the difference between the approaches exhibits a broad minimum in the domain central region and a sharp maximum in the vicinity of the heating wall. The maximum difference between respective viscosity outlet profiles takes place in the vicinity of the heating wall, where dynamical and thermal gradients are significant and difficult to be numerically represented.

Such a difference remains small: 0.033 mPa.s compared to the local viscosity (~ 5 mPa.s).

Conclusion

Looking for a general tool for representing coupled processes involving fluid flow, heat transfer and liquid food product transformation, a numerical approach was proposed by combining advantages of the Euler and Lagrange descriptions. The product transformation is evaluated along trajectories of representative fluid parcels, and such an issue provides flexibility regarding the methodology applied for evaluating the product transformation state along the thermal and dynamical histories in the processing unit.

The approach consistency was verified in studying the evolution of a starch suspension in a heat exchanger. Quite standard conditions were assumed (cylindrical geometry, laminar Newtonian flow, uniform wall heat flux, steady-state conditions), and a kinetics equation was adopted for evaluating the swelling degree of the starch suspension. Despite such an apparent simplicity, the problem is challenging. Firstly, velocity and temperature gradients are higher in the vicinity of the heating wall; care is required in meshing the domain as well as in establishing the release protocol of fluid parcels. Secondly, starch granule swelling depends on thermal histories, and the latter are progressively longer towards the heating wall; hence Lagrangian trajectories have to include a high number of positions in order to reconstruct thermal histories. Finally, granule swelling severely impacts the suspension viscosity in the vicinity of the heating wall; therefore, the accurate representation of such a result asks for a fine Eulerian mesh as well as for a dense sampling with Lagrangian trajectories. Essentially the same results at the domain outlet can be obtained from a) the Eulerian-Lagrangian approach considering an Eulerian mesh with 100x2000 rectangles and 400 Lagrangian trajectories with 10000 positions each, and b) the purely-Eulerian modeling of the whole problem in that 100x2000 mesh. We emphasize the consistency between results obtained with the Eulerian-Lagrangian approach here proposed and similar with the purely-Eulerian approach. The latter can be verified by others through any modeling software able to solve a coupled problem involving fluid flow, heat transfer and a kinetics reaction.

Laminar Newtonian flow was considered above in illustrating the Eulerian-Lagrangian approach. No reasons forbid its application to non-Newtonian flow, i.e. to most fluids encountered either in nature (gums, proteins, blood...) or in technology (polymers, emulsions, slurries...).

Further, no reasons prevent its application to flow patterns associated with the occurrence of stagnant pockets or short-circuiting features which may exist in process equipment (e.g. [START_REF] Levenspiel | Patterns of Flow in Chemical Process Vessels[END_REF]. The Eulerian-Lagrangian approach can also be applied to turbulent flows, which are often encountered in engineering. However, its application to turbulent flow could require further development, because fluid parcel trajectories can exhibit stochastic behavior; hence, additional attention should be paid in interpolating Lagrangian apparent viscosity fields into the Eulerian frame. The approach is flexible enough to be applied to more complex geometry (including 3D problems), to other strategies for meshing the computational domain and also to other boundary conditions than those assumed in the case study above. The Eulerian-Lagrangian approach was illustrated for a simple transformation process (swelling of starch granules in water, under continuous heating) whose representation can be achieved through a simple methodology (a second-order kinetics equation). These choices were deliberate. Firstly and with the consistency tests in view, we considered a methodology which could be included in both purely-Eulerian and Eulerian-Lagrangian approaches. Secondly and putting in evidence the two-way coupling between fluid flow, heat transfer and transformation process, we choose a liquid product whose progressive evolution allows significant consequences on the fluid transport properties.

The Eulerian-Lagrangian approach allows the study of transformation processes more complex than the swelling of starch granules in water; further, this approach enables the representation of transformation processes through methodologies more sophisticated than kinetics reactions. Many problems of practical interest in chemical engineering imply particles and/or droplets and/or crystals whose size distribution evolves along the process equipment. In applying the Eulerian-Lagrangian approach, the discretization of the size distribution into a great number of classes does not add convection-diffusion-reaction equations to the Eulerian algorithm for solving the fluid flow and heat transfer coupled problem; hence the population balance equation can be successfully solved along a number of representative trajectories without demanding an exorbitant increase in computational resources. It has been accomplished by [START_REF] Chantoiseau | Coupling fluid flow, heat transfer and thermal denaturation-aggregation of beta-lactoglobulin using an Eulerian/Lagrangian approach[END_REF], in studying the thermal denaturation-aggregation of whey proteins in aqueous suspension. The Eulerian-Lagrangian approach can be applied to any chemical engineering problem in which a single velocity field applies to both liquid and dispersed phases, as in oil-in-water emulsions (Almeida-Rivera and Bongers, 2010) and crystallization under batch conditions [START_REF] Falola | CrystSim: A Software environment for modelling industrial batch cooling crystallization[END_REF]. On the other hand, additional development in our approach can be required in the case of problems in which two fluid velocity phases should be retained, as in bubble columns [START_REF] Bannari | Three-dimensional mathematical modeling of dispersed two-phase flow using class method of population balance in bubble columns[END_REF] and counter-current liquid-liquid extraction columns [START_REF] Jaradat | RDC extraction column simulation using the multi-primary one secondary particle method: Coupled hydrodynamics and mass transfer[END_REF].

Finally, we argue that any conceivable methodology for evaluating the transformation state of the liquid product can be properly taken into account in applying the Eulerian-Lagrangian approach, including sophisticated tools like stochastic models and dynamic Bayesian networks (e.g., [START_REF] Baudrit | Towards a global modelling of the Camembert-type cheese ripening process by coupling heterogeneous knowledge with dynamic Bayesian networks[END_REF].

Future work includes the implementation of the Eulerian-Lagrangian approach to real liquid food transformation, in order to evaluate its performance against experimental data.

Figure Captions

Figure 5

 5 Figure 5 displays outlet profiles of axial velocity, temperature, starch swelling degree, and suspension viscosity at four stages of the iterating algorithm. Clear gray lines (k = 0) present the conditions prevailing in the exchanger before the product transformation onset, as obtained with the Eulerian modeling of fluid flow and heat transfer only. Such transformation-free conditions are characterized by i) axial velocity values remembering the fully developed parabolic flow that was assumed at the inlet; ii) temperature values associated with a developing thermal boundary layer; iii) no starch swelling; and iv) viscosity values decreasing from the exchanger center towards the heating wall, coherently with the temperature-dependence of water viscosity.Thin black lines (k = 1) display results after the first iteration of the Eulerian-Lagrangian algorithm, more precisely after the first Eulerian modeling of fluid flow and heat transfer that was performed by considering a suspension viscosity field already including the effects of the first product transformation estimation. An inflexion point appears in the axial velocity profile: relatively to the former one, fluid parcels move faster towards the exchanger center and slower towards the heating wall. Such a major change reflects the increase of the suspension viscosity in the vicinity of the heating wall, which in turn is a consequence of the starch swelling. Such a viscosity increase is due to the volume fraction of solid particles in the suspension, whose effects more than counterbalance those associated with the temperature-dependence of water viscosity.

Figure 1 :

 1 Figure 1: Main steps of the Eulerian-Lagrangian approach.

Figure 2 :

 2 Figure 2: Bulk temperature (left) and swelling degree (right) values at the domain outlet as functions of the total number of mesh cells under consideration (top) and of the number of cells in the radial direction (bottom), as obtained with the purely-Eulerian modelling of the whole problem. Three sets of results are shown in this figure, corresponding to selected values of the aspect ratio dz/dr, where dz (dr) is the dimension along the domain's length (radius) of every rectangular cell constituting the mesh considered. Dotted lines were included for illustration purposes only.

Figure 3 :

 3 Figure 3: Main results provided by the purely-Eulerian modeling of the coupled problem developed in this study. Left section presents the temperature (top), swelling degree (center)and suspension viscosity (bottom) fields in the domain under consideration; inlet is on the left, heating wall is on the bottom. Right section presents the time series of selected properties corresponding to the trajectories followed by two hypothetical fluid parcels released at the domain inlet at 0.5 and 1 mm from the heating wall. Dark blue and red colors indicate minimum and maximum respective values: about 20 and 91 °C for temperature, 0 and 1 for swelling degree, and 0.7 and 7.1 mPa.s for suspension viscosity.

Figure 4 :

 4 Figure 4: Summary of the algorithm implementing the Eulerian-Lagrangian approach.

Figure 5 :

 5 Figure 5: Profiles of axial velocity, temperature, swelling degree and suspension viscosity field at the domain outlet as obtained with the Eulerian-Lagrangian approach for solving the whole coupled problem, at four stages of the algorithm: before the first (k = 0), after the first (k = 1), after the second (k = 2), and after the last (here, k = 22) evaluation of product transformation along Lagrangian trajectories.

Figure 6 :

 6 Figure 6: Left Profiles of temperature, swelling degree and suspension viscosity field at the domain outlet as obtained with the purely-Eulerian modeling of the whole problem, by assuming the 100x2000 mesh (test 1A). Right Impact of replacing the purely-Eulerian approach by the Eulerian-Lagrangian one for solving the whole problem, with the help of

Table 2 :

 2 Summary of sensitivity tests regarding the influence of the number of trajectories on the Eulerian-Lagrangian approach for solving the whole problem. All these results correspond to the domain outlet. All these tests consider 10000 positions per trajectory. the additional trajectory is considered, as well as all the release positions at the domain inlet (**) relative to the difference between corresponding outlet and inlet bulk values

	test 2B test 2C	400 / 2 = 200 400 / 4 = 100	2010201 1010101	Difference (and relative difference**) with respect to test 2A	+4.6 10 -4 (+0.002 %) +4.8 10 -4 (+0.002 %)	1.4 10 -6 ( 0.001 %) 9.1 10 -5 ( 0.03 %)	maximum (and maximum relative**) difference with respect to test 2A
	test 2A	400	4010401	outlet bulk values	47.16	0.2796	range of outlet values	27.25 to 90.93	0.03288 to 0.9999
		number of trajectories	total number of positions (*)		temperature (°C)	swelling degree		temperature (°C)	swelling degree	(*)

Table 3 :

 3 Summary of sensitivity tests regarding the influence of the number of positions per trajectory on the Eulerian-Lagrangian approach for solving the whole problem. All these results correspond to the domain outlet. All these tests consider 400 trajectories. the additional trajectory is considered, as well as all the release positions at the domain inlet (**) relative to the difference between corresponding outlet and inlet bulk values

	test 3B test 3C	10000 / 2 = 5000 10000 / 4 = 2500	2005401 1002901	difference (and relative difference**) with respect to test 3A	+4.5 10 -4 (+0.002 %) +4.6 10 -4 (+0.002 %)	10 -5 (+0.02 %) 2.1 10 -4 (+0.08 %)	maximum (and maximum relative**) difference with respect to test 3A
	test 3A = test 2A	10000	4010401	outlet bulk values	47.16	0.2796	range of outlet values	27.25 to 90.93	0.03288 to 0.9999
		number of positions	total number of positions (*)		temperature (°C)	swelling degree		temperature (°C)	swelling degree	(*)

Table 4 :

 4 Summary of main results obtained with the purely-Eulerian (test 1A) and the Eulerian-Lagrangian (test 2A) approaches for solving the whole problem. All these results correspond to the domain outlet. the additional trajectory is considered, as well as all the release positions at the domain inlet (**) relative to the difference between corresponding outlet and inlet bulk values

	test 1A

Acknowledgement

The research leading to these results has received funding from the European Community's Seventh Framework Program (FP7/2007-2013) under the grant agreement number FP7-222 654-DREAM.