Simulation of fluid dynamics through complicated networks of channels with cellular automata

There are different options or tools available when modeling fluid dynamics: solving the analytic equations, or, if that is not possible, considering a numeric solution (finite difference, finite volume, or finite element methods) [1]. In general these methods work with non-complex borders. With complex structures, fluid dynamics numeric simulation takes a heavy burden on computation time to find solutions to solve Navier–Stokes equations. Another option is to solve multiple border problems in statistical terms or with lattice Boltzmann models to simplify the dynamics at a macroscopic level [1–3]. Another approach is to handle a mixture of particles and a continuous medium [4–6]. Recently, a method known as smoothed particles hydrodynamics has been proposed to simulate fluids in terms of the wave-particle concept [7–9]. The easiest way we found to take account of multiple borders and less computation time is the cellular automata (CA) framework [10]. CA is founded on simple rules to represent complex behaviors [10]. CA has been used to simulate gas behavior in a container [11, 12] which uses a rectangular grid; there are other similar works that instead use an hexagonal mesh [13–17]. More recent studies have begun to simulate liquids [18].

CA are the simplest discrete representation of a complex dynamic system [19, 20] and have been applied in many areas of science [21–24], since they are an alternative to study connections between the microscopic and the macroscopic world. Due to their computational speed, CA allow one to study a wide range of values of the parameters involved in the problem that otherwise would involve excessive computation time, and are a simple tool to qualitatively test predictions about how a local mechanism can generate a certain response [17, 18]. However, CA models can almost never make accurate quantitative predictions, but they are suitable for modeling processes in which the basic laws are difficult to identify. CA are simple enough to allow a detailed mathematical analysis but as complex as necessary to exhibit a wide variety of complicated phenomena [25].

CA are commonly used in 1-D or 2-D meshes. The position in the mesh is labeled with the subscripts i, j. Each position can take a number of k different states. The states can only take whole values: x_i (if you work in 1-D) or x_i,j (if you work in 2-D). For example, in 1-D with k = 2, the highway model can take only two possible values [26],

$$\begin{eqnarray*}{x}_{i} & = & 0\quad {if}\quad {is}\quad {empty}\\ {x}_{i} & = & 1\quad {if}\quad {is}\quad {occupied}.\end{eqnarray*}$$

An example in 2-D with k = 3 could be Kermack-McKendrick's model (1927), as follows representing an infectious disease [26]:

$$\begin{eqnarray*}{x}_{i,j} & = & 0\quad {if}\quad {is}\quad {suceptible}\\ {x}_{i,j} & = & 1\quad {if}\quad {is}\quad {infected}\\ {x}_{i,j} & = & 2\quad {if}\quad {is}\quad {inmune}.\end{eqnarray*}$$

All CA cells are updated simultaneously in a discrete time. The state of each position for a later time has a deterministic function in terms of the states of its closest neighbors at the present time. The update rule is homogeneous, that is, all positions use the same rule. In CA, different types of borders can be used: these can be solid or can represent a toroidal space. The CA can be deterministic or stochastic, depending on whether they are considered random events or not; that is, they are considered deterministic if they can be reconstructed in both time directions, and stochastic if they cannot be reconstructed backwards in time. The mesh used in the present work is known as a first-order Moore neighborhood, in which each individual has eight close neighbors that influence the state that the CA would take in the next time, except those close to the borders [20]; see figure 1.

Zoom In Zoom Out Reset image size

The modeling strategy using CA is to define the initial condition of the system, the system update rules and the definition of the parameters used and then represent the entire system in a programming language, and let the simulation run for the desired time.

In this paper we will follow the methodology proposed by [26], where we must first define what possible values the cells of the domain can have and then indicate the rules to follow, because we consider two fluids in some simulations. The CA proposed in this work, to solve gravitational fluid flow, has the following characteristics.

A 2-D CA, with k = 4, is considered with these possible values:

$$\begin{eqnarray*}{x}_{i,j} & = & 0\quad {if}\quad {the}\quad {cell}\quad {is}\quad {empty}\\ {x}_{i,j} & = & 1\quad {if}\quad {the}\quad {cell}\quad {has}\quad a\quad {fluid}\quad 1\quad {particle}\\ {x}_{i,j} & = & 2\quad {if}\quad {the}\quad {cell}\quad {has}\quad a\quad {fluid}\quad 2\quad {particle}\\ {x}_{i,j} & = & 3\quad {if}\quad {it}\quad {is}\quad a\quad {solid}\quad {border}.\end{eqnarray*}$$

The control rules governing the model dynamics for each cell in relation to its neighborhood for the next time step are as follows:

• Each particle can only move to one of the five spaces available (horizontal and downward).
• Movement is set in a probabilistic environment with higher probability to move downward.
• Because we want to simulate the fall of a fluid, we have considered only limiting the movement of the particles in the vertical direction, but the particles cannot move upward, except to balance hydrostatic pressure.
• Particles cannot pass solid borders (cells occupied by solid material).
• If there are no empty spaces in its neighborhood, the particle stays in its cell.
• Horizontal movement is chosen randomly to avoid any biased movement.
• When two fluids are modeled, if a heavier particle is over a lighter particle, particles interchange position.
• The particles that enter the CA system have to remain in the system.

The probability values reflect some of the characteristics that a particular phenomenon can have. For example, in this study gravity plays an important role so more importance is given to the vertical direction, since the force of gravity points in that direction, and smaller values to the lateral directions (see figure 1) and null in the vertical direction. A variant that could be programmed by this same CA is the behavior of a gas: for this case an equal probability to all the cells around the particle should be considered.

The first simulation performed was intended to see the behavior of the liquid particles in an environment where barriers prevent their passage. A constant flow from a pipe supplying fluid particles, which is positioned at coordinates (2, 4) and (2, 5), was considered. Figure 2 shows some scenes at different simulation time intervals. The fluid fills the first container before moving to the second dike, and since lateral movement is allowed, the fluid first fills the container with higher barriers (figures 2 (b) and (c)). The second example (figure 3) consists of a simulation with a finite number of liquid particles that are at the top and left at the start of the simulation, and due to the effect of gravity tend to go towards the bottom. Particles find some vertical channels that are filled with the liquid on its way down. The third simulation (figure 4) is a container with two fluids of different density, found also at the top at the start of the simulation. The less dense fluid initially floats at the top and center of a space containing both fluids, but as the simulation progresses, it moves to the left and right sides, and then falls by gravity. Another thing to keep in mind is the formation of air bubbles that appear when the liquid fills the free spaces of the container, therefore moving the air particles upward (figures 4(b)–(e)). It should also be noted that one of the containers (top) was first filled by the denser fluid so the less dense fluid had no opportunity to fill, but the less dense fluid remained above the denser fluid, as was expected. The built-in CA only work with downward and lateral movements, but do not allow upward movement of liquid particles. But there is a case where CA with these conditions would not work: when we perform a simulation on a U-tube, filling the liquid through one of the openings. It is expected that the liquid will upsurge on the other side of the tube due to the hydrostatic balance. However, this does not happen because the dynamics are based in local conditions of CA, while the hydrostatic balance can act at a much wider distance, so we decided to make an adjustment to our automaton and, only under this new scenario, to allow the movement of the fluid upwards. The last simulation was intended to make the cellular automaton work at times when its own rules prevented this from happening. Figure 5 shows the CA behavior improvement with that programming modification. Because we have the problem of upward vertical movement, we wanted to take into account the hydrostatic balance; for this we have made a first attempt, in a numerical experiment representing a U-tube, wherein the falling liquid rises from the other side for hydrostatic balance. For this we have proposed the use of the contour function, which helps us to delimit the fluid, taking into account that, once past the process of moving all the particles of the fluid, the contour function delimits all the fluid and allows us to know if there is fluid above that can be moved by hydrostatic pressure. It is important to take into account how many plots of fluids have been detected by this process, in order to adjust the hydrostatic balance one by one. The functions used only in the U-tube experiment are known as contour properties, and are as follows:

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

[C, h] = contour (input variables).

One has to be careful with the variable C, since is it necessary to distinguish which line is the right one to be used as a delimiter of the fluid.

As we have shown, complex liquid-fluid systems can be properly simulated with simple rules in a CA framework. This is an advance in the simulation of fluids in porous media, where the complexities of the cavities are a challenge during the codification process of boundary conditions in models, based on the governing equations of their movement. Different embodiments of the same experiment were performed: these experiments showed that all the realizations at the end were the same. However, we know that there is much work still to be done in this way to be able to say that CA are the solution to all fluid problems. For example, it is necessary to include the conservation of momentum when particles with different speeds collide. However the simple use of CA makes it attractive to explore more applications in fluid dynamics systems. The adjustment of the probability values can make the fluids have a better performance in these models when simulating different fluids. Although these are initial results and more complex conditions should be tested, we consider that the modeled fluid behavior was as expected in the conditions where they were tested. We believe there is a potential to model complicated scenarios with CA in very simple ways, as well as it being a tool that can be used for teaching purposes.

The authors wishes to thank the support of the UABC (PREDEPA) and SEP (PRODEP), with its academic mobility program. We also want to thank Ocean. Martha E Betancourt-Aguirre for work on the review of this paper.