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

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

λ1∂T1∂nГ1=q1,E5

λ2∂T2∂nГ2=q2,E6

and boundary conditions on the phase transition boundary

T1Гp.t.=Tp.t.,E7

T2Гp.t.=Tp.t.,E8

λ1∂T1∂nГp.t.−λ2∂T2∂nГ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=T1−Tp.t.Tn−Tp.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

∂θ1∂Fo1=Δ∗θ1,E13

∂θ2∂Fo2=Δ∗θ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

∂θ1∂n−fλ∂θ2∂n=∂n∗∂τst,E24

where

fλ=λ2/λ1,E25

and

τst=τλ1Tn−Tp.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=τλ1Tn−Tp.t.L2σρτa1L2=λ1Tn−Tp.t.a1σρ=CTn−Tp.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

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

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:

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:

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

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

or the second kind boundary condition

or the third boundary condition

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

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

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

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

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

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

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

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:

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

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

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

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

and function χis determined by the observation point position:

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πlnr−1.

7. Results of biological tissue growth calculations

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

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

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.