From a cell model with active motion to a Hele-Shaw-like system. A numerical approach
Francisco Guillén-González, Juan Vicente Gutiérrez-Santacreu
FFROM A CELL MODEL WITH ACTIVE MOTION TO A HELE-SHAW-LIKESYSTEM. A NUMERICAL APPROACH
FRANCISCO GUILL´EN-GONZ ´ALEZ † AND JUAN VICENTE GUTI´ERREZ-SANTACREU ‡ Abstract.
In this paper we deal with the numerical solution of a Hele–Shaw-like system via acell model with active motion. Convergence of approximations is established for well-posed initialdata. These data are chosen in such a way the time derivate is positive at the initial time.The numerical method is constructed by means of a finite element procedure together withthe use of a closed-nodal integration. This gives rise to an algorithm which preserves positivitywhenever a right-angled triangulation is considered. As a result, uniform-in-time a priori estimatesare proven which allows us to pass to limit towards a solution to the Hele–Shaw problem.
Keywords.
Finite-element approximation; nonlinear diffusion; free boundary problems; Hele-Shaw flows.
Contents
1. Introduction 11.1. The models 11.2. Notation 31.3. Outline 32. Spatial discretization 42.1. Finite-element approximation 42.2. The numerical scheme 72.3. Main result 83. Proof of Theorem 2.2 83.1. A priori energy estimates 83.2. Passing to the limit 124. An algorithm on unstructured meshes 175. Numerical simulation 195.1. Temporal integration 195.2. Computational experiments 19References 221.
Introduction
The models.
Tumour cells are active mechanical systems that are able to produce forceswhich cause random migration [3, 8, 14]. This movement is due to rather complicate mechanismswhich occur inside cells and give rise to changes in cell shape. Another important mechanismunder which cells move is pressure [5, 8, 13] as a consequence of space competition generated by
Date : November 9, 2018.This work was partially supported by Ministerio de Econom´ıa y Competitividad under Spanish grant MTM2015-69875-P with the participation of FEDER. a r X i v : . [ m a t h . NA ] M a y F. GUILL´EN-GONZ ´ALEZ AND J. V. GUTI´ERREZ-SANTACREU cell proliferation itself. In the setting up we take into consideration a very simplified model whichincorporates the two spatial effects for describing tumour growth.Let Ω be a connected, open, bounded set of R d , with d = 2 or 3, and [0 , T ] a time interval.Consider the cell model with active motion [11] which consists in finding a tumour cell populationdensity n : Ω × [0 , T ] → R + satisfying(1) ∂ t n − ∇ · ( n ∇ p ( n )) − ν ∆ n = n G ( p ( n )) in Ω × (0 , T ) , subject to the (natural) boundary condition(2) ∇ n · n = 0 on ∂ Ω × (0 , T ) , with n being the outwards unit normal vector on the boundary ∂ Ω, and the initial condition(3) n | t =0 = n in Ω . Here p : [0 , + ∞ ) → [0 , + ∞ ) is defined by(4) p = p ( n ) := kk − n k − ∀ n ≥ , ( k ∈ N , k ≥ , and G = G ( p ) is a truncated decreasing function such that there exists P max > G (0) > , G ( p ) = 0 ∀ p ≥ P max > , and G (cid:48) ( p ) < ∀ p ∈ (0 , P max ) . In the above, G stands for the decrease in the tumuor cell growth rate when space is limited;the lack of space is governed by the local pressure p , the parameter P max is the maximum pressurethreshold that tumour cells can exceed before entering a quiescent state, and the parameter ν > p ( n ) given in (4) is invertible for n ≥ n ( p ) := (cid:18) k − k p (cid:19) / ( k − ∀ p ≥ . In this work we assume that { n k } k ∈ N is a sequence of initial data (3) for (1) such that(7) 0 ≤ p ( n k ) ≤ P max in Ω , and that there exists a limit function n ∞ such that(8) n k → n ∞ in L p (Ω)-strongly for any p < ∞ as k → ∞ .Consequently, defining N max ( k ) := n ( P max ) with n ( · ) being given in (6), we have(9) 0 ≤ n k ≤ N max ( k ) in Ω . from which we infer that there must exist N > N max ( k ) ≤ N . Under the aboveassumptions, equation (1) generates a sequence of solutions { n k } k ∈ N which lead to a solutiondescribing the dynamics of tumour growth as a free-boundary problem. To be more precise, theconvergence of the solutions { n k } k ∈ N of the active motion cell model problem (1)-(3) towards aweak solution to a Hele–Shaw-like system, as the parameter k goes to infinity, was proven in [11].This limit system reads as follows. Find n ∞ : Ω × [0 , T ] → R + and p ∞ : Ω × [0 , T ] → R + such that(10) ∂ t n ∞ − ∆ p ∞ − ν ∆ n ∞ = n ∞ G ( p ∞ ) in Ω × (0 , T ) , subject to(11) n ∞ | t =0 = n ∞ in Ω , (12) ∇ n ∞ · n = 0 and ∇ p ∞ · n = 0 on ∂ Ω × (0 , T ) , PPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 3 jointly to the complementary relation(13) p ∞ (∆ p ∞ + G ( p ∞ )) = 0 in Ω × (0 , T ) . The key point in establishing convergence is imposing that ∂ t n k (0) ≥
0. Moreover, equation (10)is equivalent to solving(14) ∂ t n ∞ − ∇ · ( n ∞ ∇ p ∞ ) − ν ∆ n ∞ = n ∞ G ( p ∞ ) in Ω × (0 , T ) . This equivalence will be accomplished due to the equality ∇ p ∞ = n ∞ ∇ p ∞ , which comes from theequalities p ∞ ∇ n ∞ = 0 and p ∞ n ∞ = p ∞ .In this paper, we shall be concerned with the convergence of a finite element scheme, the timevariable being continuous, for the active motion cell model problem (1)-(3) towards the Hele-Shawsystem (10)-(13) as the space discrete parameter h goes to zero and k goes to infinity.1.2. Notation.
We will assume the following notation throughout this paper. Let
O ⊂ R M , with M ≥
1, be a Lebesgue-measurable set and let 1 ≤ p ≤ ∞ . We denote by L p ( O ) the space ofall Lesbegue-measurable real-valued functions, f : O → R , being p th-summable in O for p < ∞ or essentially bounded for p = ∞ , and by (cid:107) f (cid:107) L p ( O ) its norm. When p = 2, the L ( O ) space is aHilbert space whose inner product is denoted by ( · , · ). To shorten the notation, the norm (cid:107) · (cid:107) L (Ω) is abbreviated by (cid:107) · (cid:107) .Let α = ( α , α , ..., α M ) ∈ N M be a multi-index with | α | = α + α + ... + α M , and let ∂ α be thedifferential operator such that ∂ α = (cid:16) ∂∂x (cid:17) α ... (cid:16) ∂∂x d (cid:17) α M . For m ≥ ≤ p ≤ ∞ , we define W m,p ( O ) to be the Sobolev space of all functions whose m derivatives are in L p ( O ), with the norm (cid:107) f (cid:107) W m,p ( O ) = (cid:88) | α |≤ m (cid:107) ∂ α f (cid:107) pL p ( O ) /p for 1 ≤ p < ∞ , (cid:107) f (cid:107) W m,p ( O ) = max | α |≤ m (cid:107) ∂ α f (cid:107) L ∞ (Ω) , for p = ∞ , where ∂ α is understood in the distributional sense. For p = 2, W m, ( O ) will be denoted by H m ( O ).We also consider C ∞ ( O ) to be the space of functions continuously differentiable any number oftimes, and C ∞ c ( O ) to be the subspace of C ∞ ( O ) with compact support in O .Spaces of Bochner-measurable functions from a time interval [0 , T ] to a Banach space X willbe denoted as L p (0 , T ; X ) with (cid:107) f (cid:107) L (0 ,T ; X ) = (cid:82) T (cid:107) f ( s ) (cid:107) pX d s if 1 ≤ p < ∞ or (cid:107) f (cid:107) L ∞ (0 ,T,X ) =ess sup s ∈ (0 ,T ) (cid:107) f ( s ) (cid:107) X < ∞ if p = ∞ .1.3. Outline.
Next we sketch the remaining content of this work. In section 2 we present ourfinite-element spaces and some preliminary result mainly concerning interpolation operators. Fur-thermore, we set out our finite element numerical method, where the time variable remains contin-uous, and the main result of this paper. Next is section 3 which is devoted to demonstrating themain result. Firstly, a discrete maximum principle for finite-element approximations is achievedby assuming a partition of the computational domain being made up of right-angled simplexes,and a priori estimates are also established independent of ( h, k ) with h being the space parameterassociated to our finite-element space. As a result, we are able to prove positivity for the time de-rivative of finite-element approximations. Then better a priori energy estimates lead to obtainingcompactness for passing to the limit as ( h, k ) → (0 , + ∞ ). In section 4, we propose a variant of ournumerical algorithm for nonobtuse triangulations which keeps with a discrete maximum principle F. GUILL´EN-GONZ ´ALEZ AND J. V. GUTI´ERREZ-SANTACREU and positive for the discrete time but whose convergence is not clear. Finally, in section 4, somenumerical experiments are presented for studying the behavior of several parameters.2.
Spatial discretization
Finite-element approximation.
Herein we introduce the hypotheses that will be requiredalong this work.(H1) Let Ω be a bounded domain of R d ( d = 2 or 3) with a polygonal or polyhedral Lipschitz-continuous boundary.(H2) Let {T h } h> be a family of shape-regular, quasi-uniform triangulations of Ω made up ofright-angled simplexes being triangles in two dimensions and tetrahedra in three dimen-sions, so that Ω = ∪ K ∈T h K , where h = max K ∈T h h K , with h K being the diameter of K .Further, let N h = { a i } i ∈ I denote the set of all the nodes of T h .(H3) Conforming piecewise linear, finite element spaces associated to T h are assumed for approx-imating H (Ω). Let P ( K ) be the set of linear polynomials on K ; the space of continuous,piecewise P ( K ) polynomial functions on T h is then denoted as N h = (cid:8) n h ∈ C (Ω) : n h | K ∈ P ( K ) ∀ K ∈ T h (cid:9) , whose Lagrange basis is denoted by { ϕ a } a ∈N h .We now give some auxiliary results for later use. We begin by an inverse inequality whose proofcan be found in [4, Lem. 4.5.3] or [9, Lem. 1.138]. Proposition 2.1.
Under hypotheses (H1) – (H3) , it follows that, (15) (cid:107)∇ n h (cid:107) L ( K ) ≤ C inv h − K (cid:107) n h (cid:107) L ( K ) ∀ K ∈ T h , ∀ n h ∈ N h , where C inv > is a constant independent of h . Let I h be the nodal interpolation operator from C (Ω) to N h and consider the discrete innerproduct ( n h , n h ) h = (cid:90) Ω I h ( n h n h ) = (cid:88) a ∈N h n h ( a ) n h ( a ) (cid:90) Ω ϕ a f o ∀ n h , n h ∈ N h , which induces the norm (cid:107) n h (cid:107) h = (cid:112) ( n h , n h ) h defined on N h . We recall the following local errorestimate. See [4, Thm. 4.4.4] or [9, Thm. 1.103] for a proof. Proposition 2.2.
Under hypotheses (H1) – (H3) , it follows that, (16) (cid:107) ϕ − I h ϕ (cid:107) L ∞ ( K ) ≤ C app h K (cid:107)∇ ϕ (cid:107) L ∞ ( K ) ∀ K ∈ T h , ∀ ϕ ∈ W , ∞ ( K ) , where C app > is independent of h . We next state the equivalence between the norms (cid:107) · (cid:107) h and (cid:107) · (cid:107) in N h and a discrete commuterapproximation property for I h . Proposition 2.3.
Under hypotheses (H1) – (H3) , it follows that, for all n h , n h ∈ N h , (17) (cid:107) n h (cid:107) ≤ (cid:107) n h (cid:107) h ≤ / (cid:107) n h (cid:107) and (18) (cid:107) n h n h − I h ( n h n h ) (cid:107) L (Ω) ≤ C app h (cid:107) n h (cid:107) (cid:107)∇ n h (cid:107) , where C app > is independent of h . PPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 5
Proof.
We have (cid:107) n h (cid:107) = (cid:88) a ∈N h n h ( a ) (cid:90) Ω ϕ a + (cid:88) a (cid:54) = (cid:101) a ∈N h n h ( a ) n h ( (cid:101) a ) (cid:90) Ω ϕ a ϕ (cid:101) a and (cid:107) n h (cid:107) h = (cid:88) a ∈N h n h ( a ) (cid:90) Ω ϕ a . Since 1 = (cid:80) (cid:101) a ∈N h ϕ (cid:101) a , we write (cid:107) n h (cid:107) h = (cid:88) a , (cid:101) a ∈N h n h ( a ) (cid:90) Ω ϕ a ϕ (cid:101) a = (cid:88) a ∈N h n h ( a ) (cid:90) Ω ϕ a + (cid:88) a (cid:54) = (cid:101) a ∈N h n h ( a ) (cid:90) Ω ϕ a ϕ (cid:101) a . Then (cid:107) n h (cid:107) h − (cid:107) n h (cid:107) = (cid:88) a > (cid:101) a ∈N h ( n h ( a ) + n h ( (cid:101) a ) − n h ( a ) n h ( (cid:101) a )) (cid:90) Ω ϕ a ϕ (cid:101) a = (cid:88) a > (cid:101) a ∈N h ( n h ( a ) − n h ( (cid:101) a )) (cid:90) Ω ϕ a ϕ (cid:101) a ≥ . From the above equality and Young’s inequality, we have (cid:107) n h (cid:107) h = (cid:107) n h (cid:107) + (cid:88) a > (cid:101) a ∈N h ( n h ( a ) − n h ( (cid:101) a )) (cid:90) Ω ϕ a ϕ (cid:101) a ≤ (cid:107) n h (cid:107) + 2 (cid:88) a > (cid:101) a ∈N h ( n h ( a ) + n h ( (cid:101) a )) (cid:90) Ω ϕ a ϕ (cid:101) a = (cid:107) n h (cid:107) + 2 (cid:88) a ∈N h n h ( a ) (cid:90) Ω ϕ a (cid:88) (cid:101) a (cid:101) a ϕ a ≤ (cid:107) n h (cid:107) + 4 (cid:107) n h (cid:107) ≤ (cid:107) n h (cid:107) . We now prove (18). By using (16), we obtain (cid:107)I h ( n h n h ) − n h n h (cid:107) L (Ω) = (cid:88) K ∈T h (cid:107)I h ( n h n h ) − n h n h (cid:107) L ∞ ( K ) (cid:90) K ≤ C app (cid:88) K ∈T h h K (cid:107)∇ ( n h n h ) (cid:107) L ∞ ( K ) (cid:90) K . Since n h , n h ∈ P ( K ) on K ∈ T h , we write ∇ ( n h n h ) = 2 d (cid:88) i,j =1 ∂ i n h ∂ j n h . Then, from (15) and on noting that ∇ n h , ∇ n h are piecewise constant on each K ∈ T h , we deducethat (cid:107)I h ( n h n h ) − n h n h (cid:107) L (Ω) ≤ C app (cid:88) K ∈T h h K (cid:107)∇ n h (cid:107) L ∞ ( K ) (cid:107)∇ n h (cid:107) L ∞ ( K ) (cid:90) K ≤ C app (cid:88) K ∈T h h K (cid:90) K |∇ n h | |∇ n h |≤ C app C inv (cid:88) K ∈T h h K (cid:107) n h (cid:107) L ( K ) (cid:107)∇ n h (cid:107) L ( K ) ≤ C app C inv h (cid:107) n h (cid:107) (cid:107)∇ n h (cid:107) , from which we conclude that (18) holds. (cid:3) F. GUILL´EN-GONZ ´ALEZ AND J. V. GUTI´ERREZ-SANTACREU
We will need to use an (average) interpolation operator into N h with the following properties.In particular we use an extension of the Scott-Zhang interpolation operator to L (Ω) function. Werefer to [15, 10] and [2]. Proposition 2.4.
Under hypotheses (H1) – (H3) , there exists an (average) interpolation operator Q h from L (Ω) to N h such that (19) (cid:107)Q h ψ (cid:107) W s,p (Ω) ≤ C sta (cid:107) ψ (cid:107) W s,p (Ω) for s = 0 , and ≤ p ≤ ∞ , (20) (cid:107)Q h ( ψ ) − ψ (cid:107) W s,p (Ω) ≤ C app h m − s (cid:107) ψ (cid:107) W m +1 ,p (Ω) for ≤ s ≤ m ≤ , and, for all ψ ∈ C ∞ (Ω) and n h ∈ N h , (21) (cid:107)Q h ( n h ψ ) − n h ψ (cid:107) W s,p (Ω) ≤ C app h m − s (cid:107) n h (cid:107) W m,p (Ω) (cid:107) ψ (cid:107) W m +1 , ∞ for ≤ s ≤ m ≤ . The key point in proving a discrete maximum principle is the following property which is ac-complished for right-angled simplexes assumed in (H2).
Proposition 2.5.
Under hypotheses (H1) – (H3) , it follows that, for any diagonal nonnegativematrix D = diag( d i ) di =1 (with d i ≥ ), (22) D ∇ ϕ a · ∇ ϕ (cid:101) a ≤ a.e. in Ω if a (cid:54) = (cid:101) a with a , (cid:101) a ∈ N h .Proof. For every right-angled d -simplex K ∈ T h of vertices { a i } i =0 ,...,d with a being the vertexsupporting the right angle, we denote by F a i the opposite face to a i and by n a i the exterior (tothe d -simplex K ) unit normal vector to the face F a i . Let (cid:98) K be the reference unit d -simplex withvertices (cid:98) a = and (cid:98) a i = e i , i = 1 , · · · , d , where { e i } i =1 , ··· ,d is the canonical basis of R d . Let F K bethe invertible affine mapping that maps (cid:98) K onto K defined by F K (cid:98) x = a + B K (cid:98) x , where B K ∈ R d × d is orthogonal.Let (cid:98) ϕ (cid:98) a i ( (cid:98) x ) = ϕ a i ( F K (cid:98) x ). Then we have (cid:98) ∇ (cid:98) ϕ (cid:98) a i = − d | (cid:98) F (cid:98) a i || (cid:98) K | n (cid:98) a i . In particular, n (cid:98) a i = − e i if i (cid:54) = 0 and n (cid:98) a = [1 , · · · , T . Thus, we obtain (cid:98) ∇ (cid:98) ϕ (cid:98) a i · (cid:98) ∇ (cid:98) ϕ (cid:98) a j = 1 d | (cid:98) F (cid:98) a i || (cid:98) F (cid:98) a j || (cid:98) K | n (cid:98) a i · n (cid:98) a j ≤ i (cid:54) = j. Therefore, by means of the change of variable x = a + B K (cid:98) x , it follows that ∇ ϕ a i = B K (cid:98) ∇ (cid:98) ϕ (cid:98) a i andhence D ∇ ϕ a i · ∇ ϕ a j = DB K (cid:98) ∇ (cid:98) ϕ (cid:98) a i · B K (cid:98) ∇ (cid:98) ϕ (cid:98) a j = 1 d | (cid:98) F (cid:98) a i || (cid:98) F (cid:98) a j || (cid:98) K | n T (cid:98) a i B TK DB K n (cid:98) a j ≤ i (cid:54) = j because, since B K is a orthogonal matrix, the inner products defined by D and B TK DB K preservesangles. (cid:3) Remark 2.1.
When D = I d with I d being the d × d identity matrix, property (22) can be proved fornonobtuse triangulations [7] . Then property (22) can be somewhat seen a generalization restrictedfor right-angled triangulations. PPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 7
Let us now introduce the discrete Laplacian associated to the mass-lumping scalar product( · , · ) h . For any Σ h ∈ N h , let − (cid:101) ∆ h Σ h ∈ N h solve(23) − ( (cid:101) ∆ h Σ h , n h ) h = ( ∇ Σ h , ∇ n h ) ∀ n h ∈ N h . We end up with a compactness result [1, Lm. 2.4] needed in proving the equivalence betweenproblems (10) and (14).
Theorem 2.1.
Assume that (H1) - (H3) holds. Let dd +2 < (cid:96) < ∞ . Suppose that { ρ h,k } h,k ≥ ⊂ L (0 , T ; L (Ω)) is such that ρ h,k ( t, · ) ∈ N h for all t ∈ [0 , T ] and satisfies (cid:107) ρ h,k (cid:107) H (0 ,T ; L (cid:96) (Ω)) + (cid:107) ρ h,k (cid:107) L ∞ (0 ,T ; L (Ω)) ∩ L (0 ,T ; H (Ω)) + (cid:107) (cid:101) ∆ h ρ h,k (cid:107) L (0 ,T ; L (Ω)) ≤ C dat . Then there exist a subsequence { ρ h,k } h,k> (not relabeled) and a limit function ρ , such that ρ h,k → ρ in L (0 , T, H (Ω)) -strongly as ( h, k ) → (0 , + ∞ ) . Hereafter C will denote a generic constant whose value may change at each occurrence. Thisconstant may depend on the data problem and the constants C inv , C app , C com and C dat .2.2. The numerical scheme.
In order to avoid dense technical calculations, we assume for sim-plicity that each element K ∈ T h has its edges lined up with the axes.The numerical scheme relies on a finite-element method combined with a closed-nodal integrationapplied to the time-derivative and pressure-migration terms. Thus our numerical method whichconsists in finding n h,k ∈ C ([0 , T ]; N h ) such that(24) (cid:26) ( ∂ t n h,k , n h ) h + ( ∇I h (( n h,k ) k ) , ∇ n h ) + ν ( ∇ n h,k , ∇ n h ) = ( G ( p ( n h,k )) n h,k , n h ) h ∀ n ∈ N h n h,k (0) = n h,k , with p ( n h,k ) = kk − n h,k ) k − .Equivalently, we may write (24) as(25) ( ∂ t n h,k , n h ) h + ( D ( n h,k ) ∇ n h,k , ∇ n h ) + ν ( ∇ n h,k , ∇ n h ) = ( G ( p ( n h,k )) n h,k , n h ) h , where D ( n h,k ) is a piecewise constant, d × d diagonal matrix function with respect to T h definedas follows. Let K ∈ T h with vertices { a i } i =0 , ··· ,d where a corresponds to the right angle. Then(26) [ D ( n h,k ) | K ] ii = ( n h,k ) k ( a i ) − ( n h,k ) k ( a ) n h,k ( a i ) − n h,k ( a ) if n h,k ( a i ) − n h,k ( a ) (cid:54) = 0 , n h,k ( a i ) − n h,k ( a ) = 0 . By the mean value theorem, one can write(27) [ D ( n h,k ) | K ] ii = k ( n h,k ) k − ( ξ i ) , where ξ i = α a i + (1 − α ) a for a certain α ∈ (0 , { n h,k } h,k> is as follows. Let { n k } k ∈ N ⊂ H (Ω) ∩ L ∞ (Ω)satisfy (7) and (9). Then we select n h,k = Q h ( n k ) so that(28) 0 ≤ n h,k ( a ) ≤ N max ( k ) ∀ a ∈ N h , (cid:107)∇ n h,k (cid:107) ≤ C stab (cid:107)∇ n k (cid:107) , (29) n h,k → n k in H (Ω)-strongly as h → { n h,k } h,k> to be such that(30) − ( ∇I h ( n h,k ) k , ∇ n h ) − ν ( ∇ n h,k , ∇ n h ) + ( G ( p ( n h,k )) n h,k , n h ) h ≥ ∀ n h ∈ N h with n h ≥ . F. GUILL´EN-GONZ ´ALEZ AND J. V. GUTI´ERREZ-SANTACREU
Remark 2.2.
This last condition is related to imposing ∂ t n h,k (0) ≥ which is crucial to prove the k → + ∞ limit. The existence and uniqueness of a solution to scheme (24) may be readily justified by Picard’stheorem. To be more precise, one may prove that there exists a time interval [0 , T h ) for whichproblem (24) is uniquely solvable. As a consequence of a priori energy estimates, which we shallprove in the next section, one deduces that T h = T for all h > Main result.
We now are ready to state our main result of this paper. We shall prove thatscheme (24) produces a sequence of discrete solutions which satifies a priori energy bounds uniformwith respect to ( h, k ) allowing us to pass to the limit as ( h, k ) → (0 , + ∞ ) towards weak solutionsof the Hele–Shaw-like system (10)-(13). Theorem 2.2.
Assume that (H1)-(H3) hold. Then the discrete solution { ( n h,k , p h,k ) } h,k of (24) satisfies the following estimates, for all a ∈ N h and t ∈ [0 , T ] : ≤ n h,k ( a , t ) ≤ N max ( k ) , ≤ p ( n h,k ( a , t )) ≤ P max ,∂ t n h,k ( a , t ) ≥ , ∂ t p ( n h,k ( a , t )) ≥ . Furthermore, { n h,k , I h (( n h,k ) k ) } h,k converges towards weak solutions ( n ∞ , p ∞ ) of problem (10) - (13) in the sense that n h,k → n ∞ in L ∞ (0 , T ; H (Ω)) -weakly- (cid:63) and in L p ((0 , T ) × Ω) -strongly , and I h (( n h,k ) k ) → p ∞ in L ∞ (0 , T ; H (Ω)) -weakly- (cid:63) and in L p ((0 , T ) × Ω) -strongly , for any < p < ∞ provided that (H5) k h → as ( h, k ) → (0 , + ∞ ) . Proof of Theorem 2.2
A priori energy estimates.
Our goal is to prove a priori energy estimates for the discretesolution n h,k of (24) independent of ( h, k ).This first lemma will be focused on proving a discrete maximum principle for n h,k based on thehypothesis of right-angled triangulations. Moreover, some a priori energy estimates are obtained. Lemma 3.1.
Assume that (H1)-(H3) hold. Then the solution n h,k of scheme (24) satisfies (31) 0 ≤ n h,k ( a , t ) ≤ N max ( k ) ∀ a ∈ N h and ∀ t ≥ , and (32) (cid:107) n h,k (cid:107) L ∞ (0 ,T ; L (Ω)) + (cid:107) n h,k (cid:107) L (0 ,T ; H (Ω)) ≤ C, where C > is independent of ( h, k ) .Proof. We first proceed to verify (31). In doing so, we introduce a modification to scheme (25)which truncates the nonlinear diffusion term as follows:(33) ( ∂ t n h,k , n h ) h + ( D ([ n h,k ] T ) ∇ n h,k , ∇ n h ) + ν ( ∇ n h,k , ∇ n h ) = ( G ( p ([ n h,k ] T )) n h,k , n h ) h , where [ n h,k ] T is the usual truncation of n h,k from below by 0 and from above by N max ( k ). Again,by means of Picard’s theorem, one has the existence and uniqueness of a solution n h,k to (33).Let n min h,k = I h ( n − h,k ) ∈ N h be defined as n min h,k = (cid:88) a ∈N h n − h,k ( a ) ϕ a , PPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 9 where n − h,k ( a ) = min { , n h,k ( a ) } . Analogously, one defines n max h,k = I h ( n + h,k ) ∈ N h as n max h,k = (cid:88) a ∈N h n + h,k ( a ) ϕ a , where n + h,k ( a ) = max { , n h,k ( a ) } . Notice that n h,k = n min h,k + n max h,k .On choosing n h = n min h,k in (33), it follows that(34) 12 ddt (cid:107) n min h,k (cid:107) h + ( D ([ n h,k ] T ) ∇ n h,k , ∇ n min h,k ) + ν ( ∇ n h,k , ∇ n min h,k ) = (cid:107) G ( p ([ n h,k ] T )) / n min h,k (cid:107) h ≤ G (0) (cid:107) n min h,k (cid:107) h . Next observe that( D ([ n h,k ] T ) ∇ n h,k , ∇ n min h,k ) = ( D ([ n h,k ] T ) ∇ n min h,k , ∇ n min h,k ) + ( D ([ n h,k ] T ) ∇ n max h,k , ∇ n min h,k )= (cid:107)D ([ n h,k ] T ) / ∇ n min h,k (cid:107) + (cid:88) a (cid:54) =˜ a ∈N h n − h,k ( a ) n + h,k ( (cid:101) a )( D ([ n h,k ] T ) ∇ ϕ a , ∇ ϕ (cid:101) a ) . Then, using the fact that n − h,k ( a ) n + h,k ( (cid:101) a ) ≤ a (cid:54) = ˜ a and that D ([ n h,k ] T ) is a nonnegative diagonalmatrix function, one deduces, from (22), that D ([ n h,k ] T ) ∇ ϕ a · ∇ ϕ (cid:101) a ≤ ∀ a (cid:54) = (cid:101) a ∈ N h and thereby(35) ( D ([ n h,k ] T ) ∇ n h,k , ∇ n min h,k ) ≥ (cid:107)D ([ n h,k ] T ) / ∇ n min h,k (cid:107) . Analogously, one obtains(36) ν ( ∇ n h,k , ∇ n min h,k ) ≥ ν (cid:107)∇ n min h,k (cid:107) , where we have used again (22) but now for D = I d , with I d being the d × d unit matrix. Inserting(35) and (36) into (34) yields12 ddt (cid:107) n min h,k (cid:107) h + (cid:107)D ([ n h,k ] T ) / ∇ n min h,k (cid:107) + ν (cid:107)∇ n min h,k (cid:107) ≤ G (0) (cid:107) n min h,k (cid:107) h . By Gr¨onwall’s lemma, we have n min h,k ( t ) ≡ t ≥
0, since n min h,k (0) ≡ ≤ n h,k in (31). For the other inequality n h,k ≤ N max ( k ) in (31), we proceed in a similarfashion. In this case, one chooses n h = ( n h,k − N max ( k )) max in (33) and takes into account that G ( p ([ n h,k ] T )) n h,k ( n h,k − N max ( k )) max ≡ p ([ n h,k ] T ) = P max if n h,k ≥ N max ( k ).It should be noted that any solution n h,k of the modified scheme (33) satisfies the discretemaximum principle (31), and consequently [ n h,k ] T ≡ n h,k ; hence n h,k satisfies the non-truncatedscheme (24) as well. Finally, by uniqueness of solutions for scheme (24), the solution of (24) takesvalues between 0 and N max ( k ); that is (31).Now selecting n h = n h,k in (25) and invoking Gr¨onwall’s lemma, the following energy estimateholds, for all t ∈ [0 , T ]:(37) 12 (cid:107) n h,k ( t ) (cid:107) h + (cid:90) T (cid:107)D ( n h,k ) / ∇ n h,k (cid:107) + ν (cid:90) T (cid:107)∇ n h,k (cid:107) ≤ exp(2 G (0) T ) 12 (cid:107) n h,k (cid:107) h . Then the weak estimates (32) are deduced from (37) and (17). (cid:3)
A discrete maximum principle for ( n h,k ) k − and ( n h,k ) k follows as a direct consequence of (31). Corollary 3.1.
There holds (38) 0 ≤ ( n h,k ) k − ( a , t ) ≤ P max ∀ a ∈ N h and ∀ t ≥ . and (39) 0 ≤ ( n h,k ) k ( a , t ) ≤ P max N max ( k ) ∀ a ∈ N h and ∀ t ≥ . Proof.
Assertions (38) and (39) are satisfied in view of (31) and the bounds n k − h,k ( a , t ) ≤ N max ( k ) k − = k − k P max ≤ P max and n kh,k ( a , t ) ≤ N max ( k ) k = N max ( k ) k − N max ( k ) ≤ P max N max ( k ) . (cid:3) The following lemma provides the positivity and some a priori estimates for the time derivativeof n h,k and ( n h,k ) k . Lemma 3.2.
Suppose that (H1)-(H4) hold. Then it follows that (40) ∂ t n h,k ( a , t ) ≥ , ∂ t ( n h,k ( a , t )) k ≥ ∀ a ∈ N h and ∀ t ∈ [0 , T ] , and the a priori estimates (41) (cid:107) ∂ t n h,k (cid:107) L ∞ (0 ,T ; L (Ω)) ≤ C, (42) (cid:107) ∂ t ( n h,k ) k (cid:107) L (0 ,T ; L (Ω)) ≤ C, where C > is a constant independent of ( h, k ) .Proof. Let us define Σ( n h,k ) ∈ N h such thatΣ( n h,k ) = I h (( n h,k ) k ) + ν n h,k = I h (( n h,k ) k + ν n h,k ) . Moreover, let Σ (cid:48) ( n h,k ) ∈ N h and Σ (cid:48)(cid:48) ( n h,k ) ∈ N h be defined asΣ (cid:48) ( n h,k ) = k I h (( n h,k ) k − ) + ν and Σ (cid:48)(cid:48) ( n h,k ) = k ( k − I h (( n h,k ) k − ) . Then scheme (24) can be rewritten as( ∂ t n h,k , n h ) h + ( ∇ Σ( n h,k ) , ∇ n h ) = ( G ( p ( n h,k )) n h,k , n h ) h , and equivalently, from (23), as(43) ( ∂ t n h,k , n h ) h − ( (cid:101) ∆ h Σ( n h,k ) , n h ) h = ( G ( p ( n h,k )) n h,k , n h ) h . Now take n h = I h (Σ (cid:48) ( n h,k ) w h ), for any w h ∈ N h to get( ∂ t Σ( n h,k ) , w h ) h − (Σ (cid:48) ( n h,k ) (cid:101) ∆ h Σ( n h,k ) , w h ) h = (Σ (cid:48) ( n h,k ) G ( p ( n h,k )) n h,k , w h ) h . Differentiating with respect to time and defining w h,k ∈ N h such that, for each a ∈ N h and t ∈ [0 , T ], w h,k ( a , t ) := ∂ t Σ( n h,k )( a , t ) = Σ (cid:48) ( n h,k )( a , t ) ∂ t n h,k ( a , t ) , one arrives at( ∂ t w h,k , w h ) h − (Σ (cid:48) ( n h,k ) (cid:101) ∆ h w h,k , w h ) h = (Σ (cid:48)(cid:48) ( n h,k ) ∂ t n h,k (cid:101) ∆ h Σ( n h,k ) , w h ) h +(Σ (cid:48)(cid:48) ( n h,k ) ∂ t n h,k G ( p ( n h,k )) n h,k , w h ) h + k (Σ (cid:48) ( n h,k ) G (cid:48) ( p ( n h,k ))( n h,k ) k − ∂ t n h,k , w h ) h +(Σ (cid:48) ( n h,k ) G ( p ( n h,k )) ∂ t n h,k , w h ) h , PPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 11 for any w h ∈ N h . Since w h,k ( a , t ) = Σ (cid:48) ( n h,k )( a , t ) ∂ t n h,k ( a , t ) and Σ (cid:48) ( n h,k )( a , t ) ≥ ν >
0, we have ∂ t n h,k ( a , t ) = w h,k ( a , t )Σ (cid:48) ( n h,k )( a , t ) ∀ a ∈ N h ∀ t ∈ [0 , T ] . Both previous equalities yield( ∂ t w h,k , w h ) h − (Σ (cid:48) ( n h,k ) (cid:101) ∆ h w h,k , w h ) h = ( F ( n h,k ) w h,k , w h ) h , for any w h ∈ N h , where F ( n h,k ) := Σ (cid:48)(cid:48) ( n h,k )Σ (cid:48) ( n h,k ) (cid:110) (cid:101) ∆ h Σ( n h,k ) + n h,k G ( p ( n h,k )) (cid:111) + k ( n h,k ) k − G (cid:48) ( p ( n h,k )) + G ( p ( n h,k )) . Taking w h = w min h,k = I h ( w − h,k ) in the above variational formulation, we get(44) 12 ddt (cid:107) w min h,k (cid:107) h − (Σ (cid:48) ( n h,k ) (cid:101) ∆ h w h,k , w min h,k ) h ≤ (cid:107) F ( n h,k ) (cid:107) L ∞ (cid:107) w min h,k (cid:107) h . Since n h,k ∈ C ([0 , T ]; N h ) and N h is a finite dimensional space, we have that (cid:107) F ( n h,k )( t ) (cid:107) L ∞ (Ω) ≤ C h,k for all t ∈ [0 , T ], where C h,k > h and k . It should also be noted that − (Σ (cid:48) ( n h,k ) (cid:101) ∆ h w h,k , w min h,k ) h ≥
0. Indeed, choose n h = ϕ a in (23) to obtain − ( (cid:101) ∆ h w h,k )( a ) (cid:90) Ω ϕ a = ( ∇ w h,k , ∇ ϕ a ) . Then − (Σ (cid:48) ( n h,k ) (cid:101) ∆ h w h,k , w min h,k ) h = − (cid:88) a ∈N h Σ (cid:48) ( n h,k ( a ))( (cid:101) ∆ h w h,k )( a ) w min h,k ( a ) (cid:90) Ω ϕ a = (cid:88) a ∈N h Σ (cid:48) ( n h,k ( a ))( ∇ w h,k , ∇ ϕ a ) w min h,k ( a )= (cid:88) a ∈N h Σ (cid:48) ( n h,k ( a ))( ∇ w max h,k , ∇ ϕ a ) w min h,k ( a )+ (cid:88) a ∈N h Σ (cid:48) ( n h,k ( a ))( ∇ w min h,k , ∇ ϕ a ) w min h,k ( a ) . Therefore, using the fact that Σ (cid:48) ( n h,k ) ≥ ν >
0, we obtain (cid:88) a ∈N h Σ (cid:48) ( n h,k ( a ))( ∇ w max h,k , ∇ ϕ a ) w min h,k ( a ) = (cid:88) a (cid:54) =˜ a ∈N h Σ (cid:48) ( n h,k ( a )) w max h,k (˜ a ) w min h,k ( a )( ∇ ϕ ˜ a , ∇ ϕ a ) ≥ (cid:88) a ∈N h Σ (cid:48) ( n h,k ( a ))( ∇ w min h,k , ∇ ϕ a ) w min h,k ( a ) ≥ ν (cid:107)∇ w min h,k (cid:107) ≥ . Thus, (44) leads to 12 ddt (cid:107) w min h,k (cid:107) h ≤ (cid:107) F ( n h,k ) (cid:107) L ∞ (cid:107) w min h,k (cid:107) h and hence, by Gr¨onwall’s lemma, (cid:107) w min h,k ( t ) (cid:107) h ≤ exp(2 T (cid:107) F ( n h,k ) (cid:107) L ∞ ) (cid:107) w min h,k (0) (cid:107) ∀ t ∈ [0 , T ] . From (30) in (H4), we deduce that w h,k ( a ,
0) = Σ (cid:48) ( n h,k )( a , ∂ t n h,k ( a , ≥ w min h,k ( t ) ≡ w min h,k (0) ≡
0. As a result, we have that ∂ t n h,k ≥ ∂ t ( n h,k ) k = k ( n h,k ) k − ∂ t n h,k ≥
0. Thus, (40) is true.
Now we are going to obtain bounds (41) and (42). For this, we take n h = 1 in (25) and use (40)to have (cid:107) ∂ t n h,k (cid:107) L (Ω) = ( ∂ t n h,k ,
1) = ( ∂ t n h,k , h ≤ G (0) (cid:107) n h,k (cid:107) L (Ω) ≤ G (0) | Ω | N max ( k );hence estimate (41) holds. Furthermore, we have, by (39) and (40), that (cid:107) ∂ t ( n h,k ) k (cid:107) L (0 ,T ; L (Ω)) = (cid:90) T ddt (( n h,k ) k , dt = (( n h,k ) k ( T ) − ( n h,k ) k (0) , ≤ | Ω | N max ( k ) P max ;hence estimate (42) holds. (cid:3) We are now concerned with an a priori estimate for the gradient of n h,k and I h (( n h,k ) k ). Theseestimates will play an important role in obtaining compactness results which allow us to pass to thelimit as ( h, k ) → (0 , + ∞ ) from scheme (24) towards weak solutions ( n ∞ , p ∞ ) of problem (10)-(13). Lemma 3.3.
Suppose that (H1) - (H4) are satisfied. Then there exists a constant C > , indepen-dent of h and k , such that (45) (cid:107)D ( n h,k ) / ∇ n h,k (cid:107) L ∞ (0 ,T ; L (Ω)) + (cid:107)∇ n h,k (cid:107) L ∞ (0 ,T ; L (Ω)) ≤ C and (46) (cid:107)∇I h (( n h,k ) k ) (cid:107) L ∞ (0 ,T ; L (Ω)) ≤ C. Proof.
Select n h = n h,k ∈ N h in (24) to obtain( ∂ t n h,k , n h,k ) h + ( D ( n h,k ) ∇ n h,k , ∇ n h,k ) + ν (cid:107)∇ n h,k (cid:107) = (cid:107) G ( p ( n h,k )) / n h,k (cid:107) h ≤ G (0) (cid:107) n h,k (cid:107) h . From (31) and (40), we deduce that ( ∂ t n h,k , n h,k ) h ≥
0. Therefore, (cid:107)D ( n h,k ) / ∇ n h,k (cid:107) + ν (cid:107)∇ n h,k (cid:107) ≤ G (0) (cid:107) n h,k (cid:107) h . This last expression combined with (32) gives (45).Take n h = I h (( n h,k ) k ) in (24) to have( ∂ t n h,k , I h (( n h,k ) k )) h + (cid:107)∇I h (( n h,k ) k ) (cid:107) + ν ( D (( n h,k ) k ) ∇ n h,k , ∇ n h,k )= ( G ( p ( n h,k )) n h,k , I h (( n h,k ) k )) h ≤ G (0) (cid:107) ( n h,k ) k − (cid:107) L ∞ (Ω) (cid:107) n h,k (cid:107) h ≤ G (0) P max (cid:107) n h,k (cid:107) h . From this, it follows that (46) holds from (32), (40) and from noting that ( D (( n h,k ) k ) ∇ n h,k , ∇ n h,k ) ≥ (cid:3) Passing to the limit.
From estimates (31) and (45) jointly with (39) and (46), we have thatthere exist two limit functions ( n ∞ , p ∞ ) ∈ L ∞ (0 , T ; H (Ω)) and a subsequence of { ( n h,k , I h (( n h,k ) k ) } h,k ,which we still denote in the same way, such that the following convergences hold, as ( h, k ) → (0 , ∞ ):(47) n h,k → n ∞ in L ∞ (0 , T ; H (Ω) ∩ L ∞ (Ω))-weakly- (cid:63), and(48) I h (( n h,k ) k ) → p ∞ in L ∞ (0 , T ; H (Ω) ∩ L ∞ (Ω))-weakly- (cid:63). Before proceeding to pass to the limit, we need to obtain some strong convergences via an Aubin-Lions campactness lemma [16]. From (31), (41) and (45), we have that there exists a subsequence(not relabeled) such that, as ( h, k ) → (0 , ∞ ),(49) n h,k → n ∞ in L p (Ω × (0 , T ))-strongly, ∀ p < ∞ , and(50) n h,k → n ∞ in C ([0 , T ]; L q (Ω))-strongly, ∀ q < ∗ , PPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 13 where 2 ∗ stands for the conjugate exponent of 2 defined by 1 / ∗ = 1 / − /d . Analogously, from(39), (42), and (46), we have(51) I h (( n h,k ) k ) → p ∞ in L p (Ω × (0 , T ))-strongly, ∀ p < ∞ . As a result, we also have the strong convergence of p ( n h,k ) towards p ∞ , but under hypothesis(H5) in Theorem 2.2. Lemma 3.4.
Assuming hypotheses (H1) - (H5) , it follows that, as ( h, k ) → (0 , ∞ ) , (52) p ( n h,k ) → p ∞ in L p ((0 , T ) × Ω) -strongly for any p < ∞ . Moreover, (53) p ∞ n ∞ ≡ p ∞ a.e. in (0 , T ) × Ω .Proof. For each element K ∈ T h with vertices { a , · · · a d } , we associate once and for all a vertex a K of K. Thus we define a piecewise constant function P h ( n kh,k )( x ) = n kh,k ( a K ) for all x ∈ K ,which satisfies P h ( n kh,k )( x ) − n kh,k ( x ) = ∇ ( n kh,k ( ξ a K )) · ( a K − x ) = k n k − h,k ( ξ a K ) ∇ n h,k | K · ( a K − x )where ξ a K = λ a K + (1 − λ ) x with λ ∈ (0 , (cid:107)P h ( n kh,k ) − n kh,k (cid:107) L ∞ (0 ,T ; L (Ω)) ≤ C k h (cid:107) n k − h,k (cid:107) L ∞ (0 ,T ; L ∞ (Ω)) (cid:107)∇ n h,k (cid:107) L ∞ (0 ,T ; L (Ω)) ≤ C k h P max . The above argument also shows by replacing n kh,k by I h ( n kh,k ) and using (46) that (cid:107)I h ( n kh,k ) − P h ( n kh,k ) (cid:107) L ∞ (0 ,T ; L (Ω)) ≤ C h (cid:107)∇I h ( n kh,k ) (cid:107) L ∞ (0 ,T ; L (Ω)) ≤ C h.
Thus, by (51) and (H5), we deduce, the following convergence, as ( h, k ) → (0 , ∞ ):(54) n kh,k → p ∞ in L p ((0 , T ) × Ω)-strongly ∀ p < ∞ . In view of (49) and (54), there is a subsequence (not relabeled) of { ( n h,k , n kh,k ) } h,k such that, as( h, k ) → (0 , ∞ ):( n h,k ( x , t ) , n kh,k ( x , t )) → ( n ∞ ( x , t ) , p ∞ ( x , t )) a.e. ( x , t ) ∈ Ω × (0 , T ).Thus, defining (cid:101) p ∞ ( x , t ) = p ∞ ( x , t ) n ∞ ( x , t ) if n ∞ ( x , t ) (cid:54) = 0,0 otherwise,it follows that, as ( h, k ) → (0 , ∞ ), p ( n h,k ( x , t )) = kk − n kh,k ( x , t ) n h,k ( x , t ) → (cid:101) p ∞ ( x , t ) a.e. ( x , t ) ∈ Ω × (0 , T );furthermore, p ∞ ( x , t ) ← kk − n kh,k ( x , t ) = (cid:16) − k (cid:17) k − p ( n h,k ( x , t )) kk − → (cid:101) p ∞ ( x , t ) . Thus, p ∞ ≡ (cid:101) p ∞ a.e. ( x , t ) ∈ Ω × (0 , T ) and, in particular, one has equality (53) and the pointwiseconvergence p ( n h,k ( x , t )) → p ∞ ( x , t ) a.e. ( x , t ) ∈ Ω × (0 , T ).Finally, (52) is deduced from the dominated convergence theorem since p ( n h,k ) is bounded in L ∞ (Ω × (0 , T )). (cid:3) Convergence towards (10) . We are now ready to pass to the limit in scheme (24) as ( h, k ) → (0 , ∞ ). Let n ∈ C ∞ c (Ω) and φ ∈ C ∞ c (0 , T ). Consider n h = Q h ( n ) in (24), multiply by φ andintegrate on (0,T) to get − (cid:90) T ( n h,k , Q h ( n )) h φ (cid:48) ( t )dt + (cid:90) T ( ∇I h ( n kh,k ) , ∇Q h ( n )) φ ( t )dt+ ν (cid:90) T ( ∇ n h,k , ∇Q h ( n )) φ ( t )dt = (cid:90) T ( G ( p ( n h,k )) n h,k , Q h ( n )) h φ ( t )dt . We briefly outline the main steps of the passage to the limit since the arguments are quite classical.We write (cid:90) T ( n h,k , Q h ( n )) h φ (cid:48) ( t )dt = (cid:90) T ( n h,k , Q h ( n )) φ (cid:48) ( t )dt + (cid:90) T (cid:2) ( n h,k , Q h ( n )) h − ( n h,k , Q h ( n )) (cid:3) φ (cid:48) ( t )dt . It is an easy matter to show, from (20) and (49), that (cid:90) T ( n h,k , Q h ( n )) φ (cid:48) ( t )dt → (cid:90) T ( n ∞ , n ) φ (cid:48) ( t )dt , and, from (18) and (19), that (cid:90) T (cid:2) ( n h,k , Q h ( n )) h − ( n h,k , Q h ( n )) (cid:3) φ (cid:48) ( t )dt → . Therefore, (cid:90) T ( n h,k , Q h ( n )) h φ (cid:48) ( t )dt → (cid:90) T ( n ∞ , n ) φ (cid:48) ( t )dt . Analogously, we obtain (cid:90) T ( G ( p ( n h,k )) n h,k , Q h ( n )) h φ ( t )dt → (cid:90) T ( G ( p ∞ ) n ∞ , n ) φ ( t )dtfrom (20), (49) and (52). The diffusion terms are treated as follows. In view of (20), (47) and(48), it is easy to check that (cid:90) T ( ∇I h ( n kh,k ) , ∇Q h ( n )) φ ( t )dt → (cid:90) T ( ∇ p ∞ , ∇ n ) φ ( t )dtand ν (cid:90) T ( ∇ n h,k , ∇Q h ( n )) φ ( t )dt → ν (cid:90) T ( ∇ n ∞ , ∇ n ) φ ( t )dt . We have thus proved that (10) holds in the distributional sense.3.2.2.
Initial condition (11) . The initial condition (11) can be recovered from (50), which gives n h,k | t =0 → n ∞ | t =0 in L q (Ω), for 1 ≤ q < ∗ , and from (8) and (29), which give n k,h → n ∞ in L p (Ω),for 1 ≤ p < ∞ as ( h, k ) → (0 , + ∞ ).3.2.3. Equivalence between (10) and (14) . In order to see the equivalence between (10) and (14)we must prove that ∇ p ∞ ≡ n ∞ ∇ p ∞ which will be obtained by proving p ∞ ∇ n ∞ ≡ x ∈ K , we decompose p ( n h,k ( x )) ∂ x i n h,k ( x ) by using theintermediate vector ξ i given in (27) into p ( n h,k ( x )) ∂ x i n h,k ( x ) = kk − n k − h,k ( ξ i ) ∂ x i n h,k ( x ) + kk − n k − h,k ( x ) − n k − h,k ( ξ i )) ∂ x i n h,k ( x )= √ kk − n k − h,k ( ξ i ) √ k n k − h,k ( ξ i ) ∂ x i n h,k ( x ) + k ( x − ξ i ) n k − h,k ( η i )( ∂ x i n h,k ( x )) , PPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 15 where we have utilized the mean value theorem in the last term for η i = α ξ i + (1 − α ) x with α ∈ (0 ,
1) and that ∂ x i n h,k ( x ) is constant on K . Thus, by virtue of (27), we find (cid:107) p ( n h,k ) ∂ x i n h,k (cid:107) L ( K ) ≤ √ kk − (cid:107) n k − h,k ( ξ i ) √ k n k − h,k ( ξ i ) ∂ x i n h,k (cid:107) L ( K ) + k h (cid:107) n k − h,k ( η i )( ∂ x i n h,k ( x )) (cid:107) L ( K ) ≤ √ kk − √ P max (cid:107)D ( n kh,k ) / ∇ n h,k (cid:107) L ( K ) + Ck h P max (cid:107)∇ n h,k (cid:107) L ( K ) , where we have used n k − h,k ( η i ) ≤ N max ( k ) k − = ( kk − P max ) k − k − → P max as k → + ∞ in the last line.Summing over K ∈ T h , noting (45) and recalling the constraint h k → p ( n h,k ) ∇ n h,k → in L ∞ (0 , T ; L (Ω))-strongly as ( h, k ) → (0 , ∞ ) . We further know, by (47) and (52), that p ( n h,k ) ∇ n h,k → p ∞ ∇ n ∞ as ( h, k ) → (0 , ∞ ) , and hence p ∞ ∇ n ∞ ≡ × (0 , T ).3.2.4. Convergence towards the complementary relation (13) . To finish the proof of Theorem 2.2,it remains to prove that (13) holds in the distributional sense. In doing so, we will start by provingthat(55) 0 ≤ (cid:90) T ( G ( p ∞ ) n ∞ , p ∞ ψ ) − ( ∇ ( p ∞ + νn ∞ ) , ∇ ( p ∞ ψ ))d s and(56) 0 ≥ (cid:90) T ( G ( p ∞ ) n ∞ , p ∞ ψ ) − ( ∇ ( p ∞ + νn ∞ ) , ∇ ( p ∞ ψ ))d s hold for all ψ ∈ C ∞ c (Ω × [0 , T ]) with ψ ≥ • To begin with, we prove that (55) is true. We use (43) to write ∂ t n h,k − (cid:101) ∆ h Σ( n h,k ) = I h ( G ( p ( n h,k )) n h,k ) . Let ρ ε = ρ ε ( t ) be a time regularizing kernel with compact support of length ε >
0. Then, extending n h,k by zero outside [0 , T ], we have(57) ∂ t n h,k ∗ ρ ε − (cid:101) ∆ h (Σ( n h,k ) ∗ ρ ε ) = I h (( G ( p ( n h,k )) n h,k ) ∗ ρ ε ) , where we have used the equalities (cid:101) ∆ h (Σ( n h,k ) ∗ ρ ε ) = (cid:101) ∆ h (Σ( n h,k )) ∗ ρ ε and I h (( G ( p ( n h,k )) n h,k ) ∗ ρ ε ) = I h ( G ( p ( n h,k )) n h,k ) ∗ ρ ε owing to the separation between spatial and temporal variables.Since ∂ t n h,k ∗ ρ ε and ( G ( p ( n h,k )) n h,k ) ∗ ρ ε are uniformly bounded in L p (Ω × (0 , T )) for 1 ≤ p ≤ ∞ with respect to ( h, k ) for each fixed ε , we also have that − (cid:101) ∆ h (Σ( n h,k ) ∗ ρ ε ) is bounded in L p (Ω × (0 , T )).as well. In virtue of Theorem 2.1 and the above bounds combined with (49) and (51), we infer thefollowing convergence, as ( h, k ) → (0 , ∞ ):(58) ∇ (Σ( n h,k ) ∗ ρ ε ) → ∇ (( p ∞ + νn ∞ ) ∗ ρ ε ) in L (Ω × (0 , T ))-strongly. On testing (57) against Q h ( I h ( n kh,k ) ψ ) with ψ ∈ C ∞ c (Ω × [0 , T ]) such that ψ ≥
0, it follows that(59) (cid:90) T ( ∂ t n h,k ∗ ρ ε , Q h ( I h ( n kh,k ) ψ )) h = (cid:90) T (( G ( p ( n h,k )) n h,k ) ∗ ρ ε , Q h ( I h ( n kh,k ) ψ )) h − (cid:90) T ( ∇ (Σ( n h,k ) ∗ ρ ε ) , ∇Q h ( I h ( n kh,k ) ψ )) . Since ( ∂ t n h,k ∗ ρ ε , Q h ( I h ( n kh,k ) ψ )) h ≥
0, we obtain(60) 0 ≤ (cid:90) T (( G ( p ( n h,k )) n h,k ) ∗ ρ ε , Q h ( I h ( n kh,k ) ψ )) h − (cid:90) T ( ∇ (Σ( n h,k ) ∗ ρ ε ) , ∇Q h ( I h ( n kh,k ) ψ )) . Taking the limit as ( h, k ) → (0 , ∞ ) yields(61) (cid:90) T (( G ( p ( n h,k )) n h,k ) ∗ ρ ε , Q h ( I h ( n kh,k ) ψ )) h d t → (cid:90) T (( G ( p ∞ ) n ∞ ) ∗ ρ ε , p ∞ ψ )d t and(62) (cid:90) T ( ∇ (Σ( n h,k ) ∗ ρ ε ) , ∇Q h ( I h ( n kh,k ) ψ ))d t → (cid:90) T ( ∇ (( p ∞ + νn ∞ ) ∗ ρ ε ) , ∇ ( p ∞ ψ ))d t. In order to prove (61), we use the decomposition ( u h , v h ) h = ( u h , v h ) + ( I h ( u h v h ) − u h v h ,
1) for u h = ( G ( p ( n h,k )) n h,k ) ∗ ρ ε and v h = Q h ( I h ( n kh,k ) ψ ) to write (cid:90) T (( G ( p ( n h,k )) n h,k ) ∗ ρ ε , Q h ( I h ( n kh,k ) ψ )) h = (cid:90) T (( G ( p ( n h,k )) n h,k ) ∗ ρ ε , Q h ( I h ( n kh,k ) ψ ))+ (cid:90) T ( I h (( G ( p ( n h,k ))) n h,k ) ∗ ρ ε Q h ( I h ( n kh,k ) ψ ) − ( G ( p ( n h,k )) n h,k ) ∗ ρ ε Q h ( I h ( n kh,k ) ψ ) , . Then, it follows from (21), (49) and (52) that the first term converges to (cid:82) T (( G ( p ∞ ) n ∞ ) ∗ ρ ε , p ∞ ψ )d t, and, on noting that (cid:107)∇Q h ( I h ( n kh,k ) ψ ) (cid:107) ≤ C (cid:107)∇I h ( n kh,k ) (cid:107)(cid:107) ψ (cid:107) L ∞ + C (cid:107)I h ( n kh,k ) (cid:107) L ∞ (cid:107)∇ ψ (cid:107) from (19), and on recalling (18) and (46), the second term converges to zero; thereby (61) holds.In order to prove (3.2.4), we write (cid:90) T ( ∇ (Σ( n h,k ) ∗ ρ ε ) , ∇Q h ( I h ( n kh,k ) ψ )) = (cid:90) T ( ∇ (Σ( n h,k ) ∗ ρ ε ) , ∇ ( I h ( n kh,k ) ψ )) − (cid:90) T ( ∇ (Σ( n h,k ) ∗ ρ ε ) , ∇ ( I h ( n kh,k ) ψ − Q h ( I h ( n kh,k ) ψ ))) . Then, it follows from (58), (48) and (51) that the first term converges to (cid:82) T ( ∇ (( p ∞ + νn ∞ ) ∗ ρ ε ) , ∇ ( p ∞ ψ ))d t , and on noting that (cid:107)∇ ( I h ( n kh,k ) ψ − Q h ( I h ( n kh,k ) ψ )) (cid:107) ≤ h (cid:107)∇I h ( n kh,k ) (cid:107)(cid:107) ψ (cid:107) W , ∞ (Ω) from (21) and on recalling (46), the second term converges to zero; thereby (3.2.4) holds.Thus, by applying the previous convergences (61) and to (60), we arrive at0 ≤ (cid:90) T ( G ( p ∞ ) n ∞ ∗ ρ ε , p ∞ ψ ) − ( ∇ (( p ∞ + νn ∞ ) ∗ ρ ε ) , ∇ ( p ∞ ψ ))d t, and finally (55) holds by taking the limit as ε → PPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 17 • We proceed to prove (56). Write the first term on the right-hand side of (59) as(63) (cid:90) T ( ∂ t n h,k ∗ ρ ε , Q h ( I h ( n kh,k ) ψ )) h = (cid:90) T ( ∂ t n h,k ∗ ρ ε , I h ( n kh,k ) ψ ) h + (cid:90) T ( ∂ t n h,k ∗ ρ ε , Q h ( I h ( n kh,k ) ψ ) − I h ( n kh,k ) ψ ) h . These two terms are handled as follows. For the second term of (63), we have, by (17), (21) and(41), that (cid:90) T ( ∂ t n h,k ∗ ρ ε , Q h ( I h ( n kh,k ) ψ ) − I h ( n kh,k ) ψ ) h d s → h, k ) → (0 , ∞ ) . For the first term of (63), we have that, for each a ∈ N h ,( ∂ t n h,k ( a , t ) ∗ ρ ε ) n kh,k ( a , t ) = n kh,k ( a , t ) (cid:90) R ∂ t n h,k ( a , s ) ρ ε ( t − s )d s = (cid:90) R n kh,k ( a , s ) ∂ t n h,k ( a , s ) ρ ε ( t − s )d s + (cid:90) R ( n kh,k ( a , t ) − n kh,k ( a , s )) ∂ t n h,k ( a , s ) ρ ε ( t − s )d s. On integrating by parts in time and using (31) and (39), we obtain (cid:90) R n kh,k ( a , s ) ∂ t n h,k ( a , s ) ρ ε ( t − s )d s = 1 k + 1 (cid:90) R n k +1 h,k ( a , s ) ∂ t ρ ε ( t − s )d s → h, k ) → (0 , ∞ ). Furthermore, for s > t , we have that n kh,k ( a , t ) − n kh,k ( a , s ) ≤ ρ ε ) ⊂ ( − ε, (cid:90) R ( n kh,k ( a , t ) − n kh,k ( a , s )) ∂ t n h,k ( a , s ) ρ ε ( t − s )d s ≤ . Letting first ( h, k ) → (0 , ∞ ) in (63) and then ε →
0, we obtain (56) by repeating the argumentsthat led to (55).As a result of (55) and (56), we note that(64) (cid:90) T ( G ( p ∞ ) n ∞ , p ∞ ψ ) − ( ∇ ( p ∞ + νn ∞ ) , ∇ ( p ∞ ψ ))d s = 0is satisfied for all ψ ∈ C ∞ c (Ω × [0 , T ]) with ψ ≥
0, and therefore it also holds for all ψ ∈ C ∞ c (Ω × [0 , T ]).From the fact that p ∞ ∇ n ∞ = 0 and p ∞ ≥ × (0 , T ), we also deduce that ∇ p ∞ ·∇ n ∞ = 0a.e. in Ω × (0 , T ). As a consequence, the above variational equation (64) is equivalent to (cid:90) T ( G ( p ∞ ) n ∞ , p ∞ ψ ) − ( ∇ p ∞ , ∇ ( p ∞ ψ ))d s = 0which, taking into account (53), implies (13) in the distributional sense.4. An algorithm on unstructured meshes
In order to avoid using structured meshes, we propose the following scheme. Find n h,k ∈ C ([0 , T ]; N h ) such that(65) (cid:26) ( ∂ t n h,k , n h ) h + k (( n h,k ) k − ∇ n h,k , ∇ n h ) + ν ( ∇ n h,k , ∇ n h ) = ( G ( p ( n h,k )) n h,k , n h ) h ∀ n h ∈ N h ,n h,k (0) = n h,k . Equivalently, we may write (65) as( ∂ t n h,k , n h ) h + ( n h,k ∇ p ( n h,k ) , ∇ n h ) + ν ( ∇ n h,k , ∇ n h ) = ( G ( p ( n h,k )) n h,k , n h ) h . Here the finite-element space N h is constructed over a family of triangulations {T h } h> of Ωbeing shape-regular, quasi-uniform and with acute angles. This acuteness property implies (22)for the particular case where D is the d × d identity matrix [7]. We summarize the properties ofscheme (65) in the following theorem. Theorem 4.1.
Suppose that (H1) - (H4) are satisfied. Then scheme (65) satisfies the followingproperties. For all a ∈ N h and t ≥ , we have: ≤ n h,k ( a , t ) ≤ N max ( k )0 ≤ n kh,k ( a , t ) ≤ P max N max ( k ) ,∂ t n h,k ( a , t ) ≥ , ∂ t n kh,k ( a , t ) ≥ , and the a priori estimates: (cid:107) n h,k (cid:107) L ∞ (0 ,T ; L (Ω)) ∩ L (0 ,T ; H (Ω)) ≤ C, (cid:107) ∂ t n h,k (cid:107) L ∞ (0 ,T ; L (Ω)) + (cid:107) ∂ t n kh,k (cid:107) L (0 ,T ; L (Ω)) ≤ C, with C > being a constant independent of ( h, k ) .Proof. Full details of the proof are left to the interested reader since it follows mutatis mutandis the same arguments as for scheme (24). (cid:3)
Corollary 4.1.
Under hypotheses (H1) - (H4) , it follows that (66) (cid:88) K ∈T h (cid:18)(cid:90) K > | ∂ x i n kh,k ( x ) | + (cid:90) K < | ∂ x i I h n kh,k ( x ) | (cid:19) d x ≤ C, where K > = (cid:40) x ∈ K : n k − k,h ( ξ i ) n k − h,k ( x ) > (cid:41) and K < = (cid:40) x ∈ K : n k − k,h ( ξ i ) n k − h,k , ( x ) < (cid:41) , with C > being a constant independent of ( h, k ) .Proof. Choose ¯ n h = I h ( n kh,k ) to get(67)( ∂ t n h,k , I h ( n kh,k )) h +(( n h,k ) k − ∇ n h,k , ∇I h ( n kh,k ))+ ν ( ∇ n h,k , ∇I h ( n kh,k )) = ( G ( p ( n h,k )) n h,k , I h ( n kh,k )) h . It follows immediately from (39) and (41) that(68) ( ∂ t n h,k , I h ( n kh,k )) h ≥ , and from (26) that(69) ν ( ∇ n h,k , ∇I h ( n kh,k )) = ν ( D ( n h,k ) ∇ n h,k , ∇ n h,k ) ≥ . Combining (67)-(69) yields on noting (31) and (39) that(( n h,k ) k − ∇ n h,k , ∇I h ( n kh,k )) ≤ G (0) | Ω | N max ( k ) P max . PPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 19
Finally, we invoke again (26) and recall (27) to set k (( n h,k ) k − ∇ n h,k , ∇I h ( n kh,k )) = ( ∇ n kh,k , ∇I h ( n kh,k ))= (cid:88) K ∈T h (cid:90) K (cid:32) n k − k,h ( ξ i ) n k − h,k ( x ) | ∂ x i n kh,k ( x ) | + n k − h,k ( x ) n k − k,h ( ξ i ) | ∂ x i I h n kh,k ( x ) | (cid:33) d x . This completes the proof via the definitions of K < and K > . (cid:3) Remark 4.1.
Unfortunately, convergence for scheme (65) is not clear because estimate (66) doesnot provide enough control over the gradient of { n kh,k } h,k or {I h n kh,k } h,k in order to obtain compact-ness and therefore to pass to the limit as ( k, h ) → (0 , + ∞ ) . Numerical simulation
Temporal integration.
It is assumed here for simplicity that we have a uniform partition of[0 , T ] into M pieces, with time step size τ = T /M and the time values ( t m = mτ ) Mm =0 . To simplifythe notation let us denote δ t n m +1 = n m +1 − n m τ .First we present a first-order time integration for scheme (65). Algorithm 1 : Linear semi-implicit time-stepping scheme
Step ( m + 1): Given n mh,k ∈ N h , find n m +1 h,k ∈ N h solving the algebraic linearsystem(70) (cid:26) ( δ t n m +1 h,k , n h ) h + k (( n mh,k ) k − ∇ n m +1 h,k , ∇ n h )+ ν ( ∇ n m +1 h,k , ∇ n h ) = ( G ( p ( n mh,k )) n mh,k , n ) h , for all n h ∈ N h .5.2. Computational experiments.
In this section, we present several numerical experiments totest the algorithm presented herein. To do this, we consider the evolution of problem (1)-(3) with n ( x, y ) = α e − ( x + y ) on the computational domain Ω = ( − , × ( − ,
10) with α > h = 0 . τ = 10 − . The choice of the time step τ is such that it helps to mitigate the possibly numericaldeviation of the d -simplexes K ∈ T h from the right-angled structure. The resulting matrix isstrictly diagonally dominant.Our intention is to illustrate the behavior exhibited by the solution to problem (1)-(3) when thediffusion coefficient ν , the parameter k and the homeostatic pressure P max vary.We will set α and P max to be 1 and k to be 100 if not stated otherwise. Moreover, we consider G ( p ) = 200 π arctan(4( P max − p ) + ) , Analysis of the effect of α (contraction/dilation coefficient of the initial datum). In this testwe choose α = 0 . n and thepressure p ( n ) when the maximum of the initial density takes different values. In particular, we havefor α = 0 . n occurs only at the point (0 ,
0) and is 0 .
5, hence N max ( k )remains below of 1. We thus observe that the maximum increases without modifying essentiallythe exponential shape of the initial datum n until reaches N max ( k ) = 1. Once the density takesthe value 1 at t = 0 . ,
0) due to the fact that the pressure starts increasing and pushes forward the tumor cells. Then the exponential structure of the initial datum n becomes a traveling wave shapewhich moves outwards as t increases. This behavior causes that the evolution of the interface isdelayed concerning the case α = 1 as shown in Figures 1 and 2 since the maximum value 1 isreached from the beginning.Figure 3 represents the difference between the density and the pressure at times t = 0 .
1, 0 .
2, 0 . .
4, and indicates that the pressure is responsible for the advance of the tumor cells which isdeduced from the annulus shape of the difference.
Figure 1.
Evolution of the density at times t = 0 . , . , . , . α = 0 . Figure 2.
Evolution of the pressure at times t = 0 . , . , . , . α = 0 . α = 1 (bottom).5.2.2. Analysis of the effect of ν (active motion coefficient). Now we set P max = 1 and take differentvalues of ν = 0, 0 . n h,k is shown in Figure 4 where we seethat the velocity of propagation of the tumor cells increases with respect to ν as noted for times t = 0 .
1, 0 .
2, 0 . .
4. Moreover, no particular differences have been observed in the width ofthe interface between the tumor and pre-tumor cells for the different values of ν .5.2.3. Analysis of the effect of k . In this simulation we select k = 10 and 1000. The first thing wehave noted is that there is a dependence between k and τ which has been taken 0 . · − . As canbe seen in Figure 5, there are no particular differences for k = 10 and 1000 at times t = 0 .
1, 0 . . . PPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 21
Figure 3.
Evolution of the difference between the density and pressure at times t = 0 . , . , . , . α = 0 . α = 1 (bottom). Figure 4.
Comparison of the density at times t = 0 . , . , . , . ν = 0(top), 0 . Figure 5.
Comparison of the density at times t = 0 . , . , . , . k = 10 (top) and 1000 (bottom). Analysis of the effect of P max . Let us take P max = 10 and 30. Figure 6 shows that thedynamics is sensitive to the different values for the homeostatic pressure. We highlight that, for P max = 30, the evolution of the interphase is faster than the one for P max = 10. Moreover, theshape of the interphase seems different as depicted in Figure 6 for times t = 0 .
1, 0 .
2, 0 . . Figure 6.
Comparison of the density at times t = 0 . , . , . , . P max = 10 (top) and 30 (bottom). References [1]
Becker, R.; Feng, X.; Prohl, A.
Finite element approximations of the Ericksen-Leslie model for nematicliquid crystal flow.
SIAM J. Numer. Anal. 46 (2008), no. 4, 1704–1731.[2]
Bertoluzza, S
The discrete commutator property of approximation spaces.
C. R. Acad. Sci. Paris S´er. I Math.329 (1999), no. 12, 1097–1102.[3]
Betteridge, R.; Owen, M. R.; Byrne, H. M.; Alarc´on, T.; Maini, P. K.
The impact of cell crowdingand active cell movement on vascular tumour growth.
Netw. Heterog. Media 1 (2006), no. 4, 515?535.[4]
Brenner, S. C.; Scott, L. R. , The mathematical theory of finite element methods , Third edition. Texts inApplied Mathematics, 15. Springer, New York, 2008.[5]
Br´u, A.; Albertos, S.; Subiza, J. L.; Asenjo, J. A.; Broe, I.
The universal dynamics of tumor growth.
Biophys. J. 85 (2003), no. 5, 2948–2961.[6]
Byrne, H. M.; Drasdo, D. , Individual-based and continuum models of growing cell populations: a compari-son. , J. Math. Biol. 58 (2009), no. 4-5, 657–687.[7]
Ciarlet, P.G.; Raviart, P.-A.
Maximum principle and uniform convergence for the finite element method ,Comput. Methods Appl. Mech. Engrg. 2 (1973) 17–31.[8]
Drasdo, D.; Hoehme, S.
Modeling the impact of granular embedding media, and pulling versus pushing cellson growing cell clones.
New J. Phys. 14 (2012) 055025.[9]
Ern, A; Guermond, J.-L. , Theory and practice of finite elements , Applied Mathematical Sciences, 159.Springer-Verlag, New York, 2004.[10]
Girault, V.; Lions, J.-L. , Two-grid finite-element schemes for the transient Navier-Stokes problem . M2ANMath. Model. Numer. Anal. 35 (2001), no. 5, 945–980.[11]
Perthame, B.; Quir´os, F.; Tang, M.; Vauchelet, N. , Derivation of a Hele-Shaw type system from a cellmodel with active motion , Interfaces Free Bound. 16 (2014), no. 4, 489–508.[12]
Perthame, B.; Quir´os, F.; V´azquez, J. L. , The Hele-Shaw asymptotics for mechanical models of tumorgrowth , Arch. Ration. Mech. Anal. 212 (2014), no. 1, 93–127.[13]
J. Ranft, M. Basan, J. Elgeti, J.-F. Joanny, J. Prost and F. J¨ulicher , Fluidization of tissues by celldivision and apoptosis , Proceedings of the National Academy of Sciences, 107 (2010), no. 49, 20863–20868.[14]
Saut, O.; Lagaert, J.-B.; Colin, T.; Fathallah-Shaykh, H. M.
A multilayer grow-or-go model forGBM: effects of invasive cells and anti-angiogenesis on growth.
Bull. Math. Biol. 76 (2014), no. 9, 2306–2333.[15]
Scott, L.R.; Zhang, S.
Finite element interpolation of non-smooth functions satisfying boundary conditions.
Math. Comp. 54 (1990) 483–493.[16]
Simon, J.