A mapping tool using anisotropic unstructured meshes to study mixing

0 downloads 0 Views 961KB Size Report
Process. 16 (2) (2001b) 151) and Galaktionov et al. (Comput. Fluids 30 (3) (2001) 271) ... nificant to carefully control this stage of mixing in the pro- ... In 3D open ows, it requires the construction and the in- ..... between two reference meshes, the initial mesh and the de- .... tire ow mapping TP, the notion called the Lagrangian.
Chemical Engineering Science 59 (2004) 1459 – 1472

www.elsevier.com/locate/ces

A mapping tool using anisotropic unstructured meshes to study mixing in periodic %ows * Schall Yves Le Guer∗ , Eric  Laboratoire de Thermique Energ etique et Procedes (LaTEP), EA 1932, Universite de Pau et des Pays de l’Adour (UPPA), Avenue Jules Ferry, 64000 PAU, France Received 15 October 2002; received in revised form 3 July 2003; accepted 1 December 2003

Abstract The objective of the present work is to describe a new mapping tool using anisotropic unstructured meshes to study mixing within a spatially periodic pipe %ow. Instead of tracking the boundaries of elementary cell %ow domains as it was done in the original mapping method established by Kruijt et al. (A.I.Ch.E. J. 47 (5) (2001a) 1005; Int. Polym. Process. 16 (2) (2001b) 151) and Galaktionov et al. (Comput. Fluids 30 (3) (2001) 271), the deformation of elementary triangles (only three nodal points) between the inlet and exit pipe cross-sections is followed. It is however necessary to adapt the initial mesh according to criteria which takes into account the spatial stretching and folding of %uid elements. The method developed is applied to the twisted curved pipes (TCP) three-dimensional (3D) %ow. We show the evolution of concentration distributions along the TCP mixer for chaotic advection %ow regimes. This method allows the emphasis of isolated unmixed regions (KAM islands). The %exibility of the method allows also the possibility of studying multiple stirring protocols, thus contributing to a better comprehension of the physical phenomena involved in chaotic mixing. The method developed is also applicable to 2D temporally periodic %ows. ? 2004 Elsevier Ltd. All rights reserved. Keywords: Mapping method; Triangular meshes; Chaotic mixing; Mixing optimization; 3D pipe %ow

1. Introduction The problem of mixing highly viscous %uid is very common in many industrial processes including polymer processing, chemical and biochemical engineering, food engineering or material processing in civil engineering for examples. Mixing in laminar %ows is often required for typical applications that involve: • high viscosity %uids such as in polymerization reactions or the blending of polymers, • high shear %ows that provoke the segregation of shear thinning %uids, • shear sensitive biological materials. Currently, there is no general mixing theory leading to the individual study of mixing problems. During the elaboration of a material for example, the quality of the mixing will ∗

Corresponding author. E-mail addresses: [email protected] (Y. Le Guer), [email protected] (E. Schall). 0009-2509/$ - see front matter ? 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2003.12.028

condition the structure of the material and it is this structure that will determine the Bnal properties. It is thus very signiBcant to carefully control this stage of mixing in the process of development of a material. Recently Zumbrunnen and Inamdar (2001) have shown how sub-micron highly multi-layered polymer Blms could be formed by a continuous chaotic mixing %ow featuring stirring rods. For overly viscous %uids the only phenomenon which can facilitate the mixing is chaotic advection. For chaotic laminar %ows, one need eDcient tools for the analysis of their mixing. Indeed, the partially isolated mixing structures must be identiBed in order to destroy them. During the chaotic mixing process, one can observe very complex %ow structures (stretched and folded) even if the %ow is relatively simple [a two-dimensional (2D) %ow with temporal variation or a 3D %ow]. The periodic reorientation of the %uid elements is necessary to produce chaotic %ows (Ottino, 1989). For the Twisted Curved Pipe (TCP) %ow in question, this reorientation is obtained by the Dean cells found in adjacent curved pipe elements. We are also interested in the homogenization of stretching and folding distributions over a considered %uid domain. In order to optimize the mixing process, a fast

1460

Y. Le Guer, E. Schall / Chemical Engineering Science 59 (2004) 1459 – 1472

method to evaluate the mixing potentialities of various stirring protocols is needed. The presented method is a help tool for the optimization of mixing within periodic %ows. The interest of developing such a method appears very clear in the case of open 3D %ows, for which the research of the most eIective conBgurations for the mixing will be very diDcult to undertake by 3D Navier–Stokes simulations (due to multiple grids and signiBcant computing times). Moreover, this is very complex and expensive to accomplish from an experimental standpoint. In the case of the 2D %ows, whether they are treated numerically or experimentally, the optimization tests to determine the best stirring protocol are easier to do because the geometry is invariable. The testing of diIerent stirring protocols is possible using the same experimental apparatus. In 3D open %ows, it requires the construction and the introduction of new elements, both costly in time and money, so often impractical. Therefore, having powerful tools for mixing analysis becomes crucial. Knowledge of the velocity Beld is required before analyzing the mixing problem. For multi-phase %ow problems, two main classes of CFD methods have been developed to track the dynamics of the interface separating the two %uids. They are referred to as Front Capturing methods and Front Tracking methods. The most frequently employed Front Capturing methods are the Volume-of-Fluid (VOF) method (Scardovelli and Zaleski, 1999) and the Level-Set (LS) method (Osher and Sethian, 1988). These two Front Capturing methods are written for a one-%uid formalism, both phases are treated by a single set of conservation laws for the whole %ow domain, but the diIerences in the physical phase properties and surface tension are taken into account explicitly. For the Front Tracking method (Glimm et al., 2000) as well as in the Arbitrary Lagrangian–Eulerian (ALE) approach, among others, a periodic modiBcation of the computational grid is necessary. This provides sharp resolution in the deBnition of the interface. The Front Capturing method presents some diDculties in restoring the exact shape of the interface in situations where small scales are involved (typically for chaotic %ows). Moreover, Front Tracking methods have trouble handling changes in the interface topology (like droplet pinch-oI or droplet merger for example). All these highly eIective methods are diDcult to implement when studying mixing in chaotic %ows which have repetitive stretching and folding mechanisms. In addition, they also require enormous storage capacities because the data and computation times are very great. For all these reasons the mapping approach used here to study mixing in a 3D TCP mixer should present some advantages. Nevertheless, the mapping approach is limited to the %ow of similar %uids with no interfacial tension. Furthermore, it neither allows the prediction of heat and mass transfer nor the reactive process which takes place at the interface level. Beginning with the original idea of Spencer and Wiley (1951), the general principle of the mapping method was recently developed by the Materials Technology group

research team at Eindhoven University of Technology. The mapping method is based on the spatial discretization of a locally averaged concentration of %uid components within the mixture. A distribution matrix, which describes the changes in component concentration, is deBned for the %ow map. It allows the rapid prediction of mixing eDciency and the possibility of comparing multiple stirring protocols. The Eindhoven group developed the mapping method for: • idealized time periodic %ows: a 2D lid-driven cavity, a 2D journal bearing %ow (Kruijt et al., 2001a; Galaktionov et al., 2002c) and a 3D cavity with two opposite moving walls (Anderson et al., 1999; Galaktionov et al., 2001); • industrial spatially periodic %ows: a multi%ux static mixer, a corotating twin screw extruder (Kruijt et al., 2001b) and a Kenics static mixer (Galaktionov et al., 2002b). The objective of this work is to show the eDciency of the method for 3D spatially periodic pipe %ows, which are commonly used in chemical processing. The Anisotropic Unstructured Meshes Mapping Method (AUM3) developed here can yield the distribution of a passive scalar quantity (concentration or temperature) advected by the %ow after multiple periodic elements organized to follow a particular sequence of stirring protocols. The sequence of stirring protocols will be chosen considering geometric and %ow symmetries (Franjione and Ottino, 1992). Thus it is possible to optimize the distributive mixing for space or time periodic %ows. The powerful tools of grid reBnement, developed in the context of Finite Element Methods in CFD computations, can then be used. The method is exposed in detail and we explain how to rebuild an unstructured mesh, starting from the speciBc criteria related to the stretching properties of the %ow. The paper is organized as follows. In Section 2, the Eulerian velocity Beld as a basis for the generation of chaotic particle paths is presented. Then, the diIerent stirring protocols used to study the distributive mixing in the TCP sequence are emphasised. In Section 3, the implementation of the mapping method with the use of unstructured meshes for the TCP %ow is detailed. Finally in Section 4, results obtained from this method are presented and discussed while extended work is envisaged in the conclusion Section 5.

2. The 3D spatially periodic velocity eld 2.1. The Eulerian velocity 5eld The developed method will be outlined starting from a concrete example of a chaotic mixer for a 3D tubular %ow. This %ow produces very eIective mixing in the pipe cross-section. It has been subject to much study since 1989, numerical as well as experimental (Jones et al., 1989; Le Guer and Peerhossaini, 1991; Castelain et al., 2001). It falls into the class of continuous chaotic %ows such as the well-known partitioned pipe mixer (Khakhar et al., 1987;

Y. Le Guer, E. Schall / Chemical Engineering Science 59 (2004) 1459 – 1472

Kusch and Ottino, 1992). Its eIectiveness has been proved experimentally, as much for its capacity for increasing heat transfer (Mokrani et al., 1997; Acharya et al., 2001) as for its capacity for improving mass transfer (Boesinger, 2002). A similar design is also used in microchannels for which turbulent %ow cannot be considered due to the microscale of the pipe cross-section (Liu et al., 2000). In the Twisted Curved Pipe (TCP) %ow, the cross-sectional mixing is induced by centrifugal forces in the form of two recirculating cells that rotate in opposite directions (the secondary Dean cells). When diIerent curved pipe elements are Btted together, the relative position of the Dean cells in the adjacent elements changes, resulting in the possibility of obtaining chaotic particle trajectories and mixing enhancement. One may make the assumption of a spatially periodic %ow between the inlet and exit sections of the curved pipe element (or bend). The chaotic mixing is thus obtained by shifting the plane of curvature of the two adjacent curved pipe elements. When the shift angle is equal to 90◦ , the maximum perturbation between the adjacent streamlines is reached. Indeed, the destruction of the separatrix is responsible for chaotic advection. In the absence of streamline perturbation (where the shift angle is equal to zero) the stable and unstable manifolds emanating from the two hyperbolic points join linearly and form the separatrix between the Dean cells. Consequently, the particles are trapped forever between the two closed streamlines. As mentioned by Stroock et al. (2002), the TCP mixer is not appropriate for very low Reynolds numbers. This is because the transverse transport accomplished by the Dean cells is too small, so the %ow will resemble the one found in a straight tube. The %ow model chosen here is analytical. It corresponds to the solution of the Navier–Stokes equations established by Dean (1927) for a fully developed steady %ow in a curved pipe with a large radius of curvature, at a low Reynolds number. For a given geometry of the curved pipe, the motion of a passive %uid particle is determined by the advection equations, in dimensionless form: dx  = (1 − x2 − y2 )(4 − 5x2 − 23y2 dt 1152 + 8x2 y2 + x4 + 7y4 );

(1a)

dy  (1 − x2 − y2 )xy(3 − x2 − y2 ); = dt 192

(1b)

1 dz  = dt 4 (1 − x2 − y2 )

(1c)

or to a reduced form (if d z=dt = 0): dx = (4 − 5x2 − 23y2 + 8x2 y2 + x4 + 7y4 ); dz 1152

(2a)

1461

dy = xy(3 − x2 − y2 ); dz 192

(2b)

4 1 dt = ; d z  (1 − x2 − y2 )

(2c)

where x and y are the Cartesian coordinates of the advected particle in the plane of the pipe cross-section, z represents the longitudinal coordinate (i.e. the longitudinal angle of the transverse plane with respect to the center of the radius of curvature RC ). For these sets of equations:  = 64Dn2 ; =8 =

Dn2 ; Re

 = 8Re; 

where Re is the Reynolds number and Dn the Dean number: a Re = W0 :  W0 is the maximum velocity at the pipe’s center for a Poiseuille velocity proBle,  is the kinematic viscosity of the %uid: Dn2 = Re2

G 2 a7 a = : R 162 2 RC

P deBnes the pressure and G is the longitudinal pressure gradient which drives the %ow:   1 @P : G= R @z System (2) has the same form than system (1) with the independent variable being z instead of t. Therefore the motion deBned by the equations (2a) and (2b) could be analyzed like a 2D non-autonomous system. The %ow could be mapped from z to z pipe cross-sections. This is very useful for the implementation of the mapping method as the %ow considered is steady. The third equation (2c) gives the residence time of the followed particle according to z. The particle trajectories given by Eqs. (1) or (2) are toroidal and deBne closed streamlines. In the simulations, an instantaneous readjustment of the new streamlines at the bend junction will be considered. This assumption is not realistic but the main features of the chaotic %ow should be captured. For the hypothesis of spatially periodic %ow in a real device, it is suitable to consider that after a suDcient number of geometric periodic elements, the %ow itself should asymptotically tend to a periodic solution. In the experimental TCP reactor (Boesinger, 2002), a reactor with 80 bends has been studied. For some laminar %ow Reynolds numbers, complete mixing is achieved before the 10th bend element. The study of distributive mixing starts with the advective follow-up of particles within the %ow deBned by the dynamical system (2). The integration of system (2) is performed

1462

Y. Le Guer, E. Schall / Chemical Engineering Science 59 (2004) 1459 – 1472

with a four-stage accurate Runge–Kutta scheme. It is of importance to recall that the mixing is obtained within the pipe cross-section, where the secondary %ow operates %uid element reorientation. For an industrial device, when the steady %ow regime is reached, the mixing patterns observed in the pipe cross-section will never change. So, as the %ow considered is spatially periodic, it is possible to study its kinematics for diIerent stroboscopic cuts made at multiple integers of the spatial period (a bend for example). The analogue for a given spatial period SP of Eqs. (1) or (2) will be the map: Xn+1 = TP (Xn );

(3)

where the transformation TP deBnes the iterated position Xn+1 of the particle, initially located at the position Xn , and induced by the spatially periodic %ow. 2.2. The stirring protocols The %ow in each bend corresponds to the mapping TP . In this paper only four mappings TP (T1 –T4 ) are considered. This will be explained in Section 3. The reorientation of the %uid elements in the %ow is due to the  = 90◦ rotation R (or R−1 for  = −90◦ ) imposed at each junction for the boundaries of the next curved pipe element. Then a particular TCP sequence will be composed of a succession of TP R or TP R−1 mappings which corresponds to the two basic units. These mappings will be referenced as T and A protocols. Hereafter, T stands for a trigonometric rotation and A stands for an antitrigonometric one. Therefore as the Brst bend mapping (T1 –T4 ) is given, a stirring protocol can be deBne by a succession of mappings T and A. For a sequence of 100 bends (99 transitions), only the two series of protocols TA and TTTAA are here considered. One should take into account the bisector between the plane of curvature of two successive bends for a basic unit (two bends). For a TA stirring protocol sequence, the bisectrices are always oriented at the same angle with respect to the pipe cross-section (135◦ , in our conBguration). On the contrary, for a TTTAA stirring protocol sequence, the bisectors deBned for successive basic units follow a very complex sequence with angles equal to 45◦ , 135◦ , 225◦ or 315◦ . One could hope that the TTTAA stirring protocol, which does not preserve geometric symmetries, will lead to a best Bnal mixing state.

direction. The construction principle of the distribution matrix which constitutes the heart of the AUM3 is as follows. Once the inlet and exit sections, which limit the spatially periodic %ow domain, have been identiBed, it is possible to deBne the following steps: E1: the inlet section is covered with an optimized unstructured mesh (triangles). The procedure of optimization will be described below. E2: the three nodal points of each triangle are advected by the %ow (mapping TP ) for a periodic element (a bend). Then the exit deformed mesh is obtained. E3: each triangle of the deformed mesh is superimposed over the initial mesh and each matrix component ij is deBned as the area of overlap between the deformed triangle j and the one of the initial grid i, divided by the area of the triangle j,  deBnes the boundary of a triangle:  dA  | ∩ | : (4) ij = j n+1 i n dA j | n E4: the distribution matrix M is constructed from the calculated ij . E5: the concentration vector C1 is calculated according to a spatial period (a mapping) from the initial concentration distribution C0 and the distribution matrix M : C1 = MC0 :

E6: the procedure is reiterated for the same mapping or a new mapping (a sequence of stirring protocols) until the end of the mixer: Cn = M1 (M2 (M1 (M2 · · · (M1 C0 )))):

The 3D %ow inside the twisted coiled tubular mixer associated with the stirring protocols described above is considered. The mixing is studied from the mapping TP deBned between two reference meshes, the initial mesh and the deformed mesh located respectively at the inlet and at the exit of the considered spatial period (i.e. inlet and exit bend). These reference meshes are positioned in the plane of the curved pipe cross-section, perpendicular to the axial %ow

(6)

If we want to yield the smallest scales generated by the chaotic mixing, we must deBne a Bne mesh that implies a very large matrix M . The matrix M is sparse (only a few triangles of the deformed mesh intersect those of the initial mesh) but the multiplication of all the matrices Mi for a given sequence of stirring protocols will considerably reduce the sparsity of the resulting matrix. Instead of using the direct multiplication of all the distribution matrices Mi it may be better to compute the Bnal concentration vector Cn step by step: Cn = Mi Cn−1 :

3. Implementation of the mapping method

(5)

(7)

Kruijt et al. (2001a,b) and Galaktionov et al. (2001) have already successfully used this procedure for their mapping method. Here the matrix is calculated from an analytical velocity Beld, so it can be done rapidly, whereas it is a hard task if full Navier–Stokes equations are solved. Although the %ow is the same in each bend, four distributive matrices shall have to be builted since the initial unstructured mesh (triangles) is not exactly symmetric under the Ox and Oy directions, with respect to the one bend axis reference frame. In the plane of the pipe cross-section, the secondary %ow can

Y. Le Guer, E. Schall / Chemical Engineering Science 59 (2004) 1459 – 1472

N

1463

N

y

y E

x E W

W

x T1

T2

S N

(a)

S N

Fig. 2. Non-optimized coarse meshes for a quarter domain of the cross-section: (a) original mesh advected from the bend inlet and (b) deformed mesh at the bend exit.

x

W

x

E

W

y

E

y T3

S

T4

(b)

S

Fig. 1. Sketches representing the four %ow mappings studied. They are refered with respect to the %ow direction along the separatrix. T1 mapping: West to East, T2 mapping: North to South, T3 mapping: East to West, T4 mapping: South to North. The system of coordinates (x; y) is local.

induce a deformation of the initial mesh from west to east (T1 map), north to south (T2 map), east to west (T3 map), and south to north (T4 map) as shown in Fig. 1. In order to maintain the maximum of perturbation for the separatrix, an odd mapping must be followed by an even one and vice versa. Each case corresponds to the deBnition of a particular distributive matrix, M1 for T1 , M2 for T2 and so on. The original method developed by the authors quoted in the introduction is generally used for rectangular subdomains and the boundaries of each domain are composed of a great number of particles for which the displacement is registered between the periodic inlet and exit planes of the %ow. This method presents some drawbacks. Because of the stretching and folding it may be necessary to add new points (particles) in strongly elongated or convoluted regions in order to keep the continuity of the followed interface subdomain. Consequently, it is very time consuming. Franjione and Ottino (1987) have shown that this phenomenon can quickly become impossible to manage in numerical simulations for chaotic %ows. For the method using unstructured meshes, the number of particles tracked for a particular elementary domain (i.e. a triangle) does not depend on the deformation of the elementary domain (three in all cases). On the other hand, it is necessary to add triangles in the highly stretched and folded zones of the %ow. The main advantage of the presented method is to avoid tracking procedure for the creation of the mapping matrix. Furthermore, the problem of intersections between two

triangles, induced to build the mapping matrix, is less complicated than the problem of intersections between two polygons, which can present non-connected parts when the entire boundary of a subdomain is tracked. In the latter case, speciBc algorithms must be used. As an example, for a 63608 triangles mesh, the computation of the mapping matrix requires only 27 min of CPU time on a Linux Athlon (AMD) 1:5 GHz PC with 1 Go of RAM. The use of unstructured meshes also allows the possibility of complex shapes for the mixer design. On the other hand, the AUM3 requires the optimization of the initial mesh for the following reasons. A deformed triangle will not be able to account for the folding zones. So in order to avoid this skew, it is necessary to reBne the mesh in the strongly stretched and folded regions. One must also ensure that the deformed triangles will not overlap and that the sum of their areas corresponds to the total area of the exit cross-section. In order to avoid the problem of points which have zero velocity at the pipe wall, the mesh is constructed in the pipe cross-section for a radius equal to 0.999 instead of 1. 3.1. Mesh re5nement strategy The objective of the procedure is to build a mesh with a prescribed number of nodes which satisfy the following conditions: • the restitution of the exact exit area pipe cross-section, this condition cannot be met for each triangle because we are in three-dimensions, only volumes are conserved during the %ow; • the non-overlapping of deformed triangles in the exit pipe cross-section. The mesh generation procedure is decomposed as follows: • an initial coarse mesh for a quadrant of the inlet pipe cross-section is created. This quadrant corresponds to the region where the maximum deformation of the mesh was observed, it is reproduced in Fig. 2(a) for the T1 map;

1464

Y. Le Guer, E. Schall / Chemical Engineering Science 59 (2004) 1459 – 1472

(a)

where nlt is the number of edges associated with the nodal point, li being the initial edge length and lf being the stretched edge length; • the Lagrangian stretching which is obtained from local stretching. From the Eulerian velocity Beld, the local stretching #ij of a %uid element can be evaluated by the square root of the highest value of the right Cauchy– Green strain tensor Eij (Ottino, 1989; Voth et al., 2002): @Tk @Tk Eij = ; (10) @xi @xj

(b)

Fig. 3. Zones which produce triangles overlapping for the mesh represented in Fig. 2: (a) on original mesh and (b) on deformed mesh (T1 mapping).

• this mesh is projected from the map on to the outlet, its distortion is assessed for a chosen sensor; • the inlet quadrant is then re-meshed; • this procedure is repeated until the Bnal mesh satisBes the non-overlapping of the deformed triangles; • the full initial mesh is created by re%ections (this is necessary for the nodes coincidence since the mesh is anisotropic); • the full mesh is projected on to the outlet pipe cross-section for each of the four mappings considered (T1 –T4 ) in order to verify the suitability of the inlet mesh; • if the deformed mesh produces overlapping(s) (this is rare), a slight modiBcation is operated until entire satisfaction is obtained. It is the choice of the four maps T1 –T4 which imposes the 90◦ symmetry, and with this construction a single full inlet mesh is suitable for the four maps. The Bnal mixing results are always projected on to this inlet mesh, thus ensuring the possibility of studying all the stirring protocol sequences. The Brst mapping %ow, which corresponds to the deformation of the inlet mesh to the exit mesh, is shown in Fig. 2(b). The %ow domains, which produce overlapping triangles, are visualized in Figs. 3(a) for the initial mesh and 3(b) for the deformed exit mesh. From this Brst mapping, a particular scalar quantity is extracted, deBned per each node of the mesh: the sensor. We have tested three sensors based on: • the surface deformation Cnt of triangles associated with a nodal point: Cnt =

nt Sf − Si 1  ; Si nt

(8)

1

where nt is the number of triangles associated with the nodal point, Si being the initial triangle and Sf being the deformed triangle; • the deformation of edges Clt associated with a nodal point: nlt lf − l i 1  ; Clt = nlt li 1

(9)

where T is the %ow mapping from one cross-section to another and the summation rule is implied over the repeated index k =1; 2. Typically one hundred cross-section is considered to map the %ow in each bend. At each point in the cross section, #ij gives a scalar which indicates the magnitude of the stretching. In order to obtain a quantity, which characterizes the stretching experienced by a particle in the plane of the pipe cross-section for the entire %ow mapping TP , the notion called the Lagrangian stretching is introduced: #ij

lag

ds =

= 

N 1  #ij ds N

with

1

d x2 + dy2 + d z 2 ;

(11)

where N is the number of the cross-sections considered for calculation trajectories in a bend and ds is the travelling distance of a particle between two consecutive cross-sections. These criteria are reported in Fig. 4(a) for Cnt , 4(b) for Clt and 4(c) for #ij lag . The local stretching #ij is represented in Fig. 4(d). For a given sensor, the absolute value of the Hessian is used as anisotropic metric. When the metrics are obtained, they are used in a mesh regenerator, rebuilding a new mesh with about the same prescribed number of vertex (nodes). This is done by using Delaunay reconnection and vertex addition. This task was carried out with the aid of existing software, the Bidimensional Anisotropic Mesh Generator (BAMG), created by George and Hecht (1997). The results of the reBnement procedure consist in designing an unstructured mesh with respect to the above-deBned criteria. The Bnal optimized meshes are shown in Fig. 5(a) for the original mesh and 5(b) for the deformed mesh. It is of some interest to observe the improved restitution of the pipe boundary (Fig. 5(b)), in comparison to the deformed mesh obtained from the non optimized case (Fig. 2(b)). This Bnal optimized mesh is adapted to the four %ow mappings considered in T1 –T4 . The Bnal reconstructed mesh for the whole %ow domain of the pipe cross-section is shown in Fig. 6(a) for the original mesh and 6(b) for the deformed mesh (T1 mapping). As the unstructured mesh is reBned by conserving approximately the same number of nodal points, some triangles have increased in area. This can induce a worse

Y. Le Guer, E. Schall / Chemical Engineering Science 59 (2004) 1459 – 1472

(a)

1465

(b)

1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

(c)

(d)

Fig. 4. Plots showing the distribution of the scalar used to reBne the mesh for the T1 map: (a) area deformation of triangles sensor, (b) edges deformation of triangles sensor, (c) Lagrangian stretching sensor and (d) local stretching deBne for the whole cross-section.

(a)

(b)

Fig. 5. Coarse optimized unstructured meshes for a quarter domain of the cross-section: (a) original mesh advected from the bend inlet and (b) deformed mesh at the bend exit (T1 mapping).

deBnition of the striation thickness found in these zones during the mixing process, which in turn will be propagated over other %ow domains at the end of the TCP sequence. The number of iterations in the procedure needed to build the “ideal” mesh characterized the “quality” of the sensor used. The sensor based on the Lagrangian stretching of %uid particles gives the best result for the construction of the mesh.

Fig. 6. Recomposed coarse mesh (4003 nodes) for the whole pipe cross-section: (a) original mesh advected from the bend inlet and (b) deformed mesh at the bend exit (T1 mapping).

4. Results and discussion The present results are given for a parameter =150 which corresponds to the maximum admissible Dean number for the analytical Dean’s solution. Fig. 7 clearly shows the mixing of the two species along a TCP sequence of 84 elements for a TA stirring protocol. Results for TA stirring protocol are Brst presented because it produces partially mixed structures. The mesh is the Bnest one we have considered and

1466

Y. Le Guer, E. Schall / Chemical Engineering Science 59 (2004) 1459 – 1472

Fig. 7. Mixture plots for the stirring protocol TA. At the inlet of the mixer, the upper half-cross-section corresponds to concentration C = 0 and the lower one to C = 1. From the top left to the bottom right, the mixture plots corresponds to the exit of: bends 1, 4, 6, 12, 24, 36, 52, 68 and 84.

contains 63608 triangles. The Brst plot corresponds to the exit of the Brst bend where no mixing is produced. The distribution of species is then the same as the one imposed at the entrance of the mixer. Indeed, each Dean cell carries its own %uid specie, and mixing occurs only by molecular diffusion along the interface. For all the mixture plots shown, mixing is always produced at the interface between the two species, it is characterized by the apparition of a gray color C ≈ 0:5. As the %uid travels along the TCP mixer, more and more interface length is created, resulting in the increasing presence of the color gray on the plot. As the TA stirring protocol is highly periodic, similar recurrent patterns appear with more and more details inside the structures. These structures resemble those observed in time periodic %ows by Galaktionov et al. (2002a) for the case of the 2D lid-driven cavity, or else by Rothstein et al. (1999) for a forced magnetic chaotic %ow. The two black and white structures correspond to the periodic unmixed regions. Their existence persists during the whole process. These structures coincide with the KAM regular islands immersed in the chaotic sea shown in the particles distribution plot reported in Fig. 8 for the same %ow. As a direct consequence, a global mixing for the whole cross-section is never reached. From the

Fig. 8. Particle distributions (49900 points) for the TA stirring protocol: (a) at the inlet of the Brst bend and (b) at the end of the 100th bend.

observation of the mixture plots (only a few of them are reproduced in Fig. 7) and the cut-line concentration plots, we can deduce that outside the two unmixed regions, the striations (or lamellas) start to be “eaten” by numerical diIusion precisely between the 22th and 24th bends. For these striations a rapid decrease in the concentration is observed from the value C = 1 (respectively an increase in concentration from the value C = 0). Inside the unmixed regions, the concentration will normally keep its original value (C = 0

Y. Le Guer, E. Schall / Chemical Engineering Science 59 (2004) 1459 – 1472

1467

Table 1 Inlet and outlet (deformed) mesh characteristics (the pipe radius in the cross-section is equal to 1) Mesh

1

2

3

4

Number of nodes Number of triangles Minimum edge length

4003 7812 0:922 × 10−2 0:110 × 10−2 0:946 × 10−1 0.294 0:318 × 10−1 0:533 × 10−1 0:212 × 10−3 0:134 × 10−2 3.932 119.472

8113 15956 0:562 × 10−2 0:759 × 10−3 0:727 × 10−1 0.228 0:222 × 10−1 0:371 × 10−1 0:111 × 10−3 0:677 × 10−3 4.91 110.095

16257 32144 0:404 × 10−2 0:594 × 10−3 0:575 × 10−1 0.15 0:155 × 10−1 0:258 × 10−1 0:563 × 10−4 0:331 × 10−3 4.923 142.086

32043 63608 0:318 × 10−2 0:507 × 10−3 0:444 × 10−1 0.105 0:108 × 10−1 0:180 × 10−1 0:240 × 10−3 0:170 × 10−3 3.932 119.472

Maximum edge length Mean edge length Variance on edge lengths Maximum aspect ratio (Max. edge length/Min. edge length)

Inlet Outlet Inlet Outlet Inlet Outlet Inlet Outlet Inlet (I) Outlet (O)

or 1) since the KAM curves act as a barrier for the surrounding %uid elements. The specie can only penetrate them by molecular diIusion. However, along the TCP mixer, the concentration gradients inside the unmixed regions decrease slowly due to numerical diIusion. This phenomenon will be explained later. To evaluate the in%uence of the mesh-sizes on the mixing quality, four mesh sizes have been tested. Table 1 illustrates the diIerent characteristics of the inlet and deformed meshes. All the edge lengths must be compared to the pipe diameter length which is equal to 2. The dimensions considered for the inlet meshes are important because they determine the accuracy of the concentration solution. For the Bnest mesh (mesh 4), the minimum resolution, based on the triangle edge lengths should make it possible to visualize between 100 and 200 striations along a single pipe diameter. For all the mesh sizes considered, the mesh stretching is characterized by a maximum aspect ratio at the inlet and outlet pipe cross-sections. The maximum aspect ratio is deBned as the maximum ratio of triangle edge lengths found for all the mesh triangles. The I and O ratios (indicated in Table 1) are not necessarily obtained for the same triangle of the meshes. These aspect ratios show clearly the anisotropy of the meshes which is very large at the outlet pipe cross-section where %uid deformations are important. The mixture plots at the exit of the 36th bend for the TA stirring protocol are shown in Figs. 9(a)–(d) for the diIerent mesh-sizes. As the number of triangles increases, details within the mixing structures appear more and more clearly. More and more striations are also increasingly visible on the concentration plots. A cut-line done for x = 0 (vertical line) on concentration plots (Fig. 9) illustrates clearly the smoothing eIect on striations induced by the mapping method (Fig. 10). More than 23 striations can be counted for the Bnest mesh, whereas only seven can be observed for the coarsest one. As the line x=0 crosses the two regular islands, we can note that the Bnest mesh is able to restore values of concentration close to zero or one inside these unmixed regions. The blurred structures observed in Fig. 9, when a

small number of triangles is involved, are a direct consequence of the diIusive character of the mapping method. In fact, after one step of the mapping method, all the triangles will contain a portion of the black color C = 0 and another one of the white color C = 1. The resulting color of a triangle will correspond to an average concentration. This is precisely the mixing process to describe. However if the mesh is too large, a rapid average in concentration distribution will be reached, corresponding to a non-physical mechanism. This phenomenon is well explained by Galaktionov et al. (2001) and Kruijt et al. (2001a). Consequently, the mesh-size should be adapted to the Bnest structure of the %ow that must be highlighted. The restitution of a high concentration gradient cannot be obtained precisely on a small scale if the triangles are too large. As stated previously, this artiBcial diIusion is also observed within unmixed regions, where the concentration should always be zero or one for all the unmixed regions. This is not the case, especially if the mesh is too large, the homogenization of concentration is produced within the unmixed regions (Figs. 9(a), (b) and 10). An indicator of the mixing quality, according to the mesh-size considered, is given by the relative error of concentration variation *: *% =

UCth − UCAUM3 UCth

(12)

with UCth = 1 and UCAUM3 = Cmax − Cmin . This indicator is interesting to evaluate when large poor-mixing regions are present in mixture plots (Fig. 11). The relative error of concentration variation evolves from 32:4% for mesh 1, to 0:66% for mesh 4, following an exponential decreasing law: *% = 55:208 exp (−6:8738 × 10−5 Nt );

(13)

where Nt the number of triangles for a given mesh. As for turbulence, with chaotic mixing the scalar dissipation operates from large to small scales. A quantity that well characterizes the scalar dissipation due to mixing is the Discrete Intensity of Segregation or coarse-grain density

1468

Y. Le Guer, E. Schall / Chemical Engineering Science 59 (2004) 1459 – 1472

Fig. 9. Mixture plots at the exit of the 36th bend for the same stirring protocol TA and diIerent mesh-sizes. Upper row: left 7812 triangles, right 15956 triangles. Lower row: left 32144 triangles and right 63608 triangles.

1 Mesh 4 Mesh 3 Mesh 2 Mesh 1

Concentration

0.8

0.6

0.4

0.2

0 -1

-0.5

0

0.5

1

y location on pipe diameter (x = 0) Fig. 10. ProBles of concentration along the line x = 0 (vertical line) as function of the mesh-sizes for the diIerent mixture plots shown in Fig. 9.

Relative Error on Concentration Variation (%)

Y. Le Guer, E. Schall / Chemical Engineering Science 59 (2004) 1459 – 1472

1469

10

1

0.1

0

10000

20000

30000

40000

50000

60000

70000

Number of triangles Fig. 11. Relative error on concentration variation *%=(UCth −UCAUM3 )=UCth as function of the mesh-sizes for the diIerent mixture plots shown in Fig. 9.

1 Mesh 4 Mesh 3 Mesh 2 Mesh 1

Discrete Intensity of Segregation

0.8

0.6

0.4

0.2

0

0

10

20

30

40

50

60

70

80

90

100

Bend index Fig. 12. Evolution of the intensity of segregation for a sequence of 100 bends and the stirring protocol TA as function of the mesh-sizes.

(Galaktionov et al., 2001; Kruijt et al., 2001a): N

IS =

1 1 V 2 ai (Ci − C) V V C(1 − C) A i=1

(14)

with N

1 Ci a i CV = A i=1

and

A=

N 

ai ;

i=1

where Ci is the concentration associated to the elementary triangle ai . A is the pipe cross-section area and CV is the mean concentration over the entire pipe cross-section. IS = 1

represents an unmixed state (complete segregation) whereas IS = 0 characterizes the ideal mixture (uniform concentration). Fig. 12 shows the evolution of the Discrete Intensity of Segregation for a TCP sequence of 100 bends for a given mesh size. A coarse mesh produces a faster drop of the mixing criterion IS , so one can consider that the %ow mixes better, whereas the artifact is related to the resolution of the mesh. It can be observed that IS never reaches the value of zero. This indicates that the presence of regular islands prevent perfect mixing. The %attening of the curves of IS versus the number of bends is also due to the presence of islands generated by the TA stirring protocol. The error made when

1470

Y. Le Guer, E. Schall / Chemical Engineering Science 59 (2004) 1459 – 1472

Fig. 13. Comparison of mixing plots for two diIerent stirring protocols. Upper row: TA stirring protocol, exit of bends 50 and 100. Lower row: TTTAA stirring protocol, exit of bends 50 and 100.

a coarse mesh is considered is important and increase when the %ow produces numerous thin striations. Two physical behaviors from these curves can be deduced. Firstly, from the Brst to the 20th bend concentration gradients are created, IS falls abruptly. Secondly, starting with 20th bend IS decreases very slowly, which corresponds to the smoothing of the concentration gradients by the numerical diIusion (at the triangles mesh scale). This phenomenon mimics the natural diIusion that occurs in experiments and one could imagine exploiting this fact in order to introduce controlled diIusion by the adaptation of the initial size of the mesh. Nevertheless, additional work is necessary to characterize this diffusion within the diIerent regions of the %ow, with respect to the sizes of the associated mesh triangles. In Fig. 13, the mixture’s plots for the TTTAA and TA stirring protocols are compared for a 100 TCP sequence at the end of the 50th and 100th bend. The TTTAA stirring protocol mixes better than the TA one, and does not generate unmixed KAM regions. The eIect of mixing on concentration variations can be shown in the results reported in Table 2. For the TA stirring protocol the high concentration variation is due to the existence of two unmixed KAM islands. Except for these, the mixing is relatively good. For the TTTAA stirring proto-

Table 2 Concentration variations for the four mixtures plots of Fig. 13 Stirring protocol

Bend exit

Cmax

Cmin

UC

TTTAA TTTAA TA TA

50 100 50 100

0.594 0.502 0.967 0.908

0.418 0.498 0.019 0.092

0.176 0.004 0.948 0.816

col, the concentration variation is already small at the end of the 50th bend. Fig. 14 shows the evolution of the Discrete Intensity of Segregation IS for a sequence of 100 bends and for the two stirring protocols studied. The fastest decrease of IS is obtained from the TTTAA stirring protocol. An exponential decay rate is clearly shown for this stirring protocol from the 12th bend. The saturation observed from the 20th bend on the TA curve is due to the presence of KAM islands as previously mentioned. It seems that the best mixing is associated with the low periodicity observed for the recurrent mixing patterns (20 bends for the TTTAA stirring protocol instead of two bends for the TA stirring protocol).

Y. Le Guer, E. Schall / Chemical Engineering Science 59 (2004) 1459 – 1472

1471

Discrete Intensity of Segregation

1

0.01

0.0001

TTTAA stirring protocol TA stirring protocol Exponential regression on TTTAA curve

1e-06

0

10

20

30

40

50

60

70

80

90

100

Bend index Fig. 14. Evolution of the intensity of segregation for a sequence of 100 bends and the stirring protocols TA and TTTAA for the Bnest mesh-size.

5. Conclusion The main goal of this work is to present a new development of the mapping method that uses anisotropic unstructured triangular meshes and its possible applications for chaotic mixing studies. No boundaries tracking is needed to follow the deformation of the %uid elements of the %ow. However, a mesh reBnement strategy is necessary to follow the stretching and folding of the %uid. A sensor based on Lagrangian stretching gives the best result for the unstructured mesh adaptation. The AUM3, even with the use of a coarse mesh, is able to indicate which one of the stirring protocols give the best mixing quality. The %exibility and the rapidity of the AUM3 authorizes a comprehensive analysis of mixing mechanisms for diIerent stirring protocols, especially those met for chaotic 3D periodic %ows. The creation and the control of isolated mixing regions (Bresler et al., 1997; Muzzio et al., 2000) or the destruction of KAM torus by exploiting the diIerent symmetries of the %ow (Franjione and Ottino, 1992) are two examples. Here only two stirring protocols were studied, the objective of the work was to demonstrate the feasibility of the new method. Extended work concerning mixing optimizations with the combination of the use of AUM3 and Genetic Algorithms is currently in progress. We have also deBned another mixing measure, “the average scalar dissipation” which is more suitable because it takes into account the two phases of the mixing process: the Brst one being when the scalar gradients increase due to eDcient stirring, and a second one which smoothes the gradients of the scalar when the scale of striations becomes suDciently small. Self-similarity of persistent spatial patterns, and the rate of decay of length scales are also investigated for diIerent chaotic %ows in the TCP. All these studies were made possible by the use of the AUM3 and this will be developed in a future publication.

The method could also quite easily incorporate the Extended Mapping Method (Galaktionov et al., 2002a,b) for calculating sub-cell mixing, since the deformation gradient is already known for each cell, and is constant within the cell. Acknowledgements The authors particularly thank G. Carte and F. Olier for their assistance in the construction of the intersections algorithm and Dr. A. Dervieux from INRIA Sophia-Antipolis for his valuable comments. We would also like to thank the reviewers of the manuscript for their valuable suggestions. References Acharya, N., Sen, M., Chang, H.C., 2001. Analysis of heat transfer enhancement in coiled-tube heat exchangers. International Journal of Heat and Mass Transfer 44 (17), 3189–3199. Anderson, P.D., Galaktionov, O.S., Peters, G.W.M., van de Vosse, F.N., Meijer, H.E.H., 1999. Analysis of mixing in three-dimensional time-periodic cavity %ows. Journal of Fluid Mechanics 386, 149–166. Boesinger, C., 2002. M*elange diIusif et r*eactif dans des r*eacteurs tubulaires aW trajectoires complexes, Ph.D. Thesis, Universit*e de Pau et des Pays de l’Adour, France. Bresler, L., Shinbrot, T., Metcalfe, G., Ottino, J.M., 1997. Isolated mixing region: origin, robustness and control. Chemical Engineering Science 52 (10), 1623–1636. Castelain, C., Mokrani, A., Le Guer, Y., Peerhossaini, H., 2001. Experimental study of chaotic advection regime in a twisted duct %ow. European Journal of Mechanics B/Fluids 20 (2), 205–232. Dean, W.R., 1927. Note on the motion of %uid in a curved pipe. Philosophical Magazine 4, 208–223. Franjione, J.G., Ottino, J.M., 1987. Feasibility of numerical tracking of material lines and surfaces in chaotic %ows. Physics of Fluids 30 (12), 3641–3643.

1472

Y. Le Guer, E. Schall / Chemical Engineering Science 59 (2004) 1459 – 1472

Franjione, J.G., Ottino, J.M., 1992. Symmetry concepts for the geometric analysis of mixing %ows. Philosophical Transactions of Royal Society of London A 338, 301–323. Galaktionov, O.S., Anderson, P.D., Kruijt, P.G.M., Peters, G.W.M., Meijer, H.E.H., 2001. A mapping approach for 3D distributive mixing analysis. Computers and Fluids 30 (3), 271–289. Galaktionov, O.S., Anderson, P.D., Peters, G.W.M., Tucker III, C.L., 2002a. A global multi-scale simulation of laminar %uid mixing: the extended mapping. International Journal of Multiphase Flow 28 (3), 497–523. Galaktionov, O.S., Anderson, P.D., Peters, G.W.M., Meijer, H.E.H., 2002b. Morphology development in Kenics static mixers (application to the extended mapping method). Canadian Journal of Chemical Engineering 80 (6), 604–613. Galaktionov, O.S., Anderson, P.D., Peters, G.W.M., 2002c. Structure development during chaotic mixing in the journal bearing %ow. Physics of Fluids 14 (9), 3009–3017. Glimm, J., Grove, J.W., Li, X.L., Tan, D.C., 2000. Robust computational algorithms for dynamic interface tracking in three-dimensions. SIAM Journal on ScientiBc Computing 21 (6), 2240–2256. George, P.L., Hecht, F., 1997. Non-isotropic grids. In: Thompson, J.F., Soni, B.K., Weatherhill, N.P. (Eds.), Handbook of Grid Generation. CRC Press, Boca Raton, FL (Chapter 20). Jones, S.W., Thomas, O.M., Aref, H., 1989. Chaotic advection by laminar %ow in a twisted pipe. Journal of Fluid Mechanics 209, 335–357. Khakhar, D.V., Franjione, J.G., Ottino, J.M., 1987. A case study of chaotic mixing in deterministic %ows: the partitioned-pipe mixer. Chemical Engineering Science 42, 2909–2926. Kruijt, P.G.M., Galaktionov, O.S., Anderson, P.D., Peters, G.W.M., Meijer, H.E.H., 2001a. Analyzing mixing in periodic %ows by distribution matrices: mapping method. A.I.Ch.E. Journal 47 (5), 1005–1015. Kruijt, P.G.M., Galaktionov, O.S., Peters, G.W.M., Meijer, H.E.H., 2001b. The mapping method for mixing optimization. Part I: the multi%ux static mixer. Part II: transport in corotating twin

screw extruder. International Polymer Processing 16 (2), 151–171. Kusch, H.A., Ottino, J.M., 1992. Experiments on mixing in continuous chaotic %ows. Journal of Fluid Mechanics 236, 319–348. Le Guer, Y., Peerhossaini, H., 1991. Order breaking in Dean %ow. Physics of Fluids A3 (5), 1029–1032. Liu, R.H., Stremler, M.A., Sharp, K.V., Olsen, M.G., Santiago, J.G., Adrian, R.J., Aref, H., Beebe, D.J., 2000. Passive mixing in a three-dimensional serpentine microchannel. Journal of Microelectromechanical Systems 9 (2), 190–197. Mokrani, A., Castelain, C., Peerhossaini, H., 1997. The eIects of chaotic advection on heat transfer. International Journal of Heat and Mass Transfer 40 (13), 3089–3104. Muzzio, F.J., Alvarez, M.M., Cerbelli, S., Giona, M., Adrover, A., 2000. The intermaterial area density generated by time and spatial periodic 2D chaotic %ows. Chemical Engineering Journal 55 (8), 1497–1508. Ottino, J.M., 1989. The Kinematics of Mixing: Stretching, Chaos and Transport. Cambridge University Press, Cambridge. Osher, S., Sethian, J.A., 1988. Front propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations. Journal of Computational Physics 79, 12–49. Rothstein, D., Henry, E., Gollub, J.P., 1999. Persistent patterns in transient chaotic %uid mixing. Nature 401, 770–772. Scardovelli, R., Zaleski, S., 1999. Direct numerical simulation of free-surface and interfacial %ow. Annual Review of Fluid Mechanics 31, 567–603. Spencer, R.S., Wiley, R.H., 1951. The mixing of very viscous liquids. Journal of Colloid Science 6, 133–145. Stroock, A.D., Dertinger, S.K.W., Ajdari, A., Mezic, I., Stone, H.A., Whitesides, G.M., 2002. Chaotic mixer for microchannels. Science 295, 647–651. Voth, G.A., Haller, G., Gollub, J.P., 2002. Experimental measurements of stretching Belds in %uid mixing. Physics Review Letters 88, 254501. Zumbrunnen, D.A., Inamdar, S., 2001. Novel sub-micron highly multi-layered polymer Blms formed by continuous %ow chaotic mixing. Chemical Engineering Science 56 (12), 3893–3897.