Numerical method, existence and uniqueness for thermoelasticity system with moving boundary

Numerical method, existence and uniqueness for thermoelasticity system with moving boundary

M.A. Rincon^I; B.S. Santos^I; J. Límaco^II

^IInstituto de Matemática, Universidade Federal do Rio de Janeiro. E-mails: [e-mail]/ [e-mail]

^IIInstituto de Matemática, Universidade Federal Fluminense Rio de Janeiro, Brazil. E-mail: [e-mail]

ABSTRACT

In this work, we are interested in obtaining existence, uniqueness of the solution and an approximate numerical solution for the model of linear thermoelasticity with moving boundary. We apply finite element method with finite difference for evolution in time to obtain an approximate numerical solution. Some numerical experiments were presented to show the moving boundary's effects for problems in linear thermoelasticity.

Mathematical subject classification: 35A05, 35A40, 65M60, 65M06.

Key words: thermoelasticity system; moving boundary; finite element method; finite difference method.

1 Introduction

Let Q_t = {(x,t) Î ²; a(t) < x < b(t), 0 < t < T} be the non-cylindrical domain with boundary

and consider the following problem:

Existence and uniqueness of linear and nonlinear elasticity in a bounded or an unbounded cylindrical domain, has been studied by several authors, among them, [4] and [5].

In this work, we will investigate existence, uniqueness and approximate solution of the problem (I). We will also show the influence of moving boundary employing numerical examples. For this we consider the following hypotheses:

H1: a, b Î C²([0, T);),

H2: $k₁Î , such that,

H3:k > 0, and h1.h2 > 0.

We will now consider a change of variables to transform the domain Q_t into a cylindrical domain Q. Observe that, when (x,t) varies in Q_t the point (y,t) of ², with y = (x-a(t))/g(t) varies in the cylinder Q = (0,1)×(0,T). Thus, we define the application

The application belongs to C² and its inverse ^-1 is also C². The transformation of a moving boundary domain to a domain with fixed boundary has been employed elsewhere (see [2, 10, 11]).

Doing the change of variable v(y,t) = u(a(t) + g(t)y,t) and f(y,t) = q(a(t) + g(t) y,t) and applying to the problema (I), we obtain the following equivalent problem defined in a fixed cylindrical domain:

where

Let (( , )), || · || and ( , ), | · |, be respectively the scalar product and the norms in (0,1) and L²(0,1). We denote by a₁(t,v,w) and b₁(t,v,w) the bilinear forms, continuous, symmetric and coercive, defined in (0,1) by

2 Existence and uniqueness

We shall first establish the existence and uniqueness of problem (II) in Theorem 2 as auxiliary and then prove the following Theorem 1 of the original problem (I).

Theorem 1. Under the hypotheses (H1), (H2) and (H3) and given the initial data

there exist functions {u; q}: Q_t ® , solution of Problem (I) in Q_t, satisfying the following conditions:

Theorem 2. Under the hypotheses (H1), (H2) and (H3) and given the initial data

there exists functions {v; f}: Q ® , solution of Problem (II) in Q, satisfying the following conditions:

Proof of Theorem 2. To prove the theorem, we introduce the approximate solutions. Let T > 0 and denote by V_m the subspace spanned by {w₁,w₂,...,w_m}, where {w_n, l_n; n = 1,...m} are solutions of the spectral problem ((w_i,v)) = m(w_i,v), "v Î (0,1). If {v_m; f_m} Î V_m then it can be represented by

Let us consider {v_m; f_m} solutions of the system of ordinary differential equations,

where w Î V_m. The system (III) has local solution in the interval (0,T_m). To extend the local solution to the interval (0,T) independent of m, the following estimates are necessary:

A priori estimate

Taking w = and w = f_m in the equation (III)₁ and (III)₂, respectively, we get

Note that, we have the following relations:

Multiplying (4) by (h₂/h₁), adding it to (5) and using (6) we have

Knowing that a₁(t,v,w) and b₁(t,v,w) are coercive forms, by integrating (7) and applying the Gronwall's inequality, we get

Second estimate

Taking the derivative with respect to t, of approximate system (III)_1,2, and also w = , w = , respectively, we obtain

and

We also have the following relations:

Multiplying (9) by (h₁/h₂), adding it to (10) and using (11), we obtain

From (III)_1,5, |(0)|² and |(0)|² are bounded. Hence, by integrating (12) with respect t and applying the Gronwall's inequality, we get

Third estimate

Taking w = - ¶²v_m/¶y² and w = - ¶²f_m/¶y², in the approximate system (III)_1,2, we have

and

Note that, we have the following equalities:

From (14), (15) and (16) and since a₁(t,v,w) and b₁(t,v,w) are coercive forms, we obtain

The estimates obtained in (8), (13), (17) and (18), permit us to pass the limits in the approximate system (III)_1,2 in the Galerkin method and hence, we have proved the existence of solutions {v,f} in the sense defined in Theorem 2.

Uniqueness of solution

Let {,} and {,} be two solutions of Problem (II). Then v = - and f = - are also solutions of Problem (II), with null initial conditions. Then, multiplying the equation (II)_1,2, respectively by (h₁/h₂)v and f, we obtain

From Gronwall Lemma, we have |v'|² + ||v||² +|f|² = 0 and therefore, we conclude that v = f = 0 for all 0 < t < T. This completes the proof of Theorem 2.

The original problem (I)

Now let us restate the previous results for the original problem (I) in order to prove Theorem 1.

Proof of Theorem 1. Let {v,f} be a solution of Problem (II), with initial data given by

Consider the functions u(x,t) = v(y,t) and q(x,t) = f(y,t), where x = a(t) + g(t)y. To verify that u(x,t) and q(x,t), under the hypotheses of Theorem 1, are a solution of problem (I), it is sufficient to observe that the mapping: (x,t) ® ((x-a)/g , t) of the domain Q_t into Q = (0,1)×(0,T) is of class C². Since that

and from problem (II) we also have that {u,q} satisfies the problem (I).

The regularity of {v(y,t),f(y,t)} given by Theorem 2 implies that {u(x,t), q(x,t)} is a solution of problem (I) and the uniqueness of the solution of problem (I) is a direct consequence of the uniqueness of problem (II).

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 thermoelasticity system.

For convenience, our numerical analysis using finite element method approximation will be based on the equivalent problem (II) in the rectangular domain, instead of the problem (I), for which the domain depends on time. We will consider, in numerical simulations, the case in which the following change in the boundary functions, a(t) = -K(t) and b(t) = K(t), is assumed.

Note that, now we have

being the non-cylindrical domain with boundary

and consequently we have the fixed cylindrical domain Q = (-1,1)×(0, T). In this way we obtain the following relation between the functions:

3.1 Variational form of the problem

Let us consider the following variational form, given by (III)_1,2,

and

where now, using (20), the functions b_i and a_i are given by

Galerkin method and approximation

Consider the functions {v_m; f_m} Î V_m defined in (3). Taking w = j_j(y) and substituting in (22) and (23), we obtain the system of ordinary equations, given by

where

In (25) we have introduced

For numerical reasons, we can rewrite matrix B(t) in the form B = B¹ + B² by using (24), where

3.2 Finite element approximation

We now present a semi-discrete formulation for problem (25) using the Galerkin finite element method to discretize the spatial variable. We first applied the method to find the approximate solution of the exact solution v(y,t) of the Problem (II) and later, using the transformation (21) we can obtain the approximate solution of u(x,t) for the Problem (I) in the domain Q_t.

First, we divide the domain W = (0,1) in local domain W_i = (y_i, y_i+1). Then, W = int and W_iÇW_j = Æ, if i ¹ j. In finite element method, the j_i are piecewise polynomials of some degree in W and vanish on ¶W. More specifically, in this work, we have used the basis function

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 -1 = y₁ < y₂ < ...< y_m+1 = 1. Note that, if |i - j|> 2, then (j _i,j_j) = 0, and (¶j_i/¶y, ¶j_j/¶y) = 0. Hence all the matrices of system are tridiagonal.

Matrix calculation

For each W_i, we have to calculate each integral defined in (26), using the functions (24), (27) and its derivatives. Doing the calculus, we obtain, respectively the following elements, for each tridiagonal matrix A, B¹, B², C, D, E, F, G and R:

3.3 Finite difference method

The equation (25) represent a system of ordinary differential equations of second order and due to matrices characteristics (dependent on the variables y and t) of system, obtaining the solution is not always possible. So, we will apply a numerical method to obtain the approximated solution for the system (25), using the approximate Newmark's Method (see, for instance, Hugles [7], pp 493).

Let dⁿ = d(t_n) and gⁿ = g(t_n) be the approximate solution of the exact solution d(t) and g(t) of (25)_1,2, respectively, where we denote the discrete times in the interval [0,T] by t_n = nDt, n = 0,1...N.

For d> 1/4, with d Î , consider the following approximation

and for the first and second derivative, we take the difference operator in the following form

which, for this approximation the discrete error can be showed to be of order (Dt²).

Coupled system

For the system (25)_1,2 at the discrete mesh points t_n = nDt, using (29) and (30), we obtain the following coupled system:

where,

To determine the solution {dⁿ,gⁿ}, the coupled system of algebraic equations (31) may be solved by iteration, as follows: To start the iteration, we first take n = 0 in (31) and rewrite the system as

where the right-hand side is determined by the (starting) values, since the exact solutions {v(y,t); f(y,t)} are known at time t = 0 and {v⁰; f⁰} are just the initial values, i.e, d⁰ = v⁰(.) = v(.,0), g⁰ = f⁰(.) = f(.,0), where we have used (3) and (27).

We can calculate an approximation for {d^-1; g^-1} by the second order Taylor extrapolation of {v(.,t); f(.,t)} from t⁰ = 0, viz,

in which the values of d'(0), d"0) and g'(0) are given by

calculated from the equation (II)₁ and (II)₂, at t⁰ = 0 and the initial values v⁰(.) = v(.,0) and f⁰(.) = f(.,0).

The system may be solved uniquely for {d¹, g¹}, since its coefficient matrix is non singular. Having determined the values {d¹, g¹}, then for n = 1,2,...N, we obtain the approximate solution {dⁿ⁺¹, gⁿ⁺¹} for the coupled system of algebraic equations (31), which can be rewritten in the form of block matrix,

where,

The system (34) may be solved uniquely, since the matrix is non-singular. In order to solve the system we can use de Gauss Elimination, LU factorization, (see [6, 9]) or Uzwa Method (see [6]).

Note that each square matrix of linear system is of (m-1)th order, since every matrix defined by (32) is of order (m-1). So the linear system is of order 2(m-1)×2(m-1), with the coefficient block matrix and the right-hand side known from the previous iteration.

Uncoupled system

Since {dⁿ, gⁿ} must be solved simultaneously at each time step, the preceding numerical scheme is computationally coupled. From the numerical standpoint the coupled system is larger and hence harder to solve than an uncoupled system involving only dⁿ⁺¹ or only gⁿ at each time step t_n. In order to get uncoupled system, (see [1] and [12]), we replace the central difference by the backward extrapolation for the first derivative,

then by substituting it in the system (25)_1,2 together with (29) and (30), we obtain, after some simple calculation,

where

We start the iteration by taking n = 0 in (36)₂ and then taking n = 0 in (36)₁, to obtain

The terms on the right-hand side of (38)₁, now involve the values {d⁰, d^-1, d^-2, g⁰, g^-1} which are known by (33) and the term d^-2 calculated by taking n = 0 in (35),

Therefore, the value of g¹ is calculated in (38)₁ and the value of d¹ can be determined from (38)₂. Now we can go to iteration steps n = 1,2..., N similarly, by first solving (36)₂ then (36)₁ alternatively in each step. In this manner the numerical system is uncoupled in this computational scheme and we obtain the values {gⁿ,dⁿ} for n = 1,2,...N. These values together with the starting values, constitute the finite element approximate solutions to the initial boundary value problem of Problem (II).

4 Numerical simulation

A numerical example will be given to illustrate some features of the present model, using the method developed for the uncoupled system that is more efficient. In the example, we need the constants h₁, h₂ and k, which give rise to the coupling of the parabolic and hyperbolic equation in the thermoelastic system (I). The constants are given by the following formulas

where c is the specific heat; is the coefficient of thermal expansion; is the thermal conductivity; l = 2K(0) is the length of the string; r is the density of the string; q₀ is the initial temperature; l and m are the coefficients of Lamé,

where n is the coefficient of Poisson and E is the Young's modulus.

For the numerical example, these values will be calculated from the physical properties of aluminum. In this case, we have m = 26.24×10⁹ and l = 58.41×10⁹. Using the thermal and mechanical properties of aluminum, we obtain the approximate values h₁ = 0.164, h₂ = 0.161 and k = 0.177.

Let us consider in (29) the weight d = 0.5 and let W = (-K(t),K(t)) be divided into m subintervals, i.e, h = 2/m and Dt = T/N, for different values of N and T for the discrete time. To calculate the coefficients defined in (24) in each step, the function K(t) that defines the time dependence of the boundary for the non-cylindrical domain Q_t in (20) must be given. In this example it is given by K(t) = 1 - 1/exp^(t+1). Note that in this case, Q_ttends to Q rapidly as t increases. This particular function is taken in order to satisfy the hypothesis H2, i.e, K'(t) » 1. From the physical point of view, we require that the speed of the end points be less than the ''characteristic'' speed of the system.

Note that when only wave equations for small vibrations of elastic string or beam equation, both with moving boundaries, the monotonicity of those functions is not required (see [3], [8], [13]).

We consider the initial temperature, the initial position and velocity given by

In all the figures the space variable is the axis-x, by the change of variable y = (x - a(t))/g(t)).

For Fig. 1-Fig. 4 we have used Dt = 0.03 and h = 0.02, with N = m = 100 and T = 3. Fig. 1 and Fig. 2, respectively, shows the temperature q(x,t) and the displacement u(x,t) in the midpoint x = 0.

Fig. 3 and Fig. 4, show the approximate solutions q(x,t*) and u(x,t*), in the interval [0,T] = [0,3] for different values of t*, t* = 0, 0.25, 0.75, 1.5, 2.25, 3.0.

Note that the interval of the boundary has varied from [-0.63, 0.63] to [-0.98, 0.98].

To obtain Fig. 5-Fig. 8, we have used Dt = 0.1 and h = 0.04. In Fig. 5 and Fig. 6 the evolution of the displacement function u(x,t) is plotted, showing the profile of the displacement, where time varies from 0 to 5 and 0 to 10, at 0.1 interval respectively.

In Fig. 7 and Fig. 8 the evolution of the temperature function q(x,t) are plotted, showing the profile of the temperature, where time varies from 0 to 5 and 0 to 10, at 0.1 interval respectively.

Received: 03/III/04. Accepted: 06/I/05.

#625/04.

• [1] S.I. Chou and C.C. Wang, Estimates of error in finite element approximate solutions to problems in linear thermoelsticity, Part II. Computationally uncoupled numerical schemes, Arch. Rational Mech. Anal., 77 (1981), 263-299.
• [2] J. Clank, The Mathematics of Diffusion, 2rd ed., Oxford University Press (1975).
• [3] H.R. Clark, M.A. Rincon and R.D. Rodrigues, Beam Equation with Weak-Internal Damping in Domain with Moving Boundary, Applied Numerical Mathematics, Vol. 47, Fasc. 2, (2003), 139-157.
• [4] C.M. Dafermos, On the existence and the asymptotic stability of solution to the equations of linear thermoelasticity, Arch. Rational Mech. Anal., 29 (1968), 241-271.
• [5] G. Dassios and M. Grillakis, Dissipation rates and partition of energy in thermoelasticity,  Arch. Rational Mech. Anal., 87-1 (1984), 49-91.
• [6] G.H. Golub and C.F Van Loan, Matrix Computations, 3 rd ed., Johns Hopkins U. Press, Bltimore, (1996).
• [7] T.J.R. Hugles, The finite Element Method Linear Static and Dynamic Finite Element Analysis, Prentice Hall (1987).
• [8] Liu, I-Shih and M.A. Rincon, Effect of Moving Boundaries on the Vibrating Elastic String, Applied Numerical Mathematics, Vol. 47, Fasc. 2, (2003), 159-172.
• [9] Lloyd N. Trefethen and David Bau III, Numerical Linear Algebra, 1 rd ed., SIAM, Philadelphia, (1997)
• [10] L.A. Medeiros, J. Límaco and S.B. Menezes, Vibrations of Elastic Strings: Mathematical Aspects, Part One,  Journal of Computational Analysis and Applications, vol. 4 (2002), 91-127.
• [11] L.A. Medeiros, J. Límaco and S.B. Menezes, Vibrations of Elastic Strings: Mathematical Aspects, Part two,  Journal of Computational Analysis and applications, vol. 4 (2002), 211-263.
• [12] C.A. de Moura, A linear uncoupling numerical scheme for the nonlinear coupled thermoelastodynamics equations, Numerical Methods, V. Pereyra and A. Reinoza (eds.). Lecture Notes in Mathematics n. 1005, Springer-Verlag, 1983.
• [13] M.A. Rincon, Liu; I-Shih, Existence and Uniqueness of Solutions of Elastic String with Moving Ends, Mathematical Methods in the Applied Sciences, Vol. 27 (2004), 1641-1655.

Publication Dates

• Publication in this collection
  20 Apr 2006
• Date of issue
  Dec 2005

History

• Received
  03 Mar 2004
• Accepted
  06 Jan 2005