MODELING OF MOLECULAR DIFFUSION IN NON-HOMOGENEOUS MEDIA

Background. Due to combined effects of medium inhomogeneity and action of external forces, i.e. Earth’s gravitation, electro-magnetic forces, global rotation, etc., a number of specific fluid motions appear in the environmental and life systems even in absence of pure mechanical reasons. Among them are so called diffusion-induced flows which always exist around obstacles with arbitrary geometry due to breaking of natural molecular flux of a stratifying agent on an impermeable surface. Objective. The aim of the paper is to analyze numerically a diffusion induced flow structure and dynamics around motionless obstacles immersed into a stably stratified medium, including a sloping plate, a wedge-shaped obstacle, a disc, and a circular cylinder. The numerical results obtained are compared with the available analytical and experimental data. Methods. The problem is solved numerically using two different algorithms based on the finite difference and finite volume methods. The first approach is implemented in the specially developed Fortran program codes, and the second one is based on the open source package OpenFOAM with the use of C programming language for developing special own solvers, libraries, and utilities, which enable solving the problems under consideration. Results. The numerical simulation reveals a complex multi-level vortex system of compensatory circulating flows around a motionless obstacle, which structure is strongly dependent on its position relative to the horizon. The most intensive and extended high-gradient horizontal interfaces attached to sharp edges or poles of obstacles are clearly observed in the numerical computations and laboratory experiments. Diffusion-induced flows form intensive pressure deficiency zones around an obstacle, which may lead to generation of propulsion mechanism resulting in self-movement of neutral buoyancy bodies in a continuously stratified fluid, e.g. horizontal movement of a wedge, rotation of a sloping plate, etc. Conclusions. Diffusion-induced flows are a wide-spread phenomenon in biology, medicine, and environmental systems, since such flows inevitably occur in any inhomogeneous media, including different solutions, suspensions, mixtures, etc. A complex multilevel vortex structure of diffusion-induced flows on an obstacle becomes even more complicated due to self-motion of the obstacle itself and displacement of various admixtures, suspended particles, additives, etc., which are always present in the real systems.


Introduction
The environmental and life systems are mostly inhomogeneous in space and variable in time due to non-uniform distributions of dissolved matters, suspended particles, gas bubbles, temperature, medium compressibility, etc.As a result of the combined influence of medium inhomogeneity and the effects of external forces, i.e.Earth's gravitation, magnetic field, global rotation, etc., a stable stratification is formed, which gives rise to a number of effects and fluid motions absent in homogeneous media, including specific types of waves, flows, fine structure flow elements, even in the absence of purely mechanical reasons.Among such phenomena are convective flows driven by spatial variations in fluid density or the so-called diffusion-induced flows on geometrical irregularities of topography.Investigations of such flows have received much attention in laboratory studies and numerical and analytical modeling because there are abundant instances of the phenomena in environmental and life systems.
A diffusion-induced flow occurs when molecular flux of the substance existing naturally in a stably stratified fluid encounters a sloping boundary that creates a physical basis for the formation of compensating fluid motions even in the absence of mechanical reasons.Molecular diffusion requires that surfaces of constant density (isopycnals) approach a sloping boundary at an appropriate angle to ensure the condition of continuity of diffusive flux.For a barrier impermeable to a stratifying agent, which is the common salt, NaCl, in a particular case under consideration, isopycnals must be locally perpendicular to the slope in order to satisfy the no-flux boundary condition, / 0,    where  is the wall-normal coordinate.A fluid adjacent to the sloping boundary therefore differs in density from a fluid at the same level away from the boundary, and experiences a buoyancy force that drives the flow up-or down-slope for the underlying and overlying walls, respectively, the kinetic energy of motion being drawn from the internal energy of the fluid through diffusion.
Diffusion-induced flows form the basis for a variety of physical processes, including mineral transport in rocks [1], the melting of icebergs [2], the migration of tectonic plates [3], the processes of transport and mixing of passive substances, the formation of intensive valley and mountain winds in a stably stratified atmosphere [4] and density flows in the ocean [5].
A permanent interest in the problem exists due to the enlarging spectrum of practical applications including ecological and biological ones.It was found experimentally that a diffusion-induced flow may trigger a propulsion mechanism leading to a self-movement of neutral buoyancy solids with a special shape in a continuously stratified ocean medium.The laboratory experiments [6] show that when a neutrally buoyant wedge-shaped object floats in a density-stratified fluid, the diffusion-induced flow generated by molecular diffusion at its sloping boundaries produces a macroscopic sideways thrust.Computer simulations reveal that thrust results from the diffusioninduced flow creating a region of low pressure at the front, relative to the rear of an object.This discovery has implications for transport processes in regions of varying fluid density, such as marine snow aggregation at ocean pycnoclines, and wherever there is a temperature difference between immersed objects and the surrounding fluid, such as particles in volcanic clouds.An essential influence of diffusion-induced flows upon the displacement of separate elements of plankton in ocean media is maintained not only by the breaking of background diffusion flux of the stratifying agent but also by additional concentration gradients generated by the products of natural metabolism of the elementary organisms.A variety of different applications encouraged the development of new approaches for the problem under study and stimulated the application of advanced data analysis techniques for a unified characterization of experimental and numerical data sets, for the estimation of their quality and information capacity.
The first mathematical model of diffusion-induced flows was created in the middle of the past century with application to the problem on mountain and valley winds [4].The marked universality mechanism of hillside flows formation both in the stratified atmosphere and the ocean [1,7] restored the interest in the problem after a long pause.In the pioneering theoretical investigations of diffusion-induced flows the stationary solutions were obtained in the infinite plate approximation on the basis of a linearized system of equations [4].For constant slope, and assuming along-slope flow that is a function only of the wall-normal coordinate (that is, unidirectional flow), the velocity profile for a steady diffusion-induced flow is ( sin /4 ) is the buoyancy frequency, 0  is the characteristic fluid density, S  is the molecular diffusivity and  is the kinematic viscosity.The non-stationary problem was analyzed first asymptotically on an infinite slope in the work by Linden and Weber [8] and was then solved exactly by Kistovich and Chashechkin [9].

Problem statement
The goal of the present work consists in construction of a numerical algorithm for computing diffusion induced flow structure and dynamics around motionless obstacles immersed into a stably stratified medium without any restrictions on their geometrical form and angular position in space.
As a validation of the developed numerical approach the computation results on diffusion induced flows around a sloping plate, a wedge-shaped obstacle, a disc, and a circular cylinder are compared with the available analytical data and experimental images obtained by the high-resolution Schlieren technique [10].
A particular goal of the investigations consists in finding physical mechanisms for generation of propulsion force resulting in self-motion of neutralbuoyancy bodies in a stratified medium, e.g.horizontal movement of a wedge, rotation of a sloping plate, etc.For this purpose all the kinematic and dynamic characteristics of diffusion-induced flow around a wedge and a plate are analyzed in details for various geometrical parameters of the problem.
A special emphasis is made in the paper to the universality mechanism of diffusion-induced flows formation regardless of a way of forming medium inhomogeneity and a type of acting external force, which both together form a stable stratification.This enables applying the problem under consideration not only to environmental systems but, as well to a wide spectrum of problems in biological and medical sciences.

Mathematical modeling
As a mathematical model for the problem under consideration the fundamental system of equation for multicomponent inhomogeneous incom-pressible fluid in the Boussinesq approximation is considered taking into account the buoyancy and diffusion effects of a stratifying agent (salinity) [11].In the study of slow, as compared to the speed of sound, flows of fluids characterized by high thermal conductivity, one can account in calculations only for variations in density associated with concentration of the stratifying agent neglecting temperature variations.Thus, the governing equations take the following form [12] 00 00 (exp( / ) ), div 0, Here, s is the salinity perturbation including the salt compression ratio, (v , v , v ) is the vector of the induced velocity, P is the pressure except for the hydrostatic one, An obstacle with impermeable boundaries is submerged with minimum disturbances into a quiescent stratified environment.Physically valid initial and boundary conditions in the associated coordinate system are no-slip and no-flux on the surface of the obstacle for velocity components and total salinity respectively, and vanishing of all the perturbations at infinity, where, n is external normal unit vector to the surface, , of an obstacle which can be either a plate or a wedge with length, L, and thickness, , izing speed of fluid motion due to the molecular processes in a medium [13].

Numerical simulation
Numerical solution of the system (1) with the boundary conditions (2) was constructed using the finite volume method within the frame of the open source computational package OpenFOAM.New developed solver of the package, which is based on the standard one, icoFoam, solving Navier -Stockes equations for homogeneous viscous fluid, was supplemented by new equations for calculating fluid density and salinity, and specific program codes implementing no-flux boundary conditions for salinity perturbation on the impermeable surface of an obstacle.
To interpolate the convective terms in the governing equations a limited TVD-scheme was used, which ensures minimal numerical diffusion and absence of non-physical oscillations of the solution.For discretization of the time derivative a second-order implicit asymmetric three-point scheme with backward differencing was used, which ensures a good time resolution of the physical process.For calculating the diffusion terms, based on the Gauss theorem within orthogonal grid sections, a surface normal gradient was evaluated at a cell face using a second order normal-to-face interpolation of the vector connecting centers of two neighboring cells.In nonorthogonal grid regions, an iterative procedure with a user specified number of cycles was used for nonorthogonal error correction due to a grid skewness.
For solving the resulting system of linear equations conjugate (PCG) and biconjugate (PBiCG) gradient methods were used with DIC and DILU preconditioning for symmetric and asymmetric matrices respectively, which are based on simplified procedures of incomplete Cholesky and LU factorization.For coupling equations for momentum and mass conservation an unsteady well-convergent algorithm PISO (Pressure-Implicit Split-Operator) was used, which works in the most effective way for transient problems.
For fast and efficient regeneration of orthogonal computational grid when changing the geometric parameters of the problem, such as length, thickness, and angle of inclination of an obstacle, number of cells and grading level of the grid in each direction, etc., automated program was developed.The program enables according to a prescribed algorithm recreating the standard reference files used by utilities of the OpenFOAM package for spatial discretization of the problem.The algorithm for constructing a computational grid around the obstacle, oriented at an arbitrary angle to the horizon, consists in creation of two separate computational grids.The inner cylindrical grid around the plate can be ro-tated together with the obstacle when angle of its inclination is changed.The external fixed grid is merged with the rotating part at the interface using the standard utilities of the package.
The computations were terminated when the integral characteristics or their statistical evaluations take values of steady-state regime.The spatial dimensions of the computational cells were chosen from the condition of an adequate resolution of the finest flow components associated with the stratification and diffusion effects, which impose significant restrictions on the minimum spatial step.In highgradient regions of the flow, at least several computational cells must fit the minimal linear scale of the problem.Calculation time step, , t  was defined by where r  is the minimal size of grid cells and v is the lo- cal flow velocity.

Computation results
An impenetrable obstacle immersed into a stratified fluid breaks the natural diffusion transport of a stratifying agent and creates a complex system of fluid flows including basic thin jets along sloping boundaries and secondary compensatory back flows.In this section, numerical results are presented on steady-state structure and dynamics of diffusion-induced flows around an impermeable motionless sloping plate and a wedge with finite length, which are immersed into a stably stratified medium [10,14,15].Such investigations are of a great scientific and practical interest since both the stationary and non-stationary analytical solutions of the problem are constructed only in the approximation of an infinite plane.However, the most interesting features of diffusion-induced flows are located directly in the vicinity of sharp edges of an obstacle, where a complex multi-level vortex system of compensatory fluid motions and extended high-gradient horizontal interfaces are formed, which are clearly observed in the laboratory experiments [16].
Diffusion-induced flow structures around a sloping plate and a wedge have a lot of common structural elements but there are, as well, many distinctive features due to peculiarities of the geometry.
Diffusion-induced flow on a sloping plate.The clearest representation of diffusion-induced flows on a horizontal plate can be given by stream lines, which are shown in Fig. 1.The general structure of the flow consists of a four-level sequence of horizontally oriented circulating cells, which are located symmetrically relative to the plate's horizon and its central vertical plane [14,15].The contacting cells circu-late in opposite directions and are mutually coordinated in contrast to the case of thermal convection where circulations in adjacent cells have the same signs.The diffusion-induced flow around a long plate has the simplest structure with the vertical scales of the circulating cells virtually unchanged along the horizon (Fig. 1, a).The detailed consideration of fine flow structure near the plate's edges reveals that the centers of the basic circulating cells are moved away from the edges in both directions to distances comparable with the typical scale of fine structure components, .
Shortening the length of the plate has practically no influence on far flow field but leads to a significant complication of the primary structure of circulating cells together with a decrease in their intensity and thickness around the barrier.A very small plate's length represents a particular case of infinitesimal geometrical non-uniformity settled in a quiescent, stably stratified fluid (Fig. 1, b).In this extreme case, the compensating flows have the smallest intensity but nevertheless do exist, and from the general flow structure, which looks as if the stream lines of the flow converge to a point when approaching the obstacle.Even a smallest deviation from the horizontal position of the plate leads to flow symmetry breakdown and the formation of new circulation systems including unidirectional jet flows along both the sides of the sloping plate and compensating circulating flows (Fig. 2, a).A single circulation system is formed around the plate surface, where fluid flows down along the lower side of the plate, turns around in the large cell adjoining the right edge, flows up along the upper side of the plate, turns left and closes in the down-slope jet flow.At small slope angles two types of flows co-exist near the obstacle when some fluid remains closed in the vortex flow directly near the plate's edge and another part accomplishes cyclonic circulating motion around the whole surface of the obstacle.
A further increase in slope angle leads to the formation of a uniform distribution of basic jet flows along virtually the whole length of the plate excluding small regions near its edges.The area of blocked fluid decreases and the basic cyclonic circulating cell enlarges due to splitting of the main jet flows into two opposite directions while keeping the cyclonic direction of fluid circulation (Fig. 2, b).
At slope angles greater than 20°, new complicated vortex systems start arising near the plate's edges within the circulating cells directly adjacent to the basic jet flows (Fig. 2, c).The large circulating cells located at some vertical distance from the plate's edges are oriented mostly horizontally within the whole length of the computational domain.
Further increase in slope angle leads to elongation of these vortex structures to a wing-shaped form, and the general flow structure takes virtually a symmetrical form relative to the plate's plane.At the external border of the adjacent circulating cells, additional non-uniformities are revealed, which regularity correlates directly with the value of buoyancy period T b (Fig. 2, d ).
Calculated fields of density gradient  for diffusion-induced flows on a sloping plate, where both large-scale components and fine structure streaks are manifested, are in a good agreement with the highresolution schlieren images of refraction ratio gradient (the density and the refraction ratio of salt solutions are in a linear relation [16]).The calculated pattern and the schlieren image shown in Fig. 3 demonstrate clearly distinguished sharp horizontal interfaces, which are adjoined directly to the plate's edges.These typical horizontal streaks, which are lengthened with increasing sensitivity of a registration method, always exist in stratified media at an arbitrary orientation of an obstacle and are absent in homogeneous fluid.[15].The parameters of stratified fluid are chosen to be typical for the laboratory experiments [16], which significantly exceed typical values of stratification in the natural systems, i.e. atmosphere and hydrosphere of the Earth.
Fields of the basic physical variables of steady diffusion-induced flow on a wedge are shown in Figs. 4 and 5, which are constructed using the technique of polychromic map of isolines.Within such an approach brightness of green and blue colors of the images, which correspond, respectively, to positive and negative values of a field, are varied from light to dark tones between adjacent isolines.Such a technique of field construction allows extracting both qualitative and quantitative information from the patterns and analyzing flow fine structure, as well.Structure of the fields in Fig. 4 contains a number of common features manifested near both the front vertex of the wedge, where slope jet flows are formed, and the back edges, where they separate from the obstacle.Fig. 4, a demonstrates a layered structure of the salinity field due to braking of molecular diffusion flux on the impermeable surface of the obstacle that, in its turn, triggers formation of circulating fluid motions illustrated by Figs. 4, b and c.
Mechanism for formation of the diffusion-induced flow is illustrated by the pattern of salinity perturbation field in Fig. 4, a. Thin layers of salinity deficiency and excess are attached directly to the upper and lower surfaces of the wedge, respectively, that is manifested in form of blue and green streaks in the pattern.The near-wall layers are adjoined by less intensive regions of salinity excess (thin green streak above the obstacle) and deficiency (thin blue streak beneath it), which do not reach the vertices of the wedge.Horizontally extended regions of salinity deficiency and excess are also located in front of the wedge above and beneath its horizon, respectively.An extended area with virtually unperturbed distribution of salinity perturbation is located close to the base of the wedge, which supplements the obstacle to a symmetrical form relative the base ("a reflected wedge with curvilinear boundaries").Near the edges of the base fine structural asymmetrical perturbations are formed, which penetrate into the bottom area behind the wedge like in case of a sloping plate [14,15].
Pattern of horizontal component of velocity field, which is shown in Fig. 4, b, has a symmetrical form relative to the central horizontal plane of the wedge, z  0. Thickness of the ascending and descending jet flows takes values of about 0.32 cm and 0.19 cm near the tip and the base of the wedge respectively.When separating from the surface at the base and spreading further into the surrounding fluid the adjoined secondary jet flows are thickened almost 1.4 times.Thin layers of reverse flows are adjoined to internal boundaries of the free jets, which are manifested the best around the edges of the wedge.
The most complicated flow structure is observed in the field of vertical component of velocity presented in Fig. 4, c.Fluid flows away from the surface of the wedge in a layer of almost uniform thickness of 0.28 cm, but at distance of about 2 cm from the bottom section it is sharply thickened to value of 0.65 cm and in the direct vicinity of the base it is thinned again to value of about 0.3 cm.Almost through the whole length of the wedge the layers of outgoing fluid are adjoined by compensatory reverse flows.Their thickness is the smaller the closer to the area of flow bifurcation where the reverse flows are absent at all and appear only around the base of the wedge.The flows with ascending and descending velocity components interact with each other at the tip of the obstacle that results in a complex structure of vertical component of velocity field.The most complicated pattern is observed in the vicinity of the sharp edges of the base where a typical 'rosette' of dissipative gravity waves is formed [9].Fig. 5 shows patterns of pressure field around a wedge with straight, concave and convex side boundaries.By analyzing this field one can study dynamical effects of the surrounding fluid on the obstacle and, particularly, determine dimensions of flow regions with negative pressure (coloured blue) in front of the obstacle which turn to be proportional to propulsion force resulting in self-movement of a free wedge-shaped obstacle.
Pressure takes negative values in a large expanding area in front of the sharp vertex of a wedge with straight boundaries, two tapering regions being located at some distance from the surface of the wedge and thin streaks adjoining directly to the surface, as well (Fig. 5, a).The regions with negative pressure are generated due to formation of the ascending and descending flows along the sloping boundaries of the wedge as a result of breaking of salinity diffusion flux on the impermeable surface.
Structure of the pressure field around a concave wedge looks very similar to that of a wedge with straight boundaries except for some minor distinctive details (Fig. 5, b).The triangular-like regions and thin streaks over and beneath the surface of the concave wedge are slightly longer and thicker that speaks for more intensive propulsion mechanism in this case.
Thickness of the pressure deficiency zone near the sharp vertex of a convex wedge is about twice smaller than that for a concave one (Fig. 5, c).This result can be explained by deceleration of the ascending and descending flows and, as a result, decrease of propulsion intensity due to wedge opening at the apex when curvature of its lateral boundaries is increased.
With start of motion of an obstacle due to the action of the propulsion mechanisms (that is, rotation of a plate or straight line movement of a wedge) diffusion-induced flow is changed dramatically due to "inclusion" of new flow components such as upstream and attached internal waves, non-stationary vortices in the wake past an obstacle and fine structural interfaces [13].In this case, diffusion-induced flow around an obstacle becomes a complex, multiscale, and transient physical process accompanied by interactions of large and fine scale components be-tween each other and with free stream, which are characterized by their own geometries, spatial and temporal scales, manifestation level, and dissipation rate.

Conclusions
Diffusion induced flow structure and dynamics around a motionless sloping plate and a wedge-shaped obstacle immersed into a stably stratified medium were analyzed numerically based on the fundamental equations set using original Fortran program codes and solvers of own development within the framework of the open source package Open-FOAM.
Diffusion-induced flows always exist around obstacles with an arbitrary geometry as a result of breaking of natural molecular flux of a stratifying agent due to combined influence of medium inhomogeneity and effects of external forces, i.e.Earth's gravitation, magnetic field, global rotation, etc.Even a very weak stratification essentially affects flow fine structure and dynamics, especially near the sharp ed-ges and poles of an impermeable obstacle, where a complex multi-level vortex system of compensatory fluid motions and extended high-gradient horizontal interfaces are formed, which are clearly observed in the both calculations and laboratory experiments.
Spatial density variations may lead not only to the formation of relatively slow compensating convective flows but also to the generation of a propulsion mechanism resulting in a self-movement of neutral buoyancy objects with a special shape in a continuously stratified fluid.
The numerical simulation of diffusion induced flow around a wedge reveals an intensive pressure deficiency zone in front of the leading vertex of the wedge, which is responsible for the propulsion mechanism.Size of the pressure deficiency area near the sharp vertex of a concave wedge is about twice greater than that for a convex one that results in generation of much more intensive propulsion force and a faster self-movement of the wedge with concave boundaries.

Рекомендована
viscosity and salt diffusion coefficients, t is time,  and  are the Hamilton and Laplace operators respectively.
structural flow components caused by the combined dissipation and stratification effects, and typical velocity scales, ,

Fig. 1 .
Fig. 1.Stream lines of the diffusion-induced flow around a horizontal plate of different lengths in a continuously stratified fluid (N  1.26 s 1 ,   0): a, b -L  10; 0.1 cm

Fig. 2 .
Fig. 2. Stream lines of the diffusion-induced flow around a sloping plate in a continuously stratified fluid for different slope angles (L  5 cm, N  1.26 s 1 ): a -d -  1; 10;