Influence of Simulation Parameters on the Excitable Media Behaviour – The Case of Turbulent Mixing

1. Introduction

Although the turbulence is often associated with fluid dynamics, it is in fact a basic feature for most systems with few or infinity freedom degrees. It can be defined as chaotic behavior of the systems with few freedom degrees and which are far from the thermodynamic equilibrium. In this area two important zones are distinguished:

• The theory of transition from laminar smooth motions to irregular motions, characteristic to turbulence;

• Characteristic studies of completely developed turbulent systems.

In hydrodynamics, the transition problem lays back to the end of last century, at the pioneering works of Osborne Reynolds and Lord Rayleigh. Since the beginnings, it was pointed out the fruitful investigation method of considering the linear stability of basic laminar flow until infinitesimal turbulences. Nonlinearity can act in the sense of stabilizing the flow and therefore the primary state is replaced with another stable motion which is considered as secondary flow. This one can be further replaced with a tertiary stable flow, and so on. It is in fact about a bifurcations sequence, and Couette-Taylor flow is maybe the most widespread example in this sense.

The situation becomes hard to approach if the non-linearity is acting in the sense of increasing the rate of growing the unstable linear modes. Although it was anticipated that the flows can be stable according to the linear theory, in experiments it was concluded that they are unstable. It must be noticed that Reynolds himself understood this possibility, and suggested that for the transition from laminar to turbulence for a pipe flow, “the condition must be of instability at certain size perturbations and stability at smaller perturbations”. The issue of transition in flows such Poiseuille flow until Reynolds numbers under the critical value, must be due to instabilities at finite amplitude perturbations. For these flows there was nothing relevant found concerning the eventual secondary stable motions, moreover it seems that turbulence issue directly from primary flow, at a fixed Reynolds number. These strong turbulence problems are quite difficult to approach and this gives evidence that still is to be added more substance to the original Reynolds’ suggestion, after one hundred years of stability studies. For this quite new theory, unifying with the classic turbulence is far from being solving, but the recent challenges could open new research directions.

Generally, the statistical idea of a flow is represented by a map:

$$x={\Phi }_t\left(X\right),X={\Phi }_{t=0}\left(X\right)$$E1

That means X ismapped in x after a time t. In the continuum mechanics the relation (1) is named flow, and it is a diffeomorphism of class C^k.Moreover, the relation (1) must satisfy the following equation:

$$J=\det \left(D\left({\Phi }_t\left(X\right)\right)\right)=\det \left(\frac{\partial x_i}{\partial X_j}\right)$$E2

where D denotes the derivation with respect to the reference configuration, in this case X. The equation (2) implies two particles, X₁ and X₂, which occupy the same position x at a moment. Non-topological behavior (like break up, for example) is not allowed.

With respect to X there is definedthe basic measure of deformation, thedeformation gradient, F, namely the equation (3):

$$F={\left({\nabla }_X{\Phi }_t\left(X\right)\right)}^T,F_{ij}=\left(\frac{\partial x_i}{\partial X_j}\right)$$E3

where $${\nabla }_X$$ denotes differentiation with respect to X. According to equation (2), F is non singular. The basic measure for the deformation with respect to x is the velocity gradient.

After defining the basic deformation of a material filament and the corresponding relation for the area of an infinitesimal material surface, there can be defined the basic deformation measures: the length deformationλ and surface deformationη, with the following relations (Ottino, 1989):

$$\lambda ={\left(C:MM\right)}^{1/2},\eta =\left(\det F\right)\cdot {\left(C^{-1}:NN\right)}^{1/2}$$E4

where C (=F^T F) is the Cauchy-Green deformation tensor, and the vectors M,N - the orientation versors in length and surface respectively, aredefined by:

$$M=\frac{dX}{\left|dX\right|},N=\frac{dA}{\left|dA\right|}$$E5

Very often, in practice is used the scalar form of (4), namely the following relation:

$${\lambda }^2=C_{ij}\cdot M_i\cdot N_j,{\eta }^2=\left(\det F\right)\cdot C_{ij}{}^{-1}\cdot M_i\cdot N_j$$E6

where D is the deformation tensor, obtained by decomposing the velocity gradient in its symmetric and non-symmetric part.

The flow x=Φ_t(X) has a good mixing if the mean values D(lnλ)/Dt and D(lnη)/Dt are not decreasing to zero, for any initial position P and any initial orientations M and N. As the above two quantities are bounded, the deformation efficiency can be naturally quantified. Thus, there is defined (Otino, 1989) the deformation efficiency in length, e_λ = e_λ (X,M,t) of the material element dX, as follows:

$$e_{\lambda }=\frac{D\left(\ln \lambda \right)/Dt}{{\left(D:D\right)}^{1/2}}\le 1$$E7

Similarly, there is defined the deformation efficiency in surface,, e_η= e_η (X,N,t) of the area element dA: in the case of an isochoric flow (the jacobian equal 1), there exists the equation (8):

$$e_{\eta }=\frac{D\left(\ln \eta \right)/Dt}{{\left(D:D\right)}^{1/2}}\le 1$$E8

The deformation tensor F and the associated tensors C, C^-1, form the fundamental quantities for the analysis of deformation of infinitesimal elements. In most cases the flow $$x={\Phi }_t\left(X\right)$$ is unknown and has to be obtained by integration from the Eulerian velocity field. If this can be done analytically, then F can be obtained by differentiation of the flow with respect to the material coordinates X. The flows of interest belong to two classes: i) flows with a special form of $$\nabla v$$ and ii) flows with a special form of F. The second class is what we are interested for, as it contains the so-called Constant Stretch History Motion –CSHM.

2.Issues on turbulent mixing. Results and challenges

The central problem exhibited in this chapter is the challenge of unifying the theory of turbulent mixing. This implies few levels: analytical, computational and experimental. The specific literature is rich in works both on analytical and experimental models. There were realized few comparative analysis, both for analytical and computational standpoint, for 2D and 3D models of different flow types.

2.1. Recent issues in the literature

The recent discussions and works (Dimotakis, 1983) on the empirical evidence and theoretical ansatz on turbulence support the notion that fully–developed turbulence requires a minimum Reynolds of order of 10⁴ to be sustained. This value must be viewed as a necessary, but not sufficient condition for the flow to be fully developed. Presently available evidence suggests that both the fact that the phenomenon occures and the range of values of Reynolds number where it occures are universal, i.e. independent of the flow geometry.

On the other hand, how sharp this transition is doesappear to depend on the details of the flow. In particular, it is considerable sharp, as a function of Reynolds number, in the (Couette- Taylor) flow between concentric rotative cylinders. It is less well – defined for a shear layer and, among the flows considered, the least well – defined for turbulent jets. Perhaps an explanation for this variation lies in the definition of the Reynolds number itself and the manner in which the various factorsthat enter are specified for each flow. In the case of the Couette-Taylor flow, for example, both the the velocity$$U_{CT}=\Omega \cdot a$$ and the spatial scale$${\delta }_{CT}=b-a$$., where Ω is the differential rotation rate, with a and b - the inner and outer cylinder radii, are well-defined by the flow-boundary conditions.

In the case of a zero streamwise pressure gradient shear layer, the velocity $$U_{S_1}=\Delta U=U_1-U_2$$ is a constant, reasonably well specified by the flow boundary conditions at a particular station. The length scale $${\delta }_{s_1}={\delta }_{s_1}\left(x\right)$$ must be regarded as a stochastic variable in a given flow with a relatively large variance. The Reynolds number for the shear layer is then the product of a well-defined variable and a less well-defined, stochastic, variable.

In the case of a turbulent jet, both the local velocity U_j and the scale δ_jmust be regarded as stochastic flow variables, each with its own large variance. The Reynolds number for the jet is then the product of two stochastic variables, and, as a consequence, its local, instantaneous value is the least well-defined of the three.

As regards fully-developed turbulent flow, the presently available evidence does not support the notionof Reynolds – number – independent mixing dynamics, at least in the case of gas – phase shear layers for which the investigations span a large enough range. In the case of gas – phase turbulent jets, presently available evidence admits a flame length stoichiometric coefficient tending to a Reynolds number – independent behaviour. It mustbe noticed, however, that the range of Reynolds numbers spanned by experiments may not be large enough to provide us with a definitive statement, at least as evidenced by the range required in the case of shear layers.

In comparing shear layer with turbulent - jet mixing behavior, the more important conclusion may be that they appear to respond in the opposite way to Schmidt number effects, i.e. gas- vs liquid – phase behavior. Specifically, there are high Schmidt number (liquid – phase) shear layers that exhibit a low Reynolds number dependence in chemical product formation, if any. In contrast, there are gas- phase turbulent jets that exhibit an almostReynolds – number independent normalized variance of the jet-fluid concentration on the jet axis, with a strong Reynolds – number dependence found in liquid – phase jets, in the same Reynolds – number range.

To summarize, recent data on turbulent mixing suppport the notion that the fully -developed turbulent flow requires a minimum Reynolds number of 10⁴, or a Taylor Reynolds number of Re_T≈10² to be sustained. Conversely, turbulent flow below this Reynolds number cannot be regarded as fully – developed and can be expected to be qualitatively different (Dimotakis 1983).

The manifestation of the transition to this state may depend on the particular flow geometry, e.g. the appearance of streamwise vortices and three – dimensionality in shear layers. Neverless, the fact that such a transition occurs, as well as the approximate Reynolds number where it isexpected, appears to be a universal property of turbulence. It is observed in a wide variety of flows and turbulent flow phenomena. (Dimotakis, 1983)

2.2. Recent results from experimental and computational / analytical standpoint

Recently, it was realized (Ionescu, 2002, 2010) the analysis of the length and surface deformations efficiency for a mathematical 3D model associated to a vortexation phenomena. The biological material used is the aquatic algae Spirulina Platensis.

The mathematical study was done in association with the experiments realized in a special vortexation tube, closed at one end. Locally, there is produced a high intensity annular vorticity zone, which is acting like a tornado. The small scale at which the turbulence issues allows retaining the solid particles, mixing the textile fibbers or breaking up the multi-cellular filaments of aquatic algae.

The special vortex tube used for achieving the breaking up of filaments is a modified version at low pressure of a Ranque-Hilsch tube (Savulescu, 1998). Completely closing an end of the tube, there is obtained in this region a high intensity swirl. The flow in the tube is generally like a swirl, with a rate tangential velocity / axial velocity maxim near the closed end, where is also created the annular vortex structure.

The approaching aerodynamic circuit is made from: a pressure source, the box of tangential inputs, the diffusion zone, the tube where there is produced the swirl with additional pollutants inputs, and the closing end with a rotating end.

It must be noticed that this torrential flow is concentrating and intensifying the vorticity, in contrast with the usual cyclone-type flow or other flows generators. If in the installation is introduced a pollutant, there issues a turbulent mixing in the annular vorticity zone. The spatial and temporal scales revealed the existence of different domains, starting with laboratory ones and until dissipative domains or others - corresponding to fine –structure wave numbers. Thus, the applications area is very large, including collecting, aggregating, separating and fragmenting the various pollutants.

From physical standpoint (Savulescu, 1998), the vortex produced in the installation implies four mechanisms:

• Convection (scalar transport), due to streamlines which are directed towards the ending lid near the tube walls, with the pressure source near the tube centre line;

• Turbulent diffusion due to the pressure and velocity variation;

• Stratification effects due to pressure and temperature gradients;

• Turbulent mixing due to the velocity concentration in annular structures.

The convection keeps in the transport and deposit of the powder pollutant when the graining diameter is greater than 5μm. However, the graining spectra domain under 5μm undergoes a turbulent diffusion and stratification effects which are generally out of control. For revealing the local concentration of the streamlines, 2D simplified models were tested. The multiphase 3D flows simulation is still in study.

The turbulent mixing in multiphase flow reveals the following experiment components:

• The topological limit of the multiphase flow when the pollutant enters in the aggregation state or remains fragmented, in the host fluid;

• The rheological behaviour of non-Newtonian compositions, under macroscopic effects of the air velocity.

The installation was realized in two versions: a small scale vortex tube (10-20mm diameter) and a large scale one (100-300mm). They correspond generally at two particles processing classes, although in most cases the classes can be superposed.

The first category refers to collecting and separating various powders, from gaseous emissions to ceramic powders. The parameters which must be taken into account are the graining spectra, the atmosphere nature and concentration. The vorticity concentration can be used for processing the deposit of different particles in powder or cement form by an adequate closing lid.

The second category contains the processing on small scale, including deformation and breaking up mechanisms for various particles in a host fluid. It is studied the vorticity concentration near the closing lid.

During such an experiment (Ionescu, 2002), it was processed the aquatic algae Spirulina Platensis in the host fluid water. After the processing, the long chains of cellular filaments were fragmented, producing isolated cell units or – sometimes – there was recorded the breaking up of some cellular membranes (with less than 100 Angstrom thickness). The initial and final observations (after the vortexation) are exhibited in figure 1 below. It has been used the non-dimensional parameter$${\tau }_a=\frac{t\cdot Q}{D^3}$$, where t represents the time (in seconds), Q the installation debit (m³/sec) and D the diameter (m³). As it can be seen from the picture, the fragmentation degree starts to increase as τ_a grows. There can also be observed the algae form before and after the fragmentation.

Following the opinion of the specialists in cellular biology, this new technological method of processing the flows is more efficient than the classical centrifugation method.

Figure 1.

The fragmentation degree variation of Spirulina Platensis

Three specific applications were performed as fluid waste management:

• the agglomeration of short fibers (aerodynamic spinning);

• the retention of particles under 5µm without any material filter;

• the breakout of cell membranes of the phytoplankton from polluted waters and the providing of a cell content solution with important bio-stimulating features.

The mathematical modeling and analytically testing of the above experiment confirmed the experimental study. Concerning the phenomenon scale, there were taken into account medium helicoidally streamlines with approximating 10μm width. There were not gone further to molecular level. Also, it is important to notice that for this moment of the analytical testing of the model, there has not been studied the influence of strictly biological parameters (such as pH, the temperature, the humidity, etc).

The complex multiphase flow necessarily implies a theoretical approach for discovering the ways of optimization, developing and control of the installation. Numeric simulation of 3D multiphase flows is currently in study. In the mathematical framework, the flow complexity implies the following three stages:

• modeling the globalswirling streamlines;

• local modeling of the concentrated vorticity structure;

• introducing the elements of chaotic turbulence.

The mathematical model associated to the vortex phenomena is, basically, the 3D version of the widespread isochoric two-dimensional flow (Ottino, 1989):

$$\{\begin{array}{l}\stackrel{.}{x_1}=G\cdot x_2\\ \stackrel{.}{x_2}=K\cdot G\cdot x_1,-1⟨K⟨1,G⟩0\end{array}$$E9

Namely, the following vortex flow model is in study:

$$\{\begin{array}{l}\stackrel{.}{x_1}=G\cdot x_2\\ \stackrel{.}{x_2}=K\cdot G\cdot x_1\\ \stackrel{.}{x_3}=c\end{array},-1⟨K⟨1,c=const.$$E10

In the first stage, the flow model (10) was studied from the analytical standpoint. Namely, the solution of the Cauchy problem associated to (10) was found. In order to analyze the length and surface deformations of algae filaments in certain vortex conditions, the Cauchy-Green deformation tensor was calculated. There were obtained quite complex formulas for $$e_{\lambda },e_{\eta }$$(Ionescu 2002, 2010).

The second stage has been a computational / simulation standpoint. It has been realized a computational analysis of the length and surface deformation efficiency, in some specific vortexation conditions. With the numeric soft MAPLE11, the analysis has been two parts. Firstly, the following Cauchy problems has been solved:

$$\begin{array}{l}{e}_{\lambda }\left(t\right)=0,x\left(0\right)=0\\ e_{\eta }\left(t\right)=0,x\left(0\right)=0\end{array}$$E11

using a discrete numeric calculus procedure (Abell, 2005). The output of the procedure is a listing of the form$$\left(t_i,x\left(t_i\right)\right),i=\mathrm{0..25}$$. In the second part there have been realized discrete time plots, in 25 time units, following the listings obtained in the first part. The plots represent in fact theimage of the length and surface deformations, in the established scale time.

Very few irrational value sets were chosen for the length, respectively surface versors, taking into account the versor condition:

$${\displaystyle \sum M_i=1,{\displaystyle \sum N_j=1}}$$E12

The studied cases are in fact represented by the events associated to different values of length and surface orientation versors: $$\left(M_1,M_2,M_3\right)$$and $$\left(N_1,N_2,N_3\right)$$ for length and surface respectively. The events were very few, about 60. Their statistical interpretation is synthesized in (Ionescu 2010), including also the two-dimensional case.

The graphical events are illustrating the analysis of the deformations for Spirulina Platensis, in 25 time units’ vortexation. According to (Ionescu 2001, Ottino 1989), the algae filaments represent Lagrange markers. The special spiral form of the algae gives the answer to the computational context that the surface deformation e_η is significant, since it brings more cases of filaments’ break up. It was called rare event the event of very sudden breaking up of the algae filaments, this corresponding from mathematical standpoint to the break out of the program (because of impossibility of maintaining the required accuracy). Thus, a very important fact happened: the mathematical simulations matched the experimental analysis.

Sumarising these basic stages in the behavior analysis of the turbulent mixing flow, is important to notice that the filaments breaking up is due to alternating loadings that the filament undergoes, in a space-time with random events, available to the break phenomena. The modeling has been the following basic stages:

• Associating known streamlines to the medium flow (that means helical flow in a cylinder);

• Determining the linear and surface deformations from continuum fluid mechanics;

• Associating a sequence of random values to the (vectorial) length and surface orientations.

In fact, there were found four types of processes, all of which were matched by the simulations

1. Processes with relative linear behavior;

2. Linear-negative processes,whichcorrespond to alternate tasks of stretching and folding of filaments and are the most;

3. Mixing phenomena where there issue smaller or larger deviations or strong discontinuities; these concern the situations when some pieces are suddenly coming off from the whole filament, followed by the restarting of vortexation for the rest of algae;

4. Rare events – these correspond to the turbulent mixing and represent the sudden break up of spirulina filaments.

The validity of the model is confirmed by two aspects:

• the statistic increment in time of the breaking cases, according to the experiments;

• the relative singularity of the events which could produce the breaking, fact which is confirmed by the quite long duration of the experiments which led to the filaments break up; in very rare cases, the cellular membrane could be broken, and the cellular content collected.

Crossing over from 2D to 3D case, it is easy to deduce the requirement of a special analysis of the influence of parameters on the behavior of this complex mixing flow. This would include more paramater analysis types for some perturbation models in 2D and 3D case, but also another mathematical analysis types, for example spectral analysis.

2.3. New influence of simulation parameters

The analysis recently has been continued withmore computational simulations, for 2D model, both in periodic and – non – periodic case, and for 3D model, too. A lot of comparisons between periodic and non- periodic case, 2D and 3D case were realized (Ionescu 2008, 2009). In the same time, the computational appliances were varied. If innitialy, the model was studied from the standpoint of mixing efficiency, in the works that come after, new appliances of the MAPLE11 soft were tested (Abell 2005, Ionescu 2011), in order to collect more statistical datafor the turbulent mixing theory. In what follows, few of these appliances are described.

2.3.1. Few MAPLE11 apliances used in simulations

For the purpose of this chapter, there are presented only the recent used appliances of MAPLE11 soft. There have been used two graphical appliances of this soft, namely „DePlot“ and „Phaseportrait“ tools, from „DETools“ package. Both of them are based on numeric methods in producing trajectories of differential equations systems. The most used numeric method is Fehlberg fourth-fifth order Runge-Kutta method- the so-called „rkf45“ method - with degree four interpolant (Abell 2005).

DETools[DePlot] – this tool plots solutions to a system of differential equations. The calling sequences are as follows:

DEplot(deqns, vars, trange, options)

DEplot(deqns, vars, trange, inits, options)

DEplot(deqns, vars, trange, xrange, yrange, options)

DEplot(deqns, vars, trange, inits, xrange, yrange, options)

DEplot(dproc, vars, trange, number, xrange, yrange, options)

Parameters:

deqns - list or set of first order ordinary differential equations, or a single differential equation of any order;

dproc - a Maple procedure representation for first order ordinary differential equations, or a single differential equation of any order;

vars- dependent variable, or list or set of dependent variables;

trange- range of the independent variable;

number - equation of the form 'number'=integer indicating the number of differential equations when deqns is given as a function (dproc) instead of expressions;

inits - set or list of lists; initial conditions for solution curves;

xrange- range of the first dependent variable;

yrange- range of the second dependent variable;

options - (optional) equations of the form keyword=value;

Description:

• Given a set or list of initial conditions (see below), and a system of first order differential equations or a single higher order differential equation, DEplot plots solution curves, by numerical methods. This means that the initial conditions of the problem must be given in standard form, that is, the function values and all derivatives up to one less than the differential order of the differential equation at the same point.

• A system of two first order differential equations produces a direction field plot, provided the system is determined to be autonomous. In addition, a single first order differential equation produces a direction field (as it can always be mapped to a system of two first order autonomous differential equations). A system is determined to be autonomous when all terms and factors, other than the differential, are free of the independent variable. For systems not meeting these criteria, no direction field is produced (only solution curves are possible in such instances). There can be only one independent variable

• The default method of integration is method=rkf45. Other methods can be specified in the optional equations. Note that because numerical methods are used to generate plots, the output is subject to the characteristics of the numerical method in use. In particular, unusual output may occur when dealing with asymptotes.

• The direction field presented consists of either a grid of arrows or a set of randomly generated arrows. In either case, the arrows are tangential to solution curves. For each grid point, the arrow centered at (x,y) has slope dy/dx. This slope is computed using (dy/dt)/(dx/dt), where these two derivatives are specified in the first argument to DEplot. The curved arrow types (curves and comet) require additional data for the curvature of the direction field, which is computed by moving an epsilon in the direction of the slope dy/dx, and computing dy/dx, then moving an epsilon in the direction opposite the slope, and computing dy/dx. This data is then sufficient to draw a small portion of the direction field lines local to the point, which is then used to draw the curved arrows.

• By default, the two dependent variables are plotted, unless otherwise specified in the scene option.

• The deqns parameter can be given as a procedure, but must conform to the specification as given in dsolve/numeric, and the number option must be included before the initial conditions. In this instance, deqns must be of the form:

proc( N, ivar, Y, YP )

...

YP[1] := f1(ivar,Y);

YP[2] := f2(ivar,Y);

...

end proc

where N represents the number of first order equations, ivar is the independent variable, Y is a vector of length N, and YP is a vector of derivatives which is updated by the procedure (for the equivalent first order system), also of length N.

• The inits parameter must be specified as:

[ [x(t0)=x0,y(t0)=y0], [x(t1)=x1,y(t1)=y1],... ]

[ [y(t0)=y0], [y(t1)=y1],... ]

[ y(t0)=y0, y(t1)=y1,... ]

where, in the above, sets can be used in place of lists, or

[ [t0,x0,y0], [t1,x1,y1],... ]

{ [t0,x0,y0], [t1,x1,y1],... }

[ [t0,x0], [t1,x1],... ]

where the above is a list or set of lists, each sublist specifying one group of initial conditions.

• The xrange and yrange parameters must be specified as follows.

x(t) = x1..x2,y(t) = y1..y2or

x = x1..x2,y = y1..y2

More details about the parameters can be found in (Abell 2005).

DETools[Phaseportrait]. This tool has the following parameters:

deqns- list or set of first order ordinary differential equations, or a single differential equation of any order;

vars- dependent variable, or list or set of dependent variables;

trange- range of the independent variable;

inits- set or list of lists; initial conditions for solution curves;

options - (optional) equations of the form keyword=value;

Description

• Given a list (or set) of initial conditions (see below), and a system of first order differential equations or a single higher order differential equation, phaseportrait plots solution curves, by numerical methods. This means that the initial conditions of the problem must be given in standard form, that is, the function values and all derivatives up to one less than the differential order of the differential equation at the same point.

• A system of two first order differential equations also produces a direction field plot, provided the system is determined to be autonomous. In addition, a single first order differential equation also produces a direction field (as it can always be mapped to a system of two first order autonomous differential equations). For systems not meeting these criteria, no direction field is produced (only solution curves are possible in such instances). There can be ONLY one independent variable.

• All of the properties and options available in phaseportrait are also found in Deplot.

 inits should be specified as

[ [x(t0)=x0,y(t0)=y0], [x(t1)=x1,y(t1)=y1],... ]

where the above is a list (or set) of lists, each sublist specifying one group of initial conditions.

2.3.2. Comparative computational analysis

In what follows there is presented the comparative analysis for 2D and 3D mixing flow models. There are followed two main leves of comparison:

1. The 2D (non-periodic) mixing flow, namely the differential system (9):

$$\{\begin{array}{l}\stackrel{.}{x_1}=G\cdot x_2\\ \stackrel{.}{x_2}=K\cdot G\cdot x_1,-1⟨K⟨1,G⟩0\end{array}$$E13

in comparison with a „perturbed“ version of it, the system (13):

$$\{\begin{array}{l}\stackrel{.}{x_1}=G\cdot x_2+x_1\\ \stackrel{.}{x_2}=K\cdot G\cdot x_1+G\cdot \left(x_2-x_1\right)\end{array}$$E14

1. The 3D mixing flow model (10)

$$\{\begin{array}{l}\stackrel{.}{x_1}=G\cdot x_2\\ \stackrel{.}{x_2}=K\cdot G\cdot x_1\\ \stackrel{.}{x_3}=c\end{array},-1⟨K⟨1,c=const.$$E15

For allmodels it is used the sameset of parameter values (G, KG) as used in (Ionescu 2008), containing three simulation cases:

• case1:

• case2:

• case3:

Each simulation case is labelled on the figure. The time units number was succesively increased, in order to better analyse the solution behavior for each model and corresponding case. The „stepsize“ optionfor the „rkf45“ numeric method, on which the „Phaseportrait“ procedure is based, is implicitely assigned to 0.05. Also, it must be noticed that the above choice of parameters (one positive andanother negative) is optimal for analyzing the direction field associated to the models’ solutions.

The „scene“ parameter of the graphic procedure was set to {x(t), y(t)], the same in 2D and in 3D case.

Figure 2.

Case1 of simulation for the model (9). Time units number: t= 0..40

Figure 3.

Case1 ofsimulation for the model (9) with t=0..100. The movement is realized on the same trajectory

Figure 4.

Case2 of simulation for the model (9). The time units number t=0..40

Figure 5.

Case2 of simulation for the model (9) with t=0..100. The trajectory tends to “aggregate” towards the origin.

Figure 6.

Case3 of simulation for the model (9). Time units number t=0.40

Figure 7.

Case3 of simulation for the model (9) with t=0..100. The trajectory multiplies but does not draw near the origin

Figure 8.

Case1 of simulation for the model (13). The trajectory changes, it seems to be no longer periodic, but tends to infinity

Figure 9.

Case2 of simulation for the model (13). The time units number t=0..40. The trajectory changes again and becomes negative

Figure 10.

Case2 of simulation for the model (13), with t=0..100.

Figure 11.

Case2 of simulation for the model (13) with t=0..200

Figure 12.

Case3 of simulation for the model (13), t=0..40. The trajectory becomes positive

Figure 13.

Case3 of simulation for the model (13) with t=0..100

Figure 14.

Case3 of simulation for the model (13) with t=0..200. The program does not show the field arrows anymore.

Figure 15.

Case1 of simulation for the model (10). Time units number t=0..25

Figure 16.

Case1 of simulation for the model (10) with t=0..50

Figure 17.

Case2 of simulation for the model (10), t=0..25. A periodic behaviour of the trajectory

Figure 18.

Case2 of simulation for the model (10), t=0..75

Figure 19.

Case3 of simulation for the model (10), time t=0..25

Figure 20.

Case3 of simulation for the model (10), time t=0..75. The “circular loops” are fewer than in the case2, for the same simulation time.

Figure 21.

Case3 of simulation for the model (10) with t=0..100.

3. Remarks on the simulations and Maple appliances importance

Looking above to the figures, it is important to notice some special features.

• It must be pointed out the advance of the “phaseportrait” appliance of Maple soft, and generally of the “plots” package. The possibility of increasing step by step thesimulation time units offers important features of the analysis and consequently of the model behaviour. For example in fig.3, when time increases, the trajectory “moves” on the same line. Comparing to this, in fig.4 the trajectory has some “tracks” which overcrowd towards the origin. In case 3, for the same simulation time, the trajectory only multiplies and does not came toward the origin. Also, figures 8-10 illustrate in a good manner the influence of time increasing on the trajectory behaviour: this can be seen from the values variation of x(t) and y(t) on the axes.

• All the simulations for 2D case model have been realized in the same initial conditions: (x(t), y(t)) = (1,1). This gives more accuracy for the graphical comparisons. Applying this procedure step by step gives the possibility to have a total parameter control at every stage of simulation. In 3D case the initial conditions were a little modified: (x(t, y(t), z(t))=(1,0,0). It was taken into account the fact that on the z axis is represented the velocity and is natural that the mixing process starts from a zero velocity

• Also, the “scene” parameter was the same: [x(t), y(t)] both in 2D and 3D case. That means in 3D case, the trajectory is studied by watching it “from the z axis”. The field arrows are an important appliance of this plot tool, as they are very suggestive concerning the trajectory features. In this context it is extremely important to notice the figure 13. In this model (simulation case3), when increasing the time, the field arrows disappear, the program doesn’t show them anymore. This can be glossed like a “lose of equilibrium” for the 2D model, in certain parameter context, and this fact can be easily observed from the great values for x(t) and y(t) on the axis

• In figures 8-13 it must be noticed a similarity with the pictures obtained in (Ionescu, 2008), where the origin becomes an unstable focus. It is observed here that when increasing the time, e.g. in Fig. 10, the trajectory does not go periodic, it changes only the positivity.

• A special attention must be paid for the figure 17, 3D simulation case. For t=0..75 units, and looking to the trajectory from above, from the z axis,it is obvious that the vorticity of the flow appears. It is very important to notice that if in (Ionescu 2002, 2010), when searching for the mixing efficiency, there were necessary small time units for get the special events of losing the equilibrium of the model, in this case, searching for the phaseportrait needs some more time units, in order to get visible vorticity of the flow. The vorticity allure is also visible in the figures 18-20, with some changes in the loops number and their nearness to the origin.

4. Conclusion

Looking on the above computational analysis, some conclusions issue: regarding the above comparative analysis, on one hand, and, concerning the general context of random events that go with the mixing flow phenomena, on the other hand.

Looking at the above pictures, we see that each simulation case brings significant differences:

• For the modified (perturbed) 2D model, just from the case1 of simulation we see that the trajectory has a great change. It is not periodic anymore like in the initial 2D model, but is more evident its trend to infinity, just from the beginning of the simulation. The parameters have a great influence in this sense

• Concerning the 3D model, it is obvious that the field arrows could not appear, since there is a 3D simulation. Instead, the “scene” parameter was chosen to [x(t), y(t)], like in 2D case model. This is optimal for studying the trajectory behaviour, since it is like we “look” at the trajectory from above.The scene parameter is important in analysis, and a further aim would be to change it in the simulations and make some other comparative graphical analysis

Concerning the general context of random events that go with the mixing flow phenomena, the above analysis just brings at that some important conclusions:

• It is obvious that the above analysis confirms the fact that the mixing flow model is a far from equilibrium model. This is confirmed by the number of parameter / time units analysis. Using different appliances – “mixing efficiency” and computational / graphical appliances, the conclusion is the same: the model becomes far from equilibrium, in certain simulation conditions. It must be noticed that the “units” of time can be of any type, sufficiently small or large. It depends to the target of the analysis. Enlarging the set of parameter values would bring another important data for completing the panel of repetitive events. This would be a next aim. All this space-time context present in all these analysis, consolidates the basic statement that the turbulent mixing flows must be approached as chaotic systems. This is in fact regaining the idea of a system / model high sensitive to initial conditions (Ionescu 2009).

• Testing step by step each of the three models above was possible due to the flexible structure of the MAPLE11 graphical / computational appliances. In the same time this shows that these repetitive simulations are relatively easy to perform. Maybe useful comments would produce if sharing some of these files.

• Also the “step by step” analysis guides the reader by successive steps, to the so- called “far from equilibrium model”, by the possibility of sensitive modifying the model parameters and the time at each step of simulation.

• The issues of repetitive phenomena, both in 2D and 3D case,give rises to achieve some appliances of chaotic dynamical systems, whose numeric models would give new research directions on the behavior in excitable media. This would be a next aim

Acknowledgement

The basic part of my work in turbulent mixing was done in the period of my Ph. D. thesis at the Polytechnic University of Bucarest Romania. Fruitful experiments and special discussions were done togheter with the scientific team of the Institute of Applied Ecology, Bucarest. I am grateful to them.

But my entire acknowledgement points to Ph. D. Prof. Eng. Stefan Savulescu, from Politechnic University of Bucarest. Without the special vortexation installation constructed by him, the mixing flow experiments could not be possible. Due to his great experience in the turbulence field, and his fundamental observations, advices and suggestions, my work in turbulent mixing have an important sense and the great challenge of unifying the turbulence theory takes shape.