Open access peer-reviewed chapter - ONLINE FIRST

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

By Dmytro V. Yevdokymov and Yuri L. Menshikov

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

DOI: 10.5772/intechopen.93788

## 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 D1and the second phase occupies the domain D2. The boundary Гp.t.separates mentioned domains, it is a phase transition boundary and has a constant temperature Tp.t.. Let domains D1and D2are bounded by finite or infinite curves Γ1and Γ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 Γ1and Γ2. Thus, assuming thermophysical properties of materials as constant, we obtain the following mathematical model [1, 2].

T1τ=a1ΔT1,E1
T2τ=a2ΔT2,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 Γ1and Γ2:

In particular, first kind

T1Г1=T1s,E3
T2Г2=T2s,E4

or second kind

λ1T1nГ1=q1,E5
λ2T2nГ2=q2,E6

and boundary conditions on the phase transition boundary

T1Гp.t.=Tp.t.,E7
T2Гp.t.=Tp.t.,E8
λ1T1nГp.t.λ2T2nГp.t.=ρσVp.t.,E9

where a1,a2are thermal diffusivity coefficients of phases, T1s,T2s,q1,q2are temperatures and thermal fluxes on the outside parts of curves Γ1and Γ2, λ1,λ2are 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:

T10x=T10x,E10
T20x=T20x.E11

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

θ1=T1Tp.t.TnTp.t.,E12

where value Tnis chosen as one from following values maxT1sT2sor minT1sT2s(but sometimes it is possible to use maxT10T20or minT10T20) 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

θ1Fo1=Δθ1,E13
θ2Fo2=Δθ2,E14

where

Fo1=τa1L2.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=θ1s,E17
θ2Г2=θ2s,E18
θ1Гp.t.=0,E19
θ2Гp.t.=0,E20
θ10x=θ10,E21
θ20x=θ20.E22

The dimensionless Stefan condition must be considered in details. Let

Vp.t.=nτ,E23

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

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

where

fλ=λ2/λ1,E25

and

τst=τλ1TnTp.t.L2σρ,E26

τstis dimensionless time, connected with the phase transformation process.

Thus there three dimensionless time Fo1,Fo2and τstin 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 τstas the main time scale. Beside of that, as a rule, τstis the “most slow” dimensionless time. Then Eq. (13) must be replaced by

θ1τstSt=Δθ1,E27

where Stis the Stefan number determined as

St=τstFo1=τλ1TnTp.t.L2σρτa1L2=λ1TnTp.t.a1σρ=CTnTp.t.σ,E28

Eq. (14) is transformed by following way:

θ2τstfaSt=Δθ2,E29

where

fa=Fo1Fo2=a1a2.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

θτstSt=Δθ,E31
θГ1=θs,E32
θГp.t.=0,E33
θ0x=θ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<<1it must converge quickly.

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

θ1=θ10xτ+k=1Stkθ1kxτ,E36
θ2=θ20xτ+k=1Stkθ2kxτ.E37

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

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

Stτθ10xτ+Stτk=1Stkθ1kxτ=Δθ10xτ+Δk=1Stkθ1kxτ,E38
faStτθ20xτ+faStτk=1Stkθ2kxτ=Δθ20xτ+Δk=1Stkθ2kxτ.

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

Stθ10τ+k=1Stk+1θ1kτ=Δθ10+k=1StkΔθ1k,E39
faStθ20τ+fak=1Stk+1θ2kτ=Δθ20+k=1StkΔθ2k.

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

Δθ10=0,E40
Δθ11=θ10τ,E41
Δθ1i=θ1i1τ,E42
Δθ20=0,E43
Δθ21=faθ20τ,E44
Δθ1i=faθ1i1τ,E45

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

θ10Г1=θ1s,E46
θ10Гp.t.=0,E47
θ1iГ1=0,E48
θ1iГp.t.=0,E49
θ20Г2=θ2s,E50
θ20Гp.t.=0,E51
θ2iГ2=0,E52
θ2iГ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.

θ10xτn+k=1Stkθ1kxτn=faθ20xτn+k=1Stkθ2kxτn==η0xττ+k=1Stkηkxττ.E54

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

θ10nfλθ20n=η0τ,E55
θ1infλθ2in=ηiτ,E56

where functions ηkare 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 D1and D2are 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 OStj+1Mj+1,where Mj+1=maxj+1m<θ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 θiby 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:

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

with boundary conditions

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

Stefan condition

θ0n=η0τE64
θin=η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

Cx0θ10x0=Γ1Γp.t.ϕ0xx0θ10ndsΓ1Γp.t.θ10ϕ0xx0nds,E66
Cx0θ20x0=Γ2Γp.t.ϕ0xx0θ20ndsΓ2Γp.t.θ20ϕ0xx0nds,E67
Cx0θ1ix0=Γ1Γp.t.ϕ0xx0θ1indsΓ1Γp.t.θ1iϕ0xx0nds+D1ϕ0xx0θ1i1τdx,E68
Cx0θ2ix0=Γ2Γp.t.ϕ0xx0θ2indsΓ2Γp.t.θ2iϕ0xx0nds+D2ϕ0xx0θ2i1τdx,E69

where ϕ0is fundamental solution of Laplace equation, Cx0is special function reflecting control point x0position 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,1to Δτ=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 D1filled 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 D1is 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 D1is 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 D2and filtration flow in the domain D1. Thus, a general mathematical model of heat and mass transfer processes is considered system is following:

T1τ+VfT1=a1ΔT1+qT1,E70
Ci1τ+VfCi1=di1ΔCi1+qi1,i=1,N,¯E71
T2τ+VCT2=a2ΔT2,E72
Ci2τ+VCCi2=di2ΔCi2,i=1,N,¯E73

where T1is temperature in the domain D1(the one-temperature model, assuming the temperatures of frame and solution in pores are equal, is used here), Vfis filtration velocity, a1is thermal diffusivity of porous media, qT1is heat source, concerning the metabolism of cells, ci1is concentration of the i-th component in porous media, di1is diffusion coefficient of i-th component in the porous media, qi1is source (sink) of the i-th component in porous media, concerning the metabolism of cells, T2is temperature in the domain D2,V2is flow velocity in the domain D2,a2is thermal diffusivity of solution, ci2is concentration of i-th component in the domain D2, di2is 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:

Vf=0,E74
Vc=0,E75

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

T1τ=a1ΔT1+qT1,E76
Ci1τ=di1ΔCi1+qi1,i=1,N,¯E77
T2τ=a2ΔT2,E78
Ci2τ=di2ΔCi2,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 D1and D2as Γand reminder part as Γ1and Γ2correspondingly. The first kind boundary conditions can be prescribed on the boundaries Γ1and Γ2

T1Г1=T1e,E80
Ci1Г1=Сi1e,E81
T2Г2=T2e,E82
Ci2Г2=Ci2e,E83

or the second kind boundary condition

λ1T1nГ1=q1e,E84
d1iCi1nГ1=qi1e,E85
λ2T2nГ2=q2e,E86
di2Ci2nГ2=qi2e,E87

or the third boundary condition

λ1T1nГ1+α1T1Г1T1e=0,E88
di1Ci1nГ1+αi1Ci1Г1Сi1e=0,E89
λ2T2nГ2+α2T2Г2T2e=0,E90
di2Ci2nГ2+αi2Ci2Г2Сi2e=0,E91

where T1e,Ci1e,T2e,Ci2e,q1e,qi1e,q2e,qi2eare 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

T1Г=T2Г,E92
Сi1Г=Ci2Г.E93

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

λ1T1nГ=λ2T2nГ,E94
di1Ci1nГ=di2Ci2nГ.E95

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

di1Ci1nГdi2Ci2nГ=χinτ,E96

here nτis velocity of the boundary Γpropagation (velocity of biological structure growth), χiis “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 D1is 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 ωT1ci1is defined, then correspondent source terms are following

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

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

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

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

qs=χ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

Vnx=nD1qsxϕxxdx,E101

where xiis 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:

δΩD1=D1qsxdx.E102

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

δΓ=δΩD1S,E103

here Sis 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

ΔT10=qT1a1,E104
ΔCi10=qi1di1,i=1,N,¯E105
ΔT20=0,E106
ΔCi20=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).

χx0T10x0=Γ1ϕ0xx0T10ndsΓ1T10ϕ0xx0nds++Dϕ0xx0qT1a1dxdy,E108
χx0Ci10x0=Γ1ϕ0xx0Ci10ndsΓ1Ci10ϕ0xx0nds++Dϕ0xx0qi1di1dxdy,E109
χx0T20x0=Γ2ϕ0xx0T20ndsΓ2T20ϕ0xx0nds,E110
χx0Ci20x0=Γ2ϕ0xx0Ci20ndsΓ2Ci20ϕ0xx0nds.E111

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

φ0xx0=12πln1xx02+yy02,

and function χis determined by the observation point position:

χx0=0,x0D,x0Γ1/2,x0Γ1,x0D.

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

where Δϕ1=ϕ0, that is ϕ1=r28πlnr1.

## 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
00.28378030
10.28792960
20.29541710
30.30700770
40.32328381
50.34431560
60.37013320
70.40161840

### Table 1.

Mass of growing biological structure shown in Figure 5.

Time (h)Mass of biological structure
00.28751980
10.30294060
20.32969960
30.36981910
40.42709790
50.50863440
60.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,1to Δτ=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.

## How to cite and reference

### Cite this chapter Copy to clipboard

Dmytro V. Yevdokymov and Yuri L. Menshikov (October 6th 2020). Mathematical Modelling and Numerical Simulation of Diffusive Processes in Slow Changing Domains [Online First], IntechOpen, DOI: 10.5772/intechopen.93788. Available from: