Monash Workshop on Numerical Di ff erential Equations and Applications Melbourne, Australia, 10-14 February 2020 Semi-analytical solutions for transport PDEs in heterogeneous media Dr Elliot Carr � elliot.carr@qut.edu.au � @ElliotJCarr � https: // elliotcarr.github.io /
Transport equations Heterogeneous media ◮ Generic scalar transport equation: R ( x ) ∂ c Ω ⊂ R d . ∂ t = ∇ · ( D ( x ) ∇ c − v ( x ) c ) + S ( c , x ) , ◮ Heterogeneous media: coe ffi cients vary spatially. ◮ This talk is comprised of two parts: • Part 1: Semi-analytical solutions to the advection-di ff usion-reaction equation in heteroge- neous (layered) media. • Part 2: Semi-analytical solutions to the homogenization boundary value problem for di ff u- sion in 2D heterogeneous media. Dr Elliot Carr https: // elliotcarr.github.io / 1 / 31
Advection-di ff usion-reaction in layered media Problem description − Inlet Outlet c 1 ( x , t ) c 2 ( x , t ) c m − 1 ( x , t ) c m ( x , t ) x 0 ℓ 1 ℓ 2 ℓ m − 2 ℓ m − 1 L Layer m − 1 Layer 1 Layer 2 Layer m ( ) R ( x ) ∂ c ∂ t = ∂ D ( x ) ∂ c ∂ x − v ( x ) c − µ ( x ) c + γ ( x ) . ∂ x R 1 , D 1 , v 1 , µ 1 , γ 1 , 0 < x < ℓ 1 , R 2 , D 2 , v 2 , µ 2 , γ 2 , ℓ 1 < x < ℓ 2 , R ( x ) , D ( x ) , v ( x ) , µ ( x ) , γ ( x ) = . . . . . . R m , D m , v m , µ m , γ m , ℓ m − 1 < x < L . Dr Elliot Carr https: // elliotcarr.github.io / 2 / 31
Advection-di ff usion-reaction in layered media Governing equations ◮ Governing equations (Guerrero et al., 2013; van Genuchten and Alves, 1982): ∂ 2 c i ∂ c i ∂ c i R i ∂ t = D i ∂ x 2 − v i ∂ x − µ i c i + γ i , i = 1 , . . . , m , c i ( x , 0) = f i , c i ( ℓ i , t ) = c i + 1 ( ℓ i , t ) , ∂ c i ∂ c i + 1 θ i D i ∂ x ( ℓ i , t ) = θ i + 1 D i + 1 ∂ x ( ℓ i , t ) , where v i θ i = v i + 1 θ i + 1 . ◮ Nomenclature: • c i ( x , t ): solute concentration [ML − 3 ] in the i th layer • R i : retardation factor [ − ] • D i : dispersion coe ffi cient [L 2 T − 1 ] • v i : pore-water velocity [LT − 1 ] • µ i : rate constant for first-order decay [T − 1 ] • γ i : rate constant for zero-order production [T − 1 ] • θ i : volumetric water content [L 3 L − 3 ] in the i th layer Dr Elliot Carr https: // elliotcarr.github.io / 3 / 31
Advection-di ff usion-reaction in layered media Typical boundary conditions ◮ Inlet boundary condition ( x = 0): • Concentration-type: c 1 (0 , t ) = c 0 ( t ) , • Flux-type: ∂ c 1 v 1 c 1 (0 , t ) − D 1 ∂ x (0 , t ) = v 1 c 0 ( t ) , ◮ Outlet boundary condition ( x = L ): ∂ c m ∂ x ( L , t ) = 0 . ◮ General boundary conditions: ∂ c 1 Inlet: a 0 c 1 (0 , t ) − b 0 ∂ x (0 , t ) = g 0 ( t ) , ∂ c m Outlet: a L c m ( L , t ) + b L ∂ x ( L , t ) = g L ( t ) . Dr Elliot Carr https: // elliotcarr.github.io / 4 / 31
Advection-di ff usion-reaction in layered media Analytical solution via eigenfunction expansion ◮ Eigenfunction expansion solution: � ∞ c i ( x , t ) = a n T n ( λ n ; t ) X n ( λ n ; x ) . n = 1 ◮ Eigenvalues ( λ n , n ∈ N + ) are identified by substituting eigenfunctions into the boundary and interface conditions and enforcing a non-trivial solution. ◮ This yields a nonlinear transcendental equation for the eigenvalues arising from the evaluation of a 2 m × 2 m determinant f ( λ ) = 0 , A ∈ R 2 m × 2 m . where f ( λ ) : = det( A ( λ )) , ◮ For many layers (large m ) evaluating f ( λ ) is numerically unstable. ◮ Solutions tend to breakdown for m > 10 layers (Carr and Turner, 2016). ◮ Solutions for maximum of seven layers given by Liu et al. (1998) (advection-di ff usion only with µ i = γ i = 0) and Guerrero et al. (2013) (advection-di ff usion-reaction with γ i = 0). Dr Elliot Carr https: // elliotcarr.github.io / 5 / 31
Advection-di ff usion-reaction in layered media Analytical solution via Laplace transform ◮ Idea: reformulate the model into m isolated single layer problems (Carr and Turner, 2016; Rodrigo and Worthy, 2016; Zimmerman et al., 2016). ◮ Introduce unknown functions of time, g i ( t ) ( i = 1 , . . . , m − 1), at the layer interfaces (Carr and Turner, 2016; Rodrigo and Worthy, 2016): ∂ c i g i ( t ) : = θ i D i ∂ x ( ℓ i , t ) . ◮ Yields isolated single layer problems e.g. in the first layer: ∂ 2 c 1 ∂ c 1 ∂ c 1 R 1 ∂ t = D 1 ∂ x 2 − v 1 ∂ x − µ 1 c 1 + γ 1 , c 1 ( x , 0) = f 1 , ∂ c 1 a 0 c 1 (0 , t ) − b 0 ∂ t (0 , t ) = g 0 ( t ) , ∂ c 1 θ 1 D 1 ∂ x ( ℓ 1 , t ) = g 1 ( t ) . ◮ Each problem coupled together by imposing continuity of concentration at the interfaces. Dr Elliot Carr https: // elliotcarr.github.io / 6 / 31
Advection-di ff usion-reaction in layered media Analytical solution via Laplace transform ◮ Solve each layer problem expressing the solution in terms of the unknown interface functions. ◮ Taking Laplace transforms yields boundary value problems e.g. in the first layer: d 2 C 1 d C 1 d x − ( µ 1 + R 1 s ) C 1 = − R 1 f 1 − γ 1 D 1 d x 2 − v 1 s , d C 1 a 0 C 1 (0 , s ) − b 0 d x (0 , s ) = G 0 ( s ) , d C 1 θ 1 D 1 d x ( ℓ 1 , s ) = G 1 ( s ) , where C i ( x , s ) = L{ c i ( x , t ) } denotes the Laplace transform of c i ( x , t ) with transformation variable s ∈ C and G i ( s ) = L{ g i ( t ) } for i = 1 , . . . , m − 1. ◮ Laplace transforms of the boundary functions: G 0 ( s ) = L{ g 0 ( t ) } G L ( s ) = L{ g L ( t ) } are assumed to be able to be found analytically. Dr Elliot Carr https: // elliotcarr.github.io / 7 / 31
Advection-di ff usion-reaction in layered media Analytical solution via Laplace transform ◮ The boundary value problems all involve second-order constant-coe ffi cient di ff erential equations ◮ Solving using standard techniques defines the concentration in the Laplace domain: C 1 ( x , s ) = A 1 ( x , s ) G 0 ( s ) + B 1 ( x , s ) G 1 ( s ) + P 1 ( x , s ) , i = 2 , . . . , m − 1 , C i ( x , s ) = A i ( x , s ) G i − 1 ( s ) + B i ( x , s ) G i ( s ) + P i ( x , s ) , C m ( x , s ) = A m ( x , s ) G m − 1 ( s ) + B m ( x , s ) G L ( s ) + P m ( x , s ) , where the functions P i , A i and B i ( i = 1 , . . . , m ) are known functions. ◮ To determine G 1 ( s ) , . . . , G m − 1 ( s ), the Laplace transformations of the unknown interface functions g 1 ( t ) , . . . , g m − 1 ( t ), we enforce continuity of concentration at each interface in the Laplace domain: C i ( ℓ i , s ) = C i + 1 ( ℓ i , s ) , i = 1 , . . . , m − 1 . (1) ◮ Thisyieldsatridiagonal systemoflinearequations Ax = b , where x = ( G 1 ( s ) , . . . , G m − 1 ( s )) T . ◮ Summary: Concentration can be evaluated at any x and s in the Laplace domain. Dr Elliot Carr https: // elliotcarr.github.io / 8 / 31
Advection-di ff usion-reaction in layered media Analytical solution via Laplace transform ◮ Inversion of the Laplace transform is carried out numerically. ◮ Hence, our solution method is semi -analytical. ◮ Trefethen et al. (2006) defines the following approximation: � � � N c i ( x , t ) = L − 1 { C i ( x , s ) } ≈ − 2 t ℜ w k C i ( x , s k ) , k = 1 k odd where N is even, s k = z k / t and w k , z k ∈ C are the residues and poles of the best ( N , N ) rational approximation to e z on the negative real line. ◮ Summary: Concentration can be evaluated at any x and t in the time domain. ◮ Attractiveness is that the solution is completely explicit. Unlike eigenfunction expansion solutions that require a nonlinear algebraic equation to be solved for the eigenvalues: f ( λ ) = 0 , A ∈ R 2 m × 2 m . f ( λ ) : = det( A ( λ )) , where Dr Elliot Carr https: // elliotcarr.github.io / 9 / 31
Advection-di ff usion-reaction in layered media Heaviside inlet boundary condition ◮ In solute transport problems, it is common to apply a Heaviside step function at the inlet: c 0 , 0 < t < t 0 , c 0 ( t ) = c 0 H ( t 0 − t ) = 0 , t > t 0 , where c 0 is a constant and t 0 > 0 is the duration. ◮ Yields G 0 ( s ) = exp( − t 0 s ) / s and G 0 ( s ) = v 1 exp( − t 0 s ) / s for the concentration-type and flux-type boundary condition, respectively. ◮ Such exponential functions are well known to cause numerical problems in algorithms for inverting Laplace transforms (Kuhlman, 2013). ◮ To overcome this problem, we use superposition of solutions � c i ( x , t ) , 0 < t < t 0 , c i ( x , t ) = � c i ( x , t ) − � c i ( x , t − t 0 ) , t > t 0 , where � c i ( x , t ) is the solution with g 0 ( t ) = c 0 and � c i ( x , t ) is the solution with g 0 ( t ) = c 0 , f i = 0 and γ i = 0. Dr Elliot Carr https: // elliotcarr.github.io / 10 / 31
Advection-di ff usion-reaction in layered media One layer test case ∂ c 1 ∂ c 2 v 1 c 1 (0 , t ) − D 1 BCs : ∂ x (0 , t ) = v 1 c 0 , ∂ t (20 , t ) = 0 . Benchmarked against single-layer analytical solutions (van Genuchten and Alves, 1982). Absolute errors t = 10 − 3 t = 0 . 1 t = 0 . 2 t = 0 . 4 t = 0 . 6 t = 4 4 . 11 × 10 − 14 5 . 53 × 10 − 10 8 . 69 × 10 − 9 1 . 24 × 10 − 9 5 . 84 × 10 − 8 6 . 10 × 10 − 10 Dr Elliot Carr https: // elliotcarr.github.io / 11 / 31
Recommend
More recommend
