Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow Changing Domains

*Dmytro V. Yevdokymov and Yuri L. Menshikov*

## **Abstract**

Nowadays, diffusion and heat conduction processes in slow changing domains attract great attention. Slow-phase transitions and growth of biological structures can be considered as examples of such processes. The main difficulty in numerical solutions of correspondent problems is connected with the presence of two time scales. The first one is time scale describing diffusion or heat conduction. The second time scale is connected with the mentioned slow domain evolution. If there is sufficient difference in order of the listed time scale, strong computational difficulties in application of time-stepping algorithms are observed. To overcome the mentioned difficulties, it is proposed to apply a small parameter method for obtaining a new mathematical model, in which the starting parabolic initialboundary-value problem is replaced by a sequence of elliptic boundary-value problems. Application of the boundary element method for numerical solution of the obtained sequence of problems gives an opportunity to solve the whole considered problem in slow time with high accuracy specific to the mentioned algorithm. Besides that, questions about convergence of the obtained asymptotic expansion and correspondence between initial and obtained formulations of the problem are considered separately. The proposed numerical approach is illustrated by several examples of numerical calculations for relevant problems.

**Keywords:** moving boundary problem, Stefan problem, biological tissue growth, asymptotic method, heat conduction equation, diffusion equation, Laplace equation, boundary element method

### **1. Introduction**

Processes in moving boundary domains are commonly known starting from antiquity. Freezing of water with ice creation and an inverse process, when ice melts, evidently show how important this class of phenomena for environmental sciences. Similar processes of solidification and melting create physical background for most technologies of metal or other material productions. Even the mentioned two examples are enough to understand an urgency of the considered problems. However, beside of that, free surface fluid flows, shock wave propagations in gases, growing processes in biological tissues, virus dynamics can be classified as moving boundary problems too. Thus, a number of relevant problems are so many, that they require huge efforts for their mathematical modeling and numerical simulation in different fields of sciences and industries. As a result, the specific direction in numerical analysis was developed to satisfy to numerous demands of practical applications. First of all, the mentioned demands concern accuracy and effectiveness of the correspondent computational schemes.

be nonlinear too. Fortunately, asymptotic approach described above provides line-

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow…*

Taking into account all circumstances and requirements, the following conclusion can be made: only the boundary element method would be effective approach to numerical solution of the obtained asymptotic sequence of elliptic boundaryvalue problems, constructed for moving boundary problems with the governing

Two practically important problems are used as examples in the present work. The first one is well-known and mentioned above Stefan problem, in particular, case of slow phase transition will be considered below, because the condition, that Stefan number is less than 1, provides an effective applicability of asymptotic approach. The second considered problem is biological tissue growth problem, describing specific processes, which are investigated in biological, medical and

Physical theory of phase transformations on microscopic and macroscopic levels were good developed and, as a result, there are not sufficient unsolved questions, what can arise during solution of most of applied heat and mass transfer problems including phase transition [1]. However, there is not so good situation from the point of view of computational mathematics, because phase transition problems contain specific kind of nonlinearity connected with motion of phase transition boundary. As a rule, time-stepping algorithms are used for numerical solution of phase transition problems and the domain shape is fixed on the time step, that is, there is an implicit splitting of the process by the field evolution and interphase boundary motion. Thus such algorithms provide "jumping" domain shape and time step must be enough small to guarantee small domain shape "jump" and high accuracy of the field calculation. Beside of additional time step restriction, there is an additional error source, concerning the domain boundary motion. Any full review of numerical methods of Stefan problem solution requires a special investigation and cannot be included in restricted amount of the present paper. However, the following general conclusion can be made: all mentioned numerical algorithms of Stefan problem solution, based on finite element or finite difference approaches, are rather directed to fast phase transformations, because under restricted time step they require a lot of time steps for slow phase transformations. Then they are found

The quasi-stationary approximation (called Leybenzon approximation in Russian literature) is used to apply for numerical calculations of such processes [1]. However, the number of similar works was very restricted, and they were devoted to engineering design. This approach becomes popular in problems of freezing (melting) of soil, for investigation of phase transitions in solid body, in some evaporation (condensation) problems. The situation in numerical modeling of slow phase transitions was sharply changed in connection with three new problems. The first problem was simple attempt to build more accurate mathematical models for environment processes, for example, in meteorology or soil investigations. The second problem is phase transitions in microgravity conditions, which became important with starting of intensive space explorations. And finally, the third problem was connected with attempts to obtain a material with minimal residual stresses during material production or during phase transitions in solid body, what was important in material sciences. An experience of application of traditional finite difference and finite element approach to the mentioned problems was rather unsuccessful, because their numerical solution required huge computer resources and therefore their research opportunities were strongly restricted. Beside of that, the traditional methods often could not provide necessary accuracy of the numerical solution. On the other hand, the quasi-stationary approximation had difficulties

arization of obtained boundary-value problems.

*DOI: http://dx.doi.org/10.5772/intechopen.93788*

noneffective in the case of slow phase transitions.

parabolic equation.

agricultural sciences.

**49**

Let us consider the specific computational difficulties arising in moving boundary problems, in details. There exist, at least, one physical evolutionary field inside the solution domain, which determines the considered process of boundary motion. This field is described by some mathematical model with specific parameters, some of which would be grouped into dimensionless complexes like Fourier number. At least, one of such parameters must be connected with dimensionless time for evolutionary problem. Motion of solution domain boundary is completely another process, determined by different mathematical model. The mathematical model of boundary motion includes time too and allows its conversion into dimensionless time. Thus, the general problem includes, at least, two-time scales, connected with the governing physical field and with the boundary motion. Since there are not any restrictions on time scales or even connections between them, a relation of the considered time scales can have any value, including extremely small or extremely large. In last case, one of the considered processes can be defined as "slow" and other as "fast". As a rule, the boundary motion is "slow" and the field governing equation time is "fast" in the practical applications, like phase transitions in industrial technologies. To obtain an accurate solution using traditional step-by-step approximation in time, it is necessary to use extremely small-time step (to avoid sufficient approximation errors) during enough large time interval in "slow" time. Since such computational procedure requires a huge number of time steps, it is undesirable due to high-consumed computer resources and possible accumulation of computational errors. To overcome the described difficulty, an asymptotic analysis can be applied to the considered problem. A relation of dimensionless "fast" time to dimensionless "slow" time is dimensionless time-independent constant, which is enough small. It is suitable to use this value as a small parameter in correspondent asymptotic expansions. As a result, the obtained asymptotic sequence of problems gives an opportunity to solve the general problem in "slow" time instead "fast" time, what provides more effective tools of numerical or approximate analytical analysis. For example, heat conduction equation in well-known Stefan problem is replaced by sequence of Laplace equations in two-dimensional and three-dimensional in space cases and by sequence of second order ordinary differential equations in one-dimensional in space case. As a rule, the last one-dimensional case allows analytical integration of the problem.

Generally speaking, arbitrary shapes of domains and their boundaries can arise in the moving boundary problems, what is connected with rebuilding of computational grid at every time step. It means very serious computational difficulties for applications of finite difference and finite element methods, because of their computational opportunities are sufficiently determined by the grid parameters. Beside of that, transfer process of the computed solutions from nodes of old grid into the nodes of new built grid generates hard for checking errors. Boundary element method is often used for such problem because the rebuilding of boundary element grid is sufficiently simpler than similar procedure for finite element and finite difference methods.

At last, moving boundary problems are connected with specific kind of nonlinearity, which restricts a set of effective tools for the problem analysis. Beside of that, the governing equations describing temperature fields inside the phases can

### *Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow… DOI: http://dx.doi.org/10.5772/intechopen.93788*

be nonlinear too. Fortunately, asymptotic approach described above provides linearization of obtained boundary-value problems.

Taking into account all circumstances and requirements, the following conclusion can be made: only the boundary element method would be effective approach to numerical solution of the obtained asymptotic sequence of elliptic boundaryvalue problems, constructed for moving boundary problems with the governing parabolic equation.

Two practically important problems are used as examples in the present work. The first one is well-known and mentioned above Stefan problem, in particular, case of slow phase transition will be considered below, because the condition, that Stefan number is less than 1, provides an effective applicability of asymptotic approach. The second considered problem is biological tissue growth problem, describing specific processes, which are investigated in biological, medical and agricultural sciences.

Physical theory of phase transformations on microscopic and macroscopic levels were good developed and, as a result, there are not sufficient unsolved questions, what can arise during solution of most of applied heat and mass transfer problems including phase transition [1]. However, there is not so good situation from the point of view of computational mathematics, because phase transition problems contain specific kind of nonlinearity connected with motion of phase transition boundary. As a rule, time-stepping algorithms are used for numerical solution of phase transition problems and the domain shape is fixed on the time step, that is, there is an implicit splitting of the process by the field evolution and interphase boundary motion. Thus such algorithms provide "jumping" domain shape and time step must be enough small to guarantee small domain shape "jump" and high accuracy of the field calculation. Beside of additional time step restriction, there is an additional error source, concerning the domain boundary motion. Any full review of numerical methods of Stefan problem solution requires a special investigation and cannot be included in restricted amount of the present paper. However, the following general conclusion can be made: all mentioned numerical algorithms of Stefan problem solution, based on finite element or finite difference approaches, are rather directed to fast phase transformations, because under restricted time step they require a lot of time steps for slow phase transformations. Then they are found noneffective in the case of slow phase transitions.

The quasi-stationary approximation (called Leybenzon approximation in Russian literature) is used to apply for numerical calculations of such processes [1]. However, the number of similar works was very restricted, and they were devoted to engineering design. This approach becomes popular in problems of freezing (melting) of soil, for investigation of phase transitions in solid body, in some evaporation (condensation) problems. The situation in numerical modeling of slow phase transitions was sharply changed in connection with three new problems. The first problem was simple attempt to build more accurate mathematical models for environment processes, for example, in meteorology or soil investigations. The second problem is phase transitions in microgravity conditions, which became important with starting of intensive space explorations. And finally, the third problem was connected with attempts to obtain a material with minimal residual stresses during material production or during phase transitions in solid body, what was important in material sciences. An experience of application of traditional finite difference and finite element approach to the mentioned problems was rather unsuccessful, because their numerical solution required huge computer resources and therefore their research opportunities were strongly restricted. Beside of that, the traditional methods often could not provide necessary accuracy of the numerical solution. On the other hand, the quasi-stationary approximation had difficulties

growing processes in biological tissues, virus dynamics can be classified as moving boundary problems too. Thus, a number of relevant problems are so many, that they require huge efforts for their mathematical modeling and numerical simulation in different fields of sciences and industries. As a result, the specific direction in numerical analysis was developed to satisfy to numerous demands of practical applications. First of all, the mentioned demands concern accuracy and effective-

Let us consider the specific computational difficulties arising in moving boundary problems, in details. There exist, at least, one physical evolutionary field inside the solution domain, which determines the considered process of boundary motion. This field is described by some mathematical model with specific parameters, some of which would be grouped into dimensionless complexes like Fourier number. At least, one of such parameters must be connected with dimensionless time for evolutionary problem. Motion of solution domain boundary is completely another process, determined by different mathematical model. The mathematical model of boundary motion includes time too and allows its conversion into dimensionless time. Thus, the general problem includes, at least, two-time scales, connected with the governing physical field and with the boundary motion. Since there are not any restrictions on time scales or even connections between them, a relation of the considered time scales can have any value, including extremely small or extremely large. In last case, one of the considered processes can be defined as "slow" and other as "fast". As a rule, the boundary motion is "slow" and the field governing equation time is "fast" in the practical applications, like phase transitions in industrial technologies. To obtain an accurate solution using traditional step-by-step approximation in time, it is necessary to use extremely small-time step (to avoid sufficient approximation errors) during enough large time interval in "slow" time. Since such computational procedure requires a huge number of time steps, it is undesirable due to high-consumed computer resources and possible accumulation of computational errors. To overcome the described difficulty, an asymptotic analysis can be applied to the considered problem. A relation of dimensionless "fast" time to dimensionless "slow" time is dimensionless time-independent constant, which is enough small. It is suitable to use this value as a small parameter in correspondent asymptotic expansions. As a result, the obtained asymptotic sequence of problems gives an opportunity to solve the general problem in "slow" time instead "fast" time, what provides more effective tools of numerical or approximate analytical analysis. For example, heat conduction equation in well-known Stefan problem is replaced by sequence of Laplace equations in two-dimensional and three-dimensional in space cases and by sequence of second order ordinary differential equations in one-dimensional in space case. As a rule, the

last one-dimensional case allows analytical integration of the problem.

finite difference methods.

**48**

Generally speaking, arbitrary shapes of domains and their boundaries can arise in the moving boundary problems, what is connected with rebuilding of computational grid at every time step. It means very serious computational difficulties for applications of finite difference and finite element methods, because of their computational opportunities are sufficiently determined by the grid parameters. Beside of that, transfer process of the computed solutions from nodes of old grid into the nodes of new built grid generates hard for checking errors. Boundary element method is often used for such problem because the rebuilding of boundary element grid is sufficiently simpler than similar procedure for finite element and

At last, moving boundary problems are connected with specific kind of nonlinearity, which restricts a set of effective tools for the problem analysis. Beside of that, the governing equations describing temperature fields inside the phases can

ness of the correspondent computational schemes.

*Modeling and Simulation in Engineering - Selected Problems*

too, because it is related to asymptotically slow processes and doesn't take into account initial conditions. Beside of that, elliptical boundary-value problems, which must be numerically solved at every time step of quasi-stationary problem solution, are rather inconvenient for finite difference method. As a result, using of quasi-stationary approximation was very restricted last decades.

It is commonly known that any biological tissue consists of cells. Process of cell reproduction is caution of multicellular tissue growth. The growth process consists of two parts: growth of individual cells and fission of cells, that is the growth process has evidently discrete behavior. Since cells are very small and number of them is very large, consideration of each individual cell is impossible and therefore some averaging is necessary. As a rule, averaging process used in biology is similar to well-known continuous approach in mechanics. According to this approach a multicellular biological tissue is assumed as continues media with distributed sources and some diffusive properties. In fact, cells create a porous media, but pressure difference enough for filtration flow is very seldom presence in the biological tissues, therefore transport phenomena due to filtration flow can be neglected and they are provided by diffusive mechanisms. There is no single quantitative measure of biological tissue metabolism, because a lot of chemical reactions mutually interact. However, the simplest way to formulate a mathematical model for metabolism process is to introduce some numerical value called metabolism intensity and to assume, that any chemical reaction and consequently heat and mass transfer process rate is determined by (in the simplest case it is proportional to) metabolism intensity. As a rule, metabolism intensity is connected by linear relation with velocity of tissue growth. This rule is almost always right for simplest organisms, but metabolism of highest animals is more complex. All mathematical models, which will be developed in the present paper below, will be based on this assumption. Of course, it is phenomenological approach and relation function connecting metabolism intensity and consuming of nutrient substances (excrement production) must be determined experimentally. An evident advantage of such mathematical models (so-called one-parametrical models) is their flexibility and opportunity to take into account different number of

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow…*

*DOI: http://dx.doi.org/10.5772/intechopen.93788*

There are two possible mechanisms of biological tissue growth. The first one is the surface growth and the second one is the volume growth. The intensive cell fission takes place in relatively thin layer near the tissue surface in the case of surface growth. Cells situated inside the tissue have stable metabolism without intensive fission in this case, then their total volume remains constant. Reproduction of all cells takes place in the case of volume growth, although the most intensive

All considered above mathematical models are reduced to initial-boundaryvalue problems for system of diffusion equations with nonlinear sources in moving boundary domain. Motion of the domain boundary is caused by tissue growth, what is described in the work [13], which immediately preceded to the present work; however, the main attention in this paper was paid to the analytical solution of the problem. Boundary element method application for numerical solution of the obtained asymptotic problems, was, rather, pointed out briefly. Nevertheless, the problem statement from the paper [13] is reproduced in the present work.

The main aim of the present paper is to develop a universal effective numerical approach for solution of two physically different diffusive problems in slow moving boundary domains and to show a mathematical and computational similarity of the

Stefan problem is traditionally considered as a successful mathematical model of phase transition process in immovable media. Usually, the Stefan problem does not

concentration fields on different levels of consideration.

**2. Mathematical model of slow-phase transition**

fission, as a rule, takes place near the surface.

considered problems.

**51**

**2.1 Stefan problem formulation**

Velocity of phase transition is described by dimensionless parameter called Stefan number [2], which is relation of thermal energy, spent in heating (cooling) of some phase, to energy spent in phase transformation process. The term "slow" means that the Stefan number is small and therefore there are two different time scales in the problem. Asymptotic approaches give an opportunity to build a mathematical model excluding "fast" time. The first work in this direction was paper [2], where the small parameter method was applied to Stefan problem.

After the first original works of Chuang and Szekely [3, 4] a lot of papers were devoted to boundary element method application to Stefan problem. However general effectiveness of boundary element method for parabolic problems is less than similar effectiveness of finite difference method, what was shown in many papers, see, for example, [5]. Of course, using of some special boundary element method algorithms can improve the situation, but any time-stepping numerical method cannot enough effectively solve the problem of slow phase transition.

Boundary element method [6, 7] has become powerful tool for numerical solution of linear boundary-value problems. It is especially effective in comparison with traditional finite difference method and finite element method for elliptical problems in domains of complex geometrical shape. The main idea of the present paper, concerning the numerical approach, is using of boundary element method for solution of elliptical boundary-value problems, which arise for every approximation on every time step. As a result, an effective computational algorithm is developed, because of well-known advantages of boundary element method such as analysis boundary alone and high accuracy of computations. First time, the idea to use asymptotic approach to a wide class of slow phase transition problems was formulated in the paper [8], but the work [8] was mainly devoted to quasi-stationary approximation and the following terms in asymptotic expansion were not considered there. The approach was developed in the article [9], where application of boundary element method in order to solve the obtained elliptic boundary-value problems was proposed. However, the paper [9] was restricted by quasi-stationary approximation too. Full idea of boundary element application to the considered problem was briefly described in the conference paper [10]. The authors of the present work have to reproduce some results of articles [8, 9], including the initial problem statement and asymptotic formulations, because they are practically unknown for modern scientific community.

Problem of biological tissue growth [11] became very actual item at the present stage of biological science development, because a lot of processes used in agriculture and biotechnology are determined by growth of biological tissue. Beside of that, problem of tumor growth is one of the most important in medicine [11, 12]. Two kinds of circumstances determine the growth process, the first one is genetic properties of tissue and the second one is environmental conditions, for example, nutrition, temperature and so on. Most of investigations concerning a biological growth are based on phenomenological approach, considering biological tissues as a "black box" with experimentally determined properties. The growth is one of such properties of biological tissue. Full review of mathematical modeling in biological sciences requires a separate investigation, which must be sufficiently more than the present paper. Since the growth of biological tissue is the object of the present work, let us consider specific features of the mathematical models of the given processes. General simplifying assumptions must be made to formulate a mathematical model.

### *Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow… DOI: http://dx.doi.org/10.5772/intechopen.93788*

It is commonly known that any biological tissue consists of cells. Process of cell reproduction is caution of multicellular tissue growth. The growth process consists of two parts: growth of individual cells and fission of cells, that is the growth process has evidently discrete behavior. Since cells are very small and number of them is very large, consideration of each individual cell is impossible and therefore some averaging is necessary. As a rule, averaging process used in biology is similar to well-known continuous approach in mechanics. According to this approach a multicellular biological tissue is assumed as continues media with distributed sources and some diffusive properties. In fact, cells create a porous media, but pressure difference enough for filtration flow is very seldom presence in the biological tissues, therefore transport phenomena due to filtration flow can be neglected and they are provided by diffusive mechanisms. There is no single quantitative measure of biological tissue metabolism, because a lot of chemical reactions mutually interact. However, the simplest way to formulate a mathematical model for metabolism process is to introduce some numerical value called metabolism intensity and to assume, that any chemical reaction and consequently heat and mass transfer process rate is determined by (in the simplest case it is proportional to) metabolism intensity. As a rule, metabolism intensity is connected by linear relation with velocity of tissue growth. This rule is almost always right for simplest organisms, but metabolism of highest animals is more complex. All mathematical models, which will be developed in the present paper below, will be based on this assumption. Of course, it is phenomenological approach and relation function connecting metabolism intensity and consuming of nutrient substances (excrement production) must be determined experimentally. An evident advantage of such mathematical models (so-called one-parametrical models) is their flexibility and opportunity to take into account different number of concentration fields on different levels of consideration.

There are two possible mechanisms of biological tissue growth. The first one is the surface growth and the second one is the volume growth. The intensive cell fission takes place in relatively thin layer near the tissue surface in the case of surface growth. Cells situated inside the tissue have stable metabolism without intensive fission in this case, then their total volume remains constant. Reproduction of all cells takes place in the case of volume growth, although the most intensive fission, as a rule, takes place near the surface.

All considered above mathematical models are reduced to initial-boundaryvalue problems for system of diffusion equations with nonlinear sources in moving boundary domain. Motion of the domain boundary is caused by tissue growth, what is described in the work [13], which immediately preceded to the present work; however, the main attention in this paper was paid to the analytical solution of the problem. Boundary element method application for numerical solution of the obtained asymptotic problems, was, rather, pointed out briefly. Nevertheless, the problem statement from the paper [13] is reproduced in the present work.

The main aim of the present paper is to develop a universal effective numerical approach for solution of two physically different diffusive problems in slow moving boundary domains and to show a mathematical and computational similarity of the considered problems.

### **2. Mathematical model of slow-phase transition**

### **2.1 Stefan problem formulation**

Stefan problem is traditionally considered as a successful mathematical model of phase transition process in immovable media. Usually, the Stefan problem does not

too, because it is related to asymptotically slow processes and doesn't take into account initial conditions. Beside of that, elliptical boundary-value problems, which must be numerically solved at every time step of quasi-stationary problem solution,

Velocity of phase transition is described by dimensionless parameter called Stefan number [2], which is relation of thermal energy, spent in heating (cooling) of some phase, to energy spent in phase transformation process. The term "slow" means that the Stefan number is small and therefore there are two different time scales in the problem. Asymptotic approaches give an opportunity to build a mathematical model excluding "fast" time. The first work in this direction was paper [2],

After the first original works of Chuang and Szekely [3, 4] a lot of papers were devoted to boundary element method application to Stefan problem. However general effectiveness of boundary element method for parabolic problems is less than similar effectiveness of finite difference method, what was shown in many papers, see, for example, [5]. Of course, using of some special boundary element method algorithms can improve the situation, but any time-stepping numerical method cannot enough effectively solve the problem of slow phase transition.

Boundary element method [6, 7] has become powerful tool for numerical solution of linear boundary-value problems. It is especially effective in comparison with traditional finite difference method and finite element method for elliptical problems in domains of complex geometrical shape. The main idea of the present paper, concerning the numerical approach, is using of boundary element method for solution of elliptical boundary-value problems, which arise for every approximation on every time step. As a result, an effective computational algorithm is developed, because of well-known advantages of boundary element method such as analysis boundary alone and high accuracy of computations. First time, the idea to use asymptotic approach to a wide class of slow phase transition problems was formulated in the paper [8], but the work [8] was mainly devoted to quasi-stationary approximation and the following terms in asymptotic expansion were not considered there. The approach was developed in the article [9], where application of boundary element method in order to solve the obtained elliptic boundary-value problems was proposed. However, the paper [9] was restricted by quasi-stationary approximation too. Full idea of boundary element application to the considered problem was briefly described in the conference paper [10]. The authors of the present work have to reproduce some results of articles [8, 9], including the initial problem statement and asymptotic formulations, because they are practically

Problem of biological tissue growth [11] became very actual item at the present stage of biological science development, because a lot of processes used in agriculture and biotechnology are determined by growth of biological tissue. Beside of that, problem of tumor growth is one of the most important in medicine [11, 12]. Two kinds of circumstances determine the growth process, the first one is genetic properties of tissue and the second one is environmental conditions, for example, nutrition, temperature and so on. Most of investigations concerning a biological growth are based on phenomenological approach, considering biological tissues as a "black box" with experimentally determined properties. The growth is one of such properties of biological tissue. Full review of mathematical modeling in biological sciences requires a separate investigation, which must be sufficiently more than the present paper. Since the growth of biological tissue is the object of the present work, let us consider specific features of the mathematical models of the given processes. General simplifying assumptions must be made to formulate a mathematical model.

are rather inconvenient for finite difference method. As a result, using of

quasi-stationary approximation was very restricted last decades.

*Modeling and Simulation in Engineering - Selected Problems*

where the small parameter method was applied to Stefan problem.

unknown for modern scientific community.

**50**

In particular, first kind

*DOI: http://dx.doi.org/10.5772/intechopen.93788*

or second kind

initial conditions:

representation

**53**

T1j

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow…*

T2j

λ1 ∂T1 ∂n Г1

λ2 ∂T2 ∂n Г2

T1j

where a1, a2 are thermal diffusivity coefficients of phases, T1s, T2s, q1, q2 are temperatures and thermal fluxes on the outside parts of curves Γ<sup>1</sup> and Γ2, λ1, λ<sup>2</sup> are heat conduction coefficients, ρ is density of some phase (as a rule, incompressible phase, if the both phases are incompressible, the density of initial phase is chosen), Vp*:*t*:* is phase transition boundary propagation velocity. Condition (9) is called Stefan condition, and it describes the motion of the phase transition boundary in dependence on balance of thermal fluxes. Note, that the phase transformation temperature Tp*:*t*:* is constant of the material. To complete the formulation, add the

Let us become to dimensionless variables in the problem, using following

<sup>θ</sup><sup>1</sup> <sup>¼</sup> T1 � Tp*:*t*:* Tn � Tp*:*t*:*

where value Tn is chosen as one from following values max T1s f g , T2s or min T1s f g , T2s (but sometimes it is possible to use max Tf g 10, T20 or min Tf g 10, T20 ) in dependence on the particular problem. The temperature field in the second phase

> ∂θ<sup>1</sup> ∂Fo1

> ∂θ<sup>2</sup> ∂Fo2

is transformed by similar way. Thus, Eqs. (1) and (2) have forms

and boundary conditions on the phase transition boundary

λ1 ∂T1 ∂n Гp*:*t*:* � λ<sup>2</sup> ∂T2 ∂n Гp*:*t*:*

<sup>Г</sup><sup>1</sup> ¼ T1s, (3)

<sup>Г</sup><sup>2</sup> ¼ T2s, (4)

¼ q1, (5)

¼ q2, (6)

¼ ρσVp*:*t*:*, (9)

<sup>Г</sup>p*:*t*:* ¼ Tp*:*t*:*, (7)

T2j<sup>Г</sup>p*:*t*:* ¼ Tp*:*t*:*, (8)

T1ð Þ¼ 0, x T10ð Þ x , (10)

T2ð Þ¼ 0, x T20ð Þ x *:* (11)

, (12)

<sup>¼</sup> <sup>Δ</sup><sup>∗</sup> <sup>θ</sup>1, (13)

<sup>¼</sup> <sup>Δ</sup><sup>∗</sup> <sup>θ</sup>2, (14)

#### **Figure 1.**

*Two-phase Stefan problem. The first phase occupies the domain D*<sup>1</sup> *with boundary Γ*1*, and the second phaseis situated in the domain D*<sup>2</sup> *with boundary Γ*2*. The boundary part Γ<sup>p</sup>:t:* ¼ *Γ*<sup>1</sup> ∩ *Γ*<sup>2</sup> *is moving phase transition boundary.*

#### **Figure 2.**

*The most widespread one-phase Stefan problem. The first phase occupies the domain D*<sup>1</sup> *with boundary Γ*1*, and it surrounds the second phase, which is situated in the domain D*<sup>2</sup> *with boundary Γ*<sup>2</sup> ¼ *Γ<sup>p</sup>:t:* ⊂*Γ*1*. The second phase is kept under phase transition temperature T*<sup>2</sup> ¼ *Tp:t:.*

take into account a change of the substance density under phase transition, what is the main difference between Stefan problem and real phase transition processes. A density change stimulates specific flows in liquid and gaseous phases and it causes an appearing of specific additional stress fields under solid-state phase transitions. However, influence of the mentioned phenomena often can be enough small in comparison with heat conduction processes. Stefan problems can be classified by number of phases (see **Figures 1** and **2**). In particular, one-phase and two-phase Stefan problems are considered in the present work. The term "one-phase problem" does not mean that the second phase is absent, it only points out, that the second phase is kept under constant temperature of phase transition. Generally speaking, one-phase Stefan problem is particular case of two-phase problem, therefore all following formulations will be presented for two-phase problem.

Let consider a classical two-phase Stefan problem. Let the first phase occupies the domain *D*<sup>1</sup> and the second phase occupies the domain *D*2. The boundary Гp*:*t*:* separates mentioned domains, it is a phase transition boundary and has a constant temperature Tp*:*t*:*. Let domains *D*<sup>1</sup> and *D*<sup>2</sup> are bounded by finite or infinite curves Γ<sup>1</sup> and Γ2, correspondingly. For the sake of simplicity, restrict the following consideration by the cases of first and second kind boundary conditions on the outside parts of curves Γ<sup>1</sup> and Γ2. Thus, assuming thermophysical properties of materials as constant, we obtain the following mathematical model [1, 2].

$$\frac{\partial \mathbf{T}\_1}{\partial \mathbf{\tau}} = \mathbf{a}\_1 \Delta \mathbf{T}\_1,\tag{1}$$

$$\frac{\partial \mathbf{T}\_2}{\partial \mathbf{\tau}} = \mathbf{a}\_2 \Delta \mathbf{T}\_2,\tag{2}$$

Here Eqs. (1) and (2) describe the temperature fields inside the first phase and the second phase correspondingly. Boundary conditions of the first or second kinds must be prescribed on the outside parts of curves Γ<sup>1</sup> and Γ2:

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow… DOI: http://dx.doi.org/10.5772/intechopen.93788*

In particular, first kind

$$\mathbf{T\_1}|\_{\Gamma\_1} = \mathbf{T\_{1s}},\tag{3}$$

$$\mathbf{T\_2}|\_{\Gamma\_2} = \mathbf{T\_{2s}},\tag{4}$$

or second kind

$$
\lambda\_1 \frac{\partial \mathbf{T}\_1}{\partial \mathbf{n}} \Big|\_{\Gamma\_1} = \mathbf{q}\_1,\tag{5}
$$

$$
\lambda\_2 \frac{\partial \mathbf{T}\_2}{\partial \mathbf{n}} \bigg|\_{\Gamma\_2} = \mathbf{q}\_2,\tag{6}
$$

and boundary conditions on the phase transition boundary

$$\mathbf{T\_{1}}|\_{\Gamma\_{\rm p.t.}} = \mathbf{T\_{p.t.}}\tag{7}$$

$$\mathbf{T\_2}|\_{\Gamma\_{\rm p.t.}} = \mathbf{T\_{p.t.}},\tag{8}$$

$$
\lambda\_1 \frac{\partial \mathbf{T}\_1}{\partial \mathbf{n}} \bigg|\_{\Gamma\_{\rm p.t.}} - \lambda\_2 \frac{\partial \mathbf{T}\_2}{\partial \mathbf{n}} \bigg|\_{\Gamma\_{\rm p.t.}} = \rho \boldsymbol{\sigma} \mathbf{V}\_{\rm p.t.}, \tag{9}
$$

where a1, a2 are thermal diffusivity coefficients of phases, T1s, T2s, q1, q2 are temperatures and thermal fluxes on the outside parts of curves Γ<sup>1</sup> and Γ2, λ1, λ<sup>2</sup> are heat conduction coefficients, ρ is density of some phase (as a rule, incompressible phase, if the both phases are incompressible, the density of initial phase is chosen), Vp*:*t*:* is phase transition boundary propagation velocity. Condition (9) is called Stefan condition, and it describes the motion of the phase transition boundary in dependence on balance of thermal fluxes. Note, that the phase transformation temperature Tp*:*t*:* is constant of the material. To complete the formulation, add the initial conditions:

$$\mathbf{T}\_1(\mathbf{0}, \mathbf{x}) = \mathbf{T}\_{10}(\mathbf{x}),\tag{10}$$

$$\mathbf{T\_2(0,x) = T\_{20}(x)}.\tag{11}$$

Let us become to dimensionless variables in the problem, using following representation

$$\Theta\_1 = \frac{\mathbf{T\_1} - \mathbf{T\_{p.t.}}}{\mathbf{T\_n} - \mathbf{T\_{p.t.}}},\tag{12}$$

where value Tn is chosen as one from following values max T1s f g , T2s or min T1s f g , T2s (but sometimes it is possible to use max Tf g 10, T20 or min Tf g 10, T20 ) in dependence on the particular problem. The temperature field in the second phase is transformed by similar way. Thus, Eqs. (1) and (2) have forms

$$\frac{\partial \Theta\_1}{\partial \mathbf{F} \mathbf{o}\_1} = \boldsymbol{\Delta}^\* \boldsymbol{\Theta}\_1,\tag{13}$$

$$\frac{\partial \Theta\_2}{\partial \mathbf{F} \mathbf{o}\_2} = \Delta^\* \theta\_2,\tag{14}$$

take into account a change of the substance density under phase transition, what is the main difference between Stefan problem and real phase transition processes. A density change stimulates specific flows in liquid and gaseous phases and it causes an appearing of specific additional stress fields under solid-state phase transitions. However, influence of the mentioned phenomena often can be enough small in comparison with heat conduction processes. Stefan problems can be classified by number of phases (see **Figures 1** and **2**). In particular, one-phase and two-phase Stefan problems are considered in the present work. The term "one-phase problem" does not mean that the second phase is absent, it only points out, that the second phase is kept under constant temperature of phase transition. Generally speaking, one-phase Stefan problem is particular case of two-phase problem, therefore all

*The most widespread one-phase Stefan problem. The first phase occupies the domain D*<sup>1</sup> *with boundary Γ*1*, and it surrounds the second phase, which is situated in the domain D*<sup>2</sup> *with boundary Γ*<sup>2</sup> ¼ *Γ<sup>p</sup>:t:* ⊂*Γ*1*. The second*

*Two-phase Stefan problem. The first phase occupies the domain D*<sup>1</sup> *with boundary Γ*1*, and the second phaseis situated in the domain D*<sup>2</sup> *with boundary Γ*2*. The boundary part Γ<sup>p</sup>:t:* ¼ *Γ*<sup>1</sup> ∩ *Γ*<sup>2</sup> *is moving phase transition boundary.*

Let consider a classical two-phase Stefan problem. Let the first phase occupies the domain *D*<sup>1</sup> and the second phase occupies the domain *D*2. The boundary Гp*:*t*:* separates mentioned domains, it is a phase transition boundary and has a constant temperature Tp*:*t*:*. Let domains *D*<sup>1</sup> and *D*<sup>2</sup> are bounded by finite or infinite curves Γ<sup>1</sup> and Γ2, correspondingly. For the sake of simplicity, restrict the following consideration by the cases of first and second kind boundary conditions on the outside parts of curves Γ<sup>1</sup> and Γ2. Thus, assuming thermophysical properties of materials as

<sup>∂</sup><sup>τ</sup> <sup>¼</sup> a1ΔT1, (1)

<sup>∂</sup><sup>τ</sup> <sup>¼</sup> a2ΔT2, (2)

following formulations will be presented for two-phase problem.

*phase is kept under phase transition temperature T*<sup>2</sup> ¼ *Tp:t:.*

*Modeling and Simulation in Engineering - Selected Problems*

**Figure 1.**

**Figure 2.**

**52**

constant, we obtain the following mathematical model [1, 2].

must be prescribed on the outside parts of curves Γ<sup>1</sup> and Γ2:

*∂*T1

∂T2

Here Eqs. (1) and (2) describe the temperature fields inside the first phase and the second phase correspondingly. Boundary conditions of the first or second kinds where

$$\mathbf{F}\mathbf{o}\_1 = \frac{\mathbf{r}\mathbf{a}\_1}{\mathbf{L}^2}.\tag{15}$$

where St is the Stefan number determined as

Eq. (14) is transformed by following way:

phase transformation σ is quite large value, St<1.

**2.2 Asymptotic approach to Stefan problem**

transformation St< <1 it must converge quickly.

solution in a form of series

**55**

τλ1ð Þ Tn�Tp*:*t*:* L2σρ τa1 L2

> ∂θ<sup>2</sup> ∂τst

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow…*

<sup>f</sup> <sup>a</sup> <sup>¼</sup> Fo1 Fo2 ¼ a1 a2

Physical sense of the Stefan number is quite simple; it is equal to relation of heat, spent to heating (cooling) of the phase, to heat spent to phase transition. But in the same time, it can be considered as relation of two-time scales, concerning heating (cooling) and phase transition boundary motion. Since, as a rule, the latent heat of

Formulated above Stefan problem can be easy generalized for any number of phases, but one-phase Stefan problem (see **Figure 2**) has special importance for the theory. This particular case of general Stefan problem physically corresponds to situation, when one phase has temperature equal to the phase transformation temperature and it undergoes the phase transformation due to heat conduction of the another phase. The phase under constant temperature must be excluded from

> ∂θ ∂τst

> > θj

∂θ <sup>∂</sup><sup>n</sup> <sup>¼</sup> <sup>∂</sup><sup>n</sup> ∂τst

As it was noted earlier, the considered Stefan problem is nonlinear, because of complicated dependence between the temperature field and shape (motion) of the phase transition boundary. Although the question on existence and uniqueness of Stefan problem solution is still far from a complete clearness [1], the continuous dependency of Stefan problem solution on Stefan number is evident. It could be clear shown, if an integral formulation for temperature field is considered (see [6, 7]), besides of that, it is seen from the same integral formulation, that the solution can be continuously differentiated with respect to Stefan number any number of times. Then the temperature field can represented as a series of Stefan number orders, which must converge for St< 1, and in the case of slow phase

On the base of above conclusions on the solution properties we shall search

<sup>¼</sup> <sup>λ</sup><sup>1</sup> Tn � Tp*:*t*:* 

a1σρ <sup>¼</sup> C Tn � Tp*:*t*:*

 σ

faSt ¼ Δθ2, (29)

St ¼ Δθ, (31)

θj<sup>Г</sup><sup>1</sup> ¼ θs, (32)

<sup>Г</sup>p*:*t*:* ¼ 0, (33)

*:* (35)

θð Þ¼ 0, x θ0, (34)

*:* (30)

, (28)

St <sup>¼</sup> <sup>τ</sup>st Fo1 ¼

*DOI: http://dx.doi.org/10.5772/intechopen.93788*

where

the consideration. Thus

The asterisk "\*" in Laplace operator means the differentiation with respect to dimensionless coordinates

$$\mathbf{X}^\* = \mathbf{x}/\mathbf{L}.\tag{16}$$

It will be omitted in following. Dimensionless forms of the boundary and initial conditions are

$$
\Theta\_1|\_{\Gamma\_1} = \Theta\_{1s},
\tag{17}
$$

$$
\Theta\_2|\_{\Gamma\_2} = \Theta\_{2s},
\tag{18}
$$

$$\Theta\_1|\_{\Gamma\_{\rm pt.}} = \mathbf{0},\tag{19}$$

$$\left. \Theta\_{2} \right|\_{\Gamma\_{\rm pt.}} = \mathbf{0}, \tag{20}$$

$$
\theta\_1(\mathbf{0}, \mathbf{x}) = \theta\_{10}, \tag{21}
$$

$$
\Theta\_2(\mathbf{0}, \mathbf{x}) = \Theta\_{20}. \tag{22}
$$

The dimensionless Stefan condition must be considered in details. Let

$$\mathbf{V}\_{\text{p.t.}} = \frac{\partial \mathbf{n}}{\partial \mathbf{r}},\tag{23}$$

where <sup>∂</sup><sup>n</sup> <sup>∂</sup><sup>τ</sup> means the normal velocity of the phase transformation boundary motion. If n<sup>∗</sup> <sup>¼</sup> <sup>n</sup>*=*L, we have following dimensionless relation

$$\frac{\partial\Theta\_1}{\partial\mathbf{n}} - \mathbf{f}\_{\lambda}\frac{\partial\Theta\_2}{\partial\mathbf{n}} = \frac{\partial\mathbf{n}^\*}{\partial\mathbf{r}\_{\rm st}},\tag{24}$$

where

$$\mathbf{f\_{\lambda}} = \lambda\_2 / \lambda\_1,\tag{25}$$

and

$$\tau\_{\rm st} = \frac{\tau \lambda\_1 \left(\mathbf{T\_n} - \mathbf{T\_{p.t.}}\right)}{\mathbf{L^2 \sigma \rho}},\tag{26}$$

τst is dimensionless time, connected with the phase transformation process.

Thus there three dimensionless time Fo1, Fo2 and τst in the problem, however it is desirable to use only one time scale. Since the main interest on the problem is phase transition process, choose the Stefan time τst as the main time scale. Beside of that, as a rule, τst is the "most slow" dimensionless time. Then Eq. (13) must be replaced by

$$\frac{\partial \Theta\_1}{\partial \tau\_{\rm st}} \mathbf{St} = \Delta \theta\_1,\tag{27}$$

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow… DOI: http://dx.doi.org/10.5772/intechopen.93788*

where St is the Stefan number determined as

$$\text{St} = \frac{\tau\_{\text{st}}}{\text{Fo}\_1} = \frac{\frac{\pi \lambda\_1 \text{(T}\_\text{n} - \text{T}\_{\text{p.t.}})}{\text{L}^2 \sigma \rho}}{\frac{\pi \text{a}\_1}{\text{L}^2}} = \frac{\lambda\_1 \text{(T}\_\text{n} - \text{T}\_{\text{p.t.}})}{\text{a}\_1 \sigma \rho} = \frac{\mathbf{C} \left( \text{T}\_\text{n} - \text{T}\_{\text{p.t.}} \right)}{\sigma},\tag{28}$$

Eq. (14) is transformed by following way:

$$\frac{\partial \Theta\_2}{\partial \mathbf{r}\_{\rm st}} \mathbf{f}\_{\rm a} \mathbf{St} = \Delta \Theta\_2,\tag{29}$$

where

where

conditions are

where <sup>∂</sup><sup>n</sup>

where

and

replaced by

**54**

dimensionless coordinates

*Modeling and Simulation in Engineering - Selected Problems*

Fo1 <sup>¼</sup> <sup>τ</sup>a1

The asterisk "\*" in Laplace operator means the differentiation with respect to

It will be omitted in following. Dimensionless forms of the boundary and initial

θ1j

θ2j

θ1j

θ2j

The dimensionless Stefan condition must be considered in details. Let

motion. If n<sup>∗</sup> <sup>¼</sup> <sup>n</sup>*=*L, we have following dimensionless relation

∂θ<sup>1</sup> <sup>∂</sup><sup>n</sup> � <sup>f</sup> <sup>λ</sup>

Vp*:*t*:* <sup>¼</sup> <sup>∂</sup><sup>n</sup>

∂θ<sup>2</sup> <sup>∂</sup><sup>n</sup> <sup>¼</sup> <sup>∂</sup>n<sup>∗</sup> ∂τst

<sup>τ</sup>st <sup>¼</sup> τλ<sup>1</sup> Tn � Tp*:*t*:*

τst is dimensionless time, connected with the phase transformation process. Thus there three dimensionless time Fo1, Fo2 and τst in the problem, however it is desirable to use only one time scale. Since the main interest on the problem is phase transition process, choose the Stefan time τst as the main time scale. Beside of that, as a rule, τst is the "most slow" dimensionless time. Then Eq. (13) must be

> ∂θ<sup>1</sup> ∂τst

L2

<sup>∂</sup><sup>τ</sup> means the normal velocity of the phase transformation boundary

L2 *:* (15)

<sup>X</sup><sup>∗</sup> <sup>¼</sup> <sup>x</sup>*=*L*:* (16)

<sup>Г</sup><sup>1</sup> ¼ θ1s, (17)

<sup>Г</sup><sup>2</sup> ¼ θ2s, (18)

<sup>Г</sup>p*:*t*:* ¼ 0, (19)

<sup>Г</sup>p*:*t*:* ¼ 0, (20)

<sup>∂</sup><sup>τ</sup> , (23)

f<sup>λ</sup> ¼ λ2*=*λ1, (25)

σρ , (26)

St ¼ Δθ1, (27)

, (24)

θ1ð Þ¼ 0, x θ10, (21) θ2ð Þ¼ 0, x θ20*:* (22)

$$\mathbf{f\_a} = \frac{\mathbf{Fo\_1}}{\mathbf{Fo\_2}} = \frac{\mathbf{a\_1}}{\mathbf{a\_2}}.\tag{30}$$

Physical sense of the Stefan number is quite simple; it is equal to relation of heat, spent to heating (cooling) of the phase, to heat spent to phase transition. But in the same time, it can be considered as relation of two-time scales, concerning heating (cooling) and phase transition boundary motion. Since, as a rule, the latent heat of phase transformation σ is quite large value, St<1.

Formulated above Stefan problem can be easy generalized for any number of phases, but one-phase Stefan problem (see **Figure 2**) has special importance for the theory. This particular case of general Stefan problem physically corresponds to situation, when one phase has temperature equal to the phase transformation temperature and it undergoes the phase transformation due to heat conduction of the another phase. The phase under constant temperature must be excluded from the consideration. Thus

$$\frac{\partial \Theta}{\partial \mathbf{\tau}\_{\rm st}} \mathbf{S} \mathbf{t} = \Delta \theta,\tag{31}$$

$$
\left.\Theta\right|\_{\Gamma\_1} = \Theta\_\mathbf{s},
\tag{32}
$$

$$\left.\Theta\right|\_{\Gamma\_{\rm p.t.}} = \mathbf{0},\tag{33}$$

$$
\theta(\mathbf{0}, \mathbf{x}) = \theta\_0,\tag{34}
$$

$$\frac{\partial \Theta}{\partial \mathbf{n}} = \frac{\partial \mathbf{n}}{\partial \mathbf{\tau}\_{\rm st}}.\tag{35}$$

#### **2.2 Asymptotic approach to Stefan problem**

As it was noted earlier, the considered Stefan problem is nonlinear, because of complicated dependence between the temperature field and shape (motion) of the phase transition boundary. Although the question on existence and uniqueness of Stefan problem solution is still far from a complete clearness [1], the continuous dependency of Stefan problem solution on Stefan number is evident. It could be clear shown, if an integral formulation for temperature field is considered (see [6, 7]), besides of that, it is seen from the same integral formulation, that the solution can be continuously differentiated with respect to Stefan number any number of times. Then the temperature field can represented as a series of Stefan number orders, which must converge for St< 1, and in the case of slow phase transformation St< <1 it must converge quickly.

On the base of above conclusions on the solution properties we shall search solution in a form of series

*Modeling and Simulation in Engineering - Selected Problems*

$$
\boldsymbol{\Theta}\_{1} = \boldsymbol{\Theta}\_{1}^{\;0}(\mathbf{x}, \boldsymbol{\pi}) + \sum\_{\mathbf{k}=1}^{\infty} \mathbf{S} \mathbf{t}^{\mathbf{k}} \boldsymbol{\Theta}\_{1}{}^{\mathbf{k}}(\mathbf{x}, \boldsymbol{\pi}), \tag{36}
$$

θ2 0� �

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow…*

θ2 i � �

θ2 i � �

problems (40)–(45) are solved, the relation (24) can be considered as ordinary differential equation, what describes motion of phase transformation boundary. If Eq. (24) is integrated analytically (it is very seldom case, because the right hand part of relation (24) depends on phase boundary position in a complicated way), the general problem could be considered as completely solved. In general case Cauchy problem for Eq. (24) can be solved numerically by some numerical method. However, the most successful approach is method proposed in the paper [2] and based on expansion of left-hand side of Eq. (24) into a series with respect to Stefan number.

∂θ<sup>2</sup>

By the same reasons as earlier let equate factors at equal orders of Stefan number

∂θ<sup>2</sup> 0 <sup>∂</sup><sup>n</sup> <sup>¼</sup> <sup>∂</sup>η<sup>0</sup>

∂θ<sup>2</sup> i <sup>∂</sup><sup>n</sup> <sup>¼</sup> <sup>∂</sup>η<sup>i</sup>

where functions η<sup>k</sup> are subjects to determination. The relationships (56) must be complemented by some initial conditions and can be considered as series of Cauchy problems for phase transition boundary position. The first equations and boundary conditions in (40)–(53), which describes the temperature approximations indicated by index "0", coincide with well-known quasi-stationary approximation. Note that for one-dimensional (in space) case for any Stefan problems the boundary-value problems indicated by "0" and "1" can be integrated analytically but the following approximations requires some numerical method for solution of mentioned Cauchy problem. Such solutions were built for one-phase and two-phase one-dimensional (in ordinary Cartesian, polar and spherical coordinate systems) Stefan problems with different boundary conditions on the outer boundaries. There are particular cases of one-dimensional Stefan problems, analytical solutions of which are known. The mentioned analytical solutions were used to check the approach accuracy. As a

*:*

<sup>0</sup>ð Þ x, <sup>τ</sup> ∂n

þX<sup>∞</sup> k¼1

!

St<sup>k</sup> <sup>∂</sup>θ<sup>2</sup>

<sup>k</sup>ð Þ x, <sup>τ</sup> ∂n

<sup>∂</sup><sup>τ</sup> , (55)

<sup>∂</sup><sup>τ</sup> , (56)

¼

(54)

∂θ<sup>1</sup>

and obtain

**57**

<sup>0</sup>ð Þ x, <sup>τ</sup> ∂n

<sup>¼</sup> <sup>∂</sup>η<sup>0</sup>ð Þ x, <sup>τ</sup> ∂τ

þX<sup>∞</sup> k¼1

St<sup>k</sup> <sup>∂</sup>θ<sup>1</sup>

þX<sup>∞</sup> k¼1

*DOI: http://dx.doi.org/10.5772/intechopen.93788*

!

<sup>k</sup>ð Þ x, <sup>τ</sup> <sup>∂</sup><sup>n</sup> <sup>¼</sup> fa

St<sup>k</sup> <sup>∂</sup>η<sup>k</sup>ð Þ x, <sup>τ</sup> ∂τ

> ∂θ<sup>1</sup> 0 <sup>∂</sup><sup>n</sup> � <sup>f</sup><sup>λ</sup>

∂θ<sup>1</sup> i <sup>∂</sup><sup>n</sup> � <sup>f</sup> <sup>λ</sup>

result, it is shown that accuracy of test problem solutions is enough high.

The presented series of problem (40)–(53), (55), and (56) is equivalent to the initial Stefan problem (1)–(4), (7)–(11). Under restrictions imposed on functions, included into the formulation, concerning their physical sense, (for example, requirement of piecewise smooth boundary of the finite domain D, where the problem is solved, and requirement of boundedness of temperatures and thermal fluxes) the temperature functions inside domains D1 and D2 are differentiable infinite number of times [2]. Moreover, temperature derivatives with respect to time are finite too. Last assumption has physical sense, if the moment of phase

Let consider the Stefan condition (24) separately. If the series of boundary-value

<sup>Г</sup>p*:*t*:* ¼ 0, (51)

<sup>Г</sup><sup>2</sup> ¼ 0, (52)

<sup>Г</sup>p*:*t*:* ¼ 0*:* (53)

$$
\Theta\_2 = \Theta\_2^{\ 0}(\mathbf{x}, \mathbf{\tau}) + \sum\_{\mathbf{k}=1}^{\infty} \mathbf{S} \mathbf{t}^{\mathbf{k}} \theta\_2^{\ \mathbf{k}}(\mathbf{x}, \mathbf{\tau}). \tag{37}
$$

As a result of using of expansions (36) and (37), the initial problem is reduced to determination of function series θ<sup>0</sup> <sup>1</sup> , … , θ<sup>k</sup> <sup>1</sup> , … ; θ<sup>0</sup> <sup>2</sup> , … , θ<sup>k</sup> <sup>2</sup> , … .

Let substitute expansions (36) and (37) into Eqs. (27) and (29) and obtain

$$\text{St}\frac{\partial}{\partial\mathbf{\tau}}\boldsymbol{\Theta}\_{1}^{\,0}(\mathbf{x},\boldsymbol{\tau}) + \text{St}\frac{\partial}{\partial\boldsymbol{\tau}}\sum\_{\mathbf{k}=1}^{\infty}\text{St}^{\mathbf{k}}\boldsymbol{\Theta}\_{1}^{\,\mathbf{k}}(\mathbf{x},\boldsymbol{\tau}) = \boldsymbol{\Delta}\boldsymbol{\Theta}\_{1}^{\,0}(\mathbf{x},\boldsymbol{\tau}) + \boldsymbol{\Delta}\sum\_{\mathbf{k}=1}^{\infty}\text{St}^{\mathbf{k}}\boldsymbol{\Theta}\_{1}^{\,\mathbf{k}}(\mathbf{x},\boldsymbol{\tau}), \qquad \text{(38)}$$

$$\text{f}\_{\mathbf{s}}\text{St}\frac{\partial}{\partial\boldsymbol{\tau}}\boldsymbol{\Theta}\_{2}^{\,0}(\mathbf{x},\boldsymbol{\tau}) + \text{f}\_{\mathbf{s}}\text{St}\frac{\partial}{\partial\boldsymbol{\tau}}\sum\_{\mathbf{k}=1}^{\infty}\text{St}^{\mathbf{k}}\boldsymbol{\Theta}\_{2}^{\,\mathbf{k}}(\mathbf{x},\boldsymbol{\tau}) = \boldsymbol{\Delta}\boldsymbol{\Theta}\_{2}^{\,0}(\mathbf{x},\boldsymbol{\tau}) + \boldsymbol{\Delta}\sum\_{\mathbf{k}=1}^{\infty}\text{St}^{\mathbf{k}}\boldsymbol{\Theta}\_{2}^{\,\mathbf{k}}(\mathbf{x},\boldsymbol{\tau}).$$

If the functions θ<sup>k</sup> <sup>1</sup> , θ<sup>k</sup> <sup>2</sup> are bounded, the series (36) and (37) absolutely converge; therefore, correspondent series can be differentiated term by term. Thus

$$\text{St}\frac{\partial\theta\_1^{0}}{\partial\tau} + \sum\_{\mathbf{k}=1}^{\infty} \text{St}^{\mathbf{k}+1} \frac{\partial\theta\_1^{\mathbf{k}}}{\partial\tau} = \Delta\theta\_1^{0} + \sum\_{\mathbf{k}=1}^{\infty} \text{St}^{\mathbf{k}} \Delta\theta\_1^{\mathbf{k}},\tag{39}$$

$$\text{f}\_{\mathbf{a}}\text{St}\frac{\partial\theta\_2^{0}}{\partial\tau} + \text{f}\_{\mathbf{a}}\sum\_{\mathbf{k}=1}^{\infty} \text{St}^{\mathbf{k}+1} \frac{\partial\theta\_2^{\mathbf{k}}}{\partial\tau} = \Delta\theta\_2^{0} + \sum\_{\mathbf{k}=1}^{\infty} \text{St}^{\mathbf{k}} \Delta\theta\_2^{\mathbf{k}}.$$

Since the value of Stefan number is, generally speaking, arbitrary, equality (39) takes place only in case of equality of factors at equal orders of Stefan number. Hence

$$
\Delta\theta\_1^{\;0} = \mathbf{0},
\tag{40}
$$

$$
\Delta\theta\_1^{-1} = \frac{\partial\theta\_1^{-0}}{\partial\mathbf{r}}\,,\tag{41}
$$

$$
\Delta\Theta\_1^{\,i} = \frac{\partial\Theta\_1^{\,i-1}}{\partial\pi},\tag{42}
$$

$$
\Delta\Theta\_2^{\,0} = \mathbf{0},
\tag{43}
$$

$$
\Delta\Theta\_2^{-1} = \mathbf{f}\_\mathbf{a} \frac{\partial\Theta\_2^{-0}}{\partial\mathbf{\tau}},\tag{44}
$$

$$
\Delta\Theta\_1^{\ i} = \mathbf{f}\_{\text{a}} \frac{\partial\Theta\_1^{\ i-1}}{\partial\pi},
\tag{45}
$$

Let prescribe boundary conditions for Eqs. (40)–(45)

$$\left.\Theta\_1^{\;0}\right|\_{\Gamma\_1} = \Theta\_{1s},\tag{46}$$

$$\left.\Theta\_1^{\,\,0}\right|\_{\Gamma\_{\rm p.t.}} = \mathbf{0},\tag{47}$$

$$\left.\Theta\_1^{\,\,\mathrm{i}}\right|\_{\Gamma\_1} = \mathbf{0},\tag{48}$$

$$\left.\Theta\_{\mathbf{1}}\right\|\_{\Gamma\_{\mathrm{pt.}}} = \mathbf{0},\tag{49}$$

$$\left.\Theta\_2^{\,0}\right|\_{\Gamma\_2} = \Theta\_{2s},\tag{50}$$

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow… DOI: http://dx.doi.org/10.5772/intechopen.93788*

$$\left.\Theta\_{2}{}^{0}\right|\_{\Gamma\_{\text{pt.}}} = \mathbf{0},\tag{51}$$

$$\left.\Theta\_2^{\,\,i}\right|\_{\Gamma\_2} = \mathbf{0},\tag{52}$$

$$\left.\Theta\_{\mathsf{Z}}{}^{i}\right|\_{\Gamma\_{\mathsf{P}t.}} = \mathsf{0}.\tag{53}$$

Let consider the Stefan condition (24) separately. If the series of boundary-value problems (40)–(45) are solved, the relation (24) can be considered as ordinary differential equation, what describes motion of phase transformation boundary. If Eq. (24) is integrated analytically (it is very seldom case, because the right hand part of relation (24) depends on phase boundary position in a complicated way), the general problem could be considered as completely solved. In general case Cauchy problem for Eq. (24) can be solved numerically by some numerical method. However, the most successful approach is method proposed in the paper [2] and based on expansion of left-hand side of Eq. (24) into a series with respect to Stefan number.

$$\begin{split} & \frac{\partial \boldsymbol{\Theta}\_{1}^{\boldsymbol{\Phi}}(\mathbf{x},\boldsymbol{\tau})}{\partial \mathbf{n}} + \sum\_{\mathbf{k}=1}^{\boldsymbol{\Theta}} \operatorname{St}^{\mathbf{k}} \frac{\partial \boldsymbol{\Phi}\_{1}^{\boldsymbol{\Phi}}(\mathbf{x},\boldsymbol{\tau})}{\partial \mathbf{n}} = \operatorname{f}\_{\mathbf{z}} \left( \frac{\partial \boldsymbol{\Theta}\_{2}^{\boldsymbol{\Phi}}(\mathbf{x},\boldsymbol{\tau})}{\partial \mathbf{n}} + \sum\_{\mathbf{k}=1}^{\boldsymbol{\Theta}} \operatorname{St}^{\mathbf{k}} \frac{\partial \boldsymbol{\Theta}\_{2}^{\boldsymbol{\Phi}}(\mathbf{x},\boldsymbol{\tau})}{\partial \mathbf{n}} \right) = \\ & = \left( \frac{\partial \boldsymbol{\eta}^{0}(\mathbf{x},\boldsymbol{\tau})}{\partial \boldsymbol{\tau}} + \sum\_{\mathbf{k}=1}^{\boldsymbol{\Theta}} \operatorname{St}^{\mathbf{k}} \frac{\partial \boldsymbol{\eta}^{k}(\mathbf{x},\boldsymbol{\tau})}{\partial \boldsymbol{\tau}} \right). \end{split} \tag{54}$$

By the same reasons as earlier let equate factors at equal orders of Stefan number and obtain

$$\frac{\partial \boldsymbol{\Theta}\_{1}^{0}}{\partial \mathbf{n}} - \mathbf{f}\_{\boldsymbol{\lambda}} \frac{\partial \boldsymbol{\Theta}\_{2}^{0}}{\partial \mathbf{n}} = \frac{\partial \boldsymbol{\eta}^{0}}{\partial \boldsymbol{\pi}},\tag{55}$$

$$\frac{\partial \Theta\_1^{\,\,i}}{\partial \mathbf{n}} - \mathbf{f}\_{\lambda} \frac{\partial \Theta\_2^{\,\,i}}{\partial \mathbf{n}} = \frac{\partial \eta^{\,i}}{\partial \mathbf{\tau}},\tag{56}$$

where functions η<sup>k</sup> are subjects to determination. The relationships (56) must be complemented by some initial conditions and can be considered as series of Cauchy problems for phase transition boundary position. The first equations and boundary conditions in (40)–(53), which describes the temperature approximations indicated by index "0", coincide with well-known quasi-stationary approximation. Note that for one-dimensional (in space) case for any Stefan problems the boundary-value problems indicated by "0" and "1" can be integrated analytically but the following approximations requires some numerical method for solution of mentioned Cauchy problem. Such solutions were built for one-phase and two-phase one-dimensional (in ordinary Cartesian, polar and spherical coordinate systems) Stefan problems with different boundary conditions on the outer boundaries. There are particular cases of one-dimensional Stefan problems, analytical solutions of which are known. The mentioned analytical solutions were used to check the approach accuracy. As a result, it is shown that accuracy of test problem solutions is enough high.

The presented series of problem (40)–(53), (55), and (56) is equivalent to the initial Stefan problem (1)–(4), (7)–(11). Under restrictions imposed on functions, included into the formulation, concerning their physical sense, (for example, requirement of piecewise smooth boundary of the finite domain D, where the problem is solved, and requirement of boundedness of temperatures and thermal fluxes) the temperature functions inside domains D1 and D2 are differentiable infinite number of times [2]. Moreover, temperature derivatives with respect to time are finite too. Last assumption has physical sense, if the moment of phase

θ<sup>1</sup> ¼ θ<sup>1</sup>

*Modeling and Simulation in Engineering - Selected Problems*

θ<sup>2</sup> ¼ θ<sup>2</sup>

Stkθ<sup>1</sup>

∂τ X∞ k¼1

þX<sup>∞</sup> k¼1

þ fa X∞ k¼1

determination of function series θ<sup>0</sup>

<sup>0</sup>ð Þþ x, <sup>τ</sup> St <sup>∂</sup>

∂τ X∞ k¼1

<sup>0</sup>ð Þþ x, <sup>τ</sup> faSt <sup>∂</sup>

<sup>1</sup> , θ<sup>k</sup>

St <sup>∂</sup>θ<sup>1</sup> 0 ∂τ

faSt <sup>∂</sup>θ<sup>2</sup> 0 ∂τ

St <sup>∂</sup> ∂τ θ1

**56**

<sup>f</sup> aSt <sup>∂</sup> ∂τ θ2

If the functions θ<sup>k</sup>

<sup>0</sup>ð Þþ x, <sup>τ</sup> <sup>X</sup><sup>∞</sup>

<sup>0</sup>ð Þþ x, <sup>τ</sup> <sup>X</sup><sup>∞</sup>

<sup>1</sup> , … , θ<sup>k</sup>

Let substitute expansions (36) and (37) into Eqs. (27) and (29) and obtain

<sup>k</sup>ð Þ¼ x, <sup>τ</sup> Δθ<sup>1</sup>

St<sup>k</sup>θ<sup>2</sup>

therefore, correspondent series can be differentiated term by term. Thus

St<sup>k</sup>þ<sup>1</sup> <sup>∂</sup>θ<sup>1</sup>

St<sup>k</sup>þ<sup>1</sup> <sup>∂</sup>θ<sup>2</sup>

Δθ<sup>1</sup>

Δθ<sup>1</sup>

Δθ<sup>1</sup>

Δθ<sup>2</sup>

Δθ<sup>1</sup> <sup>i</sup> <sup>¼</sup> <sup>f</sup> <sup>a</sup>

> θ1 0� �

θ1 0� �

> θ1 i � �

θ1 i � �

θ2 0� �

Let prescribe boundary conditions for Eqs. (40)–(45)

Δθ<sup>2</sup>

<sup>1</sup> <sup>¼</sup> <sup>f</sup> <sup>a</sup>

k <sup>∂</sup><sup>τ</sup> <sup>¼</sup> Δθ<sup>1</sup>

> k <sup>∂</sup><sup>τ</sup> <sup>¼</sup> Δθ<sup>2</sup>

Since the value of Stefan number is, generally speaking, arbitrary, equality (39) takes place only in case of equality of factors at equal orders of Stefan number. Hence

> <sup>1</sup> <sup>¼</sup> <sup>∂</sup>θ<sup>1</sup> 0

> > i�1

∂θ<sup>2</sup> 0

∂θ<sup>1</sup> i�1

<sup>i</sup> <sup>¼</sup> <sup>∂</sup>θ<sup>1</sup>

k¼1

k¼1

As a result of using of expansions (36) and (37), the initial problem is reduced to

<sup>1</sup> , … ; θ<sup>0</sup>

<sup>k</sup>ð Þ¼ x, <sup>τ</sup> Δθ<sup>2</sup>

Stkθ<sup>1</sup>

Stkθ<sup>2</sup>

<sup>2</sup> , … , θ<sup>k</sup>

<sup>0</sup>ð Þþ x, <sup>τ</sup> <sup>Δ</sup>

<sup>2</sup> , … .

X∞ k¼1

<sup>0</sup>ð Þþ x, <sup>τ</sup> <sup>Δ</sup>

St<sup>k</sup>Δθ<sup>1</sup>

St<sup>k</sup>Δθ<sup>2</sup> k*:*

<sup>0</sup> <sup>¼</sup> 0, (40)

<sup>∂</sup><sup>τ</sup> , (41)

<sup>∂</sup><sup>τ</sup> , (42)

<sup>∂</sup><sup>τ</sup> , (44)

<sup>∂</sup><sup>τ</sup> , (45)

<sup>Г</sup><sup>1</sup> ¼ θ1s, (46)

<sup>Г</sup>p*:*t*:* ¼ 0, (47)

<sup>Г</sup><sup>1</sup> ¼ 0, (48)

<sup>Г</sup>p*:*t*:* ¼ 0, (49)

<sup>Г</sup><sup>2</sup> ¼ θ2s, (50)

<sup>0</sup> <sup>¼</sup> 0, (43)

<sup>2</sup> are bounded, the series (36) and (37) absolutely converge;

<sup>0</sup> <sup>þ</sup>X<sup>∞</sup> k¼1

> <sup>0</sup> <sup>þ</sup>X<sup>∞</sup> k¼1

Stkθ<sup>1</sup>

X∞ k¼1

St<sup>k</sup>θ<sup>2</sup>

k, (39)

<sup>k</sup>ð Þ x, <sup>τ</sup> , (38)

<sup>k</sup>ð Þ x, <sup>τ</sup> *:*

<sup>k</sup>ð Þ x, <sup>τ</sup> , (36)

<sup>k</sup>ð Þ x, <sup>τ</sup> *:* (37)

creation is excluded from the consideration, because thermal fluxes can be enough large (asymptotically infinite) in this case. Under mentioned restrictions the following theorem takes place.

Theorem 1. If the boundary and initial conditions in the boundary-value problems (27), (17), (19), (21) and (29), (18), (20), (22) are such that their solutions are differentiable infinite number of times and these derivatives are restricted, series (36) and (37) converge under any Stefan number less than 1 (St<1), and order of the remainder term is O Stjþ<sup>1</sup> Mjþ<sup>1</sup> , where Mjþ<sup>1</sup> <sup>¼</sup> max <sup>j</sup>þ1<sup>≤</sup> <sup>m</sup> <sup>&</sup>lt; <sup>∞</sup>θ<sup>m</sup> (it is assumed that the series is cut off after j-th term).

The proof of the Theorem 1 is evident and based on majorizing sequence, built by replacing of functions θ<sup>i</sup> by Mi. This sequence is geometrical progression with factor St, which is less than 1, according to the theorem conditions. Therefore, the series converge.

It is necessary to note, that the Theorem 1 proves convergence of asymptotic expansions and thus it grounds an application of asymptotic approach to the considered class of Stefan problems. More than that, an error of such approach can be estimated analytically, however such estimations, as a rule, are found useless for determination of computation error in whole, because there are other error sources in the general computational scheme.

Let consider one-phase Stefan problem, which is a particular case of above formulated problem:

$$
\Delta \boldsymbol{\theta}^0 = \mathbf{0},
\tag{57}
$$

C xð Þ<sup>0</sup> <sup>θ</sup><sup>0</sup>

*DOI: http://dx.doi.org/10.5772/intechopen.93788*

C xð Þ<sup>0</sup> <sup>θ</sup><sup>0</sup>

ð

Γ<sup>1</sup> ∪ Γp*:*t*:*

ð

θi 1

θi 2

Γ<sup>1</sup> ∪ Γp*:*t*:*

ð

Γ<sup>2</sup> ∪ Γp*:*t*:*

Γ<sup>2</sup> ∪ Γp*:*t*:*

however their applications require special investigation.

�

� ð

C xð Þ<sup>0</sup> <sup>θ</sup><sup>i</sup>

C xð Þ<sup>0</sup> <sup>θ</sup><sup>i</sup>

elements is formed.

**59**

<sup>1</sup>ð Þ¼ x0

<sup>2</sup>ð Þ¼ x0

<sup>1</sup> ð Þ¼ x0

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow…*

<sup>2</sup> ð Þ¼ x0

ϕ0ð Þ x, x0

ϕ0ð Þ x, x0

ð

ϕ0ð Þ x, x0

ϕ0ð Þ x, x0

θ0 1

θ0 2

ds þ ð

ds þ ð

where ϕ<sup>0</sup> is fundamental solution of Laplace equation, C xð Þ<sup>0</sup> is special function reflecting control point x0 position with respect to the boundary and shape of the boundary. Eqs. (66)–(69) are solved numerically by well-known boundary element

Solving mentioned systems of linear algebraic equations, corresponding to every boundary integral Eqs. (66)–(69) we can obtain the temperature distribution with required accuracy at some instant of time, that is under specific shape of interphase boundary. Then new position of interphase boundary must be determined using the Stefan condition (55) and (56). To build new interphase boundary relations (55) and (56) must be considered as an ordinary differential equation, and correspondent Cauchy problems must be solved numerically. The Euler method is used for this aim in the present work. If it is necessary, the problems (66)–(69) can be solved by boundary element method for new shape of the interphase boundary. Such timestepping process can be continued during any time interval. The simplest algorithm of the boundary element method with approximation of boundary elements by straight lines segments and approximation of known and unknown function values on boundary elements by constants [6, 7] is used in the present work. Of course, using of improved boundary element method algorithm and methods of numerical solution of Cauchy problem can make the solution procedure more effective;

D1

D2

ϕ0ð Þ x, x0

ϕ0ð Þ x, x0

∂θ<sup>0</sup> 1 ∂n ds�

∂θ<sup>0</sup> 2 ∂n ds�

ds,

ds,

∂θ<sup>i</sup>�<sup>1</sup> 1 ∂τ

∂θ<sup>i</sup>�<sup>1</sup> 2 ∂τ

dx,

dx,

(66)

(67)

(68)

(69)

∂ϕ0ð Þ x, x0 ∂n

∂ϕ0ð Þ x, x0 ∂n

Γ<sup>1</sup> ∪ Γp*:*t*:*

Γ<sup>1</sup> ∪ Γp*:*t*:*

ð

Γ<sup>2</sup> ∪ Γp*:*t*:*

Γ<sup>2</sup> ∪ Γp*:*t*:*

∂θi 1 ∂n ds�

� ð

� ð

∂ϕ0ð Þ x, x0 ∂n

∂ϕ0ð Þ x, x0 ∂n

algorithm [6, 7]. According to that algorithm, the boundaries of phases are fragmented by boundary elements; the temperatures and thermal fluxes are assumed constant on every boundary element. Thus, the system of linear algebraic equations with respect to unknown values of temperature or thermal flux on

∂θi 2 ∂n ds�

$$
\Delta\theta^1 = \frac{\partial\theta^0}{\partial\pi},
\tag{58}
$$

$$
\Delta\theta^i = \frac{\partial\theta^{i-1}}{\partial\pi} \tag{59}
$$

with boundary conditions

$$\left.\Theta^{0}\right|\_{\Gamma\_{1}} = \Theta\_{\text{s}},\tag{60}$$

$$\left.\Theta^{0}\right|\_{\Gamma\_{\rm pt}} = \mathbf{0},\tag{61}$$

$$\left.\Theta^{i}\right|\_{\Gamma\_{1}} = \mathbf{0},\tag{62}$$

$$\left.\Theta^{i}\right|\_{\Gamma\_{\rm pt.}} = \mathbf{0},\tag{63}$$

Stefan condition

$$\frac{\partial \theta^0}{\partial \mathbf{n}} = \frac{\partial \eta^0}{\partial \mathbf{\tau}}\tag{64}$$

$$\frac{\partial \theta^{\rm i}}{\partial \mathbf{n}} = \frac{\partial \eta^{\rm i}}{\partial \boldsymbol{\pi}} \tag{65}$$

### **3. Boundary element method application to slow phase transition calculations**

However, the temperature fields must be determined numerically in twodimensional and three-dimensional cases. Boundary element method [6, 7] is used in the present work. It supposes transformation of Eqs. (40)–(45) into boundary integral equations

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow… DOI: http://dx.doi.org/10.5772/intechopen.93788*

C xð Þ<sup>0</sup> <sup>θ</sup><sup>0</sup> <sup>1</sup> ð Þ¼ x0 ð Γ<sup>1</sup> ∪ Γp*:*t*:* ϕ0ð Þ x, x0 ∂θ<sup>0</sup> 1 ∂n ds� � ð Γ<sup>1</sup> ∪ Γp*:*t*:* θ0 1 ∂ϕ0ð Þ x, x0 ∂n ds, (66) C xð Þ<sup>0</sup> <sup>θ</sup><sup>0</sup> <sup>2</sup> ð Þ¼ x0 ð Γ<sup>2</sup> ∪ Γp*:*t*:* ϕ0ð Þ x, x0 ∂θ<sup>0</sup> 2 ∂n ds� � ð Γ<sup>2</sup> ∪ Γp*:*t*:* θ0 2 ∂ϕ0ð Þ x, x0 ∂n ds, (67) C xð Þ<sup>0</sup> <sup>θ</sup><sup>i</sup> <sup>1</sup>ð Þ¼ x0 ð Γ<sup>1</sup> ∪ Γp*:*t*:* ϕ0ð Þ x, x0 ∂θi 1 ∂n ds� � ð Γ<sup>1</sup> ∪ Γp*:*t*:* θi 1 ∂ϕ0ð Þ x, x0 ∂n ds þ ð D1 ϕ0ð Þ x, x0 ∂θ<sup>i</sup>�<sup>1</sup> 1 ∂τ dx, (68) C xð Þ<sup>0</sup> <sup>θ</sup><sup>i</sup> <sup>2</sup>ð Þ¼ x0 ð Γ<sup>2</sup> ∪ Γp*:*t*:* ϕ0ð Þ x, x0 ∂θi 2 ∂n ds� � ð Γ<sup>2</sup> ∪ Γp*:*t*:* θi 2 ∂ϕ0ð Þ x, x0 ∂n ds þ ð D2 ϕ0ð Þ x, x0 ∂θ<sup>i</sup>�<sup>1</sup> 2 ∂τ dx, (69)

where ϕ<sup>0</sup> is fundamental solution of Laplace equation, C xð Þ<sup>0</sup> is special function reflecting control point x0 position with respect to the boundary and shape of the boundary. Eqs. (66)–(69) are solved numerically by well-known boundary element algorithm [6, 7]. According to that algorithm, the boundaries of phases are fragmented by boundary elements; the temperatures and thermal fluxes are assumed constant on every boundary element. Thus, the system of linear algebraic equations with respect to unknown values of temperature or thermal flux on elements is formed.

Solving mentioned systems of linear algebraic equations, corresponding to every boundary integral Eqs. (66)–(69) we can obtain the temperature distribution with required accuracy at some instant of time, that is under specific shape of interphase boundary. Then new position of interphase boundary must be determined using the Stefan condition (55) and (56). To build new interphase boundary relations (55) and (56) must be considered as an ordinary differential equation, and correspondent Cauchy problems must be solved numerically. The Euler method is used for this aim in the present work. If it is necessary, the problems (66)–(69) can be solved by boundary element method for new shape of the interphase boundary. Such timestepping process can be continued during any time interval. The simplest algorithm of the boundary element method with approximation of boundary elements by straight lines segments and approximation of known and unknown function values on boundary elements by constants [6, 7] is used in the present work. Of course, using of improved boundary element method algorithm and methods of numerical solution of Cauchy problem can make the solution procedure more effective; however their applications require special investigation.

creation is excluded from the consideration, because thermal fluxes can be enough large (asymptotically infinite) in this case. Under mentioned restrictions the

Theorem 1. If the boundary and initial conditions in the boundary-value problems (27), (17), (19), (21) and (29), (18), (20), (22) are such that their solutions are differentiable infinite number of times and these derivatives are restricted, series (36) and (37) converge under any Stefan number less than 1 (St<1), and order of

The proof of the Theorem 1 is evident and based on majorizing sequence, built by replacing of functions θ<sup>i</sup> by Mi. This sequence is geometrical progression with factor St, which is less than 1, according to the theorem conditions. Therefore, the series converge. It is necessary to note, that the Theorem 1 proves convergence of asymptotic expansions and thus it grounds an application of asymptotic approach to the considered class of Stefan problems. More than that, an error of such approach can be estimated analytically, however such estimations, as a rule, are found useless for determination of computation error in whole, because there are other error sources

Let consider one-phase Stefan problem, which is a particular case of above

Δθ<sup>1</sup> <sup>¼</sup> <sup>∂</sup> <sup>θ</sup><sup>0</sup>

Δ θ<sup>i</sup> <sup>¼</sup> <sup>∂</sup> <sup>θ</sup><sup>i</sup>�<sup>1</sup> ∂ τ

θ0 

θ0 

> θi

θi 

> ∂ θ<sup>0</sup> <sup>∂</sup><sup>n</sup> <sup>¼</sup> <sup>∂</sup>η<sup>0</sup> ∂ τ

∂θ<sup>i</sup> <sup>∂</sup><sup>n</sup> <sup>¼</sup> <sup>∂</sup>η<sup>i</sup> ∂ τ

**3. Boundary element method application to slow phase transition**

However, the temperature fields must be determined numerically in twodimensional and three-dimensional cases. Boundary element method [6, 7] is used in the present work. It supposes transformation of Eqs. (40)–(45) into boundary

, where Mjþ<sup>1</sup> <sup>¼</sup> max <sup>j</sup>þ1<sup>≤</sup> <sup>m</sup> <sup>&</sup>lt; <sup>∞</sup>θ<sup>m</sup> (it is assumed

Δθ<sup>0</sup> <sup>¼</sup> 0, (57)

<sup>∂</sup> <sup>τ</sup> , (58)

<sup>Γ</sup><sup>1</sup> ¼ θs, (60)

<sup>Γ</sup>p*:*t*:* ¼ 0, (61)

<sup>Γ</sup><sup>1</sup> ¼ 0, (62)

<sup>Γ</sup>p*:*t*:* ¼ 0, (63)

(59)

(64)

(65)

Mjþ<sup>1</sup>

*Modeling and Simulation in Engineering - Selected Problems*

following theorem takes place.

the remainder term is O Stjþ<sup>1</sup>

that the series is cut off after j-th term).

in the general computational scheme.

with boundary conditions

Stefan condition

**calculations**

integral equations

**58**

formulated problem:

The main advantage of proposed time-stepping approach in comparison with time-stepping finite difference and finite element methods is time discretization with respect to slow time scale determined by interphase boundary motion unlike finite difference and finite element methods, which requires time discretization with respect to fast time scale determined by heat conduction process.

### **4. Results of slow phase transition calculations**

The restricted amount of the present paper doesn't give an opportunity to analyze a lot of calculation results, therefore the proposed approach is illustrated by only several following examples of numerical calculations. Melting process of cylinder, which had initial circular shape, is shown in **Figures 3** and **4** under different conditions of heating.

It is easy to see from comparison of **Figures 3** and **4** that the way of heating sufficiently affects shape of phase transition boundary, in particular, side heating leads to oval shape of the second phase domain, but nearly circular second phase domain takes place under multilateral heating. The times of calculations are relatively small in comparison with known finite difference and finite element calculations.

test calculations. The results of calculations can be used in investigations of a number of processes, first of all, in investigations of slow phase transformations in microgravity and in environment. They also can be used in design of space-rocket

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow…*

*DOI: http://dx.doi.org/10.5772/intechopen.93788*

Let consider D1 filled by some biological structures (in the simplest case by homogeneous or nondifferentiated cellular mass). Let restrict the following consideration by the case of homogeneous cellular structures. The tissue in the domain D1 is porous media where cells form a frame and intercellular space is porosity. Let assume that pores are filled by same liquid, which is complex solution of nutrient substances and excrements of cells. There is an intensive heat and mass transfer between the frame and the liquid in pores, what is very important specific feature of the described structure. Let the domain D1 is partially or completely surrounded by the domain D2, filled by the same solution completely. In general case there may be a convective transfer in the domain D2 and filtration flow in the domain D1. Thus, a general mathematical model of heat and mass transfer processes is considered

þ Vf ð Þ� � ∇ T1 ¼ a1ΔT1 þ qT1, (70)

þ ð Þ VC � ∇ T2 ¼ a2ΔT2, (72)

þ ð Þ VC � ∇ Ci2 ¼ di2ΔCi2, i ¼ 1, N, (73)

where T1 is temperature in the domain D1 (the one-temperature model, assuming the temperatures of frame and solution in pores are equal, is used here), Vf is filtration velocity, a1 is thermal diffusivity of porous media, qT1 is heat source, concerning the metabolism of cells, ci1 is concentration of the i-th component in

þ Vf ð Þ� � ∇ Ci1 ¼ di1ΔCi1 þ qi1, i ¼ 1, N, (71)

**5. Mathematical model of biological tissue growth**

*Sequence of phase transformation boundary positions under multilateral heating.*

∂T1 ∂τ

∂Ci2 ∂τ

∂T2 ∂τ

∂Ci1 ∂τ

technique.

**Figure 4.**

system is following:

**61**

Since there is not any known analytical solution of Stefan problem to check accuracy of the obtained numerical solutions, the authors of the present work had to change numbers of boundary elements and compare correspondent numerical results. For example, the presented above numerical calculations were made under following parameters: phase transition boundary is approximated by 40 boundary elements at every time step of calculation and outer boundaries of the domain were approximated by 4 � 40 boundary elements. Checking calculations with 80 and 4 � 80 and 160 and 4 � 160 boundary elements and dimensionless time step varying from Δτ ¼ 0, 1 to Δτ ¼ 0, 001, correspondingly, shows change of the obtained results no more than 1%. The authors made a conclusion about satisfactory accuracy of the considered numerical approach.

Some conclusions can be made here. It is evident that the above proposed approach is more general and has more opportunities than earlier existing computational methods for slow phase transitions. It gives an opportunity to effectively calculate of the processes in broader interval of Stefan numbers. Enough high accuracy of the proposed algorithms is confirmed by the above-mentioned series of

**Figure 3.** *Sequence of phase transformation boundary positions under side heating.*

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow… DOI: http://dx.doi.org/10.5772/intechopen.93788*

**Figure 4.** *Sequence of phase transformation boundary positions under multilateral heating.*

test calculations. The results of calculations can be used in investigations of a number of processes, first of all, in investigations of slow phase transformations in microgravity and in environment. They also can be used in design of space-rocket technique.

### **5. Mathematical model of biological tissue growth**

Let consider D1 filled by some biological structures (in the simplest case by homogeneous or nondifferentiated cellular mass). Let restrict the following consideration by the case of homogeneous cellular structures. The tissue in the domain D1 is porous media where cells form a frame and intercellular space is porosity. Let assume that pores are filled by same liquid, which is complex solution of nutrient substances and excrements of cells. There is an intensive heat and mass transfer between the frame and the liquid in pores, what is very important specific feature of the described structure. Let the domain D1 is partially or completely surrounded by the domain D2, filled by the same solution completely. In general case there may be a convective transfer in the domain D2 and filtration flow in the domain D1. Thus, a general mathematical model of heat and mass transfer processes is considered system is following:

$$\frac{\partial \mathbf{T}\_1}{\partial \mathbf{r}} + (\mathbf{V}\_\mathbf{f} \cdot \nabla) \cdot \mathbf{T}\_1 = \mathbf{a}\_1 \Delta \mathbf{T}\_1 + \mathbf{q}\_{\mathbf{T}1},\tag{70}$$

$$\frac{\partial \mathbf{C}\_{\mathrm{i}1}}{\partial \mathbf{r}} + (\mathbf{V}\_{\mathrm{f}} \cdot \nabla) \cdot \mathbf{C}\_{\mathrm{i}1} = \mathbf{d}\_{\mathrm{i}1} \Delta \mathbf{C}\_{\mathrm{i}1} + \mathbf{q}\_{\mathrm{i}1}, \mathrm{i} = \overline{\mathbf{1}, \mathbf{N}}, \tag{71}$$

$$\frac{\partial \mathbf{T}\_2}{\partial \mathbf{r}} + (\mathbf{V}\_\mathbf{C} \cdot \nabla)\mathbf{T}\_2 = \mathbf{a}\_2 \Delta \mathbf{T}\_2,\tag{72}$$

$$\frac{\partial \mathbf{C}\_{i2}}{\partial \mathbf{r}} + (\mathbf{V}\_{\mathbf{C}} \cdot \nabla) \mathbf{C}\_{i2} = \mathbf{d}\_{i2} \Delta \mathbf{C}\_{i2}, \mathbf{i} = \overline{\mathbf{1}, \mathbf{N}}, \tag{73}$$

where T1 is temperature in the domain D1 (the one-temperature model, assuming the temperatures of frame and solution in pores are equal, is used here), Vf is filtration velocity, a1 is thermal diffusivity of porous media, qT1 is heat source, concerning the metabolism of cells, ci1 is concentration of the i-th component in

The main advantage of proposed time-stepping approach in comparison with time-stepping finite difference and finite element methods is time discretization with respect to slow time scale determined by interphase boundary motion unlike finite difference and finite element methods, which requires time discretization

The restricted amount of the present paper doesn't give an opportunity to analyze a lot of calculation results, therefore the proposed approach is illustrated by only several following examples of numerical calculations. Melting process of cylinder, which had initial circular shape, is shown in **Figures 3** and **4** under

It is easy to see from comparison of **Figures 3** and **4** that the way of heating sufficiently affects shape of phase transition boundary, in particular, side heating leads to oval shape of the second phase domain, but nearly circular second phase domain takes place under multilateral heating. The times of calculations are relatively small in comparison with known finite difference and finite element calculations. Since there is not any known analytical solution of Stefan problem to check accuracy of the obtained numerical solutions, the authors of the present work had to change numbers of boundary elements and compare correspondent numerical results. For example, the presented above numerical calculations were made under following parameters: phase transition boundary is approximated by 40 boundary elements at every time step of calculation and outer boundaries of the domain were approximated by 4 � 40 boundary elements. Checking calculations with 80 and 4 � 80 and 160 and 4 � 160 boundary elements and dimensionless time step varying from Δτ ¼ 0, 1 to Δτ ¼ 0, 001, correspondingly, shows change of the obtained results no more than 1%. The authors made a conclusion about satisfactory

Some conclusions can be made here. It is evident that the above proposed approach is more general and has more opportunities than earlier existing computational methods for slow phase transitions. It gives an opportunity to effectively calculate of the processes in broader interval of Stefan numbers. Enough high accuracy of the proposed algorithms is confirmed by the above-mentioned series of

with respect to fast time scale determined by heat conduction process.

**4. Results of slow phase transition calculations**

*Modeling and Simulation in Engineering - Selected Problems*

accuracy of the considered numerical approach.

*Sequence of phase transformation boundary positions under side heating.*

**Figure 3.**

**60**

different conditions of heating.

porous media, di1 is diffusion coefficient of i-th component in the porous media, qi1 is source (sink) of the i-th component in porous media, concerning the metabolism of cells, T2 is temperature in the domain D2, V2 is flow velocity in the domain D2, a2 is thermal diffusivity of solution, ci2 is concentration of i-th component in the domain D2, di2 is diffusion coefficient of i-th component in the domain D2, Ni2number of components, participating of heat and mass transfer process, τ is time, Δ is Laplace operator.

Restrict the following consideration by the case:

$$\mathbf{V}\_{\mathbf{f}} = \mathbf{0},\tag{74}$$

or the third boundary condition

*DOI: http://dx.doi.org/10.5772/intechopen.93788*

λ1 ∂T1 ∂n Г1

di1 ∂Ci1 ∂n Г1

> λ2 ∂T2 ∂n Г2

di2 ∂Ci2 ∂n Г2

boundary conditions on the boundary Γ. It is evident that

þ α<sup>1</sup> T1j

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow…*

þ αi1 Ci1j

þ α<sup>2</sup> T2j

þ αi2 Ci2j

T1j

Сi1j

λ1 ∂T1 ∂n Г ¼ λ<sup>2</sup> ∂T2 ∂n Г

di1 ∂Ci1 ∂n Г ¼ di2

di1 ∂Ci1 ∂n Г � di2

considered problem as a function of metabolism intensity.

correspondent source terms are following

by following condition

here <sup>∂</sup><sup>n</sup>

**63**

where T1e, Ci1e, T2e, Ci2e, q1e, qi1e, q2e, qi2e are known functions, all coefficients in boundary conditions (80)–(91) are understood in conventional sense. Let consider

<sup>Г</sup> ¼ T2j

<sup>Г</sup> ¼ Ci2j

It is possible to formulate the second condition as a forth kind boundary condition.

Conditions (94) and (95) correspond to the case of cell fission in whole domain D1. However, it is possible the situation, when the fission of cells takes place only on the boundary Γ, then condition (94) is saved, but condition (95) must be replaced

> ∂Ci2 ∂n Г ¼ χ<sup>i</sup> ∂n

growth), χ<sup>i</sup> is "expenditure" coefficient of the i-th component during growth of biological structure. Note that condition (96) is not conventional Stefan condition (nevertheless its form coincides with Stefan condition), because right hand part of condition (96) is determined by fission process, that is by parameters determining the fission process such as the temperature, concentrations and possibly the histories, therefore right hand part of the condition (96) is prescribed. It means that the given problem is similar to phase transition problem under prescribed velocity of phase boundary motion. The moving boundary velocity is determined in the

The case, when cellular mass growth takes place in whole domain D1 is more complex than previous one. Consider a function describing metabolism intensity. As it is noted earlier, metabolism intensity is assumed proportional to a cellular mass growth (nevertheless the cell fission is very complex process with possibly enough large delay time that is with sufficient influence of previous history of the process). Let metabolism intensity function ω T1, ci1 ð Þ is defined, then

<sup>∂</sup><sup>τ</sup> is velocity of the boundary Γ propagation (velocity of biological structure

∂Ci2 ∂n Г

<sup>Г</sup><sup>1</sup> � T1e <sup>¼</sup> 0, (88)

<sup>Г</sup><sup>1</sup> � <sup>С</sup>i1e <sup>¼</sup> 0, (89)

<sup>Г</sup><sup>2</sup> � T2e <sup>¼</sup> 0, (90)

<sup>Г</sup><sup>2</sup> � <sup>С</sup>i2e <sup>¼</sup> 0, (91)

<sup>Г</sup>, (92)

<sup>Г</sup>*:* (93)

, (94)

*:* (95)

<sup>∂</sup><sup>τ</sup> , (96)

$$\mathbf{V\_c = 0,}\tag{75}$$

what corresponds to conventional multicellular structure, formed by independent cells, that is simple colony of one-cellular organisms in immovable fluid. Then:

$$\frac{\partial \mathbf{T}\_1}{\partial \mathbf{r}} = \mathbf{a}\_1 \Delta \mathbf{T}\_1 + \mathbf{q}\_{\mathbf{T}1},\tag{76}$$

$$\frac{\partial \mathbf{C}\_{\mathrm{i}1}}{\partial \mathbf{r}} = \mathbf{d}\_{\mathrm{i}1} \Delta \mathbf{C}\_{\mathrm{i}1} + \mathbf{q}\_{\mathrm{i}1}, \mathrm{i} = \overline{\mathbf{1}, \mathbf{N}}, \tag{77}$$

$$\frac{\partial \mathbf{T}\_2}{\partial \mathbf{r}} = \mathbf{a}\_2 \Delta \mathbf{T}\_2,\tag{78}$$

$$\frac{\partial \mathbf{C}\_{i2}}{\partial \mathbf{\bar{r}}} = \mathbf{d}\_{i2} \Delta \mathbf{C}\_{i2}, \mathbf{i} = \overline{\mathbf{1}, \mathbf{N}}.\tag{79}$$

If the condition (75) is not realized, it could be better to not consider the system (72) and (73), but to take into account a convective transfer using boundary conditions for Eqs. (76) and (77). This assumption is quite proved, since the system (76) and (77) describes enough slow processes.

Let prescribe boundary conditions for the systems (76)–(79). Note the common boundary of the domain D1 and D2 as Γ and reminder part as Γ<sup>1</sup> and Γ<sup>2</sup> correspondingly. The first kind boundary conditions can be prescribed on the boundaries Γ<sup>1</sup> and Γ<sup>2</sup>

$$\mathbf{T\_1}|\_{\Gamma1} = \mathbf{T\_{1e}},\tag{80}$$

$$\mathbf{C}\_{\mathbf{i}1}|\_{\Gamma1} = \mathbf{C}\_{\mathbf{i}1\mathbf{e}},\tag{81}$$

$$\mathbf{T\_2}|\_{\mathbf{r\_2}} = \mathbf{T\_{2e}},\tag{82}$$

$$\mathbf{C}\_{\mathbf{i}2}|\_{\Gamma2} = \mathbf{C}\_{\mathbf{i}2\mathbf{e}},\tag{83}$$

or the second kind boundary condition

$$
\lambda\_1 \frac{\partial \mathbf{T}\_1}{\partial \mathbf{n}}\Big|\_{\Gamma1} = \mathbf{q}\_{1\mathbf{e}},\tag{84}
$$

$$\mathbf{d}\_{\rm 1i} \frac{\partial \mathbf{C}\_{\rm i1}}{\partial \mathbf{n}} \bigg|\_{\Gamma 1} = \mathbf{q}\_{\rm i1e},\tag{85}$$

$$
\lambda\_2 \frac{\partial \mathbf{T}\_2}{\partial \mathbf{n}} \bigg|\_{\mathbf{r}2} = \mathbf{q}\_{2\mathbf{e}}, \tag{86}
$$

$$\left. \mathbf{d}\_{i2} \frac{\partial \mathbf{C}\_{i2}}{\partial \mathbf{n}} \right|\_{\Gamma2} = \mathbf{q}\_{i2\mathbf{e}}.\tag{87}$$

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow… DOI: http://dx.doi.org/10.5772/intechopen.93788*

or the third boundary condition

porous media, di1 is diffusion coefficient of i-th component in the porous media, qi1 is source (sink) of the i-th component in porous media, concerning the

metabolism of cells, T2 is temperature in the domain D2, V2 is flow velocity in the domain D2, a2 is thermal diffusivity of solution, ci2 is concentration of i-th component in the domain D2, di2 is diffusion coefficient of i-th component in the domain D2, Ni2number of components, participating of heat and mass transfer process, τ is

what corresponds to conventional multicellular structure, formed by independent cells, that is simple colony of one-cellular organisms in immovable

∂T1

∂T2

Let prescribe boundary conditions for the systems (76)–(79). Note the common boundary of the domain D1 and D2 as Γ and reminder part as Γ<sup>1</sup> and Γ<sup>2</sup> correspondingly. The first kind boundary conditions can be prescribed on the

T1j

Ci1j

T2j

Ci2j

λ1 ∂T1 ∂n Г1

d1i ∂Ci1 ∂n Г1

λ2 ∂T2 ∂n Г2

di2 ∂Ci2 ∂n Г2

If the condition (75) is not realized, it could be better to not consider the system (72) and (73), but to take into account a convective transfer using boundary conditions for Eqs. (76) and (77). This assumption is quite proved, since the system (76)

∂Ci1

and (77) describes enough slow processes.

or the second kind boundary condition

boundaries Γ<sup>1</sup> and Γ<sup>2</sup>

**62**

∂Ci2

Vf ¼ 0, (74) Vc ¼ 0, (75)

<sup>∂</sup><sup>τ</sup> <sup>¼</sup> a1ΔT1 <sup>þ</sup> qT1, (76)

<sup>∂</sup><sup>τ</sup> <sup>¼</sup> a2ΔT2, (78)

<sup>Г</sup><sup>1</sup> ¼ T1e, (80)

<sup>Г</sup><sup>1</sup> ¼ Сi1e, (81)

<sup>Г</sup><sup>2</sup> ¼ T2e, (82)

<sup>Г</sup><sup>2</sup> ¼ Ci2e, (83)

¼ q1e, (84)

¼ qi1e, (85)

¼ q2e, (86)

¼ qi2e, (87)

<sup>∂</sup><sup>τ</sup> <sup>¼</sup> di2ΔCi2, i <sup>¼</sup> 1, N*:* (79)

<sup>∂</sup><sup>τ</sup> <sup>¼</sup> di1ΔCi1 <sup>þ</sup> qi1, i <sup>¼</sup> 1, N, (77)

time, Δ is Laplace operator.

fluid. Then:

Restrict the following consideration by the case:

*Modeling and Simulation in Engineering - Selected Problems*

$$
\lambda\_1 \frac{\partial \mathbf{T}\_1}{\partial \mathbf{n}} \bigg|\_{\Gamma 1} + \mathfrak{a}\_1 \Big(\mathbf{T}\_1 |\_{\Gamma 1} - \mathbf{T}\_{1\mathbf{e}}\big) = \mathbf{0},\tag{88}
$$

$$\mathbf{d}\_{\rm i1} \frac{\partial \mathbf{C}\_{\rm i1}}{\partial \mathbf{n}} \bigg|\_{\Gamma 1} + \mathbf{a}\_{\rm i1} (\mathbf{C}\_{\rm i1}|\_{\Gamma 1} - \mathbf{C}\_{\rm i1e}) = \mathbf{0},\tag{89}$$

$$
\lambda\_2 \frac{\partial \mathbf{T}\_2}{\partial \mathbf{n}} \bigg|\_{\Gamma2} + \mathfrak{a}\_2 \left( \mathbf{T}\_2 |\_{\Gamma2} - \mathbf{T}\_{2\mathfrak{e}} \right) = \mathbf{0},\tag{90}
$$

$$\mathbf{d}\_{\rm i2} \frac{\partial \mathbf{C}\_{\rm i2}}{\partial \mathbf{n}} \bigg|\_{\rm I2} + \mathbf{a}\_{\rm i2} \left( \mathbf{C}\_{\rm i2} |\_{\rm I2} - \mathbf{C}\_{\rm i2e} \right) = \mathbf{0},\tag{91}$$

where T1e, Ci1e, T2e, Ci2e, q1e, qi1e, q2e, qi2e are known functions, all coefficients in boundary conditions (80)–(91) are understood in conventional sense. Let consider boundary conditions on the boundary Γ. It is evident that

$$\mathbf{T\_1}|\_{\Gamma} = \mathbf{T\_2}|\_{\Gamma},\tag{92}$$

$$\mathbf{C}\_{\mathbf{i}1}|\_{\Gamma} = \mathbf{C}\_{\mathbf{i}2}|\_{\Gamma}.\tag{93}$$

It is possible to formulate the second condition as a forth kind boundary condition.

$$
\lambda\_1 \frac{\partial \mathbf{T}\_1}{\partial \mathbf{n}}\Big|\_{\Gamma} = \lambda\_2 \frac{\partial \mathbf{T}\_2}{\partial \mathbf{n}}\Big|\_{\Gamma},\tag{94}
$$

$$\left. \mathbf{d}\_{\mathrm{i}1} \frac{\partial \mathbf{C}\_{\mathrm{i}1}}{\partial \mathbf{n}} \right|\_{\Gamma} = \left. \mathbf{d}\_{\mathrm{i}2} \frac{\partial \mathbf{C}\_{\mathrm{i}2}}{\partial \mathbf{n}} \right|\_{\Gamma}. \tag{95}$$

Conditions (94) and (95) correspond to the case of cell fission in whole domain D1. However, it is possible the situation, when the fission of cells takes place only on the boundary Γ, then condition (94) is saved, but condition (95) must be replaced by following condition

$$\left. \mathbf{d}\_{\mathrm{i}1} \frac{\partial \mathbf{C}\_{\mathrm{i}1}}{\partial \mathbf{n}} \right|\_{\Gamma} - \left. \mathbf{d}\_{\mathrm{i}2} \frac{\partial \mathbf{C}\_{\mathrm{i}2}}{\partial \mathbf{n}} \right|\_{\Gamma} = \chi\_{\mathrm{i}} \frac{\partial \mathbf{n}}{\partial \pi}, \tag{96}$$

here <sup>∂</sup><sup>n</sup> <sup>∂</sup><sup>τ</sup> is velocity of the boundary Γ propagation (velocity of biological structure growth), χ<sup>i</sup> is "expenditure" coefficient of the i-th component during growth of biological structure. Note that condition (96) is not conventional Stefan condition (nevertheless its form coincides with Stefan condition), because right hand part of condition (96) is determined by fission process, that is by parameters determining the fission process such as the temperature, concentrations and possibly the histories, therefore right hand part of the condition (96) is prescribed. It means that the given problem is similar to phase transition problem under prescribed velocity of phase boundary motion. The moving boundary velocity is determined in the considered problem as a function of metabolism intensity.

The case, when cellular mass growth takes place in whole domain D1 is more complex than previous one. Consider a function describing metabolism intensity. As it is noted earlier, metabolism intensity is assumed proportional to a cellular mass growth (nevertheless the cell fission is very complex process with possibly enough large delay time that is with sufficient influence of previous history of the process). Let metabolism intensity function ω T1, ci1 ð Þ is defined, then correspondent source terms are following

*Modeling and Simulation in Engineering - Selected Problems*

$$\mathbf{q}\_{\text{i1}} = \chi\_{\text{i}} \mathbf{o},\tag{97}$$

**6. Boundary element method application to biological tissue growth**

dimensional cases require some numerical method for solution of elliptic boundaryvalue problems in moving boundary domain. The most powerful tool for such problems is boundary element method [6, 7], which requires a reformulation of the

Let consider the initial boundary value problem (76)–(96). Small parameter method application to this problem is, generally speaking, similar to above onedimensional case application (see, for example, [8, 9]). Restrict the following consideration by plane case and by zero approximation of small parameter method, what corresponds to very small value of the Stefan number analogue. Thus

> <sup>1</sup> ¼ � qT1 a1

Boundary conditions for the system (104)–(107) coincide with boundary condi-

Γ1 T0 1

Γ1 C0 i1

Γ2 T0 2

> Γ2 C0 i2

ds � ð

dxdy,

ds � ð

dxdy,

, (104)

, i ¼ 1, N, (105)

dsþ

dsþ

(108)

(109)

ds, (110)

ds*:* (111)

<sup>2</sup> ¼ 0, (106)

i2 ¼ 0, i ¼ 1, N*:* (107)

∂ϕ0ð Þ x, x0 ∂n

∂ϕ0ð Þ x, x0 ∂n

∂ϕ0ð Þ x, x0 ∂n

∂ϕ0ð Þ x, x0 ∂n

ΔT<sup>0</sup>

i1 ¼ � qi1 di1

ΔT0

tions for the initial system. Let apply methods of potential theory to the system

∂T<sup>0</sup> 1 ∂n ds � ð

qT1 a1

∂C<sup>0</sup> i1 ∂n

> qi1 di1

∂T<sup>0</sup> 2 ∂n ds � ð

∂C<sup>0</sup> i2 ∂n

Here the function ϕ0ð Þ x, x0 is well-known fundamental solution of Laplace

ΔC0

ΔC<sup>0</sup>

The above-developed algorithm cannot be directly applied to the twodimensional and three-dimensional problems because boundary-value problems for partial differential equations arise in the mentioned cases instead boundary-value problems for ordinary differential equations as above. Thus two- and three-

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow…*

considered problems as boundary integral equations.

*DOI: http://dx.doi.org/10.5772/intechopen.93788*

(104)–(107).

<sup>χ</sup>ð Þ x0 T0

<sup>χ</sup>ð Þ x0 C0

<sup>χ</sup>ð Þ x0 <sup>T</sup><sup>0</sup>

<sup>χ</sup>ð Þ x0 C0

equation, which is in plane case

**65**

<sup>1</sup> ð Þ¼ x0

i1ð Þ¼ x0

<sup>2</sup> ð Þ¼ x0

i2ð Þ¼ x0

ð

ϕ0ð Þ x, x0

ϕ0ð Þ x, x0

ϕ0ð Þ x, x0

ϕ0ð Þ x, x0

ϕ0ð Þ x, x0

ϕ0ð Þ x, x0

Γ1

þ ð

D

ð

Γ1

þ ð

D

ð

Γ2

ð

Γ2

$$\mathbf{q}^{\text{L1}} = \mathbb{X}^{\text{L}} \mathbf{o} \,\text{.}\tag{98}$$

The function ω is determined experimentally. Let ω<sup>i</sup> is function of influence of the i-th parameter on the metabolism velocity function. It is evident that

$$\alpha = \prod\_{\mathbf{i}=1}^{N+1} \alpha\_{\mathbf{i}}.\tag{99}$$

As it is noted earlier the growth of cellular mass is proportional to metabolism intensity

$$\mathbf{d}^{\mathrm{v}} = \chi^{\mathrm{v}} \mathbf{o} - \chi^{\mathrm{0}} \mathbf{o} \mathbf{o} \,\mathrm{.}\tag{100}$$

Terms indicated by }0} in last relationship correspond to regular metabolism, which is specific for tissues of highest animals.

Let consider a problem about motion of the boundary Γ again, in particular, let consider the case, when local volume change is determined by relation (100). The velocity (deformation) field depends on mechanical links between cells. If cells are "free" in intercellular solution the model of distributed sources in incompressible fluid can be applied, according to which the velocity of the boundary Γ is determined as

$$\mathbf{V\_n(x\_\*) = \frac{\partial}{\partial \mathbf{n}} \int\_{\mathbf{D\_1}} \mathbf{q\_s(x)} \phi\_\*(\mathbf{x}, \mathbf{x\_\*}) d\mathbf{x},\tag{101}$$

where xi is arbitrary point of the curve-line Γ.

If the cells are linked mechanically between themselves, to determine the motion of the boundary Γ it is necessary to solve an elastoplastic problem, as a rule, under large strains. Consideration of such problems requires especial investigation and will not be made in the present work.

However, the another case is possible in biological structures; a biological structure grows saving its shape in this case. Thus, change of the structure volume can be referred to the boundary Γ uniformly:

$$
\delta \Omega\_{\mathbf{D}\_1} = \int\_{\mathbf{D}\_1} \mathbf{q}\_s(\mathbf{x}) d\mathbf{x}.\tag{102}
$$

The replacement of the boundary Γ is determined by the following relation

$$
\delta \Gamma = \frac{\delta \Omega\_{\mathcal{D}\_1}}{\mathcal{S}},
\tag{103}
$$

here S is square of surface Γ (length of curve line Γ in the plane case).

If for simplest organism's linear dependence of growth velocity on metabolism intensity defined by relationship (100) is intrinsic, a tissue existence during an enough long time without growth, but under nonzero metabolism intensity, restricted by some limits, is possible for more complex multicellular organisms.

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow… DOI: http://dx.doi.org/10.5772/intechopen.93788*

### **6. Boundary element method application to biological tissue growth**

The above-developed algorithm cannot be directly applied to the twodimensional and three-dimensional problems because boundary-value problems for partial differential equations arise in the mentioned cases instead boundary-value problems for ordinary differential equations as above. Thus two- and threedimensional cases require some numerical method for solution of elliptic boundaryvalue problems in moving boundary domain. The most powerful tool for such problems is boundary element method [6, 7], which requires a reformulation of the considered problems as boundary integral equations.

Let consider the initial boundary value problem (76)–(96). Small parameter method application to this problem is, generally speaking, similar to above onedimensional case application (see, for example, [8, 9]). Restrict the following consideration by plane case and by zero approximation of small parameter method, what corresponds to very small value of the Stefan number analogue. Thus

$$
\Delta \mathbf{T}\_1^0 = -\frac{\mathbf{q}\_{\mathbf{T}1}}{\mathbf{a}\_1},
\tag{104}
$$

$$
\Delta \mathbf{C}\_{\rm i1}^{0} = -\frac{\mathbf{q}\_{\rm i1}}{\mathbf{d}\_{\rm i1}}, \mathbf{i} = \overline{\mathbf{1}, \mathbf{N}}, \tag{105}
$$

$$
\Delta \mathbf{T}\_2^0 = \mathbf{0},
\tag{106}
$$

$$
\Delta \mathbf{C}\_{\mathbf{i}2}^{0} = \mathbf{0}, \mathbf{i} = \overline{\mathbf{1}, \mathbf{N}}.\tag{107}
$$

Boundary conditions for the system (104)–(107) coincide with boundary conditions for the initial system. Let apply methods of potential theory to the system (104)–(107).

$$\begin{split} \chi(\mathbf{x}\_{0})\mathbf{T}\_{1}^{0}(\mathbf{x}\_{0}) &= \int\_{\mathbf{r}\_{1}} \phi\_{0}(\mathbf{x},\mathbf{x}\_{0}) \frac{\partial \mathbf{T}\_{1}^{0}}{\partial \mathbf{n}} \, \mathrm{d}s - \int\_{\mathbf{r}\_{1}} \mathbf{T}\_{1}^{0} \, \frac{\partial \phi\_{0}(\mathbf{x},\mathbf{x}\_{0})}{\partial \mathbf{n}} \, \mathrm{d}s + \\ &\quad + \int\_{\mathbf{D}} \phi\_{0}(\mathbf{x},\mathbf{x}\_{0}) \frac{\mathbf{q}\_{\mathbf{T}\_{1}}}{\mathbf{a}\_{1}} \, \mathrm{d}x \, \mathrm{d}y, \\ \chi(\mathbf{x}\_{0})\mathbf{C}\_{\mathrm{il}}^{0}(\mathbf{x}\_{0}) &= \int\_{\mathbf{r}\_{1}^{0}} \phi\_{0}(\mathbf{x},\mathbf{x}\_{0}) \frac{\partial \mathbf{C}\_{\mathrm{il}}^{0}}{\partial \mathbf{n}} \, \mathrm{d}s - \int\_{\mathbf{r}\_{1}^{0}} \mathbf{C}\_{\mathrm{il}}^{0} \, \frac{\partial \phi\_{0}(\mathbf{x},\mathbf{x}\_{0})}{\partial \mathbf{n}} \, \mathrm{d}s + \\ &\quad + \int\_{\mathbf{D}} \phi\_{0}(\mathbf{x},\mathbf{x}\_{0}) \frac{\mathbf{q}\_{\mathrm{il}}}{\mathbf{d}\_{\mathrm{il}}} \, \mathrm{d}x \, \mathrm{d}y, \end{split} \tag{109}$$

$$\chi(\mathbf{x}\_{0})\mathcal{T}\_{2}^{0}(\mathbf{x}\_{0}) = \int\_{\Gamma\_{2}} \phi\_{0}(\mathbf{x}, \mathbf{x}\_{0}) \frac{\partial \Gamma\_{2}^{0}}{\partial \mathbf{n}} \, \mathrm{d}\mathbf{s} - \int\_{\Gamma\_{2}} \Gamma\_{2}^{0} \frac{\partial \phi\_{0}(\mathbf{x}, \mathbf{x}\_{0})}{\partial \mathbf{n}} \, \mathrm{d}\mathbf{s},\tag{110}$$

$$\chi(\mathbf{x}\mathbf{o})\mathbf{C}\_{\mathbf{i}2}^{0}(\mathbf{x}\mathbf{o}) = \int\_{\Gamma\_{2}} \phi\_{0}(\mathbf{x}, \mathbf{x}\_{0}) \frac{\partial \mathbf{C}\_{\mathbf{i}2}^{0}}{\partial \mathbf{n}} d\mathbf{s} - \int\_{\Gamma\_{2}} \mathbf{C}\_{\mathbf{i}2}^{0} \frac{\partial \phi\_{0}(\mathbf{x}, \mathbf{x}\_{0})}{\partial \mathbf{n}} d\mathbf{s}.\tag{111}$$

Here the function ϕ0ð Þ x, x0 is well-known fundamental solution of Laplace equation, which is in plane case

qi1 ¼ χ<sup>i</sup>

the i-th parameter on the metabolism velocity function. It is evident that

intensity

determined as

ω ¼ N Y þ1

The function ω is determined experimentally. Let ω<sup>i</sup> is function of influence of

i¼1

As it is noted earlier the growth of cellular mass is proportional to metabolism

Terms indicated by }0} in last relationship correspond to regular metabolism,

Let consider a problem about motion of the boundary Γ again, in particular, let consider the case, when local volume change is determined by relation (100). The velocity (deformation) field depends on mechanical links between cells. If cells are "free" in intercellular solution the model of distributed sources in incompressible

fluid can be applied, according to which the velocity of the boundary Γ is

∂ ∂n ð

δΩD1 ¼

D1

If the cells are linked mechanically between themselves, to determine the motion of the boundary Γ it is necessary to solve an elastoplastic problem, as a rule, under large strains. Consideration of such problems requires especial investigation and

However, the another case is possible in biological structures; a biological structure grows saving its shape in this case. Thus, change of the structure volume can be

ð

D1

The replacement of the boundary Γ is determined by the following relation

δΓ <sup>¼</sup> δΩD1

here S is square of surface Γ (length of curve line Γ in the plane case). If for simplest organism's linear dependence of growth velocity on metabolism intensity defined by relationship (100) is intrinsic, a tissue existence during an enough long time without growth, but under nonzero metabolism intensity, restricted by some limits, is possible for more complex

Vnð Þ¼ x<sup>∘</sup>

where xi is arbitrary point of the curve-line Γ.

will not be made in the present work.

referred to the boundary Γ uniformly:

multicellular organisms.

**64**

which is specific for tissues of highest animals.

*Modeling and Simulation in Engineering - Selected Problems*

ω, (97)

ωi*:* (99)

qГ<sup>1</sup> ¼ χГω*:* (98)

qs ¼ χsω � χ0ω0*:* (100)

qsð Þ x ϕ∘ð Þ x, x<sup>∘</sup> dx, (101)

qsð Þ x dx*:* (102)

<sup>S</sup> , (103)

*Modeling and Simulation in Engineering - Selected Problems*

$$\varphi\_0(\mathbf{x}, \mathbf{x}\_0) = \frac{1}{2\pi} \ln \left( \frac{1}{\sqrt{\left(\mathbf{x} - \mathbf{x}\_0\right)^2 + \left(\mathbf{y} - \mathbf{y}\_0\right)^2}} \right),$$

and function χ is determined by the observation point position:

$$\chi(\mathbf{x}\_0) = \begin{cases} \mathbf{0}, (\mathbf{x}\_0) \notin \mathbf{D}, (\mathbf{x}\_0) \notin \Gamma \\ \mathbf{1}/\mathbf{2}, (\mathbf{x}\_0) \in \Gamma \\ \mathbf{1}, (\mathbf{x}\_0) \in \mathbf{D}. \end{cases}$$

The system (108)–(111) can be easy solved by conventional boundary element method. A specific feature of the problem is boundary condition on the boundary Γ ¼ Γ<sup>1</sup> ∩ Γ2, that is boundary of growth. If the forth kind boundary conditions are prescribe on the Γ (volume growth), then correspondent integral equations are simply coupled on the curve-line Γ. If correspondent fluxes on the curve-line Γ are discontinuous, then the gap value on previous time step is used.

A quite natural problem of calculation of last domain integrals in Eqs. (108) and (109) arise during the numerical solution. As a rule, it leads to serious computational difficulties, however since the time scale of growth process is enough large and the source terms in the initial Eqs. (104) and (105) are understood as averaged in time, the considered source terms are often constant with respect to space variables. The case of constant source is considered in the present work. The domain integrals can be easy transformed in this case

$$\int\_{\mathcal{D}\_1} \phi\_0(\mathbf{x}, \mathbf{x}\_0) \mathbf{q} \mathbf{dx} \mathbf{dy} = \mathbf{q} \int\_{\mathcal{D}\_1} \text{div} \mathbf{grad} \, \phi\_1 \mathbf{dx} \mathbf{dy} = \mathbf{q} \int\_{\mathcal{E}\_1} \frac{\partial \phi\_1}{\partial \mathbf{n}} \, \text{ds},\tag{112}$$

As it can be seen from **Table 2**, the growth process is accelerated with time. A caution of such acceleration is approaching of the boundary of growing biological structure to the outer boundaries of the domain, where high nutrition concentrations are prescribed according to boundary conditions (82) and (83); as a result, correspondent diffusive fluxes are increasing, what is accelerating the growth process. However, effect of nutrition consumption by more massive increased

*Growth of biological structure, which was circular at the beginning; nutritions are going into the domain from*

**Time (h) Mass of biological structure**

 0.28378030 0.28792960 0.29541710 0.30700770 0.32328381 0.34431560 0.37013320 0.40161840

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow…*

*DOI: http://dx.doi.org/10.5772/intechopen.93788*

Similarly to above considered case of Stefan problem, since there is not any known analytical solution of biological growth problem to check accuracy of the obtained numerical solutions, the authors of the present work had to change numbers of boundary elements and compare correspondent numerical results. For example, the presented above numerical calculations were made under following parameters: phase transition boundary is approximated by 40 boundary elements at every time step of calculation and outer boundaries of the domain were approximated by 4x40 boundary elements. Checking calculations with 80 and 4x80 and 160 and 4x160 boundary elements and dimensionless time step varying from Δτ ¼ 0, 1 to Δτ ¼ 0, 001, correspondingly, shows change of the obtained results no more than 1%. The authors made a conclusion about satisfactory accuracy of the

biological structure decelerates the growth process.

*Mass of growing biological structure shown in Figure 5.*

considered numerical approach.

**Figure 6.**

**Table 1.**

**67**

*left and from below [13].*

where Δϕ<sup>1</sup> <sup>¼</sup> <sup>ϕ</sup>0, that is <sup>ϕ</sup><sup>1</sup> ¼ � <sup>r</sup><sup>2</sup> <sup>8</sup><sup>π</sup> ð Þ ln r � 1 .

### **7. Results of biological tissue growth calculations**

The results of numerical calculations of model problems of growth of one-cell organism colony are shown in **Figures 5** and **6** and in **Tables 1** and **2**.

Growth in direction of maximum concentration of nutrition is evident in both cases. Note only that the structure shown in **Figures 5** and **6** initially were the same structure and only nutrition concentrations were different.

#### **Figure 5.**

*Growth of biological structure, which was circular at the beginning; nutritions are going into the domain from above and from below [13].*

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow… DOI: http://dx.doi.org/10.5772/intechopen.93788*

#### **Figure 6.**

*φ*0ð Þ¼ *x*, *x*<sup>0</sup>

*Modeling and Simulation in Engineering - Selected Problems*

1

χð Þ¼ x0

discontinuous, then the gap value on previous time step is used.

integrals can be easy transformed in this case

ϕ0ð Þ x, x0 qdxdy ¼ q

**7. Results of biological tissue growth calculations**

structure and only nutrition concentrations were different.

ð

D1

**Figure 5.**

**66**

*above and from below [13].*

where Δϕ<sup>1</sup> <sup>¼</sup> <sup>ϕ</sup>0, that is <sup>ϕ</sup><sup>1</sup> ¼ � <sup>r</sup><sup>2</sup>

0

B@

and function χ is determined by the observation point position:

8 ><

>:

<sup>2</sup>*<sup>π</sup>* ln <sup>1</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð Þ *x* � *x*<sup>0</sup>

0, xð Þ<sup>0</sup> ∉ D, xð Þ<sup>0</sup> ∉ Γ

1*=*2, xð Þ<sup>0</sup> ∈Γ 1, xð Þ<sup>0</sup> ∈ D*:*

The system (108)–(111) can be easy solved by conventional boundary element method. A specific feature of the problem is boundary condition on the boundary Γ ¼ Γ<sup>1</sup> ∩ Γ2, that is boundary of growth. If the forth kind boundary conditions are prescribe on the Γ (volume growth), then correspondent integral equations are simply coupled on the curve-line Γ. If correspondent fluxes on the curve-line Γ are

A quite natural problem of calculation of last domain integrals in Eqs. (108) and (109) arise during the numerical solution. As a rule, it leads to serious computational difficulties, however since the time scale of growth process is enough large and the source terms in the initial Eqs. (104) and (105) are understood as averaged in time, the considered source terms are often constant with respect to space variables. The case of constant source is considered in the present work. The domain

div gradϕ1dxdy ¼ q

ð

∂ϕ<sup>1</sup> ∂n

ds, (112)

Γ1

ð

D1

organism colony are shown in **Figures 5** and **6** and in **Tables 1** and **2**.

<sup>8</sup><sup>π</sup> ð Þ ln r � 1 .

The results of numerical calculations of model problems of growth of one-cell

Growth in direction of maximum concentration of nutrition is evident in both cases. Note only that the structure shown in **Figures 5** and **6** initially were the same

*Growth of biological structure, which was circular at the beginning; nutritions are going into the domain from*

<sup>2</sup> <sup>þ</sup> *<sup>y</sup>* � *<sup>y</sup>*<sup>0</sup> � �<sup>2</sup> <sup>q</sup>

1

CA,

*Growth of biological structure, which was circular at the beginning; nutritions are going into the domain from left and from below [13].*


#### **Table 1.**

*Mass of growing biological structure shown in Figure 5.*

As it can be seen from **Table 2**, the growth process is accelerated with time. A caution of such acceleration is approaching of the boundary of growing biological structure to the outer boundaries of the domain, where high nutrition concentrations are prescribed according to boundary conditions (82) and (83); as a result, correspondent diffusive fluxes are increasing, what is accelerating the growth process. However, effect of nutrition consumption by more massive increased biological structure decelerates the growth process.

Similarly to above considered case of Stefan problem, since there is not any known analytical solution of biological growth problem to check accuracy of the obtained numerical solutions, the authors of the present work had to change numbers of boundary elements and compare correspondent numerical results. For example, the presented above numerical calculations were made under following parameters: phase transition boundary is approximated by 40 boundary elements at every time step of calculation and outer boundaries of the domain were approximated by 4x40 boundary elements. Checking calculations with 80 and 4x80 and 160 and 4x160 boundary elements and dimensionless time step varying from Δτ ¼ 0, 1 to Δτ ¼ 0, 001, correspondingly, shows change of the obtained results no more than 1%. The authors made a conclusion about satisfactory accuracy of the considered numerical approach.


**Table 2.**

*Mass of growing biological structure shown in Figure 6.*

### **8. Conclusions**

Beside formulated above conclusions, some additional conclusions can be made concerning the whole investigation. As it is mentioned above, the proposed approach is more general and has more opportunities than earlier existing computational methods for slow-phase transitions. It gives an opportunity to effectively calculate of the processes in more broad interval of Stefan numbers in domains of complex geometrical shapes. Nevertheless, an accuracy of the proposed approach has not been investigated theoretically by a proper way, and enough high accuracy of the proposed algorithms is confirmed by the series of test calculations. The proposed approach can be recommended for calculations of relevant phenomena of slow phase transformations in microgravity and in environment. They also can be used in design of space-rocket techniques. The proposed approach is the only way of calculation for a lot of practically interesting cases. It is necessary to note that the phase transition process is mainly determined by zeroth approximation; another approximations starting from the first one have small influences on the process, especially for small Stefan numbers.

The main idea of the second part of the present paper is to develop a computational method for the problem of biological structure growth, based on the fact that biological growth is relatively very slow process. Considered circumstance leads to asymptotic analysis based on smallness of relation of correspondent time scales. Nevertheless the problem was formulated in quite general form, as a result of asymptotic analysis by small parameter method it is managed to build an analytical solution in one-dimensional case (what was made in previous publications, see, for example [13]) and to propose effective boundary element algorithm for numerical solution described above. Note that the above consideration is restricted only by zero-th approximation in the asymptotic expansion of the solution.

**Author details**

Ukraine

**69**

Dmytro V. Yevdokymov and Yuri L. Menshikov\*

provided the original work is properly cited.

\*Address all correspondence to: menshikov2003@mail.ru

Mechanical and Mathematical Faculty, Oles Honchar Dnipro National University,

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow…*

*DOI: http://dx.doi.org/10.5772/intechopen.93788*

© 2020 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium,

Calculations of specific biological structures did not concern the aim of the present work. However, the examples of calculations of special model problems show workability and effectiveness of the proposed method.

There is a quite natural question about applicability of the algorithm to the very important problem of tumor growth. The answer remains unclear at the moment, because it is unclear whether the used metabolism model can describe a tumor growth process or not. However there is no mathematical insuperable hindrance but only biological.

*Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow… DOI: http://dx.doi.org/10.5772/intechopen.93788*

### **Author details**

**8. Conclusions**

**Table 2.**

only biological.

**68**

especially for small Stefan numbers.

*Mass of growing biological structure shown in Figure 6.*

*Modeling and Simulation in Engineering - Selected Problems*

Beside formulated above conclusions, some additional conclusions can be made

**Time (h) Mass of biological structure**

 0.28751980 0.30294060 0.32969960 0.36981910 0.42709790 0.50863440 0.63095590

The main idea of the second part of the present paper is to develop a computational method for the problem of biological structure growth, based on the fact that biological growth is relatively very slow process. Considered circumstance leads to asymptotic analysis based on smallness of relation of correspondent time scales. Nevertheless the problem was formulated in quite general form, as a result of asymptotic analysis by small parameter method it is managed to build an analytical solution in one-dimensional case (what was made in previous publications, see, for example [13]) and to propose effective boundary element algorithm for numerical solution described above. Note that the above consideration is restricted only by

Calculations of specific biological structures did not concern the aim of the present work. However, the examples of calculations of special model problems

There is a quite natural question about applicability of the algorithm to the very important problem of tumor growth. The answer remains unclear at the moment, because it is unclear whether the used metabolism model can describe a tumor growth process or not. However there is no mathematical insuperable hindrance but

zero-th approximation in the asymptotic expansion of the solution.

show workability and effectiveness of the proposed method.

concerning the whole investigation. As it is mentioned above, the proposed approach is more general and has more opportunities than earlier existing computational methods for slow-phase transitions. It gives an opportunity to effectively calculate of the processes in more broad interval of Stefan numbers in domains of complex geometrical shapes. Nevertheless, an accuracy of the proposed approach has not been investigated theoretically by a proper way, and enough high accuracy of the proposed algorithms is confirmed by the series of test calculations. The proposed approach can be recommended for calculations of relevant phenomena of slow phase transformations in microgravity and in environment. They also can be used in design of space-rocket techniques. The proposed approach is the only way of calculation for a lot of practically interesting cases. It is necessary to note that the phase transition process is mainly determined by zeroth approximation; another approximations starting from the first one have small influences on the process,

> Dmytro V. Yevdokymov and Yuri L. Menshikov\* Mechanical and Mathematical Faculty, Oles Honchar Dnipro National University, Ukraine

\*Address all correspondence to: menshikov2003@mail.ru

© 2020 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## **References**

[1] Rubinstein L. The Stefan Problem*.* Translations of Mathematical Monographs. Vol. 27. Providence, R.L.: American Mathematical Society; 1971. pp. 347-419

[2] Lock G, Gunderson J, Quon D, Donnelly J. International Journal of Heat and Mass Transfer. 1969;**12**(11):130-145

[3] Chuang YK, Szekely J. International Journal of Heat and Mass Transfer. 1971; **14**:1285-1294

[4] Chuang Y, Szekely J. International Journal of Heat and Mass Transfer. 1972; **15**:1171-1174

[5] Khrutch V, Yevdokymov D. Bulletin of Dnipropetrovsk University. Series: Mechanics. 2002;**6**(1):9-16. (in Russian)

[6] Brebbia C, Telles J, Wrobel LC. Boundary Element Techniques. Berlin, New York: Springer-Verlag; 1984. p. 345

[7] Banerjee P, Butterfield R. Boundary Element Methods in Engineering Science. London, New York: McGraw-Hill; 1981. p. 452

[8] Bevza E, Yevdokymov D, Kochubey O. Bulletin of Donetsk National University. Series A: Natural Sciences. 2002;**12**:249-253. (in Russian)

[9] Bevza E, Yevdokymov D, Kochubey O. Bulletin of Kharkov National University. 2003;**590**:26-31. (in Russian)

[10] Brazaluk I, Yevdokymov D. Asymptotic approach and boundary element method for calculation of slow phase transitions. In: Proceedings of the 4th International Conference on Computational Methods for Thermal Problems; Georgia Tech, Atlanta, USA. 2016. pp. 286-289. Available from: http://www.thermacomp.com/public/ index.php?node=6&nm=ThermaComp

[11] Greenspan H. Journal of Theoretical Biology. 1976;**56**:229-242

**Chapter 4**

**Abstract**

surface temperature

**1. Introduction**

**71**

Environment

Modeling of the Two-Dimensional

*Nencho Deliiski, Ladislav Dzurenda and Natalia Tumbarkova*

A two-dimensional mathematical model has been created, solved, and verified for the transient nonlinear heat conduction in logs during their thawing in an air environment. For the numerical solution of the model, an explicit form of the finite-difference method in the computing medium of Visual FORTRAN Professional has been used. The chapter presents solutions of the model and its validation towards own experimental studies. During the validation of the model, the inverse task of the heat transfer has been solved for the determination of the logs' heat transfer coefficients in radial and longitudinal directions. This task has been solved also in regard to the logs'surface temperature, which depends on the mentioned coefficients. The results from the experimental and simulative investigation of 2D nonstationary temperature distribution in the longitudinal section of poplar logs with a diameter of 0.24 m, length of 0.48 m, and an initial temperature of approximately –30°C during their many hours thawing in an air environment at

**Keywords:** heat conduction, modeling, logs, thawing, heat transfer coefficients,

The duration and the energy consumption of the thermal treatment of frozen logs aimed at their thawing and plasticizing for the production of veneer in winter are very high [1–9]. For example, thawing and plasticizing of poplar and pine logs with an initial temperature of –10°C and moisture content of 0.6 kgkg<sup>1</sup> about 53

In the specialized literature, there are few reports about the temperature fields subjected to thawing in agitated water or steam frozen logs [7–21], and there is very scarce information about research of the temperature distribution in frozen logs during their thawing in an air environment given by the authors only [22].

The computation of the temperature field in logs during their thawing in water or steam is carried out using mathematical models, which solve the so-called direct task of the heat transfer. This is the task when all variables in the model are known,

The computation of the temperature field in logs during their thawing in an air environment requires solving of the so-called inverse task of the heat transfer. This

kWh<sup>m</sup><sup>3</sup> and 64 kWh<sup>m</sup><sup>3</sup> thermal energy, respectively, are needed [9].

and this allows computing the temperature field in the body [23, 24].

Thawing of Logs in an Air

room temperature are presented, visualized, and analyzed.

[12] Byrne H, Chaplin M. Mathematical Biosciences. 1995;**130**:151-181

[13] Brazaluk I, Gubin O, Yevdokymov D. Bulletin of Dnipropetrovsk University. Series: Mechanics. 2015:96-114

### **Chapter 4**

**References**

pp. 347-419

**14**:1285-1294

**15**:1171-1174

Hill; 1981. p. 452

(in Russian)

**70**

[8] Bevza E, Yevdokymov D, Kochubey O. Bulletin of Donetsk National University. Series A: Natural Sciences. 2002;**12**:249-253. (in Russian)

[9] Bevza E, Yevdokymov D, Kochubey O. Bulletin of Kharkov National University. 2003;**590**:26-31.

[10] Brazaluk I, Yevdokymov D. Asymptotic approach and boundary element method for calculation of slow phase transitions. In: Proceedings of the

4th International Conference on Computational Methods for Thermal Problems; Georgia Tech, Atlanta, USA. 2016. pp. 286-289. Available from: http://www.thermacomp.com/public/ index.php?node=6&nm=ThermaComp

[1] Rubinstein L. The Stefan Problem*.*

*Modeling and Simulation in Engineering - Selected Problems*

[11] Greenspan H. Journal of Theoretical

[12] Byrne H, Chaplin M. Mathematical

Biology. 1976;**56**:229-242

[13] Brazaluk I, Gubin O, Yevdokymov D. Bulletin of Dnipropetrovsk University. Series:

Mechanics. 2015:96-114

Biosciences. 1995;**130**:151-181

Monographs. Vol. 27. Providence, R.L.: American Mathematical Society; 1971.

[3] Chuang YK, Szekely J. International Journal of Heat and Mass Transfer. 1971;

[4] Chuang Y, Szekely J. International Journal of Heat and Mass Transfer. 1972;

[5] Khrutch V, Yevdokymov D. Bulletin of Dnipropetrovsk University. Series: Mechanics. 2002;**6**(1):9-16. (in Russian)

[6] Brebbia C, Telles J, Wrobel LC. Boundary Element Techniques. Berlin, New York: Springer-Verlag; 1984. p. 345

[7] Banerjee P, Butterfield R. Boundary Element Methods in Engineering Science. London, New York: McGraw-

[2] Lock G, Gunderson J, Quon D, Donnelly J. International Journal of Heat and Mass Transfer. 1969;**12**(11):130-145

Translations of Mathematical
