Incremental Static Analysis of 2D Flow by Inter-Colliding Point-Particles and Use of Incompressible Rhombic Element

A simple method is proposed, for incremental static analysis of a set of inter-colliding particles, simulating 2D flow. Within each step of proposed algorithm, the particles perform small displacements, proportional to the out-of-balance forces, acting on them. Numerical experiments show that if the liquid is confined within boundaries of a set of inter-communicating vessels, then the proposed method converges to a final equilibrium state. This incremental static analysis approximates dynamic behavior with strong damping and can provide information, as a first approximation to 2D movement of a liquid. In the initial arrangement of particles, a rhombic element is proposed, which assures satisfactory incompressibility of the fluid. Based on the proposed algorithm, a simple and short computer program (a “pocket” program) has been developed, with only about 120 Fortran instructions. This program is first applied to an amount of liquid, contained in a single vessel. A coarse and refined discretization is tried. In final equilibrium state of liquid, the distribution on hydro-static pressure on vessel boundaries, obtained by proposed computational model, is found in satisfactory approximation with corresponding theoretical data. Then, an opening is formed, at the bottom of a vertical boundary of initial vessel, and the liquid is allowed to flow gradually to an adjacent vessel. Almost whole amount of liquid is transferred, from first to second vessel, except of few drops-particles, which remain, in equilibrium, at the bottom of initial vessel. In the final equilibrium state of liquid, in the second vessel, the free surface level of the liquid confirms that the proposed rhombing element assures a satisfactory incompressibility of the fluid.


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 intercolliding 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].

Proposed Method
In 2D (two dimensional space), a single vessel (Figure 1 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 x u y u , propor- tional to out-of-balance forces x F , y F , 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.

Constant Input Data
The boundary conditions are given, that is, position, configuration and dimensions of linear boundaries of single vessel (Figure 1

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 b bu l l ≥ and p pu l l ≥ .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, 0 x F = , y F w = − .

Any Step of Algorithm
For all out-of-balance nodal forces x F , y F , acting on particles, the maximum absolute value is found, ( ) Then, every particle performs small displacements x u , y u , proportional to its out-of-balance forces x F , y F , as follows: From the new positions (x, y) of the particles, the boundary-particle distances b l and the inter-particle distances p l are found.And the boundary-particle and inter-particle repulsive forces are determined, The out-of balance force, acting on every particle, is found, as shown in Figure 3: where w weight of the particle, b R repulsive force from adjacent boundary and p R repulsive force from adjacent particle.

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 x u , y u , present position (co-ordinates x, y) and out-of-balance forces x F , y F .For particles adja- cent to boundaries, with b bu l l < , the external repulsive reactions For every couple of adjacent particles, with p pu l l < , the mutual internal repulsive reaction ( ) p p R l .Also, the maximum absolute value of all nodal forces ( ) x y F F F = can be received as output of present step.

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

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 b b R l − function + only 20 for subroutine of inter-particle p p R l − 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.

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 l , for every couple of adjacent particles, with p pu l l < ), where the limit distance is 2 0.707 , 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 pu a l > , are activated later.For small deformations 0.293 u a < 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),   ( ) ( ) 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 0.293 u a > , 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 incompressibility, as will also be confirmed in the following applications.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 50 N w = . So, the weight of total amount of liquid is 180 50 N 9.0 kN W = × = .In , so that distance, from boundary to adjacent particle, will be not reduced more than 10%.
On the other hand, the pu R 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, 0.1 u ∆ = 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 b bu l l ≥ and p pu l l ≥ ; 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, 0 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.

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 ticles, alternated with 18 rows of 19 particles, that is totally 2 × 18 × 19 = 684 particles.Weight of a particle is 12.5 N w = , thus total weight 684 12.5 8.55 kN W N = × = .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 , 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, 0 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 x F , y F , 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 100 cm in y ≈ (not yet in equilibrium) to about 75 cm eq y ≈ 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.

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.  Third application.In-complete flow from initial to adjacent vessel.Ouput.After final equilibrium state, in initial vessel, an opening is formed at bottom of its right vertical boundary, and the liquid is allowed to flow gradually to an adjacent vessel.In 100,000 steps of algorithm and 5.0 min of computing time, more than half of total amount of liquid has been transferred from first to second vessel.However, there is not yet an equilibrium state.There is a significant resultant of horizontal external reactions  However, there is not yet an equilibrium state, as significant out-of-balance nodal forces still remain, much larger, in absolute values x F , y F , 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 0 x F ∑  , that is, there exists a significant resultant of horizontal boundary reactions, directed to right.

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 x F , y F , less than 0.1 N, so the algorithm is interrupted and the state of Figure 10 is obtained.Where (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 b b R l − (Figure 1(c)) describes the variation of repulsive force b R , from a boundary to an adjacent particle, versus distance b l of particle, perpendicularly to boundary.As shown in Figure 1(c), for b bu l l ≥ , 0 b R = .Also, an inter-particle function p p R l − (Figure 1(d)) describes the variation of mutual repulsive force p R versus distance p l , for every couple of adjacent particles.As shown in Figure 1(d), for p pu l

Figure 1 .
Figure 1.In 2D, amount of a liquid is simulated by a set of inter-colliding point-particles.
(a)) or set of inter-communicating vessels (Figure 1(b)).The weights w of the particles.The boundary-particle function b b R l − (Figure 1(c)).Inter-particle function p p R l − (Figure 1(d)).Also, the steplength Δu of the algorithm is given.

Figure 2 .
Figure 2. Flow-chart of proposed algorithm, for incremental static analysis of a set of inter-colliding point-particles, simulating 2D flow.

Figure 3 .R
Figure 3. Forces acting on a particle: w particle weight.b R repulsive force from adjacent boundary, 1 p R , 2 p R repulsive force from adjacent particle.
(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 p p R l − function (mutual repulsive force p R versus distance p

Figure 4 .
Figure 4. (a) The proposed rhombic element; (b) Inter-particle p p R l − function; (c) Rhomb mechanism, for small deformation u < 0.293a, before activation of internal diagonals; (d) Reduction of initial un-deformed area of rhomb mechanism, until limit deformation u = 0.293a.

1 .Figure 5 .
Figure 5. First application.Single vessel.Coarse discretization.Input data.(a) Vessel configuration and dimensions.Initial positions (x, y co-ordinates) of point-particles; (b) Dimensions in enlarged rhombic element and distances from boundaries; (c) Boundary-particle b b R l − function; (d) Inter-particle p p R l − function.

Figure 5 (Figure 5 (
Figure 5(d), respectively.The parameters of these diagrams have been determined as follows: 5.0 cm bu l =

xFFigure 6 .
Figure 6.First application single vessel coarse discretization ouput data.(a) Positions of particles in final equilibrium state.Free surface level of liquid has descended to about 72.5 cm eq y ≈ .External repulsive reactions on boundaries, obtained by coarse discretization, compared to corresponding theoretical forces.(b) Theoretical hydrostatic pressure distribution on vessel boundaries, in final equilibrium state.

Figure 7 .
Figure 7. Second application.Single vessel.Refined discretization.Input data.(a) Configuration and dimensions of vessel.Initial positions of particles.Single vessel; (b) Dimensions in enlarged rhombic element and distances from boundaries; (c) Boundary-particle b b R l − function; (d) Inter-particle p p R l − function.

Figure 7 (
b).The boundary-particle function b b R l − and the inter-particle function p p R l − are shown in Figure 7(c) and Figure 7(d), respectively.The parameters of these functions, as well as the step-length 0.025 mm u ∆ =

Figure 8 .
Figure 8. Second application.Single vessel.Refined discretization.Ouput data.(a) Positions of particles in final equilibrium state.Free surface level of liquid has descended to about 75 cm eq y ≈ .External repulsive reactions on boundaries, obtained by refined mesh, compared to corresponding theoretical forces.(b) Theoretical hydro-static pressure distribution on boundaries, in final equilibrium state.

Figure 9 .
Figure 9. Third application.In-complete flow from initial to adjacent vessel.Ouput.After final equilibrium state, in initial vessel, an opening is formed at bottom of its right vertical boundary, and the liquid is allowed to flow gradually to an adjacent vessel.In 100,000 steps of algorithm and 5.0 min of computing time, more than half of total amount of liquid has been transferred from first to second vessel.However, there is not yet an equilibrium state.There is a significant resultant of horizontal external reactions

Figure 10 .
Figure 10.Fourth application.Complete flow from initial to adjacent vessel.Final equilibrium state.Ouput.From the final state of incomplete flow of third application, the liquid continues to flow gradually, from first to second vessel.In 300,000 steps of algorithm and 15 min of computing time, from beginning of first application, almost whole amount of liquid has been transferred to second vessel, except of fourteen drops-particles remaining, in equilibrium, at bottom of first vessel.The state, in second vessel, is few steps of algorithm before equilibrium.The free surface level of liquid, in second vessel, confirms that a satisfactory incompressibility of the fluid is assured by the proposed rhombic element.

Table 1 .
Comparison of computational data of four applications, as regards number of particles, step-length Δu of algorithm, approximate values for number of algorithm steps required and computing time consumed.The last two rows refer to applications, by refined mesh, not presented in this work.