Received 18 April 2016; accepted 30 May 2016; published 2 June 2016

1. Introduction

In Computational Fluid Mechanics, recently, a combination of grid method (in Eulerian co-ordinates) with a mesh-free particle method (in Lagrangian co-ordinates) is recommended [1] - [3] or a purely particle method [4] - [6] . However, in particle hydro-dynamics, a lot of artificial and computing time consuming, high frequency oscillations are created, which complicate the computation. So, additional techniques are devised, in order to suppress these high frequency oscillations. We can see the method SPH (Smoothed Particle Hydro-Dynamics). [7] . Whereas, in the MPS method (Moving Particles Semi-Implicit), first developed by S. Koshizuka [8] and, then, used by other Researchers, too [9] [10] some high frequency oscillations are suppressed, depending on time step-length Δt, because of the implicit nature of the method, with no additional technique for suppression.

Recently, the movement of a fluid is treated by incremental static analysis [11] , by which dynamic behavior with strong damping is approximated and the computation is significantly simplified. This incremental static analysis is adopted in present work. Also, recently, point-particles (that is with zero volume or zero area in 2D) are used in Computational Fluid Mechanics [11] , which significantly simplifies the computation, too. These point-particles are also adopted here.

Current work has two motivations: 1) To develop a simple algorithm for incremental static analysis of inter- colliding particles, which approximates dynamic behavior with strong damping and is simpler than particle dynamics, which creates artificial high frequency oscillations which have to be suppressed by an additional technique. We can see SPH [7] . 2) To use point-particles, assisted by a proposed rhombing element assuring incompressibility. These point-particles are much simpler than finite size particles [6] .

2. Proposed Method

In 2D (two dimensional space), a single vessel (Figure 1(a)) or a set of inter-communicating vessels (Figure 1(b)) is considered, containing amounts of liquids, which are simulated by sets of inter-colliding point-particles (that is, with zero area). A boundary-particle function (Figure 1(c)) describes the variation of repulsive force, from a boundary to an adjacent particle, versus distance of particle, perpendicularly to boundary. As shown in Figure 1(c), for,. Also, an inter-particle function (Figure 1(d)) describes the variation of mutual repulsive force versus distance, for every couple of adjacent particles. As shown in Figure 1(d), for,.

Instead of dynamic analysis of particles which creates artificial high frequency oscillations, a simple incremental static analysis is preferred here, which approximates dynamic behavior with strong damping.

Within each step of incremental static analysis, every particle performs small displacements , proportional to out-of-balance forces, , acting on it. Numerical experiments show that, if the liquid is confined within boundaries of a single vessel or a set of inter-communicating vessels, then, the proposed incremental static analysis converges to a final equilibrium state.

The algorithm, based on proposed incremental static analysis, for the movement of a set of inter-colliding point-particles, simulating 2D flow, is briefly shown by the flow-chart of Figure 2, and is described below, in more detail.

3. Description of Algorithm

3.1. Constant Input Data

The boundary conditions are given, that is, position, configuration and dimensions of linear boundaries of single vessel (Figure 1(a)) or set of inter-communicating vessels (Figure 1(b)). The weights w of the particles. The boundary-particle function (Figure 1(c)). Inter-particle function (Figure 1(d)). Also, the step- length Δu of the algorithm is given.

3.2. Initial Conditions

The initial conditions of the particles (x, y co-ordinates in 2D) are given, as shown in Figure 1(a) and Figure 1(b). The particles are initially, sparsely arranged, so that all and. Thus, all boundary-particle and inter-particle repulsive forces are initially zero, ,. And the only out-of-balance forces are, initially, for all particles, those due to their weights, ,.

3.3. Any Step of Algorithm

For all out-of-balance nodal forces, , acting on particles, the maximum absolute value is found, .

Then, every particle performs small displacements, , proportional to its out-of-balance forces, , as follows:

From the new positions (x, y) of the particles, the boundary-particle distances and the inter-particle distances are found. And the boundary-particle and inter-particle repulsive forces are determined, for and for, respectively (see Figure 1(c) and Figure 1(d)).

The out-of balance force, acting on every particle, is found, as shown in Figure 3:

where weight of the particle, repulsive force from adjacent boundary and repulsive force from adjacent particle.

3.4. Output in Present Step of Algorithm

At the end of present step of algorithm, the following can be received, as output data: For every particle, displacements, , present position (co-ordinates x, y) and out-of-balance forces,. For particles adja-

cent to boundaries, with, the external repulsive reactions. For every couple of adjacent particles, with, the mutual internal repulsive reaction. Also, the maximum absolute value of all nodal forces can be received as output of present step.

3.5. End of Algorithm

If maximum absolute nodal force is less than a predetermined lower bound (where final equilibrium state is assumed), or, if the steps of algorithm counter n exceed a predetermined upper bound, the algorithm is interrupted. Otherwise, we go to the next step of the algorithm.

4. Computer Program

Based on the above-described algorithm, for incremental static analysis of a set of inter-colliding point-particles, simulating 2D flow, a simple and short computer program (“pocket” program) has been developed, with totally only about 120 Fortran instructions (50 for main program of incremental static analysis + 50 for subroutine of boundary-particle function + only 20 for subroutine of inter-particle function).

The above program runs with the version Force 2.0 of Fortran, whose Compiler is free available in Google, even in Internet Cafés.

5. Proposed Incompressible Rhombic Element

The particles are initially arranged, in such a way, so that to form elementary rhombs, with particles at their four external nodes, as shown in Figure 1(a) and Figure 1(b). Such an enlarged elementary rhomb is presented in Figure 4(a), whereas, in Figure 4(b), is shown the inter-particle function (mutual repulsive force versus distance, for every couple of adjacent particles, with), where the limit distance is , that is equal to lengths of four external small sides of the rhomb. That is, the four external small sides are activated, from the beginning of the algorithm, whereas the two internal diagonals, along x, y, with lengths, are activated later. For small deformations of the rhomb, as shown in Figure 4(c), the internal diagonals are not activated, so rhombic element is a mechanism, which, however, exhibits a sufficient incompressibility. Indeed, the present area of deformed rhomb mechanism, is, according to Figure 4(a) and Figure 4(c),

where initial area of undeformed rhomb.

For the limit deformation, until which, rhomb is a mechanism, its area is

,

,

that is, the area of elementary rhomb has been reduced by only 8.6%, as shown in Figure 4(d).

After this limit deformation, that is, for, an internal diagonal, the vertical one, in Figure 4(a), is activated, so the rhombic element turns from mechanism to rigid, and further reduction of its area is prevented.

So, the proposed rhombic element, in the initial arrangement of particles, exhibits a satisfactory incompres- sibility, as will also be confirmed in the following applications.

6. Applications

6.1. First Application Single Vessel Coarse Discetization

Input data. In 2D (two dimensional space), a square vessel, with sides 100 cm, is considered, as shown in Figure 5(a), which contains an amount of liquid, simulated by a set of inter-colliding point-particles (that is, with zero area). The initial positions (x, y co-ordinates) of particles are shown in Figure 5(a). The particles are initially

arranged, in such a way, so that to form elementary rhombs, which have an in-compressibility property, as described in previous Section 5. There are 10 rows of 9 particles, alternated with 9 rows of 10 particles, that is totally 2 × 9 × 10 = 180 particles. The weight of every particle is. So, the weight of total amount of liquid is.

In Figure 5(b), an enlarged rhombic element, of the initial arrangement of particles, is shown, with, for the lengths of two long internal diagonals, along x, y, and, for the lengths of four short inclined external sides of rhomb, as well as perpendicular distances 5.0 cm from boundaries to adjacent particles.

The boundary-particle function and the inter-particle function are shown in Figure 5(c) and Figure 5(d), respectively. The parameters of these diagrams have been determined as follows: from nearest initial boundary-particle distance and, from lengths of four external short inclined sides of elementary rhomb, both as shown in Figure 5(b).

The parameter has been determined, after pre-estimation, by trials, of approximate maximum external repulsive reaction, as, so that distance, from boundary to adjacent particle, will be not reduced more than 10%.

On the other hand, the has been determined, in such a way, so that inter-particle stiffness is half of boundary-particle stiffness, which is reasonable.

The step-length of algorithm, mm, has been determined, in such a way, so that to produce, within each step of algorithm, inter-particle reaction increments not more than 1/5 of particle weight.

The particles are initially sparsely arranged, so that all and; thus, initially, boundary-particle and mutual inter-particle repulsive forces are all zero, ,. And the only out-of-balance nodal forces, acting on particles, are those due to particle weight, that is, for all particles, ,.

Output data. The present first application run by the previously described, in Sections 3, 4, simple and short computer program, for incremental static analysis of a set of inter-colliding point-particles, simulating 2D flow. And, in about 2.000 steps of algorithm and only 5.0 sec of computing time, the amount of liquid, simulated by 180 particles (coarse mesh) reached to final equilibrium state, shown in Figure 6(a), by the final positions x, y

of particles, where all out-of-balance nodal forces, in absolute values, , are less than 0.1 N.

It is observed, in Figure 6(a), that the free surface level of the liquid has descended, from the initial value (not yet in equilibrium) to about, in final equilibrium state. In same Figure 6(a),

the external repulsive reactions at boundaries, obtained by coarse discretization, are compared to corresponding forces, obtained from theoretical hydro-static pressure distribution, shown in Figure 6(b). The approximation between computational and theoretical results can be considered as satisfactory. This is mainly due to the proposed incompressible rhombic element, described in previous Section 5.

6.2. Second Application Single Vessel Refined Discetization

Input data. A more refined mesh is tried, for the same previous problem of first application (same amount of liquid, in same 2D vessel). The initial positions of particles are shown in Figure 7(a). There are 19 rows of 18 par-

ticles, alternated with 18 rows of 19 particles, that is totally 2 × 18 × 19 = 684 particles. Weight of a particle is, thus total weight. The particles are again initially arranged, in such a way, so that to form rhombs, which are now smaller than those of first application, 1/2 in lengths and 1/4 in area. An enlarged elementary rhomb, with its dimensions, as well as distances from boundaries, are shown in Figure 7(b). The boundary-particle function and the inter-particle function are shown in Figure 7(c) and Figure 7(d), respectively. The parameters of these functions, as well as the step-length, are determined in the same way, as in first application. The particles are again initially sparsely arranged, so that initially the only nodal forces are those due to particles weight, that is, for every particle, initially, ,.

Output data. The present second application, with single vessel and refined mesh, has also run by the same simple and short computer program, and, in about 8000 steps of algorithm and 1.5 min of computing time, reached to final equilibrium state, where all out-of-balance nodal forces are, in absolute values, , less than 0.1 N. The positions of particles of this state are shown, in Figure 8(a), and it is observed that the free surface level of the liquid has descended from the initial value (not yet in equilibrium) to about in final equilibrium state. In Figure 8(a), are also shown the external repulsive reactions at boundaries, obtained by the refined mesh, in comparison with corresponding forces, obtained from the theoretical hydrostatic pressure distribution on boundaries, in final equilibrium state, as shown in Figure 8(b).

A slightly improved approximation, between computational and theoretical data is achieved by refined mesh, as shown in Figure 8(a). However, in order to simplify input-output and save computing time, for the two following large applications, the coarse mesh is chosen. And, as will be confirmed after the end of applications, in Table 1, significant simplification of input-output and saving of computing time are achieved by use of the coarse discretization.

6.3. Third Application Incomplete Flow from Initial to Adjacent Vessel Not Yet in Equilibrium

After the obtained final equilibrium state of liquid, in the initial single vessel, of first application, now, an opening, with height 20 cm, is formed at the bottom of right vertical boundary, and the liquid is allowed to flow gradually, by incremental static analysis, to an adjacent vessel, communicating with the first one, as shown in Figure 9. In 100,000 steps of algorithm and about 5.0 min of computing time, from beginning of first application, more than half of total amount of the liquid has been transferred to the second vessel, as shown in Figure 9.

However, there is not yet an equilibrium state, as significant out-of-balance nodal forces still remain, much larger, in absolute values, , than 0.1 N.

The external repulsive reactions, at the boundaries, are also noted in Figure 9. These also clearly show that there is not yet a final equilibrium state, as, that is, there exists a significant resultant of horizontal boundary reactions, directed to right.

6.4. Fourth Application Complete Flow from Initial to Adjacent Vessel Final Equilibrium State

In previous third application, the flow, from initial to adjacent vessel, was in-complete and there was not an equilibrium state. Now, the flow, from first to second vessel, continues. In about 300,000 steps of algorithm and 15 min of computing time, from beginning of first application, all out-of-balance nodal forces are, in absolute values, , less than 0.1 N, so the algorithm is interrupted and the state of Figure 10 is obtained. Where

almost whole amount of liquid has been transferred to second vessel, except of few, fourteen, drops-particles, which remain, in equilibrium, at the bottom of first vessel, as shown in Figure 10; this is reasonable and was expected.

The above state is obviously a few steps of algorithm before final equilibrium, in second vessel. However, it is preferred to demonstrate, here, this state, because final equilibrium, in a vessel, has already been demonstrated in two first applications.

The initial vessel, with width 100 cm, is narrower than the adjacent one, which has double width, 200 cm. On the other hand, in final equilibrium state of first vessel (Figure 6(a)), the free surface level of liquid was. So, after subtraction of fourteen drops-particles, remaining at bottom of first vessel, it was expected that the free surface level of liquid, in final equilibrium state, in second vessel, would be . Indeed, the free surface level of liquid, in second vessel, as observed in Figure 10, in almost final equilibrium state (a few steps before equilibrium) is about , which confirms that a satisfactory in-compressibility of the liquid is assured by the proposed in-compressible rhombic element (see Section 5).

If the present fourth application had run with the refined mesh, the input-output would be complicated, with 684 particles and step-length Δu = 0.025 mm, compared to 180 particles and Δu = 0.1 mm of coarse mesh used. And a lot of algorithm steps, 1,200,000, and 4.5 hr of computing time would be required, by refined mesh, much more than corresponding values of 300,000 steps and only 15 min of computing time, of the coarse mesh used.

In Table 1, the computational data, of above four applications, are compared to each other, as regards number of particles, step-length Δu of algorithm, number of algorithm steps required, computing time consumed. The last two rows refer to applications, by refined mesh, not presented in this work.

7. Conclusions

A simple incremental static algorithm is proposed, for 2D flow, simulated by inter-colliding point-particles (with zero area). Within each step of proposed algorithm, every particle performs a small displacement, proportional to the out-of-balance force, acting on it. Numerical experiments show that, if the liquid is confined within boundaries of a set of inter-communicating vessels, the algorithm converges to a final equilibrium state.

The proposed incremental static method approximates dynamic behavior of a liquid with strong damping. So, artificial high frequency oscillations are suppressed. And there is no more need for an additional technique to suppress them (see SPH_Smoothed Particle Hydrodynamics).

A rhombic element, in the initial arrangement of particles, is proposed, with particles at its four external nodes which assure sufficient incompressibility of the liquid.

Based on the proposed simple incremental static algorithm, for a set of inter-colliding point-particles, simulating 2D flow, a simple and short computer program has been developed, with totally only about 120 Fortran instructions.

The above simple and short computer program is first applied to an amount of liquid, contained in a singe vessel. A coarse discretization and a refined one are tried. In final equilibrium state, the hydrostatic pressure distribution, on vessel boundaries, obtained by proposed computational model, is compared to corresponding theoretical data, and is found in satisfactory approximation with them.

At the bottom of a vertical boundary of initial vessel, an opening is formed, and the liquid is allowed to flow gradually, by incremental static analysis, to an adjacent vessel. Almost whole amount of liquid is transferred to second vessel, except of few drops-particles, which remain, in equilibrium, at the bottom of first vessel; this is reasonable and is expected. In the final equilibrium state, in second vessel, when all out-of-balance forces, acting on particles, are in absolute values, less than 0.1 N, the free surface level of liquid, confirms that a satisfactory fluid incompressibility is assured by the proposed rhombic element.

If the large application, of complete flow from initial to adjacent vessel, had run by the refined discretization, input-output would be complicated and excessive computing time would be consumed.