The use of a time-fractional transport model for performing computational homogenisation of 2D heterogeneous media exhibiting memory effects
TThe use of a time-fractional transport model for performing computationalhomogenisation of 2D heterogeneous media exhibiting memory effects
Libo Feng a , Ian Turner a,b, ∗ , Patrick Perr´e c , Kevin Burrage a,b,d a School of Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, QLD. 4001, Australia b Australian Research Council Centre of Excellence for Mathematical and Statistical Frontiers (ACEMS), Queensland University ofTechnology (QUT), Brisbane, Australia c Laboratory of Chemical Engineering and Materials (LGPM), CentraleSup´elec, Universit´e Paris-Saclay, France d Visiting Professor, Department of Computer Science, University of Oxford, OXI 3QD, UK
Abstract
In this work, a two-dimensional time-fractional subdiffusion model is developed to investigate the underlying trans-port phenomena evolving in a binary medium comprised of two sub-domains occupied by homogeneous material.We utilise an unstructured mesh control volume method to validate the model against a derived semi-analyticalsolution for a class of two-layered problems. This generalised transport model is then used to perform compu-tational homogenisation on various two-dimensional heterogenous porous media. A key contribution of our workis to extend the classical homogenisation theory to accommodate the new framework and show that the effectivediffusivity tensor can be computed once the cell problems reach steady state at the microscopic scale. We verifythe theory for binary media via a series of well-known test problems and then investigate media having inclusionsthat exhibit a molecular relaxation (memory) effect. Finally, we apply the generalised transport model to estimatethe bound water diffusivity tensor on cellular structures obtained from environmental scanning electron microscope(ESEM) images for Spruce wood and Australian hardwood. A highlight of our work is that the computed diffusivityfor the heterogeneous media with molecular relaxation is quite different from the classical diffusion cases, beingdominated at steady-state by the material with memory effects.
Keywords: control volume method, homogenisation theory, two-layered problems, time-fractional derivative,periodic boundary, heterogeneous medium
1. Introduction
Originating from Fourier’s and Darcy’s law, numerous real-world problems involve transport processes in het-erogeneous media, especially multi-scale transport processes, in which the medium properties differ spatially. Oneimportant multi-scale problem is the distributed microstructure model or double/dual-porosity model [1, 2, 3],where a large scale domain has a number of small isolated or disconnected inclusions embedded within it.To model this type of problem, macroscopic and microscopic equations are used to describe the global transportin the connected domain and the local transport in the isolated domains, respectively. Some common approachesfor the dual-scale problems include the heterogeneous multi-scale method [1], distributed microstructure model [2]and “equation-free” approach [4]. In [5], Carr et al. proposed an extended distributed microstructure model, whichcan effectively describe the non-equilibrium field in the inclusions for the water flow in unsaturated soils. However,a number of recent studies highlight that abnormal transport phenomena is evident in some complex heterogeneousmedia [6, 7, 8], which deviates from the classical Fickian diffusion. Generally, asymptotic long-time behaviour ofthe mean square displacement is utilised to characterise anomalous diffusion, which has the form (cid:104) x ( t ) (cid:105) ∼ K γ Γ(1 + γ ) t γ , t → ∞ , ∗ Corresponding author.
Email addresses: [email protected] (Libo Feng), [email protected] (Ian Turner), [email protected] (Patrick Perr´e), [email protected] (Kevin Burrage)
Preprint submitted to arXiv February 5, 2021 a r X i v : . [ m a t h . NA ] F e b here K γ is the generalised diffusion coefficient and γ (0 < γ <
1) is the anomalous diffusion exponent. It turnsout that mathematical models involving fractional derivatives are useful tools to treat the anomalous diffusivetransport in heterogeneous media [6, 9, 10], which motivates us to explore these generalised transport equations inthis research.In the seminal paper by Whitaker on the method of volume averaging, macroscopic equations describing thedrying process were derived with nonlinear effective parameters within their definitions. This theory also included aset of closure equations that relates these parameters to the pore morphology of the medium, which facilitates theirprediction. Homogenization theory can also be used to predict these parameters from the microstructure of theporous medium [12, 13]. The process starts with a full description of the model at the scale of the heterogeneitiesand upscales the relevant information to derive the macroscopic conservation law by replacing the actual geometrywith a periodic idealisation. The macroscopic equation for this homogenised medium is then derived in the limitas the heterogeneities tend to zero. In the classical diffusion model, the effective parameters are defined in termsof the solution of an elliptic equation on the period, subject to periodic boundary conditions (see for example,[1, 2, 3, 5, 14]). A key contribution of our work is the derivation of a time-evolutionary homogenisation theory thatis applied to predict the diffusivity of porous media that exhibit memory effects. A good example of such a materialis wood, in which the cellular structure undergoes a molecular relaxation phenomenon during drying [15, 16].To model this phenomenon we propose a generalised transport equation involving a time-fractional operator. Achallenge is then the way to deal with the fluxes at the interfaces between two different media when one exhibitsmemory and the other does not. We show how to correctly model this scenario, paying careful attention to thetreatment of the interfacial boundary between the different media.The types of fractional operators utilised in these models mainly are time-fractional derivatives, which is generallyused to describe the algebraic decaying or power-law waiting time in a L´evy process. For the time-fractionalderivative, it has a wide application in many problems [17, 18, 19, 20, 21, 22, 23, 24]. We start with the followingtwo-dimensional time-fractional subdiffusion equation with variable coefficients ∂u ( x, y, t ) ∂t = R D − γt [ ∇ · ( Q ∇ u ( x, y, t ))] , ( x, y, t ) ∈ Ω × (0 , T ] , (1)subject to the initial and Neumann boundary conditions given, respectively, as u ( x, y,
0) = φ ( x, y ) , ( x, y ) ∈ Ω , R D − γt [ Q ∇ u ( x, y, t ) · n ] = ψ ( x, y, t ) , ( x, y, t ) ∈ ∂ Ω × (0 , T ] , (2)where ∇ = [ ∂∂x , ∂∂y ] T , Q = diag[ q ( x, y ) , q ( x, y )], ψ ( x, y, t ) = diag[ ψ ( x, y ) , ψ ( x, y )], n is the unit vector normalto ∂ Ω outward to Ω and D − γt (0 < γ <
1) is the Riemann-Liouville fractional derivative defined by R D − γt u ( x, y, t ) = 1Γ( γ ) ∂∂t (cid:90) t u ( x, y, ξ )( t − ξ ) − γ dξ. As an application, we will extend the equation to a time-fractional model for simulating transport in heterogeneoussystems, which is motivated by the models in [5] and [25]. In [5], Carr et al. considered an extended distributedmicrostructure model for gradient-driven transport on a binary medium without memory. In [25], Zeng et al.investigated a coupled system of Caputo time-fractional partial diffusion equations. Different from the models in[5] and [25], this paper deals with a time-fractional system based on the Riemann-Liouville fractional derivative ondomains including irregularly shaped inclusions with memory. We derive the semi-analytical solution for a classof two-layered problems to validate the proposed numerical method and present an homogenisation theory andnumerical calculation, for two-dimensional periodic structures with memory effects. The main contributions of thispaper are summarised as follows: • A two-dimensional time-fractional subdiffusion equation with variable coefficients is considered, in which anunstructured mesh control volume method (CVM) is applied to solve the problem. Given the non-smoothnessof the solution, the modified weighted shifted Gr¨unwald-Letnikov (WSGL) formula with starting weights isutilised to deal with the non-smooth problem, which is further extended to solve the time-fractional transportmodels in heterogeneous media. • A novel time-fractional transport model on a binary medium (a porous medium comprised of two distinctphases) with regular or irregular-shaped inclusions is simulated using the CVM. To verify the accuracy of2he computational model we derive a semi-analytical solution for a class of two-layered problems with quasi-periodic boundary conditions, in which the finite Fourier and Laplace transforms together with a numericalinverse Laplace transform technique are used. In addition, the mass balance equation for the layered mediumis also developed. • Homogenisation theory is extended from a classical diffusion equation to a time-fractional diffusion equationfor carrying out computational homogenisation on some well-known test problems and then on media havinginclusions that exhibit a molecular relaxation effect. The generalised transport model is also applied toestimate the bound water diffusivity tensor on cellular structures obtained from ESEM images for Sprucewood and Australian hardwood (Eucalyptus pilularis). An important finding is that, unlike the classicalcase, the fractional-order indices have a significant effect on the mass transfer and the equivalent diffusivityis dominated by the material with memory effects.The structure of this paper is as follows. In Section 2, the modified WSGL formulae with starting weights for thetime-fractional operator is proposed and the unstructured mesh control volume method for the two-dimensionalsubdiffusion problem (1) is introduced. We show that the proposed computational scheme offers second orderconvergence in time and space. In Section 3, a time-fractional transport model on a binary medium is proposedand the semi-analytical solution for a class of two-layered problems is derived. In Section 4, homogenisation theoryfor the time-fractional diffusion equation is developed. A series of numerical examples are considered in Sections 5and 6 to verify the homogenisation theory, in which different heterogeneous media morphologies and wood cellularperiodic structures are simulated. The numerical results show that the fractional-order indices have a significanteffect on the mass transfer for the time-fractional transport equation, which is very different to the classical case.Finally, some conclusions are drawn in Section 7.
2. The control volume method for the subdiffusion equation
Firstly, we present the grid partition used to discretise equation (1) in the temporal direction. Define t n = nτ , n = 0 , , , . . . , N , where τ = TN is the uniform temporal step. The Riemann-Liouville time-fractional derivativecan be rewritten in a more general form as R D αt v ( t ) = − α ) ddt (cid:82) t v ( ξ )( t − ξ ) α dξ , 0 < α <
1. To approximate theRiemann-Liouville fractional derivative, the shifted Gr¨unwald-Letnikov difference operator was proposed [26, 27]: A ατ,p v ( t ) = τ − α (cid:80) ∞ k =0 g ( α ) k v ( t − ( k − p ) τ ), where g ( α ) k = ( − k (cid:0) αk (cid:1) and p is an integer. To improve the accuracy ofthe approximation formula, some high order schemes based on the WSGL formula have been developed [27, 28].Here we focus on the second order WSGL formula, which has the form B ατ,p,q v ( t ) = α − q p − q ) A ατ,p v ( t ) + 2 p − α p − q ) A ατ,q v ( t ) , where p and q are two non equal integers. According to [28], p and q should be chosen satisfying | p | ≤ | q | ≤ p, q ) = (0 , − α − q p − q ) = α , p − α p − q ) = − α . Then at t = t n , we obtain R D αt v ( t n ) = B α,nτ,p,q v + O ( τ ) = τ − α n (cid:88) k =0 ω ( α ) n − k v ( t k ) + O ( τ ) , where B α,nτ,p,q v = α A ατ,p v ( t n ) − α A ατ,q v ( t n ) and the convolution weights ω ( α ) k are defined as ω ( α )0 = α g ( α )0 , ω ( α ) k = α g ( α ) k − α g ( α ) k − , k ≥
1. Define D α,nτ,p,q v = ( B α,nτ,p,q v + B α,n − τ,p,q v ), then at t = t n − , we have R D αt v ( t n − ) = D α,nτ,p,q v + O ( τ ) = τ − α n (cid:88) k =0 D ( α ) n − k v ( t k ) + O ( τ ) , (3)where D ( α )0 = ω ( α )0 , D ( α ) k = ω ( α ) k + ω ( α ) k − , k ≥ Remark 2.1.
To guarantee the second-order convergence of (3), the required conditions are v ( t ) ∈ L ( R ) and R D α +2 t v ( t ) and its Fourier transform belongs to L ( R ) [28]. .2. The implementation of the control volume method To implement the control volume method, we partition the solution domain Ω into a mesh comprised of non-overlapping triangles. Denote T h as the triangulation, N e the number of triangles and h the maximum diameter ofthe triangular elements. For a group of triangles with the same vertex, a control volume can be formed by joiningthe midpoints and barycenters of each triangle. Integrating (1) over each control volume V i ( i = 1 , , . . . , N p ) andapplying Gauss’s Divergence Theorem, yields (cid:90) V i ∂u∂t dV i = R D − γt (cid:73) ∂V i ( Q ∇ u ) · n d Γ i , where d Γ i is an infinitesimal small line segment on ∂V i . Utilising a lumped mass approach to approximate the timederivative, gives (cid:52) V i ∂u i ∂t = R D − γt (cid:73) ∂V i ( Q ∇ u ) · n d Γ i , i = 1 , . . . , N p , (4)where (cid:52) V i is the volume of the control volume. We assume that the integration path is anticlockwise and the trianglevertices are numbered counter-clockwise. We choose piecewise linear polynomials on the triangles. Within anelement triangle e p , p = 1 , , . . . , N e , the shape function can be defined in the form N i ( x, y ) = ep ( a i x + b i y + c i ), i = 1 , ,
3, where a i , b i and c i are some constants and ∆ e p is the area of triangle element p . Combining theseshape functions, we can construct the basis function l k ( x, y ) and approximate u ( x, y, t ) as u ( x, y, t ) ≈ u h ( x, y, t ) = (cid:80) N P k =1 u k l k ( x, y ). The line integral in (4) can be approximated by the midpoint formula at each control surface: (cid:73) ∂V i ( Q ∇ u ) · n d Γ i = (cid:73) Γ i q ( x, y, t ) ∂u∂x dy − (cid:73) Γ i q ( x, y, t ) ∂u∂y dx = m i (cid:88) j =1 2 (cid:88) r =1 Å q ( x, y, t ) ∂u∂x ã (cid:12)(cid:12)(cid:12)(cid:12) ( x r ,y r ) ∆ y ij,r − m i (cid:88) j =1 2 (cid:88) r =1 Å q ( x, y, t ) ∂u∂y ã (cid:12)(cid:12)(cid:12)(cid:12) ( x r ,y r ) ∆ x ij,r = N P (cid:88) k =1 m i (cid:88) j =1 2 (cid:88) r =1 u k Å q ( x, y, t ) ∂l k ( x, y ) ∂x ã (cid:12)(cid:12)(cid:12)(cid:12) ( x r ,y r ) ∆ y ij,r − N P (cid:88) k =1 m i (cid:88) j =1 2 (cid:88) r =1 u k Å q ( x, y, t ) ∂l k ( x, y ) ∂y ã (cid:12)(cid:12)(cid:12)(cid:12) ( x r ,y r ) ∆ x ij,r , (5)where ( x r , y r ) is the mid-point of the control face and m i is the number of sub-control volumes associated with thenode i . For more details, the readers can refer to [29]. Then from (4), we can derive the following ODE system M d u dt = R D − γt Ku + F b , (6)where M = diag[ (cid:52) V , (cid:52) V , . . . , (cid:52) V N p ], u = [ u , u , . . . , u N p ] T , K is the stiffness matrix derived from (5) and F b is the contribution from the boundary conditions. Now we discuss the treatment of the boundary conditions. Forexample, a point i is on the boundary with a control volume and k and j are its adjacent points on the boundary(see Figure 1); k is the midpoint of k and i and j is the midpoint of j and i ; i is the midpoint of j and i and i is the midpoint of k and i . Then according to (2) and (5), we obtain F b ( i ) = R D − γt (cid:34) Å q ( x, y, t ) ∂u∂x ã (cid:12)(cid:12)(cid:12)(cid:12) ( x i ,y i ) ( y j − y i ) − Å q ( x, y, t ) ∂u∂y ã (cid:12)(cid:12)(cid:12)(cid:12) ( x i ,y i ) ( x j − x i ) (cid:35) + R D − γt (cid:34) Å q ( x, y, t ) ∂u∂x ã (cid:12)(cid:12)(cid:12)(cid:12) ( x i ,y i ) ( y i − y k ) − Å q ( x, y, t ) ∂u∂y ã (cid:12)(cid:12)(cid:12)(cid:12) ( x i ,y i ) ( x i − x k ) (cid:35) = ψ ( x i , y i , t )( y j − y i ) − ψ ( x i , y i , t )( x j − x i ) + ψ ( x i , y i , t )( y i − y k ) − ψ ( x i , y i , t )( x i − x k ) . At t = t n − , applying the Crank-Nicolson scheme to the integer time derivative and the WSGL formula to thetime-fractional derivative, we obtain M u n − u n − τ = τ γ − n (cid:88) k =1 D (1 − γ ) n − k K ( u k − u ) + F n − b + t γ − n − Ku Γ( γ ) . (7)4 igure 1: An illustration of a point with control volume on the boundary. Since the WSGL formulae require some restrictions on the solution ( R D α +2 t v ( t ) ∈ L ( R )), they may not exhibitthe expected high order accuracy when the solution is not smooth enough, especially if it involves terms of theform t σ r . To overcome this issue, the modified WSGL formulae with correction terms were proposed to deal withproblems exhibiting a non-smooth solution [33, 34]. At t = t n − , we have R D αt v ( t n − ) = τ − α n (cid:88) k =0 D ( α ) n − k v ( t k ) + τ − α m (cid:88) k =1 E ( n,α ) k v ( t k ) + O ( τ ) , (8)where the starting weights E ( n,α ) k are chosen such that (8) is exact for v ( t ) = t σ r (1 ≤ r ≤ m ), which leads to thefollowing linear system: m (cid:88) k =1 E ( n,α ) k k σ r = Γ( σ r + 1)Γ( σ r + 1 − α ) Å n − ã σ r − α − n (cid:88) k =0 D ( α ) n − k k σ r , r = 1 , . . . , m. In addition, we also need to define the starting weights for the first order time derivative at t = t n − , i.e., dv ( t n − ) dt = v ( t n ) − v ( t n − ) τ + τ − m (cid:88) k =1 P ( n,α ) k v ( t k ) + O ( τ ) , (9)where the starting weights P ( n,α ) k are chosen such that (9) is exact for v ( t ) = t σ r (1 ≤ r ≤ m ), which can beobtained from the following equations: m (cid:88) k =1 P ( n,α ) k k σ r = σ r Å n − ã σ r − − (cid:0) n σ r − ( n − σ r (cid:1) , r = 1 , . . . , m. In this paper, we choose the correction terms σ r = rα , r = 1 , , . . . , m , such that ( m + 1) α ≥
2. For more details,the readers can refer to [33]. Then the modified numerical scheme of (7) with correction terms is
M u n − u n − τ + 1 τ m (cid:88) k =1 P ( n,α ) k M ( u k − u )= τ γ − (cid:34) n (cid:88) k =1 D (1 − γ ) n − k K ( u k − u ) + m (cid:88) k =1 E ( n,α ) k K ( u k − u ) (cid:35) + F n − b + t γ − n − Ku Γ( γ ) , (10)which also can be extended to solve the time-fractional transport model (11). We complete this section by assessing the accuracy and convergence order of the CVM method with the modifiedWSGL formulae for the two-dimensional time-fractional subdiffusion model (1) on a square domain Ω = [0 , × [0 , Q = I with initial condition u ( x, y,
0) = sin x sin y and boundary condi-tions R D − γt î Q ∂v ( x,y,t ) ∂x · n ó = t γ − E γ,γ ( − t γ ) cos x sin y, ( x, y, t ) ∈ ∂ Ω × (0 , T ], where the Mittag-Leffler function is E α,β ( z ) = (cid:80) ∞ k =0 z k Γ( kα + β ) . The exact solution of this problem is u ( x, y, t ) = E γ, ( − t γ ) sin x sin y .5 able 1: The error and convergence order of numerical schemes (7) and (10) for problem (1) for different h and γ = γ = γ with τ = 1 × − at t = 1, in which the number of correction terms m for γ = 0 . γ = 0 . m = 3 and m = 2, respectively. No correction Apply correction h ( γ = 0 .
5) error Order error Order2.6170E-01 1.2688E-03 – 1.3558E-03 –1.5176E-01 1.9206E-03 -0.76 4.4375E-04 2.058.4435E-02 2.1995E-03 -0.23 1.0488E-04 2.464.1624E-02 2.2646E-03 -0.04 2.8136E-05 1.862.1733E-02 2.2821E-03 -0.01 7.3674E-06 2.06 h ( γ = 0 .
8) error Order error Order2.6170E-01 1.2909E-03 – 1.3607E-03 –1.5176E-01 3.7649E-04 2.26 4.4321E-04 2.068.4435E-02 5.3879E-05 3.32 1.0482E-04 2.464.1624E-02 5.6789E-05 -0.07 2.8160E-05 1.862.1733E-02 7.2710E-05 -0.38 7.4373E-06 2.05We apply the numerical scheme (7) and the modified numerical scheme (10) to solve (1) and present the errorand convergence order for different h and γ with τ = 1 × − at t = 1 in Table 1. It is straightforward toobserve that the regularity of the exact solution is low particularly when γ is small, which can cause difficulty forthe numerical method to attain the theoretical temporal convergence order. This finding is evident in the resultspresented in Table 1, where we can see that the numerical scheme (7) without correction terms fails to provide goodconvergence behaviour. However, when correction terms are added, the singularity of the solution can be capturedeffectively and global second-order convergence is achieved. We conclude that the modified numerical scheme (10)with correction terms is a very effective computational method, thus justifying its choice for solving the generalisedtransport models considered throughout the following sections.
3. Application to a heterogeneous medium
In this section, we extend the control volume method to solve a time-fractional transport model on a binarymedium consisting of a connected phase Ω and an inclusion (disconnected phase) Ω . The time-fractional indicesare different for the two media (see Figure 2). Here we consider the following system: ® ∂u ( x,y,t ) ∂t = R D − γ t [ ∇ · ( Q ∇ u ( x, y, t ))] , ( x, y, t ) ∈ Ω × (0 , T ] , ∂v ( x,y,t ) ∂t = R D − γ t [ ∇ · ( Q ∇ v ( x, y, t ))] , ( x, y, t ) ∈ Ω × (0 , T ] , (11)subject to u ( x, y,
0) = φ ( x, y ) , ( x, y ) ∈ Ω ∪ Γ , v ( x, y,
0) = φ ( x, y ) , ( x, y ) ∈ Ω ∪ Γ ∪ Γ , R D − γ t [ Q ∇ v ( x, y, t ) · n ] = ψ ( x, y, t ) , ( x, y, t ) ∈ Γ × (0 , T ] . In addition, we also need the following boundary conditions at the interface Γ between the two media u ( x, y, t ) = v ( x, y, t ) , R D − γ t ( Q ∇ u ( x, y, t ) · n ) = R D − γ t ( Q ∇ v ( x, y, t ) · n ) , ( x, y, t ) ∈ Γ × (0 , T ] . (12)Integrating (11) over V i ∈ Ω and V j ∈ Ω , respectively, yields ® (cid:82) V i ∂u∂t dV i = R D − γ t (cid:72) ∂V i ( Q ∇ u ) · n d Γ i , (cid:82) V j ∂v∂t dV j = R D − γ t (cid:72) ∂V j ( Q ∇ v ) · n d Γ j . igure 2: An illustration of a binary medium comprised of two sub-domains occupied by homogeneous material. The square domain(medium 2) is connected and the circular inclusion (medium 1) is isolated or disconnected. Furthermore, the matrix form can be derived as ® M d u dt = R D − γ t K u + F b , M d v dt = R D − γ t K v + F b − F b , (13)where u = [ u , u ] T , v = [ v , v ] T , u ∈ Ω , u = v ∈ Γ (which is the shared part at the interface Γ ) and v ∈ Ω , F b is the contribution from the boundary conditions at the internal interface Γ and F b is the contributionfrom the boundary conditions at the external interface Γ . Denote U = [ u , u , v ] T , then the system (13) canbe recast as M d U dt = R D − γt KU + F b , which can be solved using the same technique discussed in the previous section. In this section, we consider the semi-analytical solution for the following 2D time-fractional two-layered problem(see Figure 3). ® ∂u ( x,y,t ) ∂t = R D − α t [ D ∆ u ( x, y, t )] , ( x, y, t ) ∈ Ω × (0 , T ] , ∂v ( x,y,t ) ∂t = R D − α t [ D ∆ v ( x, y, t )] , ( x, y, t ) ∈ Ω × (0 , T ] , (14)where Ω = ( l , l ) × (0 , l y ) and Ω = ( l , l ) × (0 , l y ) \ Ω . Since the medium is homogeneous in the y -direction andwe impose no flux boundary conditions at y = 0 and y = l y , solving problem (14) can be reduced to solving thefollowing 1D two-layered problem: ∂u ( x,t ) ∂t = R D − α t (cid:16) D ∂ u ( x,t ) ∂x (cid:17) , ( x, t ) ∈ ( l , l ) × (0 , T ] , ∂v ( x,t ) ∂t = R D − α t (cid:16) D ∂ v ( x,t ) ∂x (cid:17) , ( x, t ) ∈ ( l , l ) ∪ ( l , l ) × (0 , T ] . (15) Figure 3: An illustration of a two-layered medium.
In the following, we will consider the equivalent form of (15), which can be written as ∂X i ( x, t ) ∂t = R D − γ i t Å D bi ∂ X i ( x, t ) ∂x ã , x ∈ ( l i − , l i ) , i = 1 , , , (16)with the initial conditions X i ( x,
0) = X i ( x ) and quasi-periodic boundary conditions of the following form R D − γ t Å D b ∂X ∂x ( l , t ) ã = R D − γ t Å D b ∂X ∂x ( l , t ) ã , X ( l , t ) = X ( l , t ) + q , R D − γ i t Å D b,i ∂X i ∂x ( l i , t ) ã = R D − γ i +1 t Å D b,i +1 ∂X i +1 ∂x ( l i , t ) ã , X i ( l i , t ) = X i +1 ( l i , t ) , i = 1 , , u = X , v = X ∪ X , q is a constant, D b = D , D b = D b = D , γ = α and γ = γ = α . Theseboundary conditions have been chosen to validate the time evolutionary homogenisation theory proposed in theproceeding section. The problem (16) will be solved using a combination of finite Fourier and Laplace transforms,which is an extension of the semi-analytical approach discussed in [30]. First define d i = l i − l i − , i = 1 , , ϕ i , i = 1 , , − d dx ϕ i = λ i ϕ i , dϕ i dx ( l i − ) = 0 and dϕ i dx ( l i ) = 0 , is given by ‹ X i ( λ i,m , t ) := (cid:104) X i , ϕ i,m (cid:105) = (cid:82) l i l i − X i ( x, t ) ϕ i,m ( x ) dx , where the eigenvalues λ i,m = m π d i , m ≥
0. Thecorresponding eigenfunctions are ϕ i,m ( x ) = √ d i for m = 0 and ϕ i,m ( x ) = » d i cos [ λ m ( x − l i − )] for m >
0. Weintroduce the time fractional potential for each layer as X i = R D − γ i t X i , i = 1 , , i th layer to obtain ≠ ∂X i ∂t , ϕ i,m ∑ = D bi ≠ ∂ X i ∂x , ϕ i,m ∑ . The term on the righthand side is integrated by parts to obtain D bi ≠ ∂ X i ∂x , ϕ i,m ∑ = − D bi ß≠ − d dx ϕ i,m , X i ∑ + ï − ∂ X i ∂x ( l i , t ) ϕ i,m ( l i ) + ∂ X i ∂x ( l i − , t ) ϕ i,m ( l i − ) ò™ . (17)Next, define the interfacial flux terms as D b ∂ X ∂x ( l , t ) = D b ∂ X ∂x ( l , t ) = v ( t ), D b ∂ X ∂x ( l , t ) = D b ∂ X ∂x ( l , t ) = v ( t ), D b ∂ X ∂x ( l , t ) = D b ∂ X ∂x ( l , t ) = v ( t ). Substituting the relevant boundary condition information into (17),the transformed layer equations are then given by ∂ ‹ X ∂t = − D b λ ,m ‹ X + [ v ( t ) ϕ ,m ( l ) − v ( t ) ϕ ,m ( l )] ,∂ ‹ X ∂t = − D b λ ,m ‹ X + [ v ( t ) ϕ ,m ( l ) − v ( t ) ϕ ,m ( l )] ,∂ ‹ X ∂t = − D b λ ,m ‹ X + [ v ( t ) ϕ ,m ( l ) − v ( t ) ϕ ,m ( l )] , together with the transformed initial conditions ‹ X i ( λ i,m ,
0) = (cid:104) X i, ( x ) , ϕ i,m (cid:105) , i = 1 , ,
3. We now apply the Laplacetransform in time and denote ‹ X i ( λ i,m , s ) = L ¶ ‹ X i ( λ i,m , t ) © , i = 1 , ,
3. Since L ¶ R D − γ i t ‹ X i © = s − γ i ‹ X i ( λ i,m , s ),the Laplace transformation of the three layer equations given above can be rearranged as ‹ X ( λ ,m , s ) = ‹ X ( λ ,m , η ,m ( s ) + ¯ v ( s ) ϕ ,m ( l ) − ¯ v ( s ) ϕ ,m ( l ) η ,m ( s ) , (18) ‹ X ( λ ,m , s ) = ‹ X ( λ ,m , η ,m ( s ) + ¯ v ( s ) ϕ ,m ( l ) − ¯ v ( s ) ϕ ,m ( l ) η ,m ( s ) , (19) ‹ X ( λ ,m , s ) = ‹ X ( λ ,m , η ,m ( s ) + ¯ v ( s ) ϕ ,m ( l ) − ¯ v ( s ) ϕ ,m ( l ) η ,m ( s ) , (20)where η i,m ( s ) = s + D bi s − γ i λ i,m . In order to determine the three unknown interfacial flux values ¯ v ( s ), ¯ v ( s )and ¯ v ( s ), we need the boundary conditions at the interfaces: X ( l , s ) = X ( l , s ), X ( l , s ) = X ( l , s ) + q s and X ( l , s ) = X ( l , s ). Noting that X i ( x, s ) = (cid:80) ∞ m =0 ‹ X i ( λ i,m , s ) ϕ i,m ( x ), i = 1 , , × ∞ (cid:88) m =0 ‹ X ( λ ,m , s ) ϕ ,m ( l ) = ∞ (cid:88) m =0 ‹ X ( λ ,m , s ) ϕ ,m ( l ) , ∞ (cid:88) m =0 ‹ X ( λ ,m , s ) ϕ ,m ( l ) = ∞ (cid:88) m =0 ‹ X ( λ ,m , s ) ϕ ,m ( l ) + q s , ∞ (cid:88) m =0 ‹ X ( λ ,m , s ) ϕ ,m ( l ) = ∞ (cid:88) m =0 ‹ X ( λ ,m , s ) ϕ ,m ( l ) , (21)that can be solved for ¯ v ( s ), ¯ v ( s ) and ¯ v ( s ) at a given value of s using Cramer’s rule. Finally, the solutionswithin each layer can be determined by applying the inverse Laplace transform resolved numerically: ‹ X i ( λ i,m , t ) = L − (cid:110) ‹ X i ( λ i,m , s ) (cid:111) = 12 πi (cid:90) Γ e st ‹ X i ( λ i,m , s ) ds = 12 πi (cid:90) Γ e z t ‹ X i ( λ i,m , z/t ) dz ≈ − (cid:60) Ñ K/ (cid:88) k =1 c k − ‹ X i ( λ i,m , z k − /t ) t é , where z = st , c k − and z k − are the residues and poles of the best ( K, K ) rational approximation of e z on the neg-ative real line as computed by the Carath´eodory–Fej´er method. Full details of this numerical approach can be foundin Trefethen et al. [31]. Then the solution in each layer can be computed using X i ( x, t ) = (cid:80) ∞ m =0 ‹ X ( λ i,m , t ) ϕ i,m ( x ), i = 1 , ,
3. We will use this semi-analytical solution to validate the numerical solutions we compute using the finitevolume method summarised in Section 2.
We now determine mass balance equations to confirm the accuracy of our simulation results. Define L = (cid:80) j =1 d j and the average value of the variable X as follows: (cid:104) X (cid:105) = 1 L (cid:90) l l X ( x, t ) dx = 1 L (cid:88) i =1 (cid:90) l i l i − X i ( x, t ) dx = (cid:88) i =1 Å d i L ã (cid:104) X i (cid:105) , where the average value of X in the i th layer is given by (cid:104) X i (cid:105) = d i (cid:82) l i l i − X i ( x, t ) dx . Using the continuity of thefluxes at the interfaces, we now integrate (16) across each layer to obtain: d (cid:104) X (cid:105) dt = 0 , (22)subject to (cid:104) X (cid:105) (0) = L (cid:82) l l X ( x ) dx = (cid:104) X (cid:105) = (cid:80) i =1 (cid:0) d i L (cid:1) (cid:104) X i (cid:105) , where the i th layer average initial condition is (cid:104) X i (cid:105) = d i (cid:82) l i l i − X i ( x, dx . Then we can obtain that the average value of X must remain constant at its initialaverage for the duration of the simulation, namely (cid:104) X (cid:105) = (cid:104) X (cid:105) . (23)
4. Homogenisation theory for time-fractional differential equation
In this section, we will extend the homogenisation theory from classical differential equations to fractionaldifferential equations in heterogeneous media. Let Ω be a bounded domain in R with (sufficiently smooth) boundary ∂ Ω. At first, we note that the following two equations are equivalent. ∂u∂t = R D − γt [ ∇ ( D · ∇ u )] , (24) C D γt u = [ ∇ ( D · ∇ u )] . (25)9his equivalence can be verified by imposing the operator R D γ − t on both sides of (24) and using the relationshipbetween the Riemann-Liouville fractional operator and the Caputo fractional operator [32] to obtain (25). Dueto this equivalence, for ease of mathematical exposition we focus on the homogenisation theory for the fractionaldifferential equation (25) involving the Caputo fractional derivative. We consider the following generalised fractionaltransport model: C D γt u ε ( x, t ) = ∇ · ( D ε ( x, t ) ∇ u ε ( x, t )) , in Ω , t ∈ (0 , T ] , (26) u ε ( x,
0) = f ( x ) , in Ω , u ε ( x, t ) = 0 , on ∂ Ω , t ∈ (0 , T ] , where the Caputo derivative is defined as C D γt u ( x, t ) = − γ ) (cid:82) t ( t − η ) − γ ∂u ( x,η ) ∂η dη , and 0 < γ < D ε ( x, t ) is symmetric with D εij ∈ L ∞ (Ω) and D εij ( x, t ) ξ i ξ j ≥ Cξ i a.e. in Ω, for some constant C > ξ i ∈ R . We define D ε ( x, t ) := D (cid:0) xε , tε p (cid:1) , p >
0. The choiceof p will be discussed below and we assume there is no periodicity in the time variables. We follow the discussionfrom [14] by assuming further conditions on the structure of the functions D εij , namely:1. D εij ( x, t ) := D ij (cid:0) xε , tε p (cid:1) , p > D ij ( y, τ ) is Y -periodic as a function of y , where the fast variables in space and time are y = xε and τ = tε p ,respectively.3. D ij ( y, τ ) ∈ L ∞ ( R y × R τ ).Assuming the fast variable y and slow variable x are independent as ε →
0, we seek an asymptotic expansion of thesolution u ε ( x, t ) in the form u ε ( x, t ) = (cid:80) ∞ k =0 ε k u k ( x, y, t, τ ), with u k being Y -periodic in y and we will add furtherconstraints on u k regarding the time variable τ below. Before proceeding with the two-scale asymptotic expansionof the solution u ε , we first consider a general function f ( t, τ ), with τ = tε p , p > C D γt f ( t, τ ) = 1Γ(1 − γ ) (cid:90) t ( t − η ) − γ f (cid:48) ( η, τ ) dη. (27)Next, the chain rule gives ddt = ∂∂t + ε p ∂∂τ , i.e., f (cid:48) ( η, τ ) = ∂f∂t + ε p ∂f∂τ . Substituting this expression into (27), weobtain C D γt f ( t, τ ) = 1Γ(1 − γ ) (cid:90) t ( t − η ) − γ Å ∂f ( η, ζ ) ∂t + 1 ε p ∂f ( η, ζ ) ∂ζ ã dη = C D γt f ( t, τ ) + 1Γ(1 − γ ) ε p (cid:90) t ( t − η ) − γ ∂f ( η, ζ ) ∂ζ dη. The second term on the RHS can be written as1Γ(1 − γ ) ε p (cid:90) tεp ( t − ε p ζ ) − γ ∂f ( η, ζ ) ∂ζ ε p dζ = 1Γ(1 − γ ) (cid:90) τ ∂f ( η, ζ ) ∂ζ dζ ( τ − ζ ) γ ε pγ , where we have set ζ = ηε p and t = ε p τ . Our investigation of different choices for the parameter p identified that p = γ gives C D γt f ( t, τ ) = C D γt f ( t, τ ) + ε C D γτ f ( t, τ ). We will use this relation below to derive the time-fractionalunit cell model.Next, denote A ε = − ∂∂x i î D ij ( y, τ ) ∂∂x j ó and then expand A ε as: A ε = ε − A + ε − A + ε A , where A = − ∂∂y i î D ij ( y, τ ) ∂∂y j ó , A = − ∂∂y i î D ij ( y, τ ) ∂∂x j ó − ∂∂x i î D ij ( y, τ ) ∂∂y j ó , A = − ∂∂x i î D ij ( y, τ ) ∂∂x j ó . Sub-stituting the expressions for u ε and A ε into system (24) gives:LHS: Å C D γt + 1 ε C D γτ ã ( u + εu + ε u + . . . ) = ε − C D γτ u + ε − C D γτ u + ε (cid:0) C D γt u + C D γt u (cid:1) + . . . RHS: − A ε ( u + εu + ε u + . . . ) = − ε − A u − ε − ( A u + A u ) − ε ( A u + A u + A u ) + . . . ε , we obtain: O (cid:0) ε − (cid:1) : C D γτ u = − A u , (28) O (cid:0) ε − (cid:1) : C D γτ u = − ( A u + A u ) , (29) O (cid:0) ε (cid:1) : C D γt u + C D γτ u = − ( A u + A u + A u ) . (30)Due to the complexity introduced by the time-fractional operator and the time dependent periodic boundaryconditions at the external and interfacial boundaries, we choose to solve these problems to steady-state in τ atwhich point u = u ( x, t ) is independent of y and τ and (29) can be expressed as C D γτ u + A u = ∂D ij ∂y i ∂u∂x j . Next, we introduce θ j as the steady state Y -periodic solution of the following time-fractional problem C D γτ θ j = − A ( θ j + y j ) , j = 1 , . The solution u can be defined up to an additive constant as u ( x, y, t, τ ) = θ j ∂u∂x j + u ( x, t ) . Hence, the unit cell problem becomes to solve: C D γτ θ j = ∂∂y i ï D ij ( y, τ ) ∂ ( θ j + y j ) ∂y j ò , j = 1 , , subject to periodic boundary conditions in y and initially θ j ( y,
0) is specified. Denote by I Y ( · ) = | Y | (cid:82) Y · dy theaverage of a quantity in y . At steady-state in τ (30) becomes C D γt u = − I y ( A u + A u ) , which gives the homogenised equation as C D γt u = ∂∂x i Å D ij ( y, τ ) ∂u∂x j ã , where D = 1 | Y | (cid:90) Y D ( I + J Tθ ) dy, D = ï D D D D ò , J θ = ñ ∂θ ∂y ∂θ ∂y ∂θ ∂y ∂θ ∂y ô . We can express the unit cell problem in terms of the new variable ϕ j = θ j + y j , j = 1 ,
2, as C D γτ ϕ j = ∂∂y i ï D ij ( y, τ ) ∂ϕ j ∂y j ò . (31)This variable change affects the boundary conditions: ϕ ( L , y , τ ) = ϕ (0 , y , τ ) + L , ϕ ( y , L , τ ) = ϕ ( y , , τ ) ,ϕ ( L , y , τ ) = ϕ (0 , y , τ ) , ϕ ( y , L , τ ) = ϕ ( y , , τ ) + L . Also we need the boundary conditions (12) at the interface. The initial conditions become ϕ j ( y , y ,
0) = θ j ( y , y , y j , j = 1 ,
2. In this case we obtain the effective diffusivity tensor as D = 1 | Y | (cid:90) Y DJ Tϕ dy. . Numerical examples In this section, we present the results obtained from the computational homogenisation simulations for binarymedia having the three different morphologies shown in Figure 4: Morphology 1 is a two-layered problem with arectangular inclusion; Morphology 2 is a binary medium with a circular inclusion; and Morphology 3 is a binarymedium with an L-shaped inclusion. All of the computations were carried out using MATLAB R2018a on a DELLdesktop with the configuration: Intel(R) Core(TM) i7-6700 CPU3.40GHz and RAM 16.0 GB. (a) Morphology 1: Rectangular inclusion(b) Morphology 2: Circular inclusion(c) Morphology 3: L-shaped inclusionFigure 4: Schematics of three binary media taken from [35], each with volumetric fraction for Ω taken as 0.25. The triangulation andcontrol volume partitions are also shown. (a) Morphology 1 Ω = [ , ] × [0 , = [0 , × [0 , \ Ω . (b) Morphology 2Ω = { ( x, y ) | ( x − . + ( y − . < π } , Ω = [0 , × [0 , \ Ω . (c) Morphology 3 Ω = [ , ] × [0 , ] ∪ [ , ] × [ , ] ,Ω = [0 , × [0 , \ Ω . To perform the homogenisation, we solve (31) with periodic boundary conditions and nonzero initial conditions.The effective diffusivity (32) is computed once steady-state is achieved. For the homogenisation problems consideredwe take D = D = 0 and denote D = ï D b x D b xy D b yx D b y ò . (32)In order to ensure the correct interfacial boundary conditions are invoked, we revert to the following time-fractionalmodel with the Riemann-Liouville fractional derivative: ® ∂u ( x,y,t ) ∂t = R D − γ t [ ∇ · ( D b ∇ u ( x, y, t ))] , ( x, y, t ) ∈ Ω × (0 , T ] , ∂v ( x,y,t ) ∂t = R D − γ t [ ∇ · ( D b ∇ v ( x, y, t ))] , ( x, y, t ) ∈ Ω × (0 , T ] , R D − γ t [ − D b ∇ u ( x, , t ) · n ] = R D − γ t [ D b ∇ u ( x, , t ) · n ] , R D − γ t [ − D b ∇ v (0 , y, t ) · n ] = R D − γ t [ D b ∇ v (1 , y, t ) · n ] , R D − γ t [ − D b ∇ v ( x, , t ) · n ] = R D − γ t [ D b ∇ v ( x, , t ) · n ] , and quasi-periodic boundary conditions v (1 , y, t ) = v (0 , y, t )+1 , v ( x, , t ) = v ( x, , t )+1, and with additional bound-ary conditions at the interface Γ are u ( x, y, t ) = v ( x, y, t ), R D − γ t ( D b ∇ u ( x, y, t ) · n ) = R D − γ t ( D b ∇ v ( x, y, t ) · n ).We use the time evolutionary homogenisation theory from Section 4 to estimate the effective parameters forthese three different morphologies as comparisons with the values of the diffusivity tensors can be made with thetest problems given in [35]. To calculate the equivalent diffusivity, we advance the numerical scheme in time untilthe steady-state is reached. We then calculate the equivalent diffusivity D b x in the x -direction and D b y in the y -direction using the following formulae obtained by approximating (32) in Section 4 using a midpoint quadraturerule applied to the integral over the unit cell D b x = (cid:80) i (cid:0) χD b ∂u∂x + (1 − χ ) D b ∂v∂x (cid:1) i S (cid:52) i (cid:80) i S (cid:52) i , D b y = (cid:80) i Ä χD b ∂u∂y + (1 − χ ) D b ∂v∂y ä i S (cid:52) i (cid:80) i S (cid:52) i , where S (cid:52) i is the area of the i th triangular element and the indicator function is defined as χ = ß , ( x, y ) ∈ Ω , , ( x, y ) ∈ Ω . Throughout this section we will frequently make reference to the harmonic average ( K ) and arithmetic average( K ), which can be calculated by K = (cid:16) (cid:15) D b + (cid:15) D b (cid:17) − , K = (cid:15) D b + (cid:15) D b , where (cid:15) is the volumetric fractionof Ω and (cid:15) is the volumetric fraction of Ω . Morphology 1 : Rectangular inclusion (layered composite material) .For Morphology 1, we are able to use the derived semi-analytic solution from Section 3.2 to validate the CVMmethod to simulate time-fractional transport in a layered medium with quasi-periodic boundary conditions. Thisproblem is chosen to validate the CVM method when applied to the homogenisation theory presented in Section 4.We exhibit the numerical solution at the central line y = 0 . u ( x, y,
0) = v ( x, y,
0) = u . The maximum error betweenthe numerical and semi-analytical solutions is presented in Table 2 for different γ and γ with h = 1 . × − , τ = 1 × − , D b = 10, D b = 1, u = q = 1 at t = 1. We conclude from the results in Table 2 and Figure 5 thatthe numerical and semi-analytical solutions agree very well, with the maximum error in all cases ≈ . × − . Table 2: The maximum error between the numerical solution and semi-analytical solution at the central line y = 0 . γ and γ with h = 1 . × − , τ = 1 × − , D b = 10, D b = 1, u = q = 1 at t = 1. γ =0.2 γ =0.5 γ =0.8 γ =0.2 3.2198E-04 3.2206E-04 3.2088E-04 γ =0.5 3.2002E-04 3.1897E-04 3.1818E-04 γ =0.8 3.2182E-04 3.1903E-04 3.1749E-04 γ =1.0 3.2632E-04 3.2194E-04 3.1917E-04Next, we investigate the impact of varying the fractional indices γ and γ on the solution profile at steady-state.At first, we fix γ = 1 to reduce the diffusion on Ω to normal diffusion and then change γ to observe the solutionbehaviour at t = 10 s (see Figure 5a). At this time, the diffusion for γ = 1 has reached its steady-state, whichcorresponds with the solution of the classical two-layered diffusion problem. The diffusion for γ = 0 . x = 0 .
5. 13 a) Different γ (b) Different t (c) Different γ Figure 5: A comparison between the numerical solution (symbol) and semi-analytical solution (line) for different γ , γ and t with D b = 10, D b = 1 and q = u = 1, in which other parameters are: (a) γ = 1 at t = 10 s ; (b) γ = 0 . γ = 1; (c) γ = 0 . t = 10 s . We conclude that the fractional index γ has a significant impact on the diffusion process for the two-layeredmedia with memory. Another interesting finding is that a larger value of the fractional index γ will delay thediffusion process from reaching steady-state. Furthermore, the final steady-state profiles for all the fractional casesare identical when γ = 1. To investigate this phenomenon further, we now fix γ = 1 and choose γ = 0 . t = 1 s to t = 10 s (see Figure 5b). We can see that the diffusion for γ = 0 . t = 10 s to reach its steady-state, which is much longer than required for the case γ = 0 . t = 10 s ). Weagain find that the steady-state solution profile for γ = 0 . γ = 0 . γ = 0 . γ on the diffusion, the reverse phenomenon to that observedabove is apparent (see Figure 5c). It appears that the larger the choice of γ , the quicker it reaches the steady-statesolution. (a) Different D b (b) Different q (c) Different u Figure 6: A comparison between the numerical solution (symbol) and semi-analytical solution (line) for different D b , q and u with γ = 0 . γ = 1, in which other parameters are: (a) D b = 1, q = u = 1 at t = 1 s ; (b) D b = 10, D b = 1, u = 1 at t = 10 s ;(c) D b = 10, D b = 1, q = 1 and t from 0 s to 10 s . In Figure 6, the effect of the other model parameters D b , D b , q and u on the diffusion process is exhibited.Figure 6a shows the impact of the diffusivity coefficient for the memory part ( D b ) with fixed D b = 1, γ = 0 . γ = 1, from which we see that the smaller the value of D b , the easier it becomes for the process to attain steadystate. Figures 6b and 6c show the solution profiles for different values of q and u for both the numerical andsemi-analytical solutions, which validates the accuracy of computational model and the mass balance equation (23).In particular, it is evident that as the value of q is decreased, the steady-state solution profile becomes flatter.Furthermore, for all choices of u , the average solution value remains constant at u during the simulation.We now give a comparison between the Riemann-Liouville (RL) and Caputo variants of the model discussed inSection 4: ß C D γ t u ( x, y, t ) = D b ∆ u ( x, y, t ) , ( x, y, t ) ∈ Ω × (0 , T ] , C D γ t v ( x, y, t ) = D b ∆ v ( x, y, t ) , ( x, y, t ) ∈ Ω × (0 , T ] , (33)14ith the initial and periodic boundary conditions for the Caputo model now expressed in terms of the classical flux u ( x, y,
0) = v ( x, y,
0) = u , D b ∂v (0 , y, t ) ∂x = D b ∂v (1 , y, t ) ∂x , v (1 , y, t ) = v (0 , y, t ) + q , and the boundary conditions at the interface Γ between the two media u ( x, y, t ) = v ( x, y, t ) , D b ∂u ( x, y, t ) ∂x = D b ∂v ( x, y, t ) ∂x , ( x, y, t ) ∈ Γ × (0 , T ] . (34) Figure 7: A comparison between the numerical solution (symbol) and semi-analytical solution (line) between the classical andfractional two-layered problems for the different time-fractional transport models, in which the the parameters are γ = 1, γ = 1 forthe classical case and γ = 0 . γ = 1 for the fractional case with D b = 10, D b = 1, q = u = 1 at t = 10 . As stated in Section 4, the fractional partial differential equations (24) and (25) are equivalent, however a mostimportant observation is that the interfacial boundary conditions (12) and (34) are not. Following a similar approachas outlined in Section 3.2, we have derived the semi-analytical solution for the Caputo model (33) and also exhibitedits behaviour along the central line y = 0 .
5. Figure 7 compares the steady-state solution profiles between the classicaland fractional two-layered problems for the two different time-fractional transport models. A key finding is thatfor the Caputo model (33), there is no observable difference in the steady-state profiles computed for the classicaland fractional cases whatever the choice of γ . This means that the impact of the internal layer memory effect isdiminished, which may not be a true accord of the physics associated with the anomalous transport phenomena.The main reason for this outcome is due to the incorrect treatment of the interfacial boundary conditions (34),which also explains why we need to transform the Caputo fractional derivative to the Riemann-Liouville fractionalderivative to perform the computational homogenisation. Otherwise, the interfacial boundary conditions (12) wouldneed to be imposed on the Caputo model (33), which is more complicated to calculate.We complete this section with a summary of our findings for the homogenised diffusivities computed using thetime-fractional homogenisation theory for the layered binary medium. To ease the mathematical exposition, wedenote α = 1 − γ and α = 1 − γ . At first we consider the case α = α = 0, which reduces the problem tothe classical two-layered problem. Figures 9a and 9b show the equivalent diffusivity, the harmonic average andarithmetic average curves of the two-layered medium with different diffusivity ratios D b D b . We can see that D b x and D b y coincide with the harmonic average and arithmetic average, respectively, which is a well known result. When D b D b = 10, D b x = 1 . D b y = 3 . D b D b = 0 . D b x = 0 . D b y = 0 . D b x for different α with increasing timeand α = 0 under different diffusivity ratios D b D b . An interesting observation from these figures is that for thistwo-layered composite material, the equivalent diffusivity tends to the diffusivity of the material with memory (i.e. D b x → D b = 10 when D b = 10 and D b x → D b = 0 . D b = 0 . D b x = 1 . . α has a significantimpact on the iterative process and it delays the diffusion to reach its steady-state (at least t > s ) significantly15 able 3: The relative error of the numerical results with the theoretical equivalent diffusivity in [35] for Morphology 1. D b D b D b x Theoretical Relative error D b y Theoretical Relative error10 1 1.2903 1.290 2.50E-04 3.2500 3.250 00.1 1 0.3077 0.309 4.23E-03 0.7750 0.775 0compared to the time for the classical case ( t < s ). We also notice that the final equivalent diffusivity D b x tendsto the same constant, which is independent of the fractional-order α for a fixed α = 0. (a) Morphology 1: D b D b = 10 (b) Morphology 1: D b D b = 0 . D b D b = 10 (d) Morphology 2: D b D b = 0 . D b x using different α for Morphologies 1 and 2 with increasing time underdifferent diffusivity ratios D b D b , where D b = 1 and α = 0. We believe these special properties for the composite material with memory are due to the impact of the fractionaloperators. In addition, it appears that the fractional-order α has no influence on the equivalent diffusivity D b y ,which may be due to the homogeneity in the y direction (see Figures 9a and 9b). Now, we fix D b = 1 and decreasethe value D b to observe the evolution of the equivalent diffusivity for the classical and fractional cases ( α = 0 . D b for both cases but they never meet (see Figure 11a). Figures 12a and 12b depict the steady-state homogenisedsolutions for the classical two-layered problem and the time-fractional two-layered problem. An interesting outcomeis that a very different diffusion profile is obtained for the case with memory. Again it can be concluded that thefractional-order has a significant effect on the diffusion process. Morphology 2 : Circular inclusion .At first, we consider the case α = α = 0. Figures 9c and 9d illustrate the equivalent diffusivity of the binarymedium under different diffusivity ratios D b D b . We determined that when D b D b = 10, D b = 1 . D b D b = 0 . D b = 0 . a) Morphology 1: D b D b = 10 (b) Morphology 1: D b D b = 0 . D b D b = 10 (d) Morphology 2: D b D b = 0 . D b D b = 10 (f) Morphology 3: D b D b = 0 . D b D b , in which we fix D b = 1. The equivalent diffusivity for this model is between the harmonic average and arithmetic average, which isdifferent to the two-layered problem considered on binary medium 1. In Table 4, the numerical equivalent diffusivityis compared with the theoretical value given in [35] using the initial condition u ( x, y,
0) = v ( x, y,
0) = 1. We cansee the numerical equivalent diffusivity agrees very well with the theoretical value.Next, we observe the case when Ω has memory, i.e., α (cid:54) = 0, α = 0 and the evolution of the equivalentdiffusivity D b for different α with increasing time under different diffusivity ratios D b D b is depicted in Figures 8c17nd 8d, respectively. When D b D b = 10, D b ≈ . D b D b = 0 . D b ≈ . D b value is notbounded between the harmonic average and arithmetic average when D b D b = 10 (see Figures 9c and 9d). Also, theequivalent diffusivity D b is independent of α with fixed α = 0. Table 4: The relative error of the numerical results with the theoretical equivalent diffusivity in [35] for Morphology 2. D b D b D b Theoretical Relative error10 1 1.5148 1.5200 3.43E-030.1 1 0.6602 0.6590 1.85E-03
Figure 10: The evolution of the equivalent diffusivity D b for different α for Morphology 3 with increasing time under differentmedium diffusivity ratios D b D b , where we fix the parameters D b = 1 and α = 0. The effects of the fraction-order α on the diffusion are similar to that reported for Morphology 1. Furthermore,we study the effects of the diffusivity ratio D b D b for fixed D b = 1 on the equivalent diffusivity and the related resultsare displayed in Figure 11b. It can be seen that, with a decreasing diffusivity for the memory part D b , the memoryeffect is weakened and the equivalent diffusivity D b coincides with the value for the classical case when D b is small( D b ≤ − ), which tends to a constant D b → . has a significant impact on the steady-state solution of the model.18 a) Morphology 1 (b) Morphology 2Figure 11: The effects of the diffusivity ratio D b D b for fixed D b = 1 on the equivalent diffusivity, where we use the parameter α = 0 . Morphology 3 : L-shaped inclusion .Again, we first consider the classical α = α = 0 case. Figures 9e and 9f illustrate the equivalent diffusivity forthe problem on binary medium 3 under different medium diffusivity ratios D b D b with initial condition is u ( x, y,
0) = v ( x, y,
0) = 1. The relative error between the numerical equivalent diffusivity and the theoretical diffusivity ispresented in Table 5.
Table 5: The relative error of the numerical results with the theoretical equivalent diffusivity in [35] for Morphology 3. D b D b D b x Theoretical Relative error10 1 1.4841 1.4800 2.75E-030.1 1 0.5328 0.5330 3.72E-04 D b D b D b y Theoretical Relative error10 1 1.8828 1.8800 1.49E-030.1 1 0.6758 0.6750 1.15E-03 D b D b D b xy Theoretical Relative error10 1 -7.9601E-02 -0.0796 1.80E-050.1 1 -2.8599E-02 -0.0286 4.51E-05When Ω has memory ( α (cid:54) = 0, α = 0), the evolution of the equivalent diffusivity D b for different α withincreasing time under different medium diffusivity ratios D b D b is displayed in Figure 10. We see that the equivalentdiffusivity tends to a fixed value at steady-state, which is consistent with the findings reported for Morphologies1 and 2. When D b D b = 10, D b x → . D b y → . D b xy → . D b D b = 0 . D b x → . D b y → . D b xy → − . D b x and D b y for theclassical case all lie between the harmonic average and arithmetic average, while when D b D b = 10, the values forthe fractional case are out of this bound (see Figures 9e and 9f). An interesting finding is that the sign of thecross-diffusion component D b xy has changed in the whole iterative process when D b D b = 10 (see Fig10). Finally, thecomparison of the steady-state in the x -direction and y -direction between the classical case and fractional case isshown in Figures 12d to 12h, from which some clear differences can be observed in the steady-state solution.In summary, the equivalent diffusivity for the classical two-layered medium can be bounded by the harmonicaverage and arithmetic average. However, the equivalent diffusivity for the binary medium with memory effect isnot bounded by these limits. For the fractional homogenisation problem, the equivalent diffusivity value tends toa fixed value after a long time iteration. In addition, the equivalent diffusivity deceases and recovers the classicalcase as the diffusivity coefficient of the memory part is decreased to zero.19 a) Morphology 1 (b) Morphology 1 (c) Morphology 2 (d) Morphology 2(e) Morphology 3 (f) Morphology 3 (g) Morphology 3 (h) Morphology 3Figure 12: The steady-state for classical and time-fractional case of the three different kinds of two-layered problems with D b D b = 10.
6. Application of the model to estimating the bound water diffusivity of wood
We now use the time-fractional transport model to perform computational homogenisation using different mi-croscale images of Spruce earlywood, Spruce latewood and an Australian hardwood (blackbutt) to estimate thebound water diffusivity when the cellular structure exhibits anomalous diffusion or memory effects. This type ofbehaviour has been observed experimentally in [6, 15, 16, 36]. The ESEM imaging of these wood species were per-formed at the Laboratoire G´enie des Proc´ed´es et Mat´eriaux, Centrale–Sup´elec in France and are shown in Figure13. Recent experimental work highlighted that traditional diffusion equations cannot adequately describe the ab-sorption of water in a cell wall [15, 36]. These effects are exacerbated even further for a modified cell wall occurringdue to heat treatment [16]. The reduction in diffusivity is thought to be due to the degradation of hemicelluloses.The thermal treatment leads to a change in the mechanism of attaining hygroscopic equilibrium by the modifiedwood. For these applications, traditional conservation laws prove inadequate for describing the underlying physicalprocesses, while models based on fractional operators succeed [6].Fractional-order derivatives provide excellent alternatives to their classical counterparts for such applications byinterpolating between the integer orders of differential equations to capture nonlocal relations in time using power-law memory kernels. A major attraction is that the dissipative and dispersive properties observed in experimentaldata are representative of a variety of different forms of generalised constitutive laws. Diffusion involving molecularrelaxation in the complex porous microstructure induces a second time constant in the macroscopic physical phe-nomena that gives rise to ‘sub-diffusive’ transport behaviour (see chapter 17 [37]), which motivates why we haveused a fractional-in-time derivative for the generalised conservation law. For the computational homogenisationsperformed over the given unit cell pore structures, the following modified transient diffusion equation can be usedto model the transport of water within the wood cell, which is an extension of the wood drying problem discussedin [38, 39] ∂ Ψ w ∂t + R D αt ( ∇ · Q w ) = 0 . Introducing the indicator variable χ = ß , ( x, y ) ∈ C ( g ) , , ( x, y ) ∈ C ( s ) , then we can define the fractional indexes and the mass fluxes as α = (1 − χ ) α + χα , Q w = (1 − χ ) Q ( s ) w + χ Q ( g ) w ,where Q ( s ) w = − ρ s D b ∇ X m , Q ( g ) w = − ρ g D v − ω v ∇ ω v , where X m is the cell wall bound water moisture content and ρ v is20he vapour density. The conserved quantity is defined as Ψ w = (1 − χ ) ρ s X m + χρ v . The quasi-periodic boundaryconditions are X m ( L x , y ) = X m (0 , y ) + 1 L x ∂X B ∂x , < y < L y ,X m ( x, L y ) = X m ( x,
0) + 1 L y ∂X B ∂y , < x < L x , R D αt [ Q w (0 , y ) · n ∂C ] = R D αt [ Q w ( L x , y ) · n ∂C ] , < y < L y , R D αt [ Q w ( x, · n ∂C ] = R D αt [ Q w ( x, L y ) · n ∂C ] , < x < L x , where n ∂C is the unit vector normal to ∂C outward to C , ∂X B ∂x and ∂X B ∂y are the gradients imposed over the cell,respectively. The closure conditions for the problem are P a + P v = P atm , ρ v = P v M v R ( T +273 . , ρ a = P a M a R ( T +273 . , wherethe parameters are calculated using the same formulae in [38] at T = 20 ℃ . (a) Spruce earlywood(b) Spruce latewood(c) Australian hardwoodFigure 13: The anatomical images of two different species of wood cells (100 µm × µm ) and their boundary contour andtriangulation, in which the volumetric fraction of solid part (Ω ) is approximately (a) 0.3427, (b) 0.8445, (c) 0.4566, respectively. Since wood has a periodic structure [40], a representative elementary volume (unit cell), needs to be chosen.We used ESEM images of the real wood cellular structure by extracting the boundary of the pores and then adigital representation of the cell was formed. Furthermore, we use the mesh generator
Gmsh [41] to perform thetriangulation, which can divide the cell domain [0 , L x ] × [0 , L y ] into two sub-sets: the cell lumens (Ω or C ( g ) ) andthe solid phase (Ω or C ( s ) ). Here, we refer the x -direction as the radial direction and refer the y -direction as thetangential direction. We denote the harmonic average as the series value and the arithmetic average as the parallelvalue to be consistent with [40]. For the treatment of the periodic boundary conditions, one can refer to [38]. Herewe will consider the problem on two different species of wood cells: softwood Spruce earlywood and latewood andAustralian hardwood (see Figure 13). 21 a) Spruce earlywood (b) Spruce latewood (c) Australian hardwoodFigure 14: The equivalent diffusivity for three different species of wood cells. Note that all the values have been scaled throughdivision by D v . Figure 14a shows the dimensionless equivalent diffusivity of Spruce earlywood for two different cases: the solidparts without memory (classical case) and with memory (fractional case). Since the solid cell walls have a lowdiffusivity of water, the equivalent diffusivity value is very close to the series model. As the cell walls in earlywoodare aligned with the radial direction, the value in the radial direction is larger than in the tangential direction. Forthe classical case, the radial value ( D b x ≈ . × − ) and tangential value ( D b y ≈ . × − ) are bounded bythe series and parallel values, which agrees with the results presented in [39, 40]. While for the fractional case (wherethe solid part has memory), the radial value ( D b x ≈ . × − ) and tangential value ( D b y ≈ . × − ) arenot bounded by the series and parallel values, being smaller than the classical case.In the following discussion, we consider the equivalent diffusivity of Spruce latewood cell problem. Figure14b displays the dimensionless equivalent diffusivity of Spruce latewood for two different cases: the solid partswithout memory (classical case) and with memory (fractional case). In contrast to Spruce earlywood, the tracheidsin latewood are more flattened in the tangential direction, which blocks the mass flux in the radial direction.An inverted anisotropic ratio is observed. For the classical case, the radial value is D b x ≈ . × − andtangential value is D b y ≈ . × − . For the fractional case (the solid part has memory), the radial value is D b x ≈ . × − and tangential value is D b y ≈ . × − , both of which are smaller than the classicalcase. Next, we consider the Australian hardwood case. The dimensionless equivalent diffusivity of hardwood forthe classical case and the fractional case is presented in Figure 14c, from which a similar result can be seen. Forthe classical case, the radial value ( D b x ≈ . × − ) and tangential value ( D b y ≈ . × − ) are boundedby the series and parallel values. While for the fractional case (the solid part has memory), the radial value( D b x ≈ . × − ) and tangential value ( D b y ≈ . × − ) are not bounded by the series and parallelvalues. (a) Spruce earlywood (b) Spruce latewood (c) Australian hardwoodFigure 15: The evolution of the equivalent diffusivity for three different species of wood cells. Note that all the values have been scaledthrough division by D v .
7. Conclusions
In this paper, we considered a two-dimensional time-fractional subdiffusion equation, in which the unstructuredmesh control volume method is applied. Our chosen numerical examples showed that the method was stableand effective. As an application, we successfully simulated a time-fractional transport model in a binary mediumconsisting of regular and irregular inclusions and derived a semi-analytical solution for a class of two-layered problemssubjected to quasi-periodic boundary conditions. An important contribution was the extension of the classicalhomogenisation theory to accommodate the new framework and to show that the effective diffusivity tensor can becomputed once the cell problems reach steady state. We found for all morphologies considered that the effectivediffusivity slowly converged to the diffusivity of the material in the unit cell with memory, which was a very differentfinding to that observed for homogenised diffusivity parameters for standard materials.In order to harness the main findings of our work, we postulate that the time evolutionary behaviour of the effec-tive parameters must be used in the homogenised macroscopic model until steady-state in the unit cell is achieved,at which time the homogenised effective parameter is dominated by the diffusivity of the material exhibiting memoryeffects. Furthermore, the fractional index used in the macroscopic time-fractional equation must be averaged overthe unit cell. In our future research, we plan to investigate this homogenisation strategy in more detail and verifyit using experimental data. We will also use the time-fractional multi-scale problem to simulate more complicatedheterogeneous systems.
Acknowledgment
This work was supported financially by a visiting Professorial Fellowship that enabled Turner to work at theUniversit´e Paris-Saclay, France for a period of three months in 2017-18. We acknowledge the financial support forthis research received through the Australian Research Council Discovery Grant DP150103675.
References [1] E. J. Carr, I. W. Turner, Dual-scale computational modelling of water flow in unsaturated soils containingirregular-shaped inclusions, Int. J. Numer. Methods Eng. 98(3) (2014) 157-173.[2] J. Lewandowska, A. Szymkiewicz, K. Burzynksi, M. Vauclin, Modeling of unsaturated water flow in double-porosity soils by the homogenization approach, Adv. Water Resour. 27 (2004) 283-296.[3] A. Szymkiewicz, J. Lewandowska, Micromechanical approach to unsaturated water flow in structured geo-materials by dual-scale computations, Acta Geotech. 3 (2008) 37-47.[4] G. Samaey, D. Roose, I.G. Kevrekidis, The gap-tooth scheme for homogenization problems, Multiscale Model.Simul. 4(1) (2005) 278-306. 235] E. J. Carr, P. Perr´e, I. W. Turner, The extended distributed microstructure model for gradient-driventransport: a two-scale model for bypassing effective parameters, J. Comput. Phys. 327 (2016) 810-829.[6] I. Turner, I. Ilic, and P. Perr´e, Modelling non-Fickian behavior in the cell walls of wood using a fractional-in-space diffusion equation, Drying Technology 29(16) (2011) 1932-1940.[7] A. Bueno-Orovio, D. Kay, V. Grau, B. Rodriguez, K. Burrage, Fractional diffusion models of cardiac electricalpropagation: role of structural heterogeneity in dispersion of repolarization, Journal of The Royal SocietyInterface, 11(97) (2014) 20140352.[8] P. Perr´e, Coupled heat and mass transfer in biosourced porous media without local equilibrium: A macro-scopic formulation tailored to computational simulation, International Journal of Heat and Mass Transfer140 (2019) 717-730.[9] S. A. Fomin, V. A. Chugunov, T. Hashida, Non-Fickian mass transport in fractured porous media, Advancesin Water Resources 34(2) (2011) 205-214.[10] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach,Physics Reports 339(1) (2000) 1-77.[11] S. Whitaker, Coupled transport in multiphase systems: a theory of drying, in: Y.I. Cho, J.P. Hartnett, T.F.Irvine, G.A. Greene (Eds.), Advances in Heat Transfer, vol. 31, Elsevier, 1998, pp. 1-104.[12] G. Allaire, Homogenization and two-scale convergence, SIAM Journal on Mathematical Analysis 23(6) (1992)1482-1518.[13] U. Hornung, Homogenization and porous media (Vol. 6), Springer Science & Business Media, 1996.[14] A. Bensoussan, J. L. Lions, G. Papanicolaou, Asymptotic analysis for periodic structures, American Math-ematical Society, Providence, Rhode Island, (2011).[15] W. Olek, P. Perr´e, J. Weres, Implementation of a relaxation equilibrium term in the convective boundarycondition for a better representation of the transient bound water diffusion in wood, Wood Sci. Technol. 45(2011) 677-691.[16] W. Olek, R. R´emond, J. Weres, P. Perr´e, Non-Fickian moisture diffusion in thermally modified beech woodanalyzed by the inverse method, International Journal of Thermal Sciences 109 (2016) 291-298.[17] S. B. Yuste, K. Lindenberg, Subdiffusion-limited A + A reactions, Physical Review Letters 87(11) (2001)118301.[18] R. Schumer, D. A. Benson, M. M. Meerschaert, B. Baeumer, Fractal mobile/immobile solute transport,Water Resources Research 39(10) (2003) 1296.[19] T. A. M. Langlands, Solution of a modified fractional diffusion equation, Physica A: Statistical Mechanicsand its Applications 367 (2006) 136-144.[20] B. I. Henry, T. A. M. Langlands, S. L. Wearne, Fractional cable models for spiny neuronal dendrites, PhysicalReview Letters 100(12) (2008) 128103.[21] A. Chang, H. Sun, Time-space fractional derivative models for CO2