A new numerical solution scheme for tracking sharp moving fronts

A new numerical solution scheme for tracking sharp moving fronts

19th European Symposium on Computer Aided Process Engineering – ESCAPE19 J. JeĪowski and J. Thullie (Editors) © 2009 Elsevier B.V. All rights reserved...

126KB Sizes 7 Downloads 87 Views

19th European Symposium on Computer Aided Process Engineering – ESCAPE19 J. JeĪowski and J. Thullie (Editors) © 2009 Elsevier B.V. All rights reserved.

907

A new numerical solution scheme for tracking sharp moving fronts Ala Eldin Bouaswaig a and Sebastian Engell a a

Process Dynamics and Operations Group , Technische Universität Dortmund, 44221 Dortmund, Germany, E-mail: [email protected]

Abstract The numerical solution of a hyperbolic or a convection dominated parabolic partial differential equation is challenging due to the large local gradients that are present in the solution. In this paper, a novel approach that is based on combining the high order weighted essentially non-oscillatory (WENO) scheme with a static moving grid method is presented. The proposed algorithm is tested on two case studies and enhancements in the performance are observed when compared with the conventional WENO scheme on a uniform grid making it a promising alternative when dealing with problems of this nature. Keywords: WENO scheme, Adaptive gridding, Burgers' equation, Population balance equation, Emulsion polymerization

1. Introduction The numerical solution of hyperbolic or convection-dominated parabolic partial differential equations (PDE) is challenging due to the large local gradients in the solution profile. If a uniform fixed grid is used for the discretization, it has to be sufficiently fine to be able to track the steep fronts that are associated with these large gradients. Such a distribution of the grid nodes over the discretization domain is, however, likely to be inefficient since in smooth regions less nodes are actually required and only in the vicinity of the steep front does the grid have to be dense. This stimulates the need to use a nonuniform grid and to adapt it continuously to track the evolving steep front in the solution. Grid adaptation techniques can be categorized into two groups, static and dynamic [1]. In the dynamic version, the grid is continuously moved with the solution and both subproblems, the grid adaptation and the solution of the original PDE, are tackled simultaneously. In contrast, when static grid adaptation is used, the grid is either refined (local grid refinement method) [2] or the grid cells are moved (moving mesh method) at discrete points of time [3]. Grid adaptation must be combined with an appropriate numerical method for discretizing the PDE. Recently high order discretization methods have gained increasing interest [1,2]. For example, the weighted essentially non-oscillatory (WENO) scheme has emerged as a promising alternative for handling hyperbolic PDEs on a uniform grid. Combining the WENO scheme with grid adaptation merges the benefits of having a high resolution scheme with those of adapting the grid. This has been applied to dynamic grid adaptation by Lim et al. [1] and static grid adaptation with local grid refinement by Smit et al. [2]. However, a combination of the WENO scheme with a moving mesh static grid adaptation technique has not yet been investigated. The aim of this contribution is to describe how the WENO scheme can be applied on a nonuniform grid that is composed of a fixed number of nodes and that moves at discrete

A.E. Bouaswaig and S. Engell

908

points of time. To illustrate the benefits that are gained by such an implementation, the quality of the results obtained by using this novel approach is evaluated by testing it on the viscid Burger’s equation and on a challenging hyperbolic PDE, the population balance equation (PBE) that describes the growth of polymer particles in emulsion polymerization. The remainder of this paper is divided as follows: In Section 2, the mathematical formulation of the WENO scheme on a nonuniform grid is presented and the static grid adaptation algorithm is briefly addressed. In Section 3 the performance of the method is investigated on the two case studies mentioned above, and finally, Section 4 gives a summary and conclusions.

2. The numerical approach The WENO scheme on a non-uniform grid The WENO scheme on a nonuniform grid is presented here based on the paper of Shu [4], and the used notation is similar to that in [4]. Starting from the discrete form of a general one dimensional homogeneous hyperbolic conservation law of a quantity ( u ):

dui + dt



i+

1 2

− fˆ

i−

1 2

= 0,

Δxi

cell averages of the function x

1 fi ≡ Δxi

i+

x

i = 1,2,..., N ,

(1)

f (x) are defined by Eq. (2):

1 2

³ f (ϕ )dϕ.

(2)

1 i− 2

For a WENO scheme of order k in smooth regions, k stencils are defined and on each th

of them a k order approximation of the flux is derived as follows: k −1

f ( r1) = ¦ cr( ,kj),i f i − r + j , i+

r = 0,...., k − 1,

(3)

j =0

2

k k § · ¨ k ¦ω =0,ω ≠m ∏ p =0, p ≠ m,ω ( xi + 1 − xi −r + p − 1 ) ¸ 2 2 ¸ =¨ ¦ Δxi −r + j . k ¨ m= j +1 ¸ ( x − x ) ∏ω =0,ω ≠m i−r +m− 1 i −r +ω − 1 ¸ ¨ 2 2 © ¹

cr( k, j),i

(4)

The desired numerical flux in Eq. (1) is now constructed by applying Eq. (5):



i+

k −1

1 2

= ¦ ω ( q1) f ( q1) ,

where

q =0

i+

2

i+

2

ω ( q ) are nonlinear weights that depend on the smoothness of the solution.

(5)

A New Numerical Solution Scheme for Tracking Sharp Moving Fronts

909

Grid adaptation The static grid adaptation technique is performed by firstly adapting the grid to suite the initial profile. The profile is then interpolated to obtain new initial conditions on the new nonuniform grid and the PDEs are discretized. Finally, the resulting set of ODEs is integrated for a certain time duration while keeping the grid fixed and the procedure is then repeated until t final is reached [3]. The location of the grid nodes is determined so as to equally distribute a monitor function m( x ) , which in this work is taken to be the arc-length of the solution defined by Eq. (6):

∂f ( x) m( x ) = 1 + ∂x

2

(6) 2

Thus, to perform the grid adaptation step, all points on the solution profile are connected with each other and the total length of the resulting polygon is divided into equal parts. The new position of the discretization nodes is consequently determined by projecting the points that result from this equidistribution on the x-axis.

3. Numerical experiments All simulation runs shown in this section are performed on a PC with an Intel® Core™2 Duo Desktop Processor and with 2.0 GB of RAM. gPROMS® is used as a DAE solver and for implementing grid adaptation, it is connected to a dynamic link library programmed in Compaq Visual FORTRAN. Additionally, the l2 error norm used for analyzing the results is defined as follows:

§ N · l 2 = ¨ ¦ (u i − uˆ i ) 2 ¸ , © i =1 ¹ where ui stands for the numerical solution and

(7)

uˆi is the corresponding analytical

solution. The Viscid Burgers’ equation (BE) The BE, given by Eq. (8), is frequently used to test grid adaptation [1,5]:

∂ 2u ∂ (0.5u 2 ) ∂u +μ 2 . =− ∂x ∂t ∂x

(8)

Fig. 2 compares the results of WENO35 with adaptive gridding (WENO35-AG) with the analytical solution of the BE taken from [5]. The corresponding grid is depicted in Fig. 4 at discrete points of time. The performance of the WENO35-AG is outstanding and the numerical solution can hardly be distinguished from the analytical one. This observation is also reflected in the results of the error norm shown in Fig. 1. It is obvious from the results that, using grid adaptation enhances the performance of the WENO scheme. For example, it follows from Fig. 1 that for the classical WENO35 to produce results as accurate as WENO35-AG on a grid made of 100 nodes, the uniform grid must consist of around 250 nodes.

A.E. Bouaswaig and S. Engell

910 0

10

1

Classical WENO35 WENO35-AG ] -[ -1 10 m r o N l2

0.8 0.6 x( u

0.4 0.2

-2

10

50

100

150 200 250 # Grid nodes [-]

0 0

300

0.2

0.4

0.6

0.8

1

x [-]

Figure 1. Influence of grid refinement on the l2 error norm for WENO35 and WENO35-AG scheme at t = 1.0 (-) when applied to the BE.

Figure 2. Comparison of the analytical solution of the BE with the solution of WENO35-AG. Number of grid nodes=100.

100 80

1

Classical WENO35 WENO35-AG

0.8 ] s[ e m ti U P C

60

] -[ t

40 20 0 50

0.6 0.4 0.2

100

150 200 # Grid nodes [-]

250

300

Figure 3. Comparison of the CPU time for the classical WENO35 scheme and WENO35-AG when applied to the BE.

0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x [-]

Figure 4. Location of grid nodes at discrete points of time during the simulation of the BE for theWENO35-AG scheme.

The required CPU time for the WENO35 on this uniform grid and for the WENO35-AG on a 100 nodes gird are approximately 58 [s] and 26 [s] respectively as illustrated in Fig. 3. Hence, even though the calls of the grid adaptation subroutines in the WENO35AG method entail an additional computational load, the total CPU time required to achieve a certain accuracy is considerably less than that for the classical WENO35 scheme due to the possibility of using a coarser grid with WENO35-AG. Growth of particles in emulsion polymerization Emulsion polymerization is a radical polymerization process that is used to manufacture several polymers of commercial importance. The modeling of the process is quite involved, and when it is desired to describe the particle size distribution (PSD), a PBE is usually used, and this significantly increases the complexity of the model. The reader is referred to [6] for detailed information about the process. The polymer particles are of different sizes ( r ) and are located at different positions in the reactor. Assuming a well-mixed reactor, the spatial variation can be neglected, furthermore, to simplify the task the particles are assumed to be colloidally stable (no coagulation) and only the growth of a seed of polymer particles is considered. The PBE for this case reads:

∂n(r , t ) ∂ (r(r , t )n(r , t )) , =− ∂t ∂r

(9)

A New Numerical Solution Scheme for Tracking Sharp Moving Fronts

1-

911

(a)

]

120

m

1-d

t=0.75 [-] t=1 [-]

100

l l o m [ n oi ct n uf y t si n e D

t=0.5 [-] t=0.25 [-]

80 60

t=0 [-]

40 20 0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.3 0.4 0.5 0.6 0.7 Normalized particle radius [-]

0.8

0.9

1

(b) ] 1 -[ e m ti 0.75 h ct 0.5 a b d e 0.25 zil a m 0 r o 0 N

0.1

0.2

Figure 5. (a) Comparison of the MOC solution (continuous line) with the simulation results of the WENO35-AG scheme (− · − · −) for the PBE in emulsion polymerization on a 100 nodes grid. (b) Location of grid nodes during simulation for the WENO35-AG scheme.

where n(r , t ) is the population density function and r( r , t ) is the growth rate of polymer particles and is given by Eq.(10):

r(r , t ) =

k p MWM

[M ] p n (r , t ).

4πr ρ p N A 2

(10)

In Eq. (10), kp is the propagation rate coeffecient, MWM is the molecular weight of the monomer, NA is Avogadro’s number, ρp is the polymer density, [M]p is the monomer concentration in the particle phase and n is the average number of radicals per particle. Fig. 5 illustrates the evolution of the PSD and the grid over the normalized batch time. As can be observed, the grid moves to track the movement of the PSD and the nodes are concentrated in the regions of pronounced spatial variations. 120 1-

]

m

100

l l o m [ n oi ct n uf y t si n e D

80

1-d

60 40 20 0 0.35

0.4 0.45 Normalized particle radius [-]

0.5

Figure 6. Comparison of WENO35 on a grid made of 220 nodes (filled circles) with the WENO35-AG (dotted line) and MOC (continuous line) on a grid made of 100 nodes when applied to the PBE in emulsion polymerization.

A.E. Bouaswaig and S. Engell

912

In addition, the solution profile is neither spoiled by spurious oscillations nor by numerical diffusion. This fact is more evident in Fig. 6 that depicts the PSD at the end of the batch. The results shown in this figure are for WENO35-AG on a grid consisting of 100 nodes, WENO35 on a uniform grid composed of 220 nodes and a semi-analytical solution obtained by applying the method of characteristics (MOC) [7]. The results are comparable, however, when WENO35 is used, the required computation time is 318 [s] while its counterpart for the WENO35-AG is 205 [s]. This supports the previous conclusion that, being able to use a coarser grid with WENO-AG reduces the total computation time despite the additional computational load that arises from the calls of the grid adaptation subroutine.

4. Conclusions To avoid the undesired excessive numerical diffusion when dealing with hyperbolic and convection dominated parabolic PDEs a combination of a high order numerical discretization method and grid adaptation can be applied. This will guarantee concentrating the grid nodes in the locations where they are needed the most and will enable the use of a relatively coarser grid without having to sacrifice the accuracy. It is shown here that coupling the WENO scheme with static grid adaptation provides a tool that has proven to be efficient when applied to two case studies of different nature: a non-linear parabolic PDE, the viscid Burgers’ equation, and a challenging nonlinear hyperbolic PDE, the PBE in the context of emulsion polymerization. Results in both cases reveal that the proposed algorithm provides better performance when compared with the classical WENO scheme making it a promising option for similar problems.

5. Acknowledgements The financial support of the General People’s Committee for Higher Education (Libya) and of The Graduate School of Production Engineering and Logistics (Technische Universität Dortmund) is gratefully acknowledged.

References [1] Y. I. Lim, J. M. Le Lann, and X. Joulia, Moving mesh generation for tracking a shock or steep moving front, Comp. Chem. Eng., 25(2001) 653. [2] J. Smit, M. van Sint Annaland and J. A. M. Kuipers, Grid adaptation with WENO schemes for non-uniform grids to solve convection-dominated partial differential equations, Chem. Eng. Sci., 60(2005) 2609. [3] J. M. Sanz-Serna and I. Christie, A simple adaptive technique for nonlinear wave problems, J. Comput. Phys., 67(1986) 348. [4] C.-W. Shu, Essentially non-oscilllatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, In: A. Quarteroni (Ed.), Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, Lecture Notes in Math. vol. 1697, Springer-Verlag, Berlin, 1998. [5] A. Vande Wouwer, P. Saucez and W. E. Schiesser, Simulation of distributed parameter systems using a Matlab-based method of lines toolbox: Chemical engineering Applications, Ind. Eng. Chem. Res., 43(2004) 3469. [6] R.G. Gilbert, Emulsion polymerization: A mechanistic approach, Academic press, London, 1995. [7] S. Sarra, The method of characteristics & conservation laws, Journal of Online Mathematics and Applications, 3(2003).