Open access peer-reviewed chapter

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

Written By

Dmytro V. Yevdokymov and Yuri L. Menshikov

Submitted: December 6th, 2019 Reviewed: August 28th, 2020 Published: October 6th, 2020

DOI: 10.5772/intechopen.93788

From the Edited Volume

## Modeling and Simulation in Engineering

Edited by Jan Valdman and Leszek Marcinkowski

Chapter metrics overview

View Full Metrics

## 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 initial-boundary-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.

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 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 boundary-value 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 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.

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.

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-boundary-value 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 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 1 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 T p . t . . Let domains D 1 and D 2 are bounded by finite or infinite curves Γ 1 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 Γ 1 and Γ 2 . Thus, assuming thermophysical properties of materials as constant, we obtain the following mathematical model [1, 2].

T 1 τ = a 1 Δ T 1 , E1
T 2 τ = a 2 Δ T 2 , E2

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 Γ 1 and Γ 2 :

In particular, first kind

T 1 Г 1 = T 1 s , E3
T 2 Г 2 = T 2 s , E4

or second kind

λ 1 T 1 n Г 1 = q 1 , E5
λ 2 T 2 n Г 2 = q 2 , E6

and boundary conditions on the phase transition boundary

T 1 Г p . t . = T p . t . , E7
T 2 Г p . t . = T p . t . , E8
λ 1 T 1 n Г p . t . λ 2 T 2 n Г p . t . = ρσ V p . t . , E9

where a 1 , a 2 are thermal diffusivity coefficients of phases, T 1 s , T 2 s , q 1 , q 2 are temperatures and thermal fluxes on the outside parts of curves Γ 1 and Γ 2 , λ 1 , λ 2 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), V p . 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 T p . t . is constant of the material. To complete the formulation, add the initial conditions:

T 1 0 x = T 10 x , E10
T 2 0 x = T 20 x . E11

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

θ 1 = T 1 T p . t . T n T p . t . , E12

where value T n is chosen as one from following values max T 1 s T 2 s or min T 1 s T 2 s (but sometimes it is possible to use max T 10 T 20 or min T 10 T 20 ) 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

θ 1 Fo 1 = Δ θ 1 , E13
θ 2 Fo 2 = Δ θ 2 , E14

where

Fo 1 = τ a 1 L 2 . E15

The asterisk “*” in Laplace operator means the differentiation with respect to dimensionless coordinates

X = x / L . E16

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

θ 1 Г 1 = θ 1 s , E17
θ 2 Г 2 = θ 2 s , E18
θ 1 Г p . t . = 0 , E19
θ 2 Г p . t . = 0 , E20
θ 1 0 x = θ 10 , E21
θ 2 0 x = θ 20 . E22

The dimensionless Stefan condition must be considered in details. Let

V p . t . = n τ , E23

where n τ means the normal velocity of the phase transformation boundary motion. If n = n / L , we have following dimensionless relation

θ 1 n f λ θ 2 n = n τ st , E24

where

f λ = λ 2 / λ 1 , E25

and

τ st = τλ 1 T n T p . t . L 2 σρ , E26

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

Thus there three dimensionless time Fo 1 , Fo 2 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

θ 1 τ st St = Δ θ 1 , E27

where St is the Stefan number determined as

St = τ st Fo 1 = τλ 1 T n T p . t . L 2 σρ τ a 1 L 2 = λ 1 T n T p . t . a 1 σρ = C T n T p . t . σ , E28

Eq. (14) is transformed by following way:

θ 2 τ st f a St = Δ θ 2 , E29

where

f a = Fo 1 Fo 2 = a 1 a 2 . E30

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

θ τ st St = Δ θ , E31
θ Г 1 = θ s , E32
θ Г p . t . = 0 , E33
θ 0 x = θ 0 , E34
θ n = n τ st . E35

### 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

θ 1 = θ 1 0 x τ + k = 1 St k θ 1 k x τ , E36
θ 2 = θ 2 0 x τ + k = 1 St k θ 2 k x τ . E37

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

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

St τ θ 1 0 x τ + St τ k = 1 St k θ 1 k x τ = Δ θ 1 0 x τ + Δ k = 1 St k θ 1 k x τ , E38
f a St τ θ 2 0 x τ + f a St τ k = 1 St k θ 2 k x τ = Δ θ 2 0 x τ + Δ k = 1 St k θ 2 k x τ .

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

St θ 1 0 τ + k = 1 St k + 1 θ 1 k τ = Δ θ 1 0 + k = 1 St k Δ θ 1 k , E39
f a St θ 2 0 τ + f a k = 1 St k + 1 θ 2 k τ = Δ θ 2 0 + k = 1 St k Δ θ 2 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

Δ θ 1 0 = 0 , E40
Δ θ 1 1 = θ 1 0 τ , E41
Δ θ 1 i = θ 1 i 1 τ , E42
Δ θ 2 0 = 0 , E43
Δ θ 2 1 = f a θ 2 0 τ , E44
Δ θ 1 i = f a θ 1 i 1 τ , E45

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

θ 1 0 Г 1 = θ 1 s , E46
θ 1 0 Г p . t . = 0 , E47
θ 1 i Г 1 = 0 , E48
θ 1 i Г p . t . = 0 , E49
θ 2 0 Г 2 = θ 2 s , E50
θ 2 0 Г p . t . = 0 , E51
θ 2 i Г 2 = 0 , E52
θ 2 i Г p . t . = 0 . E53

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.

θ 1 0 x τ n + k = 1 St k θ 1 k x τ n = f a θ 2 0 x τ n + k = 1 St k θ 2 k x τ n = = η 0 x τ τ + k = 1 St k η k x τ τ . E54

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

θ 1 0 n f λ θ 2 0 n = η 0 τ , E55
θ 1 i n f λ θ 2 i n = η i τ , E56

where functions η k 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 D 1 and D 2 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 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 St j + 1 M j + 1 , where M j + 1 = max j + 1 m < θ m (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 θ i by M i . 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:

Δ θ 0 = 0 , E57
Δ θ 1 = θ 0 τ , E58
Δ θ i = θ i 1 τ E59

with boundary conditions

θ 0 Γ 1 = θ s , E60
θ 0 Γ p . t . = 0 , E61
θ i Γ 1 = 0 , E62
θ i Γ p . t . = 0 , E63

Stefan condition

θ 0 n = η 0 τ E64
θ i n = η i τ E65

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

However, the temperature fields must be determined numerically in two-dimensional 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

C x 0 θ 1 0 x 0 = Γ 1 Γ p . t . ϕ 0 x x 0 θ 1 0 n ds Γ 1 Γ p . t . θ 1 0 ϕ 0 x x 0 n ds , E66
C x 0 θ 2 0 x 0 = Γ 2 Γ p . t . ϕ 0 x x 0 θ 2 0 n ds Γ 2 Γ p . t . θ 2 0 ϕ 0 x x 0 n ds , E67
C x 0 θ 1 i x 0 = Γ 1 Γ p . t . ϕ 0 x x 0 θ 1 i n ds Γ 1 Γ p . t . θ 1 i ϕ 0 x x 0 n ds + D 1 ϕ 0 x x 0 θ 1 i 1 τ dx , E68
C x 0 θ 2 i x 0 = Γ 2 Γ p . t . ϕ 0 x x 0 θ 2 i n ds Γ 2 Γ p . t . θ 2 i ϕ 0 x x 0 n ds + D 2 ϕ 0 x x 0 θ 2 i 1 τ dx , E69

where ϕ 0 is fundamental solution of Laplace equation, C x 0 is special function reflecting control point x 0 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 time-stepping 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.

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.

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 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 D 1 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 D 1 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 D 1 is partially or completely surrounded by the domain D 2 , filled by the same solution completely. In general case there may be a convective transfer in the domain D 2 and filtration flow in the domain D 1 . Thus, a general mathematical model of heat and mass transfer processes is considered system is following:

T 1 τ + V f T 1 = a 1 ΔT 1 + q T 1 , E70
C i 1 τ + V f C i 1 = d i 1 Δ C i 1 + q i 1 , i = 1 , N , ¯ E71
T 2 τ + V C T 2 = a 2 ΔT 2 , E72
C i 2 τ + V C C i 2 = d i 2 Δ C i 2 , i = 1 , N , ¯ E73

where T 1 is temperature in the domain D 1 (the one-temperature model, assuming the temperatures of frame and solution in pores are equal, is used here), V f is filtration velocity, a 1 is thermal diffusivity of porous media, q T 1 is heat source, concerning the metabolism of cells, c i 1 is concentration of the i -th component in porous media, d i 1 is diffusion coefficient of i -th component in the porous media, q i 1 is source (sink) of the i -th component in porous media, concerning the metabolism of cells, T 2 is temperature in the domain D 2 , V 2 is flow velocity in the domain D 2 , a 2 is thermal diffusivity of solution, c i 2 is concentration of i -th component in the domain D 2 , d i 2 is diffusion coefficient of i -th component in the domain D 2 , N i 2 number of components, participating of heat and mass transfer process, τ is time, Δ is Laplace operator.

Restrict the following consideration by the case:

V f = 0 , E74
V c = 0 , E75

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

T 1 τ = a 1 ΔT 1 + q T 1 , E76
C i 1 τ = d i 1 Δ C i 1 + q i 1 , i = 1 , N , ¯ E77
T 2 τ = a 2 ΔT 2 , E78
C i 2 τ = d i 2 Δ C i 2 , i = 1 , N . ¯ E79

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 D 1 and D 2 as Γ and reminder part as Γ 1 and Γ 2 correspondingly. The first kind boundary conditions can be prescribed on the boundaries Γ 1 and Γ 2

T 1 Г 1 = T 1 e , E80
C i 1 Г 1 = С i 1 e , E81
T 2 Г 2 = T 2 e , E82
C i 2 Г 2 = C i 2 e , E83

or the second kind boundary condition

λ 1 T 1 n Г 1 = q 1 e , E84
d 1 i C i 1 n Г 1 = q i 1 e , E85
λ 2 T 2 n Г 2 = q 2 e , E86
d i 2 C i 2 n Г 2 = q i 2 e , E87

or the third boundary condition

λ 1 T 1 n Г 1 + α 1 T 1 Г 1 T 1 e = 0 , E88
d i 1 C i 1 n Г 1 + α i 1 C i 1 Г 1 С i 1 e = 0 , E89
λ 2 T 2 n Г 2 + α 2 T 2 Г 2 T 2 e = 0 , E90
d i 2 C i 2 n Г 2 + α i 2 C i 2 Г 2 С i 2 e = 0 , E91

where T 1 e , C i 1 e , T 2 e , C i 2 e , q 1 e , q i 1 e , q 2 e , q i 2 e 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

T 1 Г = T 2 Г , E92
С i 1 Г = C i 2 Г . E93

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

λ 1 T 1 n Г = λ 2 T 2 n Г , E94
d i 1 C i 1 n Г = d i 2 C i 2 n Г . E95

Conditions (94) and (95) correspond to the case of cell fission in whole domain D 1 . 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

d i 1 C i 1 n Г d i 2 C i 2 n Г = χ i n τ , E96

here n τ is velocity of the boundary Γ propagation (velocity of biological structure growth), χ i 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 D 1 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 ω T 1 c i 1 is defined, then correspondent source terms are following

q i 1 = χ i ω , E97
q Г 1 = χ Г ω . E98

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

ω = i = 1 N + 1 ω i . E99

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

q s = χ s ω χ 0 ω 0 . E100

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

V n x = n D 1 q s x ϕ x x dx , E101

where x i 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:

δ Ω D 1 = D 1 q s x dx . E102

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

δ Γ = δ Ω D 1 S , E103

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.

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

The above-developed algorithm cannot be directly applied to the two-dimensional 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-dimensional cases require some numerical method for solution of elliptic boundary-value 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 one-dimensional 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

ΔT 1 0 = q T 1 a 1 , E104
Δ C i 1 0 = q i 1 d i 1 , i = 1 , N , ¯ E105
ΔT 2 0 = 0 , E106
Δ C i 2 0 = 0 , i = 1 , N . ¯ E107

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).

χ x 0 T 1 0 x 0 = Γ 1 ϕ 0 x x 0 T 1 0 n ds Γ 1 T 1 0 ϕ 0 x x 0 n ds + + D ϕ 0 x x 0 q T 1 a 1 dxdy , E108
χ x 0 C i 1 0 x 0 = Γ 1 ϕ 0 x x 0 C i 1 0 n ds Γ 1 C i 1 0 ϕ 0 x x 0 n ds + + D ϕ 0 x x 0 q i 1 d i 1 dxdy , E109
χ x 0 T 2 0 x 0 = Γ 2 ϕ 0 x x 0 T 2 0 n ds Γ 2 T 2 0 ϕ 0 x x 0 n ds , E110
χ x 0 C i 2 0 x 0 = Γ 2 ϕ 0 x x 0 C i 2 0 n ds Γ 2 C i 2 0 ϕ 0 x x 0 n ds . E111

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

φ 0 x x 0 = 1 2 π ln 1 x x 0 2 + y y 0 2 ,

and function χ is determined by the observation point position:

χ x 0 = 0 , x 0 D , x 0 Γ 1 / 2 , x 0 Γ 1 , x 0 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 Γ = Γ 1 Γ 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

D 1 ϕ 0 x x 0 qdxdy = q D 1 div grad ϕ 1 dxdy = q Γ 1 ϕ 1 n ds , E112

where Δ ϕ 1 = ϕ 0 , that is ϕ 1 = r 2 8 π 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.

Time (h) Mass of biological structure
0 0.28378030
1 0.28792960
2 0.29541710
3 0.30700770
4 0.32328381
5 0.34431560
6 0.37013320
7 0.40161840

### Table 1.

Mass of growing biological structure shown in Figure 5.

Time (h) Mass of biological structure
0 0.28751980
1 0.30294060
2 0.32969960
3 0.36981910
4 0.42709790
5 0.50863440
6 0.63095590

### Table 2.

Mass of growing biological structure shown in Figure 6.

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.

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.

## 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.

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.

## References

1. 1. Rubinstein L. The Stefan Problem. Translations of Mathematical Monographs. Vol. 27. Providence, R.L.: American Mathematical Society; 1971. pp. 347-419
2. 2. Lock G, Gunderson J, Quon D, Donnelly J. International Journal of Heat and Mass Transfer. 1969;12(11):130-145
3. 3. Chuang YK, Szekely J. International Journal of Heat and Mass Transfer. 1971;14:1285-1294
4. 4. Chuang Y, Szekely J. International Journal of Heat and Mass Transfer. 1972;15:1171-1174
5. 5. Khrutch V, Yevdokymov D. Bulletin of Dnipropetrovsk University. Series: Mechanics. 2002;6(1):9-16. (in Russian)
6. 6. Brebbia C, Telles J, Wrobel LC. Boundary Element Techniques. Berlin, New York: Springer-Verlag; 1984. p. 345
7. 7. Banerjee P, Butterfield R. Boundary Element Methods in Engineering Science. London, New York: McGraw-Hill; 1981. p. 452
8. 8. Bevza E, Yevdokymov D, Kochubey O. Bulletin of Donetsk National University. Series A: Natural Sciences. 2002;12:249-253. (in Russian)
9. 9. Bevza E, Yevdokymov D, Kochubey O. Bulletin of Kharkov National University. 2003;590:26-31. (in Russian)
10. 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. 11. Greenspan H. Journal of Theoretical Biology. 1976;56:229-242
12. 12. Byrne H, Chaplin M. Mathematical Biosciences. 1995;130:151-181
13. 13. Brazaluk I, Gubin O, Yevdokymov D. Bulletin of Dnipropetrovsk University. Series: Mechanics. 2015:96-114

Written By

Dmytro V. Yevdokymov and Yuri L. Menshikov

Submitted: December 6th, 2019 Reviewed: August 28th, 2020 Published: October 6th, 2020