## Abstract

We investigate how a plug of obstacles inside a two-dimensional channel affects the drainage of high viscous fluid (oil) when the channel is invaded by a less viscous fluid (water). The plug consists of an Apollonian packing with, at most, 17 circles of different sizes, which is intended to model an inhomogeneous porous region. The work aims to quantify the amount of retained oil in the region where the flow is influenced by the packing. The investigation, carried out with the help of the computational fluid dynamics package *ANSYS-FLUENT*, is based on the integration of the complete set of equations of motion. The study considers the effect of both the injection speed and the number and size of obstacles, which directly affects the porosity of the system. The results indicate a complex dependence in the fraction of retained oil on the velocity and geometric parameters. The regions where the oil remains trapped is very sensitive to the number of circles and their size, which influence in different ways the porosity of the system. Nevertheless, at low values of Reynolds and capillary numbers *Re*<4 and *n*_{c}≃10^{−5}, the overall expected result that the volume fraction of oil retained decreases with increasing porosity is recovered. A direct relationship between the injection speed and the fraction of oil is also obtained.

## 1. Introduction

Flows of bi-phase fluids in porous media occur in a large number of industrial applications and in natural systems. Among them, we may indicate boiling and pressurized water reactors, phase separators used in oil production, geo-meteorological phenomena, including sedimentation and soil erosion, blood flow [1–4].

As well known, the phenomenological concept ‘porous medium’ encompasses a multitude of physical conditions that share a common property: the presence of solid matter within the region where the flow is observed. The medium can be in the form of a single plug, inside which many narrow channels allow for the fluid flow, or in the form of a large number of solid disconnected obstacles, which also lead to the presence of very narrow channels as compared to the linear dimensions of region [5]. The first situation often occurs in the subsurface, where large pressures compact geological layers to form porous sedimentary rocks, which form the preferential storage system for petroleum, gas and water. The second situation can be thought of as sand and gravel systems on river bed, but also on subsurface systems where the sedimentary rocks have been fractured either by natural or anthropogenic factors [6,7].

Both situations were first described in a uniform way, by means of the permeability *k*, a single phenomenological parameter characterizing the global effect of the solid matter that forces the fluid to flow along a much more complex pattern, with much higher energy dissipation when compared with the corresponding obstacle-free flow. Here *k* depends on several geometrical and physical properties of the medium, like the porosity *ϕ* [3,4]. The investigation of the dependence of *k* on the actual distribution of solid matter largely benefited from computational fluid dynamics (CFD) methods [8,9], i.e. from the direct integration of flow equations. Another strategy in the study of fluid motion that also received much attention consists in replacing the actual medium by cellular models [10]. Macroscopic quantities are treated locally, by assuming that the displacement of the fluid among neighbouring cells must obey the proper conservation laws. In the case of porous media, this amounts to consider a discrete assembly of interconnected channels, and assuming that within each channel the flow is obstacle free [11,12].

Both possibilities, initially proposed to single-phase problems, can be conveniently extended to consider multi-phase flows. Here, the scope of the problem is enhanced, and questions beyond understanding the details of resistance of the medium to the flow arise. The most direct one is to describe how the obstacles change the regions that different fluids occupy.

In this work, we consider the flow of a two-phase Newtonian fluid (e.g. water and oil) in a two-dimensional channel that accommodates a number of circular obstacles (discs), geometrically displaced according to the Apollonian packing [13–16]. This particular assembly, which has been used in previous investigations [17,18] about the dependence of *k* with respect to packing composed of obstacles of different sizes, aims to describe the main properties of a fractured medium with grains of different sizes.

For such packing, *ϕ* and *k* depend on the number and size of discs [17,18], a fact that influences the results discussed herein. We would like to call attention that usual porous medium studies consider that the whole region where the fluid flows is somewhat homogeneous. This differs from the described model set-up, which is characterized by the presence of obstacles in a limited region of the channel, the rest being obstacle free. This limitation is due to our intention to concentrate the study to understanding the effect of a number of obstacles in retaining some fraction of a resident (defending) viscous fluid (oil) when the region is invaded by a less viscous fluid.

The main goal of this study is to understand how the presence of the packing limits the access of each phase to different regions of inside the channel. Since we consider, as initial conditions, that oil occupies the space where the packing is placed, we obtain information about the fraction of oil retained inside the channel due to the presence of the packing. Our investigation explored the dependence of the retained oil with respect to the input velocity and the two geometric aspects (size and number of discs) defining the porosity of the packing. Since no analytic solution of the Navier–Stokes for this problem is available, it can only be addressed within CFD or other numerical approach as, for instance, cellular automata and lattice Boltzmann methods.

## 2. Model

We consider a two-dimensional channel in the plane *XY* , limited by two parallel lines placed at *y*=±*h*/2, where the fluid flow takes places under a pressure gradient applied along the *X*-direction. In its interior, a set of discs, grouped according to an Apollonian packing, models a fractured porous region, composed of isolated obstacles. This packing was introduced by Apollonius of Perga in the third century BC as the solution to the problem of determining the minimum number of non-overlapping circles necessary to cover a region in the plane [19]. In our model, the packing is constructed according to the following rule. Initially, at the first packing generation *G*=1, four tangent discs of the same size, with radius *R*=*h*/4, are placed together within the channel, in such a way that their centres are located on the vertices of a square of side *h*/2. The second generation *G*=2 is obtained by adding three more discs: one in the centre of the *G*=1 packing, which touches the four *G*=1 discs, and two discs in the space limited by two discs and one of the walls, which are considered to be circles with infinite radii placed at . Higher order generations *G*>2 are obtained by adding up new tangent circles in all spaces limited by any three discs in the previous generation. The number of discs in the packing at generation *G* is indicated by *N*_{G}. We recall that, given any set of three mutually tangent circles with known radii and centre locations, it is possible to determine the radii and the coordinates of the centres of a fourth circle that is tangent to the three previous ones [19].

The touching discs in the described configuration refrain fluid movement under the pressure gradient for any value of *G*. To allow for fluid displacement, it becomes necessary to open channels between the discs, which is achieved by uniformly reducing their radii *R*_{i} by a factor *s*<1. In this work, we consider *s*=0.7, 0.8 and 0.9, and *G*=1, 2 and 3. In this way, the packing porosity, which takes account of the space between the discs within a square of area *h*^{2}, depends both on *s* and on *G*, i.e. *ϕ*=*ϕ*(*s*,*G*). For the purpose of finding the numerical solution within the CFD procedure, the otherwise infinite channel is considered to have a finite length *L*.

The model set-up is illustrated in figure 1, for *G*=1 (*a*) and *G*=3 (*b*), and *s*=0.9. For the purpose of keeping track of the presence of the defending (oil) and invading (water) fluids in different regions of the channel, it is divided into four zones *Z*_{i} *i*=1 4, indicated by different colours (grey tones). With this division, it is possible to use a facility in the CFD application to choose the area in which we want to calculate the volume of each of the phases. At the initial integration time (*t*=*t*_{0}=0), the invading fluid (water) occupies zone *Z*_{1}. The void space of zone *Z*_{2}, as well as zones *Z*_{3} and *Z*_{4} are fully filled with oil.

This model has already been used to investigate the flow of single-phase fluid under different boundary conditions. Although the particular choice of disc sizes and geometrical placement can hardly be found in natural systems, the use of the Apollonian packing to model porous geometry has some advantages, which result from this specific packing geometry. Indeed, among many other properties, it is observed that the disc sizes decay exponentially as their number in the packing grow, the distribution of discs with respect to their number of neighbours in the packing decays as a power law, the network obtained connecting the centres of touching discs, which reproduces the structure of contact forces, is topologically invariant. When it is used in any numerical modelling, the packing is limited to a finite number of sites but, due to its deterministic nature, it is easy to improve the quality of the results, just by increasing the value of *G* (and number of discs) in a controlled way. As we mentioned in the Introduction, the geometrical properties of the packing provide an interesting way to analyse the properties of flows in a region occupied by objects of different sizes, for instance, a fractured rock. The model can also be changed by randomly disturbing the geometrical setting. This amounts to changing the positions and/or sizes of discs, as well as individually selecting the reduction factor for each object. By controlling amplitude of such changes, it is possible gradually approach to more realistic description of porous media formed by isolated objects [20]. Finally, the number of discs can also be increased by adding more than one set of discs defined by the Apollonian packing, which can be subjected to random disturbances as discussed earlier. Of course, the investigation of all possible generalizations is far beyond the scope of this work. Therefore, we concentrate on the results for the pure deterministic packing. In such case, it is not necessary to take averages over discs distributions. Nevertheless, we have performed several independent integrations for the same geometry, as the CFD results may be susceptible to uncontrolled rounding-off error propagations.

The geometrical set-up considered the channel width *h*=4.0 mm, while the total length *L*=33 mm is assigned to the four zones according to the following scheme: *Z*_{1}:1 mm; *Z*_{2}:4 mm (packing zone); *Z*_{3}:8 mm; *Z*_{4}:20 mm. By properly re-scaling the linear dimensions, velocities and fluid features of the fluid used in this study, our results can also be applied to similar systems of different sizes. The fact that the permeability of similar systems of different size does not differ much has been pointed out as an advantage of numerical methods in other approaches to porous media problems [20,21].

## 3. Methods

The equations of mass and momentum conservation are integrated with the help of *ANSYS-FLUENT* CFD package, under the assumption of incompressible and isothermal regime [22]. The equation of conservation of mass, or the continuity equation, is written as
3.1This equation is valid for both compressible fluids as for incompressible fluids, and the term *S*_{m} is the flow mass added to the continuous phase from a dispersed phase. The incompressibility assumption requires that *ρ* does not depend on time. Thus, if *S*_{m}=0, the continuity equation reduces to
3.2The equation for momentum conservation, in an inertial system of reference, is written as
3.3where *p* is the pressure, is the stress tensor, *ρ*** g** is the gravitational force on the fluid mass and

**represents any external force acts in the fluid.**

*F*The stress tensor is given by
3.4where *μ* is the molecular viscosity, *I* is the unit tensor and the second term on the r.h.s. describes the effect of volume expansion. Thus, ignoring external forces, and using the non-slip condition and the continuity equation in equations (3.3) and (3.4), the momentum conservation equation is reduced to the Navier–Stokes equation
3.5

The numerical integration takes into account the usual non-slip boundary condition (** v**=0) between all solid–fluid interfaces. The CFD framework requires the definition of a mesh, which results from a division of the domain into discrete, non-overlapping cells [9]. The differential equations expressing mass and momentum conservation are transformed into algebraic equations for the values of the pertinent fields on the faces of the cells. The resulting values for each of the two fluid phases are interpolated to the centre of cells. For a model formed by circles, as the Apollonian packing, it is advisable to use a discretization based on triangular cells, which allow for a best fit of the curved walls. The mesh generator used to discretize the domain ignores the symmetry of the packing geometry, resulting in a triangulation, which does not obey mirror symmetry with respect to the centre line. In this study, the end effect is that the mesh does not exactly reproduce the reflection symmetry with respect to the horizontal line at centre of the channel. In the next section, it is possible that asymmetric oil retention patterns result as stationary state at the end of invading process. Therefore, a completely symmetric mesh would be relevant in this study in order to discern if any pattern asymmetry results from discretization or from the intrinsic instability of a more symmetric solution. Nevertheless, we were able to greatly circumvent the mesh asymmetry by imposing the same mesh boundaries at all fixed walls, which significantly reduced mesh distortions.

In the numerical solution of equations (3.2) and (3.5) for two-phase systems, the Courant–Friedrichs–Lewy (CFL) condition is taken in account to analyse the stability of the solution [9] and controlling its accuracy. For two-dimensional systems, the CFL condition is expressed by
3.6where *C* is the (dimensionless) Courant number, Δ*t* is the time step, *v*_{x} (*v*_{y}) is the velocity along the *x* (*y*)-direction, and Δ*x* (Δ*y*) is the displacement of the fluid phase travelled within the cell along *x* (*y*)-direction. The CFL condition relates the typical time interval required for the fluid parcel to traverse any given cell. If the CFL condition is not satisfied, errors associated with un-physical deformation of the two fluid interface may propagate, so that the numerical solution becomes meaningless. For two-phase systems, condition (3.6) must be satisfied for the velocities of both fluids, which leads to difficulties in adjusting space discretization and time step in order to optimize the performance of the package.

Indeed, when *ANSYS-FLUENT* performs a time-dependent multi-phase fluid flow calculation [22], the time step used for the volume fraction calculation will not be the same as the time step used for the rest of the transport equations. In the region near the fluid interface, the software divides the area of each cell by the sum of the outgoing fluxes. The result represents the time it would take for the fluid to empty out that cell (with known size). The characteristic transit time is taken as the minimum value of this time taken over all the cells. It is used to verify whether condition (3.6) holds.

## 4. Results

The results we obtained were limited to the low Reynolds number regime (*Re*<10), where
4.1For this range of *Re* values, the properties of the flow also depend on the capillary number *n*_{c}, defined by
4.2In equations (4.1) and (4.2), *h* indicates the channel height, the constant inlet speed of the invading fluid, *ρ* and *μ* represent the invading fluid density and viscosity, respectively, while *σ* is the surface tension between the two fluids [5,23]. In this investigation, the capillary number, which expresses the relationship between the viscous force and surface tension between the two fluids, takes values in the interval 3.0×10^{−6}<*n*_{c}<3.5×10^{−5}. So it covers not only a range conditions where capillary forces dominate viscous forces *n*_{c}≅10^{−5}, but also the opposite situation in which viscous forces dominate the flow.

The effect of the packing on the drainage of oil from the channel by the less viscous fluid depends on the values of the following parameters: contact angle, inlet speed, channel size, fluid viscosities and surface tension [23]. For each of the quoted values of reduction factor *s* and generation *G*, we carried on numerical integration of the system for five different values of inlet speed indicated in table 1. During this procedure, the values of all other parameters were kept fixed.

For the invading fluid, we used the properties of water at 20°C, which corresponds to *ρ*_{inv}=998.2 kg m^{−3} and *μ*_{inv}=0.001003 kg ms^{−1}, while the defending fluid is characterized by *ρ*_{def}=922 kg m^{−3} and *μ*_{def}=0.05 kg ms^{−1}. This corresponds to oil with typical density (in respect to that of water) expressed by American Petroleum Institute (API) index 22. Further parameter values are *σ*=0.0312 N m^{−1}, contact angle between fluids *θ*=*π*/2, and time integration step *τ*=2.5×10^{−4} s.

The choice of the inlet speed range was based on the analysis of the flow of the same two-phase system in a channel without any obstacle. We selected values of such that, in the packing-free condition, all oil could be expelled by the invading water. In this way, we can be sure that any oil retention was caused only by the presence of the packing. The corresponding value of *Re* were low enough to be in the laminar (table 1). There is an intermediate range of values of the inlet velocity in which the regime is still laminar but the invading fluid is not able to completely expel the oil from the channel. This regime is characterized by a particular fingering pattern, in which a single finger is opened connecting inlet to outlet, through which water flows. In the stationary regime, a certain amount of oil remains adhered to the channel walls [23].

We measure the effectiveness of the oil drainage by the quantities *B*_{2} and *B*_{2,3}. Let *O*_{i}, *i*=2,3, denote the amount of oil retained in zones *Z*_{2} and *Z*_{3}. The amount of oil initially placed in the zone *Z*_{2} corresponds to , where *r*_{i}(*s*) indicates the radius of each disc in the packing with reduction factor *s*. Then let us define
4.3so that *B*_{2} informs how much oil was expelled from *Z*_{2} by the invading fluid, and *B*_{2,3} measures the total amount of oil that could not be recovered at the end of the process.

### (a) Process of oil retention

When water is injected at a constant speed, it first pressures the amount of water initially placed in *Z*_{1} against oil in *Z*_{2}, forcing it into *Z*_{3}. After some time interval, water also reaches *Z*_{3} but leaves behind a certain amount of oil between the discs in the zone *Z*_{2}. Depending on the features of the two fluids and on the value of , the channels opened by the invading fluid (consequently, the regions where oil remains stored) may have different patterns. The fraction of retained oil depends essentially on the resulting channel pattern. Variations on the retained fraction (e.g. when is increased) are much smaller when the channel pattern is not changed when compared with the opposite condition, in which a new channel pattern arises. In all investigated situations, water flow through the central channel, horizontally placed between the large *G*=1 discs, has been observed. The central channel is the preferential path taken by the invading fluid.

Water reaches *Z*_{3} in three different patterns: (i) exclusively through the central channel; (ii) through the central channel and one of the two channels formed by a large circle and the wall; and (iii) through all three channels. In this last case, the flow through *Z*_{3} and *Z*_{4} is much like that observed in the packing-free channel at the same value of . The invading fluid completely expels oil from these regions. The amount of retained oil (if any) is trapped within the channels between the discs and, eventually, a small drop attached to the right part of the discs in the downstream direction.

In case (i), the invading fluid flowing through the central channel pressures oil towards the walls. For the chosen values of , the invading fluid causes the progressive decrease of the width of the oil layers adhered to the walls, until the rupture of the resulting thin films occurs. This happens at a certain distances (not necessarily the same for the two walls) from the point marking the beginning of *Z*_{3}. At the moment this occurs, the large surface tension causes the oil left behind to flow back in the −*x*-direction. An amount of oil remains attached to the large circle, part in zone *Z*_{2}, part in *Z*_{3}. Similarly, the region occupied by the downstream oil is also reshaped towards the usual interface pattern of the packing-free invasion process. Finally, the dynamics in case (ii) follows a combination of the processes observed in (i) and (iii).

The different scenarios described above depend on the number and size of discs through changes in the total area of the walls of the channel, which impact on the influence of the viscous forces. It is possible to identify a trend towards the enhancement in the amount of oil expelled with decreasing speed, which results from changes between capillary forces and viscous forces. The changes in oil retention profile are also related to accumulation in the downstream region after the packing, as can be seen in figures 2 and 3.

The qualitative analysis of the flow by the direct comparison of resulting patterns indicates that the presence of new circles causes the water to penetrate the porous region, as can be seen by the comparison between figures 2 and 3.

These illustrations indicate that, in a laminar flow with low *Re*, despite the larger energy dissipation, owing to the increase in the fluid–wall interface caused the increase of the number of discs, geometric aspects of the porous region can facilitate the invasion of water into thinner channels. However, this is not a guarantee of an enhancement of oil removal, because oil can be trapped in other regions. These can be better understood by combining the observation of trapping patterns with quantitative measures expressed by the parameters *B*_{2} and *B*_{2,3}.

### (b) Analysis of the volume fraction of retained oil

Figure 4 illustrates the dependence of *B*_{2} and *B*_{2,3} with respect to for *G*=1,2, and 3 and *s*=0.7,0.8 and 0.9. It shows that *B*_{2} remains almost constant for different values of when *G*=1. This is consistent with the illustrations in figure 3, showing that the invading fluid flows predominantly through the central channel, while the four secondary channels between circles and walls remain blocked. The only exception observed for the smallest value of when *s*=0.7 results from opening one channel between one of the left discs and the wall. On the other hand, for the same values of *s* and , the dependence of *B*_{2,3} is subject to larger variations. This results from the different amounts oil that remain attached to the right of the discs, most of which occupy regions in zone *Z*_{3}.

This picture changes for the higher generations, mainly because of the presence of the disc with the second largest radius in the central part of the packing. This disc alters the flow pattern in different ways for different values of *s* and . The dependence of *B*_{2} as function of is characterized by large variations for all values of *s*, some of which with opposite tendencies. Fluctuations in the behaviour of *B*_{2,3} are enhanced, since it results from oil stacked within the small channels as well as in the region behind the large circles zone *Z*_{3}.

Despite the large fluctuations in the values of *B*_{2} and *B*_{2,3} with respect to , an overall general tendency is the increase in the fraction of retained oil as the value of *s* increases. Although this aspect becomes more evident for *B*_{2,3}, it can also be detected in the case of *B*_{2}. This corresponds to the reported behaviour in other previous studies and the overall known resistance to remove a fluid from dense packing and media with small porosity. The larger differences between *B*_{2} and *B*_{2,3} when the value of *s* increases is explained by the fact that the amount of oil *O*_{3} is less sensitive to variations in the value of *s* than *O*_{2} and *A*_{2}.

It is important to note that the two different measures can guide expected results for other packing arrangements that can be built in similar way. If the whole channel is filled with discs, one should expect the value of *B*_{2} as a good approximation. However, for inhomogeneous packing, where regions with obstacles and obstacle-free alternate, *B*_{2,3} should provide a better approximation for the fraction of retained oil. This is useful when one compares the dependency of these measures as function of as well as for the analysis of the dependency with respect to the porosity, that is carried out in the following subsection.

### (c) Porosity

The results for the fraction of retained oil can also be analysed as a function of the packing porosity, i.e. as function of *ϕ*=*A*_{2}/*h*^{2}. This is illustrated in figure 5, where we show graphs for *B*_{2} and *B*_{2,3} as a function of *ϕ* for each of the five different values of . Once *ϕ*=*ϕ*(*s*,*G*), each panel comprise nine points. For the sake of a better visualization, we identify the points corresponding to the same value of *G* by different symbols and/or colours, connected by straight lines, so that the variation of *ϕ* results from different values of *s*.

The graphs for *B*_{2,3} show a decrease in the fraction of retained oil with increasing porosity, that is, it becomes easier to expel the oil from the packing for larger values of *ϕ*, when the channels between the discs and between the discs and the walls are wider.

This dependency becomes weaker in the case we compare the behaviour of *B*_{2}. For instance, when *G*=1, the fraction of oil retained within the packing region can be regarded as rather insensitive to the value of *s* for almost all values of . For *G*=2 and 3, the tendency of a decrease in *B*_{2} with respect to *ϕ* is again recognized, but the range of retained values is clearly smaller.

As already discussed, the changes in the trend of the reduction of retained oil with respect to *s* are related to different accumulation patterns indicated in figures 2 and 3. This occurs more often at low values of . The graphs also make evident that, at *G*=1 there is less variation, what can be explained by the dominating role of the pattern where only the central channel is opened by the invading fluid.

## 5. Conclusion

This work discussed the numerical investigation of a two-phase flow. The main focus was to obtain a quantitative relationship between the amount of oil retained within a two-dimensional channel with a packing of discs of different sizes in its interior when the system is invaded by a less viscous fluid. Although the flow differs from that occurring in a Hele–Shaw set-up by neglecting the viscosity effect between fluids and limiting plates, viscous fingering phenomena are also present [24]. The packing is constructed by a geometric recurrence rule, and the flow proceeds at low speed to warrant a laminar flow. The regime of low velocity was also reduced to warrant that, at all invading velocities used in the investigation, no oil would be retained inside the channel in the packing-free condition.

The total amount of retained oil depends on the invading fluid velocity as well as on several geometrical features, including the number, position and size of the discs. No simple relation could be obtained between this amount and any single specific geometric feature of the packing. Global tendencies, however, can be detected. They are consistent with empirical findings that extraction is difficult when porosity is small and the invading velocity increases.

The quantitative results were expressed in terms of the fraction of oil retained in the region where the packing was contained, as well as in terms of all remaining oil. This makes it possible to extrapolate the results to situations in which a larger number of obstacles is uniformly placed within the channel, or when they follow a sequence of regions with and without obstacles.

Different effects follow the increase of the number of discs inside the channel. On one side, it creates barriers within the porous region and increases the area of contact of fluid with the walls. At the same time, changes in the flow pattern may enhance the efficiency of the invading fluid in expelling the resident one.

Preliminary results for more complex configurations, with two or more Apollonian plugs placed inside the channel, confirm that the devised framework can lead to similar results for more realistic patterns (see §2), which will be discussed in a forthcoming contribution. The computational efforts to investigate many distinct configurations are much larger, as we have to consider all difficulties that this study made evident for the simulation of two-phase flows within CFD. A good mesh discretization becomes crucial due to the necessity to account for the fractions of two fluids in the same cell. The control of Courant number, which is not relevant to the single-phase flow in the same packing, becomes an important issue. For a fixed mesh, this control refrains the range of velocities of the invading fluid for which the condition expressed by equation (3.6) is not violated. Finer meshes would be required to reach larger velocities, which constitute a limitation to the investigation as compared to the single-phase system.

## Authors' contributions

A.S.S.P. conceived the study, prepared the CFD framework, carried out numerical work and result analyses, and helped draft the manuscript; R.S.O. conceived the study and prepared the CFD framework; R.F.S.A. conceived the study, coordinated the study and wrote the manuscript. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

The authors acknowledge the financial support of Brazilian agencies CNPq, CAPES, FINEP, and FAPESB. R.F.S. Andrade benefit from the support of the Instituto Nacional de Ciência e Tecnologia para Sistemas Complexos (INCT-SC).

## Footnotes

One contribution of 13 to a theme issue ‘Topics on non-equilibrium statistical mechanics and nonlinear physics (II)’.

- Accepted July 27, 2015.

- © 2015 The Author(s)