The Stefan problem with moving boundary

A mathematical model of the linear thermodynamic equations with moving ends, based on the Stefan Problem is considered. In this work, we are interested in obtaining existence, uniqueness and regularity using the Faedo-Galerkin method. For numerical solutions, we shall employ the finite element method together with the Crank-Nicolson method. A numerical experiment is presented to show the moving boundary for the problem.


Introduction
One wide variety of problems with moving ends has been studied by many years.The thermal equations with moving ends give origin to several physical problems, such as, Stefan Problem, who can be interpreted as flame propagation problem, study of atom movements, study of combustion theory.Flame propagation problems have been established in [2] that studied a mathematical formulation suitable in the analyze for the theory of diffusion flames.Studies the mathematical model of moving ends arising in combustion theory was developed by Schmidt-Lainé [1] and [3].The discharge or absorption of thermal energy during the change in atom arrangement have been studied in [8] and [7].We can also see the problem as a model for certain processes involving chemical reactions.The objective of this paper is to obtain existence and uniqueness of solutions and too we apply the finite element method with a finite difference method in time to obtain an approximated numerical solution.Some numerical experiments are presented to show the effect of moving boundary in the solutions of the linear thermodynamic, whose model is based on the Stefan Problem, (see [11]) given by, u(x, 0) = ϕ(x), u(0, t) = ψ(t), u(s(t), t) = ρ(t), ∀t > 0 where u(x, t) is the temperature and s(t) is the solution of equation: −λ(t)s ′ (t) + u x (s(t), t) = h(t) s(0) = s 0 , must be determined and ψ, ρ, h, λ, ϕ, s 0 are initial date of the problem defined for t ≥ 0 and ϕ(x) for 0 ≤ x ≤ s 0 .

Formulation Problem
Our objective is study the approximated solutions of the following thermodynamics system, (I) where k > 0 is the thermal conductivity.We consider the domain Q t ⊂ IR 2 defined by and the horizontal length of the interval is defined by Let us study the existence, uniqueness, and the approximated solution of the Problem (4).For this we consider the following hypothesis: where γ 0 and γ 1 are positive constants.Consider now the change of variable to transform the domain Q t in a cylindric domain Q, given by Note that τ ∈ C 2 (Q t ) and by inverse function Theorem, the application τ −1 also is C 2 (Q) (see for instance [5,9]).The change of variable v(y, t) = u(α(t) + γ(t)y, t), implies that , t .Using the application τ , we obtain the following cylindrical problem: (II) ( where For convenience, we have also used the prime ( ′ ) to denote the derivative with respect to time t and (( .)), ( .), .and | .|, the scalar product and the norms in H 1 0 (0, 1) and L 2 (0, 1), respectively.As our primary objective for this work is to present numerical results of the problem.We shall show the existence and uniqueness from the solution of the cylindrical Problem (3) who implies in the existence and uniqueness of the solution of the noncylindrical Problem (4), since they are equivalent.

Existence and Uniqueness
We will investigate global existence and uniqueness solution for Problem(3) and consequently for Problem (4).For this, consider the notation Ω = (0, 1), Ω t = (α(t), β(t)) and Ω 0 = (α(0), β(0)).Theorem 3.1 Under the hypotheses (H1) and (H2) and the initial data v 0 ∈ H 1 0 (Ω) and g ∈ L 1 (0, ∞; L 2 (Ω)) ∩ L 2 (0, T ; L 2 (Ω)), there exists only one strong solution v : Q −→ IR of the problem (II), that is, satisfying the following conditions: And consequently we obtain the following result, Theorem 3.2 Under the hypotheses (H1) and (H2) and the initial data , there exists only one strong solution u : Q t −→ IR of the problem (I), that is, satisfying the following conditions: The details of the demonstration can be found in [10].In this work we are omitting the calculations and will find just the idea of the demonstration.
Proof: We introduce the approximate solutions.Let T > 0 and denote by V m the subspace spanned by {w 1 , w 2 , ..., w m }, where {w ν ; ν = 1, • • • m} are the first m eigenvectors of the space H 1 0 (Ω), solution of the spectral problem - where v m is the solution of the system of ordinary differential equations The system (4) has local solution in the interval (0, T m ), by Carathéodory Theorem.To extend the local solution to the interval (0, T ) independent of m the following a priori estimate is needed.
Integrating the third term of (5), by parts and using the boundary conditions, we have, Using the norma equivalence of norms in H 1 0 (Ω), Hölder inequality ( 6) and substituting into (5), we have Integrating from [0, t), with t ∈ [0, T m ), we get From the hypotheses (H1) and (H2), where c is a positive constant.Let By the use of Gronwall's Lemma in the last inequality, we obtain the following Estimate II: The second term on left hand side in (8) implies that by equivalence of norms in H 1 0 (Ω).For the third term in (8) can be estimated as, Substituting ( 9), (10) in the equation ( 8), we get From elementary inequality Therefore, We observe that, from the hypothesis (H1), γ(t) < γ 1 , ∀t ∈ [0, T ].Integrating in 0 < t < T m , using the hypothesis (H1) and the definition a(t), we obtain Note that, from hypothesis (H2), where c is a positive constant.Then, there is a positive constant c, such that By Gronwall inequality we have Estimate III: 4) and after to integrate the second term by parts we have The first term on left hand side in (13) implies that The second term on left hand side in (13) give us As 0 < y < 1 and from the hypothesis (H1), we have Then Similarly we have Substituting ( 14), ( 15), ( 16) and (17) in the equation (13), we obtain Denoting C 1 by We have then where Integrating from 0 to t < T m and as g ∈ L 2 (0, T ; By Gronwall inequality we conclude that The estimates obtained in ( 7), ( 12) and (18), permit us to pass the limits in the approximate system (4) in the Faedo-Galerkin method, see for instance ( [5,9]) and hence, we have proved the existence of solutions v(y, t) in the sense defined in Theorem 3.1.

⊔ ⊓
Uniqueness: The uniqueness of solution is made by contradiction and does not prove that work.⊔⊓ Hence, now we can to prove the Theorem 3.2: Proof: Let v the solution of the Problem (3), with the following initial data Consider the function u(x, t) = v(y, t), where x = α(t) + γ(t)y.To verify that u(x, t), under the hypotheses of Theorem 3.2, is the solution of Problem 1, it is sufficient to observe that the mapping τ : Since the Problems 1 and 3 are equivalently, then u satisfy the Problem 1.⊔ ⊓ In the next theorem, we shall also prove the following regularity result, under further requirements of initial conditions and the regularity for α(t) and β(t), given by So we establish the following regularity result, whose proof can be found in [10]: Theorem 3.3 Under the hypotheses (H2) and (H3) and given the initial data )) and g ∈ H k (0, T ; H k (0, 1)), then there exists a unique solution v : Q → IR, of the problem (II), satisfying the following conditions The regularity of v(y, t) given by Theorem (3.3), implies the regularity of u(x, t) solution of the Problem (4) and the uniqueness of solutions of the Problem ( 4) is a direct consequence of uniqueness of the Problem (3).⊔ ⊓

Approximate Solution
Our goal in this section is the numerical implementation of approximate solutions.To obtain the numerical approximate solutions we will use both finite element method and finite difference method.Moreover, some numerical experiments will be presented to analyze the effect of the moving boundary in the thermodynamics system (4).For convenience, our numerical analysis using finite element method approximation will be based on the equivalent problem (3) in the rectangular domain, instead of the problem (4), for which the domain depends on time.The variational formulation of the problem (3), is given by 4.1.Faedo-Galerkin Method and Approximation.Let V m the subspace spanned by {ϕ 1 , ϕ 2 , . . ., ϕ m+1 }, where Restricting the equation ( 19) in the subspace V m and taking v h in (20) we find that We define where A, B(t) and D are m × m square matrix and G is a vector known m × 1. Substituting the matrices in ( 22), we obtain the following ordinary differential system, where c = [c 1 , c 2 , ..., c m ] t is a m × 1 c is the vector to be determined in each discrete time.

Finite Difference Method.
In order to solve the system (23) in each discrete time, we will apply the numerical method due to Crank-Nicolson (see, for instance, Hugles [6] ).For T > 0, N a given positive integer, we define the time step ∆t and let c n = c(t n ) be the approximate solution of the exact solution c(t) of ( 23), where we denote the discrete time in the interval [0, T ] by t n = n∆t, n = 0, 1 • • • N , and the values of W at the discrete time t n by W n .
Consider the following approximation for the function c(t), using difference finite method in the discrete time n = 0, 1, Setting t = t n in (23), using (24) and after multiplying the equation by ∆t, we obtain the iterative method (25) Suppose that the matrices, D = D ij , A = A ij and B = B ij are known, then the iterative method (25) can be easily implemented.From initial conditions, c 0 is known.Then taking t = 0 into (25) yields, Solving the linear system then get the vector , are determined with an iterative process in each instant of time t n = n∆t.
Note that the matrix D and A are positive definite and symmetric, and the B matrix is anti-symmetric.In addition, the matrix of coefficients (D +(∆t/2)(a n A− B n )) of the system (25) is not singular for each instant of time.

Finite Element Approximation.
To calculate the matrices of linear system (25), we need to introduce the basis function ϕ i ∈ V m .In finite element method, the bases functions are piecewise polynomials of some degree in Ω and vanish on ∂Ω (see, for instance, [4]).More specifically, in this work, we have used the hat function, i.e, piecewise linear functions, as the basis function of V m subspace, defined in the following way where we are considering the uniform mesh, h = h i = y i+1 − y i , i = 1, 2, . . ., m in the discretization in m-parts, with 0 = y 1 < y 2 < • • • < y m+1 = 1.Note that, if |i − j| > 2, then (ϕ i , ϕ j ) = 0, and (∂ϕ i /∂y, ∂ϕ j /∂y) = 0. Hence all the matrices of system (25) are tridiagonal.Note that we can interpolate the function g(y, t) in Ω, using the base function ϕ i (y), i.e, we take Thus, the function G(t) can be rewritten as follows: and for discrete time, we have G(t n ) = G n = Dg n , where D is the matrix defined before.

Numerical Simulation
In this section, attention is turned to implementation of approximate solutions, whose programs were developed in the thesis [12].A numerical experiments will be given to analyze the effect of the moving ends of Ω t .The goal is to obtain exact solutions to compare with the approximate numerical one, which are obtained by the iterative method (25).For this, we introduce an appropriate external force f (x, t) sufficiently regular on the right-hand side of the original Problem (1).Therefore, the exact solution is known for suitable choice of f (x, t), which allow us to get the numerical simulation for Problem (1) or Problem (3).The error estimates, in the semi-discrete or fully discrete, can be found [10].This will be the procedure for analysis of the convergence.Notice that, for problem (1) the equivalent external force acting in the equation (3) 1 would be a function g(y, t) defined from f via the diffeomorphism T , defined in (2).This way, we will work with the system (3).
Example -.Let v(y, t) = 1 π 2 sin(πy)cos(πt) be a exact solution of the o Problem (3) for suitable choice of g(y, t).Then, the initial and boundary conditions read, respectively Numerical Error -.In the norm L ∞ (0, T ; L 2 (Ω t )) the numerical error between the approximate u h (x, t) and exact u(x, t) solution is calculated by : The table that follow, is shown the errors obtained between exact and approximate solutions to the problem (1) when the step time is fixed ∆t = 0.01 and for several sizes of mesh h = 0.1; 0.05; 0.02; 0.01; 0.001 that represent respectively m = 10; 20; 50; 100; 1000 nodes of the interval [α(t), β(t)].For this case were made 100 iterations in the iterative method (25).Note that each error depends algebraically on h.That is, the error E L ∞ (0,T ;L 2 (Ωt)) decrease as the mesh size h decrease, as expected, albeit slowly.
The results in the table of error, encourages us that the numerical method used is good, thus, employing the same procedure to determine the numerical solution of the problem with zero external force, g(y, t) = 0 in the problem (3), i.e, f (x, t) = 0 the original problem (1).Moreover, when f (x, t) = 0 we have the exponential decay of energy, which is dependent on the thermal conductivity k as shown in the following theorem: Asymptotic behavior of the energy -.The goal in this section is to establish a rate decay for the energy of system (1).Therefore, the asymptotic behavior, as t → ∞, of the natural energy will be obtained inside of the time dependent domain Q t .Thus, we can state Theorem 5.1 Assuming the hypotheses of Theorem 3.2, then the energy E(t) associated with the global weak solution of system (1) satisfies where γ 1 is a positive real constant, given in (H1).
In order to prove the Theorem 5.1 one needs to establish the Poincarè inequality in Ω t .Thus one has.
Proof: In fact, from the fundamental theorem of calculus one has From this and Cauchy-Schwartz inequality one has Integrating in Ω t one gets Integrating from 0 to t we conclude the demonstration of Theorem 5.1 ⊔ ⊓ In the following, we show graphically the asymptotic decay of energy for different k, representing the materials, that are inserted in the table: We can see that the rate of decay of energy depends on the value of k, as expected.
The Fig. 2, depicts the solution when f ≡ 0 in the Problema 1 and taking k = 1, 71 (silver).That is, the effect of the moving boundary in the thermodynamics system (1).One can see the decline of the solution when the time t is increasing.

Figure 1 -
Figure 1 -Asymptotic decay of energy

Figure 2 -
Figure 2 -Approximated solution u(x, t) for the Problem with moving boundary.