**Optimization of Pouring Velocity for Aluminium Gravity Casting**

Y. Kuriyama, K. Yano and S. Nishido

*Gifu National College of Technology Mie University AISIN TAKAOKA CO., LTD Japan* 

#### **1. Introduction**

574 Fluid Dynamics, Computational Modeling and Applications

The high temporal and spatial resolution technologies employed for the reconstruction of the fluid-dynamic field inside the Multidune allows recognizing the flow field within the apparatus with both arrangements is characterized by three areas: the main transport current and, in each chamber, recirculation areas above and below. The fluid-dynamic behaviour is substantially similar in each chamber but in the first and the last one (C1 and C8), where the inlet and outlet nozzles prevent the formation of similar vortical structures. With Arrangement A, the characteristic velocity of the principal current appears significantly larger than the velocity within the upper and lower recirculation areas; this aspect is amplified with increasing hydraulic head at the apparatus inlet. With Arrangement B they appear more comparable. With both arrangements, the increase of the hydraulic head augments the transport effectiveness of the main current without improving the capture feasibility of both the upper and lower recirculation zones. The apparatus will then lose its effectiveness in separating plastic particles increasing both the hydraulic head and the transiting flowrate. Furthermore, with both arrangements, C3, C4, C5 and C6 present analogous recirculation areas. The same circumstance occurs at each flowrate establishing in the apparatus. Then, the presence of eight chambers assures plastic particle separation even if a chamber should be filled with the settled material allowing the following chambers to

The authors would like to acknowledge Dr. Leonardo Cherubini and Dr. Emanuela Lupo for

De Sena G., Nardi C., Cenedese A., La Marca F., Massacci P., Moroni M. (2008). The

La Marca F., Moroni M., Cherubini L., Lupo E., Cenedese A. (2011). Recycling of plastic

Moroni M., Cenedese A. (2005). Comparison among feature tracking and more consolidated

Moroni M., Kleinfelter N., Cushman J.H. (2008). Alternative Measures of Dispersion Applied to Flow in a Convoluted Channel. *Advances in Water Resources* 32(5), 737-749. Moroni M., La Marca F., Cherubini L., Cenedese A. (2011). Recovering plastics via the

PlasticsEurope (2007b). Press Release 9 May 2007, 1 p. Association of Plastics Manufacturers

Plastic Separation Feasibility. *Waste Management* 29(9), 1560-1571. Hussain F. (1986). Coherent structures and turbulence, J. Fluid Mech. 173, 303. Joeng G., Hussain F. (1995). On the identification of a vortex, J. Fluid Mech. 285, 69.

*Measurement Science and Technology* 16, 2307-2322.

Europe. PlasticsEurope, 21 pp.

in Europe (AISBL), Brussels, Belgium.

Hydraulic Separator Multidune: Preliminary Tests on Fluid-Dynamic Features and

waste via the hydraulic separator Multidune. *Waste Management*, ISSN: 0956-053X,

velocimetry image analysis techniques in a fully developed turbulent channel flow.

hydraulic separator Multidune: flow analysis and efficiency tests. *International Journal of Environmental Science and Technology*, ISSN: 1735-1472, (submitted). PlasticsEurope (2007a). An analysis of plastics production, demand and recovery for 2005 in

**6. Conclusion** 

become effective in the separation process.

their contribution during the experiments.

**7. Acknowledgments** 

(submitted).

**8. References** 

In current casting factories, tilting type automatic pouring machines are often used to pour the molten metal into the mold, with the operator relying on experience, perception and repeated testing to manually determine the pouring velocity. However, seeking an optimum multistep pouring velocity through trial and error results in an enormous number of combinations and is very difficult. For this reason, it cannot be said that suitable casting that realizes a high-quality cast is being carried out, inviting a decline in the yield rate due to product defects.

Furthermore, the extension of the production preparatory phase and increase in costs due to this kind of trial operation also become a significant problem.

Until now, much research relating to product defects in aluminum gravity molding has been conducted (Yutaka et al., 2001)(Takuya, 2004). Meanwhile, research applying casting CAE for the purpose of improving the quality of castings and production efficiency is coming to attention (Itsuo, 2006). Furthermore, in recent years, optimization of the casting problem has begun to be carried out in accompaniment with developments in computers (Takuya et al., 2007)(Ken'ichi et al., 2008). However, these all target comparatively short calculation time die-casting and adoption of optimization technology in sand mold casting and gravity casting is lagging.

Accordingly, in this research, the objective is to stabilize the fluid speed in the mold gate and derive an optimum pouring velocity to realize a mitigation of defects such as pin holes and blow holes in aluminum gravity casting through invoking a fluid behavior simulator, swiftly filling the sprue cup and controlling the liquid level at a fixed high level of liquid. For the automatic pouring machine, the multistep velocity input is designed for actual product of an intake manifold. The effectiveness of this research is shown by a fluid analysis simulation and an actual test.

#### **2. Experimental apparatus**

An overview of the automatic pouring machine is shown in Fig. 1. It is a tilting type automatic pouring machine with one degree of freedom in the forward and back direction of the ladle. In the pouring machine, the tilting angular velocity and velocity switching angle are configured by teaching pendant and the pouring velocity is determined. The setup enables

Optimization of Pouring Velocity for Aluminium Gravity Casting 577

In this pouring machine, the tilting movement to input is unknown because the machine doesn't have the output device for the angular velocity or the angle. At the simulation of the fluid of molten metal, the identification of the tilting movement to input is needed. Thus, to get the unknown parameter of *V*max and *α*, the movement to input of the actual poring machine is filmed with the motion capture system. Table 1 shows the input of the pouring

Switching angle (deg) Pouring velocity (%)

The unknown parameter of *V*max and *α*, is identified by using the data of the angular velocity from the motion capture. The angular acceleration and maximum angular velocity are *α*=200 (deg/s2), *V*max =51.9 (deg/s) respectively. Fig.4 shows the result of identification of actual pouring machine, where the solid line is the path of actual pouring machine and the broken

From the Fig.4, it can be seen the motion is accorded very well. In this result, the plant

1 22 10 2 32 30 3 42 50 4 -- 30

Fig. 3. Tilting velocity of pouring machine.

machine using analysis the motion.

line is the path of plant model.

**3. Identification the motion of the pouring machine** 

Table 1. Setting of the tilting input for identification the motion.

model can be recreated the motion of the actual movement.

four steps of tilting angle velocity and three of switching angle for a total of seven variables as shown in Fig.2. The tilting angular velocity command is assigned in the form of raised trapezoidal shapes as shown in Fig.3.

Fig. 1. Overview of the pouring machine


Fig. 2. Input setting of the automatic pouring machine.

As the tilting angular velocity setting is displayed as a percentage, each step of tilting angular velocity is displayed by the following formula using the maximum angular velocity.

$$V\_{\rm n} = \frac{\upsilon\_{\rm n}}{100} V\_{\rm max} \tag{1}$$

Here, it is necessary to carry out parameter identification, since maximum angular velocity *V*max and the angular acceleration *α* are unknown parameters.

four steps of tilting angle velocity and three of switching angle for a total of seven variables as shown in Fig.2. The tilting angular velocity command is assigned in the form of raised

trapezoidal shapes as shown in Fig.3.

Fig. 1. Overview of the pouring machine

Fig. 2. Input setting of the automatic pouring machine.

*V*max and the angular acceleration *α* are unknown parameters.

As the tilting angular velocity setting is displayed as a percentage, each step of tilting angular velocity is displayed by the following formula using the maximum angular velocity.

<sup>n</sup> <sup>n</sup> max <sup>100</sup>

Here, it is necessary to carry out parameter identification, since maximum angular velocity

*<sup>v</sup> V V* <sup>=</sup> (1)

Fig. 3. Tilting velocity of pouring machine.

#### **3. Identification the motion of the pouring machine**

In this pouring machine, the tilting movement to input is unknown because the machine doesn't have the output device for the angular velocity or the angle. At the simulation of the fluid of molten metal, the identification of the tilting movement to input is needed. Thus, to get the unknown parameter of *V*max and *α*, the movement to input of the actual poring machine is filmed with the motion capture system. Table 1 shows the input of the pouring machine using analysis the motion.


Table 1. Setting of the tilting input for identification the motion.

The unknown parameter of *V*max and *α*, is identified by using the data of the angular velocity from the motion capture. The angular acceleration and maximum angular velocity are *α*=200 (deg/s2), *V*max =51.9 (deg/s) respectively. Fig.4 shows the result of identification of actual pouring machine, where the solid line is the path of actual pouring machine and the broken line is the path of plant model.

From the Fig.4, it can be seen the motion is accorded very well. In this result, the plant model can be recreated the motion of the actual movement.

Optimization of Pouring Velocity for Aluminium Gravity Casting 579

(m) and the number of wires is 50 in the vertical and 55 in the horizontal. In this research a porous baffle (hereinafter, baffle) was used to reproduce this filter in the CFD simulator. The air porosity *b*p, the linear velocity drop coefficient *b*l and the two-dimensional velocity drop coefficient *b*q are assigned as setting parameters for the baffle. *b*l and *b*q are defined by

( ) *bu b u u*

Here, *B* denotes the baffle flow loss, *u* the flow speed within the baffle and *L* the length in which the flow loss is produced. The air porosity of the wire mesh can be calculated from the area ratio of the metal wires and opening between them. As the linear velocity drop is dominant for the baffle flow loss, the 2-D velocity drop coefficient is set at *b*q=0 and the results of searching for the linear velocity drop coefficient *b*l are shown in Table 4. Furthermore, the search range was carried out in 0.05 increments over *b*l=0.00–1.50.

*<sup>B</sup>* <sup>l</sup> <sup>5</sup> <sup>q</sup> 0. <sup>1</sup> <sup>=</sup> <sup>+</sup> (2)

*L*

Simulation results considering the baffle loss are shown in Fig. 6.

an equation for baffle flow loss shown by Eq. (2).

Fig. 5. Mesh setting of CFD

Fig. 4. Result of identification of actual pouring machine.

#### **4. Construction of fluid behavior simulator and flow evaluation**

The fluid analysis software FLOW-3D was used in this research. AC2B is taken as the subject fluid. The cast quantity is 1.863×10-3(m3) and the product part volume is 1.429×10-3 (m3). The physical values of AC2B for the analysis were set as in Table 2. An outline map of the mesh in the simulation domain is shown in Fig. 5 and the mesh parameters are shown in Table 3.


Table 2. Fluid parameters


Table 3. Mesh parameters

In an actual plant, a molten metal filter (wire mesh) is installed in the sprue runner with the purpose of removing slag as shown in Fig. 5. The thickness of the wires in the mesh is 0.5×10-3

Fig. 4. Result of identification of actual pouring machine.

Table 2. Fluid parameters

Table 3. Mesh parameters

**4. Construction of fluid behavior simulator and flow evaluation** 

The fluid analysis software FLOW-3D was used in this research. AC2B is taken as the subject fluid. The cast quantity is 1.863×10-3(m3) and the product part volume is 1.429×10-3 (m3). The physical values of AC2B for the analysis were set as in Table 2. An outline map of the mesh in the simulation domain is shown in Fig. 5 and the mesh parameters are shown in Table 3.

> Fluid parameters AC2B aluminum alloy Density 2550 kg/m2 Viscosity 0.00125 Pa・s

Specific heat 1071 J/(kg・K) Thermal conductivity 100 W/(m・K)

Mesh block Cell size (m) Number of cell X-direction 0.004~0.080 108 Y-direction 0.005 54 Z-direction 0.004~0.080 65 Total number of cell 379080

In an actual plant, a molten metal filter (wire mesh) is installed in the sprue runner with the purpose of removing slag as shown in Fig. 5. The thickness of the wires in the mesh is 0.5×10-3

Temperature of fluid 993 K

(m) and the number of wires is 50 in the vertical and 55 in the horizontal. In this research a porous baffle (hereinafter, baffle) was used to reproduce this filter in the CFD simulator.

The air porosity *b*p, the linear velocity drop coefficient *b*l and the two-dimensional velocity drop coefficient *b*q are assigned as setting parameters for the baffle. *b*l and *b*q are defined by an equation for baffle flow loss shown by Eq. (2).

Fig. 5. Mesh setting of CFD

$$B = \frac{1}{L}(b\_l u + 0.5b\_\mathbf{q} u |u|) \tag{2}$$

Here, *B* denotes the baffle flow loss, *u* the flow speed within the baffle and *L* the length in which the flow loss is produced. The air porosity of the wire mesh can be calculated from the area ratio of the metal wires and opening between them. As the linear velocity drop is dominant for the baffle flow loss, the 2-D velocity drop coefficient is set at *b*q=0 and the results of searching for the linear velocity drop coefficient *b*l are shown in Table 4. Furthermore, the search range was carried out in 0.05 increments over *b*l=0.00–1.50. Simulation results considering the baffle loss are shown in Fig. 6.

Optimization of Pouring Velocity for Aluminium Gravity Casting 581

improper velocity setting and unstable liquid levels inside the sprue cup etc. Accordingly, in this research, liquid level control is realized through optimizing the pouring input using a

GA is an algorithm that models natural selection and mutation in the processes of inheritance and evolution in biological groups in the processes of evolution and inheritance

The three steps of switching angle and four steps of pouring velocity for a total of seven set parameters are taken as variables that are the settings input for the pouring machine, and an optimum tilting velocity pattern is calculated within the limitations of the real machine. In order to swiftly fill the sprue cup and stabilize the liquid surface at a uniform level, the tilting end time is taken as a performance function and the optimization problem is expressed by Eq. (3) with the liquid level inside the sprue as a limiting condition. Through taking the tilting end time as a performance function, the filling time is reduced and through already taking the

*h* ≥ 0.025 Here, *t*p denotes the tilting end time, *J*p a penalty function denoted by Eq. (4), and h the

In Eq. (4), the penalty clause *w*1,*w*2=100 is imposed in the case the displacement from the

Optimization by GA was carried out using the calculation parameters shown in Table 5 Forty-eight hours was required for optimization with 41st generation using an Intel Core2

> Number of variable 7 Number of population 10 Number of elite preservation 1 Mutation fraction 0.01 Crossover fraction 0.80

At that time the evaluation value was *J*=4.668. The tilting angles and velocities obtained from the optimum parameters are shown in Fig. 7. And the simulation results of the liquid

As a result, it can be seen that the objective liquid level is not reached in the case of any control and a satisfactory liquid surface is not maintained. Conversely, it can be seen that a satisfactory liquid level control that swiftly fills the sprue cup is realized through the

sprue exterior to the sprue interior liquid level drops to below 0.025(m) (*h* < 0.025).

<sup>p</sup> <sup>p</sup> minimize : *<sup>J</sup>* <sup>=</sup> *<sup>t</sup>* <sup>+</sup> *<sup>J</sup>* (3)

<sup>p</sup> *w*<sup>1</sup> *w*<sup>2</sup> *J* = + (4)

for populations in the natural world and is a probabilistic optimization method.

liquid level as a limiting condition, liquid level control becomes possible.

displacement from the sprue exterior to the molten metal.

**6. Optimization with CFD simulator** 

Table 5. Parameters of genetic algorithms

optimization of pouring control input.

level control is shown in Fig. 8.

Quad CPU equipped PC.

**6.1 Optimization result and melt flow analysis** 

genetic algorithm (GA) (Thomas et al., 1993).


Table 4. Parameters of porous baffle

Fig. 6. Comparison of flow in sprue cup between simulation with molten metal filter and without molten metal filter.

Comparing these with the pouring test results, it can be seen that a satisfactory reproduction of the molten metal behavior inside the sprue is achieved.

### **5. Derivation of optimum pouring input using a genetic algorithm**

Through swiftly filling the inside of the sprue cup with molten metal and controlling the liquid level to preserve a liquid level with high uniformity, the flow velocity in the mold gate is made constant.

This is regarded as making possible the production of high-quality casting that mitigates product defects. However, in the case the pouring velocity is determined by operator trial and error, problems occur such as the overflowing of molten metal from the melt due to

Linear loss coefficient : *b*l 0.90 Quadratic loss coefficient : *b*q 0.00

Fig. 6. Comparison of flow in sprue cup between simulation with molten metal filter and

**5. Derivation of optimum pouring input using a genetic algorithm** 

Comparing these with the pouring test results, it can be seen that a satisfactory reproduction

Through swiftly filling the inside of the sprue cup with molten metal and controlling the liquid level to preserve a liquid level with high uniformity, the flow velocity in the mold

This is regarded as making possible the production of high-quality casting that mitigates product defects. However, in the case the pouring velocity is determined by operator trial and error, problems occur such as the overflowing of molten metal from the melt due to

Table 4. Parameters of porous baffle

without molten metal filter.

gate is made constant.

of the molten metal behavior inside the sprue is achieved.

Void 0.655

improper velocity setting and unstable liquid levels inside the sprue cup etc. Accordingly, in this research, liquid level control is realized through optimizing the pouring input using a genetic algorithm (GA) (Thomas et al., 1993).

GA is an algorithm that models natural selection and mutation in the processes of inheritance and evolution in biological groups in the processes of evolution and inheritance for populations in the natural world and is a probabilistic optimization method.

The three steps of switching angle and four steps of pouring velocity for a total of seven set parameters are taken as variables that are the settings input for the pouring machine, and an optimum tilting velocity pattern is calculated within the limitations of the real machine. In order to swiftly fill the sprue cup and stabilize the liquid surface at a uniform level, the tilting end time is taken as a performance function and the optimization problem is expressed by Eq. (3) with the liquid level inside the sprue as a limiting condition. Through taking the tilting end time as a performance function, the filling time is reduced and through already taking the liquid level as a limiting condition, liquid level control becomes possible.

$$\text{minimize}: J = t\_{\text{p}} + J\_{\text{p}} \tag{3}$$

*h* ≥ 0.025

Here, *t*p denotes the tilting end time, *J*p a penalty function denoted by Eq. (4), and h the displacement from the sprue exterior to the molten metal.

$$J\_p = \mathbf{w}\_1 + \mathbf{w}\_2 \tag{4}$$

In Eq. (4), the penalty clause *w*1,*w*2=100 is imposed in the case the displacement from the sprue exterior to the sprue interior liquid level drops to below 0.025(m) (*h* < 0.025).

#### **6. Optimization with CFD simulator**

#### **6.1 Optimization result and melt flow analysis**

Optimization by GA was carried out using the calculation parameters shown in Table 5 Forty-eight hours was required for optimization with 41st generation using an Intel Core2 Quad CPU equipped PC.


Table 5. Parameters of genetic algorithms

At that time the evaluation value was *J*=4.668. The tilting angles and velocities obtained from the optimum parameters are shown in Fig. 7. And the simulation results of the liquid level control is shown in Fig. 8.

As a result, it can be seen that the objective liquid level is not reached in the case of any control and a satisfactory liquid surface is not maintained. Conversely, it can be seen that a satisfactory liquid level control that swiftly fills the sprue cup is realized through the optimization of pouring control input.

Optimization of Pouring Velocity for Aluminium Gravity Casting 583

depends on whether or not the intensity of the turbulence is enough to overcome the

Turbulence transport models characterize turbulence by a specific turbulent kinetic energy *Q* and a dissipation function *D*. The characteristic size of turbulence eddies is then

0.1 *<sup>Q</sup> <sup>L</sup>*

This scale is used to characterize surface disturbances. The disturbance kinetic energy per unit volume (i.e., pressure) associated with a fluid element raised to a height *L*, and with

surface tension energy based on a curvature of *L* is given by Eq. (6).

3

*<sup>D</sup>* <sup>=</sup> (5)

surface-stabilizing forces of gravity and surface tension as shown in Fig. 10.

Fig. 9. Photograph of blow hole

Fig. 10. Model of air entrainment

given by Eq. (5).

Fig. 7. Tilting input of optimization result.

Fig. 8. Tilting input of optimization result.

#### **6.2 Evaluation of the optimum input**

Air entrainment is one of the defect origins of such as pin holes and blow holes as shown in Fig.9. Thus, using the evaluation function of air entrainment which is one of the functions of *Flow-3D*, the optimum input is evaluated.

The air entrainment at the liquid surface is based on the concept that turbulent eddies raise small liquid elements above the free surface that may trap air and carry it back into the body of the liquid. The extent to which liquid elements can be lifted above the free surface depends on whether or not the intensity of the turbulence is enough to overcome the surface-stabilizing forces of gravity and surface tension as shown in Fig. 10.

Fig. 9. Photograph of blow hole

582 Fluid Dynamics, Computational Modeling and Applications

Air entrainment is one of the defect origins of such as pin holes and blow holes as shown in Fig.9. Thus, using the evaluation function of air entrainment which is one of the functions of

The air entrainment at the liquid surface is based on the concept that turbulent eddies raise small liquid elements above the free surface that may trap air and carry it back into the body of the liquid. The extent to which liquid elements can be lifted above the free surface

Fig. 7. Tilting input of optimization result.

Optimum input

Conventional input

Fig. 8. Tilting input of optimization result.

**6.2 Evaluation of the optimum input** 

*Flow-3D*, the optimum input is evaluated.

Fig. 10. Model of air entrainment

Turbulence transport models characterize turbulence by a specific turbulent kinetic energy *Q* and a dissipation function *D*. The characteristic size of turbulence eddies is then given by Eq. (5).

$$L = 0.1 \frac{\sqrt{Q^3}}{D} \tag{5}$$

This scale is used to characterize surface disturbances. The disturbance kinetic energy per unit volume (i.e., pressure) associated with a fluid element raised to a height *L*, and with surface tension energy based on a curvature of *L* is given by Eq. (6).

Optimization of Pouring Velocity for Aluminium Gravity Casting 585

lower figure shows the optimum input, where the color of this figure is indicated the value of the air entrapment, and red color is the most air entarinment point. From the figure, the optimum input can be decrease the air entrainment, and it can be expect to realize a

mitigation of defects such as pin holes and blow holes in experiment.

Fig. 12. Flow of molten metal with the mold

$$P\_d = p \text{g}L + \frac{\sigma}{L} \tag{6}$$

where *ρ* is the liquid density, *σ* is the coefficient of surface tension, and *g* is the component of gravity normal to the free surface.

For air entrainment to occur, the turbulent kinetic energy per unit volume, *Pt*=*ρQ*, must be larger than *Pd*; i.e., the turbulent disturbances must be large enough to overcome the surfacestabilizing forces. The volume of air entrained per unit time, *Va*, is given as Eq.(7).

$$\frac{\partial V\_a}{\partial\_t} + \nabla V\_a = Rt(1 - V\_a) \tag{7}$$

where *RC P P* = − *air t d* 2( ) / ρ , *u* is fluid velocity, *t* is time, and *C*air is a coefficient of proportionality : *C*air =0.5; i.e., assume on average that air will be trapped over about half the surface area. If *Pt* is less than *Pd* then *Va* is zero. Futhermore, Fig.11 is shown the meaasurement point of out flow rate and *Va*.

Fig. 11. Measurment plane of out flow rate and *Va* 

The air entrainment is expressed in Eq.(8).

$$A = \sum\_{k=1}^{n} V\_{ak} F\_{jk} V\_{jk} V\_{ck} \tag{8}$$

where *A* is the quantity of air entrainment, *Va* is the volume of air entrained per unit time, *Ff*  is the fluid fraction, *Vf* is the volume fraction, *Vc* is the volume of the mesh cell, and *n* is the aggregate number of mesh cells.

Fig.12 shows the flow of molten metal with the mold. It can also be seen the satisfactory liquid level control that swiftly fills the sprue cup is realized. Fig.13 shows the air entrainment at the maching surface. Upper figure shows using the conventional input,

σ

= + (6)

∂ (7)

, *u* is fluid velocity, *t* is time, and *C*air is a coefficient of

*P pgL <sup>d</sup> <sup>L</sup>*

where *ρ* is the liquid density, *σ* is the coefficient of surface tension, and *g* is the component of

For air entrainment to occur, the turbulent kinetic energy per unit volume, *Pt*=*ρQ*, must be larger than *Pd*; i.e., the turbulent disturbances must be large enough to overcome the surface-

(1 ) *<sup>a</sup> a a*

*<sup>V</sup> V Rt V* <sup>∂</sup> +∇ = −

proportionality : *C*air =0.5; i.e., assume on average that air will be trapped over about half the surface area. If *Pt* is less than *Pd* then *Va* is zero. Futhermore, Fig.11 is shown the

1

where *A* is the quantity of air entrainment, *Va* is the volume of air entrained per unit time, *Ff*  is the fluid fraction, *Vf* is the volume fraction, *Vc* is the volume of the mesh cell, and *n* is the

Fig.12 shows the flow of molten metal with the mold. It can also be seen the satisfactory liquid level control that swiftly fills the sprue cup is realized. Fig.13 shows the air entrainment at the maching surface. Upper figure shows using the conventional input,

*k A V FVV* =

*ak fk fk ck*

<sup>=</sup> (8)

*n*

stabilizing forces. The volume of air entrained per unit time, *Va*, is given as Eq.(7).

*t*

ρ

meaasurement point of out flow rate and *Va*.

Fig. 11. Measurment plane of out flow rate and *Va* 

The air entrainment is expressed in Eq.(8).

aggregate number of mesh cells.

gravity normal to the free surface.

where *RC P P* = − *air t d* 2( ) /

lower figure shows the optimum input, where the color of this figure is indicated the value of the air entrapment, and red color is the most air entarinment point. From the figure, the optimum input can be decrease the air entrainment, and it can be expect to realize a mitigation of defects such as pin holes and blow holes in experiment.

Fig. 12. Flow of molten metal with the mold

Optimization of Pouring Velocity for Aluminium Gravity Casting 587

Fig. 14. Experimental results with conventional input and optimum input

Fig. 15. Defect fraction of machining surface

Fig. 13. Air entrainment at machining surface

### **7. Pouring control experiment using optimum velocity input**

The desired optimum velocity input was applied to an actual machine and a pouring experiment was conducted. The conventional inputs were used to compare the optimal solution. Eleven trial runs were carried out and the frequency of defects appearing in the machined surface was evaluated after machining the product part from which the sprue, gating system and gate riser part had been removed. Fig. 14 shows the experiment results using each input and Fig. 15 shows the defect frequency for each.

As a result, it was possible to reduce the manufacturing defect to 1/6 rate by using an optimum pouring velocity as shown in Fig. 14 and Fig. 15. From this it was seen that the optimum input parameter derived using the proposed technique was valid for mitigating the occurrence of product defects.

Fig. 13. Air entrainment at machining surface

the occurrence of product defects.

**7. Pouring control experiment using optimum velocity input** 

using each input and Fig. 15 shows the defect frequency for each.

The desired optimum velocity input was applied to an actual machine and a pouring experiment was conducted. The conventional inputs were used to compare the optimal solution. Eleven trial runs were carried out and the frequency of defects appearing in the machined surface was evaluated after machining the product part from which the sprue, gating system and gate riser part had been removed. Fig. 14 shows the experiment results

As a result, it was possible to reduce the manufacturing defect to 1/6 rate by using an optimum pouring velocity as shown in Fig. 14 and Fig. 15. From this it was seen that the optimum input parameter derived using the proposed technique was valid for mitigating

Fig. 14. Experimental results with conventional input and optimum input

Fig. 15. Defect fraction of machining surface

**0**

**27**

*Italy*

Marco Marcon

*via Ponzio 34/5, 20133 Milano*

**Fluid Dynamics Without Fluids**

*Politecnico di Milano, Dipartimento di Elettronica e Informazione,*

This chapter will discuss some interesting real applications where Fluid Dynamics equations found fruitful applications without dealing with "strictly speaking" fluids. In particular, thanks to the large set of analyses performed over different kinds of fluids in different operating and boundary conditions, a wide range of Computational Fluid Dynamics algorithms flourished tackling different aspects, from convergence rate, to stability according to the discretization, to multigrid and linearization problems. This robust and thorough background, both on theoretical and on practical aspects, made Computational Fluid Dynamics (CFD) appealing also to other sciences and applications where Fluid Dynamics equations, or similar equations very close to them, can be useful in describing complex phenomena not related to fluids. Some applications that will be discussed concern, e.g., Geometry of liquid snowflakes whose contour is growing steered by curvature, staring from a circle. Furthermore Image Restoration and Segmentation can also benefit from CFD since a set of evolutionary algorithms, based on level-set curvature flow equations, plays a fundamental role in steering active contours or snakes through the noise present in the image till the complete warping of the desired framed object. Also in this case advanced techniques like Ghost Fluids Method for two competing fluids dynamics can be used to separate different objects in images. Other interesting applications that will be described concern applicability of CFD to surface extraction from cloud of points. This is a common problem when complex clouds of points, representing 3D objects or scenes are obtained by laser scanners or multi-camera vision systems. These points represent unambiguous features from corners or sharp edges and the final 3D closed surface must fit on these points smoothly interpolating empty space between them. Also in this case CFD can provide useful tools to define the evolution of a 3D surface representing the border between two competing fluids, one representing the "inside" and the other the "outside" of the object itself. The two fluids evolution will stop when surface sticks on all the 3D points: the viscosity of the two fluids will control the smoothness of this surface that will wrap the cloud and turbulence is used to model injection into grooves or narrow holes. This chapter will also discuss another interesting application of CFD to robotic navigation in complex environments where we are looking for the best path, both in terms of length and distance from objects, through a set of obstacles, different terrains traversability or path slope. Also in this case an imaginary fluid with a predefined viscosity floods from the robot position through the whole environment, its front

**1. Introduction**

## **8. Conclusion**

In this research, an analysis technique that enables a reduction in calculation time in fluid analysis simulations was proposed and pouring control input was optimized with the purpose of reducing occurrence of defects such as blow holes and pin holes in aluminum gravity casting.

As a result of optimizing pouring control input using GA, optimum pouring control input realized liquid level control and its validity in mitigating product defects was seen through a real machine pouring experiment.

## **9. References**


## **Fluid Dynamics Without Fluids**

Marco Marcon

*Politecnico di Milano, Dipartimento di Elettronica e Informazione, via Ponzio 34/5, 20133 Milano Italy*

#### **1. Introduction**

588 Fluid Dynamics, Computational Modeling and Applications

In this research, an analysis technique that enables a reduction in calculation time in fluid analysis simulations was proposed and pouring control input was optimized with the purpose of reducing occurrence of defects such as blow holes and pin holes in aluminum

As a result of optimizing pouring control input using GA, optimum pouring control input realized liquid level control and its validity in mitigating product defects was seen through

Yutaka K., Nobuaki N. & Hideaki O., (2001). Classification of Pin Hole Defect of Iron

Takuya S.,(2004). Prediction of Gas Defects by Mold-Filling Simulation with Consideration

Itsuo O., (2006). State of the Art of Simulation of Casting. *Journal of Japan Foundry Engineering* 

Takuya S. & Ichiyo N., (2007). Optimization of Injection Speed for Reduction of Cold Shut in

Ken'ichi Y., Koutarou H., Yoshifumi K.& Seishi N., (2008). "Optimum Velocity Control of

*Society*, Vol.,78, No.,12, (2006 Dec.) pp.602-608, 1342-0429

Casting. *Journal of Japan Foundry Engineering Society*, Vol.,73, No.,4, (2001 Apr.) pp.

of Surface Tension. *Journal of Japan Foundry Engineering Society*, Vol.,76, No.,7, (2004)

Die Casting, *Journal of Japan Foundry Engineering Society*, Vol.79, No.10,(2007 Oct.),

Die Casting Plunger Accounting for Air Entrapment and Shutting," *International Journal of Automation Technology*, Vol.2, No.4,(2008 Jun.) pp. 259-265, 1881-7629 Thomas B. & Hans- Paul S.,(1993). "An Overview of Evolutionary Algorithm for Parameter

Optimization, " *Evolutionary Computation*, Vol.1, No.1,(1993 Apr.) pp. 1-23, 1063-

**8. Conclusion** 

gravity casting.

**9. References** 

a real machine pouring experiment.

258-263, 1342-0429

pp. 562-569, 1342-0429

pp.592-600, 1342-0429

6560

This chapter will discuss some interesting real applications where Fluid Dynamics equations found fruitful applications without dealing with "strictly speaking" fluids. In particular, thanks to the large set of analyses performed over different kinds of fluids in different operating and boundary conditions, a wide range of Computational Fluid Dynamics algorithms flourished tackling different aspects, from convergence rate, to stability according to the discretization, to multigrid and linearization problems. This robust and thorough background, both on theoretical and on practical aspects, made Computational Fluid Dynamics (CFD) appealing also to other sciences and applications where Fluid Dynamics equations, or similar equations very close to them, can be useful in describing complex phenomena not related to fluids. Some applications that will be discussed concern, e.g., Geometry of liquid snowflakes whose contour is growing steered by curvature, staring from a circle. Furthermore Image Restoration and Segmentation can also benefit from CFD since a set of evolutionary algorithms, based on level-set curvature flow equations, plays a fundamental role in steering active contours or snakes through the noise present in the image till the complete warping of the desired framed object. Also in this case advanced techniques like Ghost Fluids Method for two competing fluids dynamics can be used to separate different objects in images. Other interesting applications that will be described concern applicability of CFD to surface extraction from cloud of points. This is a common problem when complex clouds of points, representing 3D objects or scenes are obtained by laser scanners or multi-camera vision systems. These points represent unambiguous features from corners or sharp edges and the final 3D closed surface must fit on these points smoothly interpolating empty space between them. Also in this case CFD can provide useful tools to define the evolution of a 3D surface representing the border between two competing fluids, one representing the "inside" and the other the "outside" of the object itself. The two fluids evolution will stop when surface sticks on all the 3D points: the viscosity of the two fluids will control the smoothness of this surface that will wrap the cloud and turbulence is used to model injection into grooves or narrow holes. This chapter will also discuss another interesting application of CFD to robotic navigation in complex environments where we are looking for the best path, both in terms of length and distance from objects, through a set of obstacles, different terrains traversability or path slope. Also in this case an imaginary fluid with a predefined viscosity floods from the robot position through the whole environment, its front

Fig. 1. A typical Tyndall figure (Liquid snowflake)

variation across the interface and depends from curvature:

to (Davis (2001)) a six fold symmetrical crystal will have a:

accordingly to (Pimpinelli & Villain (1999));

line segments (Dobrushin et al. (1992)).

*γ*<sup>6</sup> = *γ*<sup>0</sup>

Where *TM* is the equilibrium temperature of a planar solid-liquid interface (273.15*K* for water), *κ* is the curvature, *γ* is the surface energy per unit area for the solid-liquid interface, *ρ* is the density of a specific phase but in this case we are assuming equal density both for ice and water and *L* is the latent heat of fusion. The steering interface force is then due to the pressure

Fluid Dynamics Without Fluids 591

Where the subscript *l* and *s* refer to liquid and solid state respectively and *s* indicates the

A further improvement in the Gibbs-Thomson model considers surface energy anisotropy due to crystal ice structure (six fold symmetry), this will add a dependence of the surface energy *γ* from the relative orientation of the surface with respect to the underlying crystal. Accordingly

*σn*

Where *σn* represents the anisotropy degree and the complete interface temperature will be,

In particular, accordingly to equations 3 and 4, whenever |*σn*| < 1 we are in a week anisotropy condition, i.e. the surface energy *γ* is always positive, and, at equilibrium, the surface will resemble a smooth line, while, for anisotropy values greater that one, as for the snowflakes water example, the thermodynamic equilibrium implies a closed polygonal made of straight

Once the Temperature of the interface is determined the shape of the boundary between the liquid and solid part and its time evolution can be recovered from the energy conservation, in particular, in order for the ice-water interface to advance inside the solid part, two energy contributions must be provided: one for solid melting and one for interface stretching. The melting energy that must be provided per unit of time for a surface *A* that is moving inside

⎛ ⎝1 − *κ*

*<sup>n</sup>*<sup>2</sup> <sup>−</sup> <sup>1</sup> cos (*nθ*)

*<sup>γ</sup>*<sup>6</sup> <sup>+</sup> *<sup>d</sup>*2*γ*<sup>6</sup> *dθ*<sup>2</sup> *ρL*

� *TI* (*κ*) *TM*

�� � � � *n*=6

⎞

⎠ (4)

− 1 �

(2)

(3)

*ps* − *pl* = *ρ* (*sl* − *ss*) (*TI* − *TM*) = *ρL*

entropy per unit mass, for further details refer to (Hills & Roberts (1993)).

� 1 +

*TI* (*κ*) = *TM*

evolution speed, accordingly to CFD, will be slower in narrow passages and, once it reaches the target, it will define the easiest way.

#### **2. Snowflakes and phase interfaces**

In many physical problems outlines of a considered fluid have different speeds at different points, and in many cases local speed is directly connected to the curvature. In the following we will show how a geometrical contour, represented by a 2D line or a 3D surface of the evolving shape changes from point to point accordingly to the local geometrical properties of the contour itself. The presented approach can be applied to a generic N-dimensional case and we can refer to the moving contour as a general hypersurface, anyway in the following, we will limit without loss of generalization our considerations to a 3D space where fluid boundary is represented by a closed surface.

In particular Sethian (1989) showed how algorithms based on direct parametrization of the surface evolution could present ambiguous and error-prone solutions due to the local error propagation while global approaches based on implicit representations can result in a much more robust solution. This kind of approaches usually relies on higher dimensional functions for which the considered surface is just represented by a level-set. The contour motion can then be described applying the proper Hamilton-Jacobi equation related to a hyperbolic conservation law to the implicit function and tracking the particular level-set representing the contour itself.

Implicit function formulation is much more robust to numerical techniques with respect to direct methods and can be implemented starting from an original closed and non-intersecting surface: two famous physical examples are the combustion model of a flame and the grow of a snowflake: Markstein (1951) assumed that thin flame fingers close to their ends are cooler than their inner part and move slower that hot nonconvex regions: the proposed evolution equation assumes that flame contour speed is inversely proportional to its curvature *κ*, (where *κ* is positive for convex regions and negative for nonconvex ones). This assumption come from the fact that in highly convex regions, flame particles collide with a high number of slower and cooler surrounding particles that slow down their motion. The symmetrical effect is present in spikes growing in snowflakes. In particular these two last examples, flames and snowflakes, represent two very good examples where simple equations can explain very well complex phenomena. For a comprehensive description of ice crystal growth we suggest Libbrecht (2005), while in the following we describe how solid-liquid interface dynamics are steered by curvature in liquid snowflakes (or Tyndall figures). We refer to some recent studies (Hennessy (2010), Hobbs (2010)) on cylindrical discs of liquid in superheated crystals of ice. In this case the snowflake is made of water inside ice and its geometrical shape evolution is determined by boundary condition between solid and liquid phase. In Fig. 1 there is a typical snapshot of a Liquid snowflake.

#### **2.1 Liquid snowflakes surface evolution**

The typical modelization of snowflakes surface growing is the Gibbs-Thomson equation where the interface temperature *TI*, between the solid and the liquid part is ruled by the equation:

$$T\_I\left(\kappa\right) = T\_M \left(1 - \kappa \frac{\gamma}{\rho L}\right) \tag{1}$$

2 Will-be-set-by-IN-TECH

evolution speed, accordingly to CFD, will be slower in narrow passages and, once it reaches

In many physical problems outlines of a considered fluid have different speeds at different points, and in many cases local speed is directly connected to the curvature. In the following we will show how a geometrical contour, represented by a 2D line or a 3D surface of the evolving shape changes from point to point accordingly to the local geometrical properties of the contour itself. The presented approach can be applied to a generic N-dimensional case and we can refer to the moving contour as a general hypersurface, anyway in the following, we will limit without loss of generalization our considerations to a 3D space where fluid

In particular Sethian (1989) showed how algorithms based on direct parametrization of the surface evolution could present ambiguous and error-prone solutions due to the local error propagation while global approaches based on implicit representations can result in a much more robust solution. This kind of approaches usually relies on higher dimensional functions for which the considered surface is just represented by a level-set. The contour motion can then be described applying the proper Hamilton-Jacobi equation related to a hyperbolic conservation law to the implicit function and tracking the particular level-set representing the

Implicit function formulation is much more robust to numerical techniques with respect to direct methods and can be implemented starting from an original closed and non-intersecting surface: two famous physical examples are the combustion model of a flame and the grow of a snowflake: Markstein (1951) assumed that thin flame fingers close to their ends are cooler than their inner part and move slower that hot nonconvex regions: the proposed evolution equation assumes that flame contour speed is inversely proportional to its curvature *κ*, (where *κ* is positive for convex regions and negative for nonconvex ones). This assumption come from the fact that in highly convex regions, flame particles collide with a high number of slower and cooler surrounding particles that slow down their motion. The symmetrical effect is present in spikes growing in snowflakes. In particular these two last examples, flames and snowflakes, represent two very good examples where simple equations can explain very well complex phenomena. For a comprehensive description of ice crystal growth we suggest Libbrecht (2005), while in the following we describe how solid-liquid interface dynamics are steered by curvature in liquid snowflakes (or Tyndall figures). We refer to some recent studies (Hennessy (2010), Hobbs (2010)) on cylindrical discs of liquid in superheated crystals of ice. In this case the snowflake is made of water inside ice and its geometrical shape evolution is determined by boundary condition between solid and liquid phase. In Fig. 1 there is a typical snapshot of

The typical modelization of snowflakes surface growing is the Gibbs-Thomson equation where the interface temperature *TI*, between the solid and the liquid part is ruled by the

<sup>1</sup> <sup>−</sup> *<sup>κ</sup> <sup>γ</sup> ρL* 

(1)

*TI* (*κ*) = *TM*

the target, it will define the easiest way.

**2. Snowflakes and phase interfaces**

boundary is represented by a closed surface.

contour itself.

a Liquid snowflake.

equation:

**2.1 Liquid snowflakes surface evolution**

Fig. 1. A typical Tyndall figure (Liquid snowflake)

Where *TM* is the equilibrium temperature of a planar solid-liquid interface (273.15*K* for water), *κ* is the curvature, *γ* is the surface energy per unit area for the solid-liquid interface, *ρ* is the density of a specific phase but in this case we are assuming equal density both for ice and water and *L* is the latent heat of fusion. The steering interface force is then due to the pressure variation across the interface and depends from curvature:

$$p\_s - p\_l = \rho \left(s\_l - s\_s\right) \left(T\_I - T\_M\right) = \rho L \left(\frac{T\_I \left(\kappa\right)}{T\_M} - 1\right) \tag{2}$$

Where the subscript *l* and *s* refer to liquid and solid state respectively and *s* indicates the entropy per unit mass, for further details refer to (Hills & Roberts (1993)).

A further improvement in the Gibbs-Thomson model considers surface energy anisotropy due to crystal ice structure (six fold symmetry), this will add a dependence of the surface energy *γ* from the relative orientation of the surface with respect to the underlying crystal. Accordingly to (Davis (2001)) a six fold symmetrical crystal will have a:

$$\gamma\_6 = \gamma\_0 \left[ 1 + \frac{\sigma\_n}{n^2 - 1} \cos \left( n\theta \right) \right] \Big|\_{n=6} \tag{3}$$

Where *σn* represents the anisotropy degree and the complete interface temperature will be, accordingly to (Pimpinelli & Villain (1999));

$$T\_I\left(\kappa\right) = T\_M \left(1 - \kappa \frac{\gamma\_\ell + \frac{d^2 \gamma\_\ell}{d\theta^2}}{\rho L}\right) \tag{4}$$

In particular, accordingly to equations 3 and 4, whenever |*σn*| < 1 we are in a week anisotropy condition, i.e. the surface energy *γ* is always positive, and, at equilibrium, the surface will resemble a smooth line, while, for anisotropy values greater that one, as for the snowflakes water example, the thermodynamic equilibrium implies a closed polygonal made of straight line segments (Dobrushin et al. (1992)).

Once the Temperature of the interface is determined the shape of the boundary between the liquid and solid part and its time evolution can be recovered from the energy conservation, in particular, in order for the ice-water interface to advance inside the solid part, two energy contributions must be provided: one for solid melting and one for interface stretching. The melting energy that must be provided per unit of time for a surface *A* that is moving inside

respect to the closed surface. The discrete formulation requires that the value of *φ* is initially evaluated over all the mesh points and the steering equations would require the update of each point at each time step. This computational complexity of order *O*(*n*3) at each time step could be unacceptable in many practical cases so many faster update solutions have been

Fluid Dynamics Without Fluids 593

• Narrow band: the updated mesh points are only those in a small stripe surrounding the interface, this will reduce computation to *O*(*n*2) but as soon as the interface is approaching the stripe boundary all values must be recomputed (more details in

• Octree: relevant mesh points and high detailed surface regions can be divided in sub-grids, e.g. a cube can be split in eight sub-cubes (the method name come from here) and the procedure can be iterated for higher precision. More details can be found in (Losasso et al.

• Sparse block grid: The whole domain is subdivided in blocks representing clusters of the original grid points (they can be overlapped or not) then the update is evaluated only on

Accordingly to the previous description the set of points belonging to the surface will satisfy the zero level-set condition:*φ* (**p**(*t*), *t*) = 0 where **p**(*t*) represents the path of a surface point.

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> <sup>∇</sup>*<sup>φ</sup>* ·

Considering the normal unitary vector to a generic point of *φ* for a generic level-set (contour),

*<sup>v</sup>* (*t*) <sup>=</sup> *<sup>d</sup>***<sup>p</sup>** (*t*)

In the cases considered in the following and in many physical cases the interface evolution is ruled by a Hamilton-Jacobi equation and the local velocity can be assumed proportional to

*<sup>y</sup>* <sup>−</sup> <sup>2</sup>*φyφxφxy* <sup>+</sup> *<sup>φ</sup>yyφ*<sup>2</sup>

*d***p**(*t*)

*dt* <sup>=</sup> 0 (8)

*dt* · **<sup>n</sup>** (9)

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>v</sup>* (*t*)|∇*φ*<sup>|</sup> <sup>=</sup> <sup>0</sup> (10)

*v* = *v*<sup>0</sup> − *�κ* (11)

3/2 (12)

*x*

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>*v*<sup>0</sup> |∇*φ*<sup>|</sup> <sup>+</sup> *�κ* |∇*φ*<sup>|</sup> (13)

blocks containing the surface narrow band. (Further details in (Bridson (2003))

*dφ* (**p**(*t*),*t*)

*dφ* (**p**(*t*),*t*)

*<sup>κ</sup>* <sup>=</sup> *<sup>φ</sup>xxφ*<sup>2</sup>

*∂φ*

 *φ*2 *<sup>x</sup>* + *φ*<sup>2</sup> *y*

*dt* <sup>=</sup> *∂φ*


*dt* <sup>=</sup> *∂φ*

proposed; the most common are:

(2004))

**n** = <sup>∇</sup>*<sup>φ</sup>*

obtaining:

.

**3.1 Level-set evolution**

(Adalsteinsson & Sethain (1995))

The first temporal derivative will become:

the local curvature of the surface itself:

where *�* is a constant and

We can rewrite *φ* as follows:

the solid at a normal speed *vn* is *lmel* = *ρLAvn* while the energy contribution per unit of time due to surface stretching is *lsur* = *γ dA dt* ; which, accordingly to the differential geometry can be written as:

$$d\_{sur} = \gamma \frac{dA}{dt} = A \gamma \kappa v\_n \tag{5}$$

The surface advancement is steered by two opposite heat diffusive fluxes: one inside the liquid part and one inside the solid part; respectively: −*kl*∇*Tl* and *ks*∇*Ts* where *k* indicates the thermal conductivity of the considered phase.

The equilibrium of energy equations per area and time unit can be written as:

$$\frac{l\_{mel}}{A} + \frac{l\_{sur}}{A} = (-k\_l \nabla T\_l + k\_s \nabla T\_s) \cdot \mathbf{n} \tag{6}$$

i.e.:

$$\left(\rho L + \gamma \kappa\right) v\_{\text{\tiny{\tiny{\Gamma}}}} = \left(-k\_{\text{I}} \nabla T\_{\text{I}} + k\_{\text{s}} \nabla T\_{\text{s}}\right) \cdot \mathbf{n} \tag{7}$$

Where **n** is surface normal oriented from the liquid part towards the solid one; this equation is called the Stefan condition. The numerical solution for the liquid-solid interface evolution, accordingly to the Stefan condition, requires some particular precaution due to accuracy that must be placed in surface tracking. Since the analyzed problem is far from thermodynamic equilibrium due to the superheated temperature of the solid, some typical approaches, e.g. the Enthalpy model (Shastri & Allen (1998)), are not accurate. In fact, Enthalpy approach, even if requires only the evaluation of a single parabolic PDE at each time step, assumes that the interface is within a mushy solid-liquid region and the evolving interface is just approximated. A possible extension of Enthalpy approach to account for superheated solid state could be phase-field models; anyway an accurate interface tracking requires a very small grid spacing resulting in a high computational cost (Biben et al. (2005)). One of the method that in recent years is getting growing attention is Level-set: in the following of this chapter we will show it versatility and accuracy in many CFD applications and in liquid snowflakes shape definition.

#### **3. Level-set methods**

For most of real CFD systems symbolic formulation that easily yield to symbolic mathematical solution is unfeasible, and the goal becomes to find the system modelization with the best trade-off between accuracy and computational complexity accounting for local, global and independent properties of the moving surface. In particular we will refer to "local" as properties associated to local surface features (e.g. curvature, local density or temperature), "global" as properties associated to overall surface shape or extension, while independent properties are not linked to the surface itself, like a flow below the surface. The Level-set method aims to track a thin propagating interface over time, it outperforms many other approaches handling topological complexities such as corners and cusps, and in handling complexities in the evolving interface such as entropy conditions and weak solutions. A detailed description of level-set can be found in (Sethian (1999)).

The general idea is to define a function *φ* (*p*, *t*) over the whole space-time domain, where *p* indicates a point in the space and *t* represents the time variable. The initial surface is the zero level-set of an implicit function for *t* = 0: i.e. the starting surface is the set of points *p* so that *φ* (*p*, *t*) = 0 the value of the *φ* function for non-interface points is, the signed distance from the surface where 'signed' means that different signs indicate internal or external position with 4 Will-be-set-by-IN-TECH

the solid at a normal speed *vn* is *lmel* = *ρLAvn* while the energy contribution per unit of time

*dA*

The surface advancement is steered by two opposite heat diffusive fluxes: one inside the liquid part and one inside the solid part; respectively: −*kl*∇*Tl* and *ks*∇*Ts* where *k* indicates the

Where **n** is surface normal oriented from the liquid part towards the solid one; this equation is called the Stefan condition. The numerical solution for the liquid-solid interface evolution, accordingly to the Stefan condition, requires some particular precaution due to accuracy that must be placed in surface tracking. Since the analyzed problem is far from thermodynamic equilibrium due to the superheated temperature of the solid, some typical approaches, e.g. the Enthalpy model (Shastri & Allen (1998)), are not accurate. In fact, Enthalpy approach, even if requires only the evaluation of a single parabolic PDE at each time step, assumes that the interface is within a mushy solid-liquid region and the evolving interface is just approximated. A possible extension of Enthalpy approach to account for superheated solid state could be phase-field models; anyway an accurate interface tracking requires a very small grid spacing resulting in a high computational cost (Biben et al. (2005)). One of the method that in recent years is getting growing attention is Level-set: in the following of this chapter we will show it versatility and accuracy in many CFD applications and in liquid snowflakes shape definition.

For most of real CFD systems symbolic formulation that easily yield to symbolic mathematical solution is unfeasible, and the goal becomes to find the system modelization with the best trade-off between accuracy and computational complexity accounting for local, global and independent properties of the moving surface. In particular we will refer to "local" as properties associated to local surface features (e.g. curvature, local density or temperature), "global" as properties associated to overall surface shape or extension, while independent properties are not linked to the surface itself, like a flow below the surface. The Level-set method aims to track a thin propagating interface over time, it outperforms many other approaches handling topological complexities such as corners and cusps, and in handling complexities in the evolving interface such as entropy conditions and weak solutions. A

The general idea is to define a function *φ* (*p*, *t*) over the whole space-time domain, where *p* indicates a point in the space and *t* represents the time variable. The initial surface is the zero level-set of an implicit function for *t* = 0: i.e. the starting surface is the set of points *p* so that *φ* (*p*, *t*) = 0 the value of the *φ* function for non-interface points is, the signed distance from the surface where 'signed' means that different signs indicate internal or external position with

*lsur* = *γ*

The equilibrium of energy equations per area and time unit can be written as:

*lsur*

*lmel <sup>A</sup>* <sup>+</sup>

detailed description of level-set can be found in (Sethian (1999)).

*dt* ; which, accordingly to the differential geometry can be

*<sup>A</sup>* <sup>=</sup> (−*kl*∇*Tl* <sup>+</sup> *ks*∇*Ts*) · **<sup>n</sup>** (6)

(*ρL* + *γκ*) *vn* = (−*kl*∇*Tl* + *ks*∇*Ts*) · **n** (7)

*dt* <sup>=</sup> *<sup>A</sup>γκvn* (5)

due to surface stretching is *lsur* = *γ dA*

thermal conductivity of the considered phase.

written as:

i.e.:

**3. Level-set methods**

respect to the closed surface. The discrete formulation requires that the value of *φ* is initially evaluated over all the mesh points and the steering equations would require the update of each point at each time step. This computational complexity of order *O*(*n*3) at each time step could be unacceptable in many practical cases so many faster update solutions have been proposed; the most common are:


#### **3.1 Level-set evolution**

Accordingly to the previous description the set of points belonging to the surface will satisfy the zero level-set condition:*φ* (**p**(*t*), *t*) = 0 where **p**(*t*) represents the path of a surface point. The first temporal derivative will become:

$$\frac{d\phi\left(\mathbf{p}(t),t\right)}{dt} = \frac{\partial\phi}{\partial t} + \nabla\phi \cdot \frac{d\mathbf{p}(t)}{dt} = 0\tag{8}$$

Considering the normal unitary vector to a generic point of *φ* for a generic level-set (contour), **n** = <sup>∇</sup>*<sup>φ</sup>* |∇*φ*<sup>|</sup> the speed function *<sup>v</sup>* along this direction is then:

$$v\left(t\right) = \frac{d\mathbf{p}\left(t\right)}{dt} \cdot \mathbf{n} \tag{9}$$

obtaining:

.

$$\frac{d\phi\left(\mathbf{p}\left(t\right),t\right)}{dt} = \frac{\partial\phi}{\partial t} + v\left(t\right)\left|\nabla\phi\right| = 0\tag{10}$$

In the cases considered in the following and in many physical cases the interface evolution is ruled by a Hamilton-Jacobi equation and the local velocity can be assumed proportional to the local curvature of the surface itself:

$$
v = v\_0 - \epsilon \kappa \tag{11}$$

where *�* is a constant and

$$\kappa = \frac{\phi\_{xx}\phi\_y^2 - 2\phi\_y\phi\_x\phi\_{xy} + \phi\_{yy}\phi\_x^2}{\left(\phi\_x^2 + \phi\_y^2\right)^{3/2}}\tag{12}$$

We can rewrite *φ* as follows:

$$\frac{\partial \phi}{\partial t} = -v\_0 \left| \nabla \phi \right| + \epsilon \kappa \left| \nabla \phi \right|\tag{13}$$

contour equal to a circle of radius *<sup>R</sup>* the curvature will be equal to *<sup>κ</sup>* <sup>=</sup> <sup>−</sup> <sup>1</sup>

analytical integration of eq. 16 refer to (Hennessy (2010)).

are presented for the larger curve of the raw above

slow down.

of the level-set accordingly to the aforemention equations is represented in Fig. 3 where a very good agreement of the zero level-set with real data is reached. For further details about

Fluid Dynamics Without Fluids 595

Fig. 3. on the top raw: different snapshots of the zero level set evolution varying the relative absorbition *αr*; instead on the bottom line different level set values (different temperatures)

Fluid dynamics and Level-set formulation can also be fruitfully adopted in Motion planning Algorithms. In particular the flexibility of the presented technique for time optimal motion planning with moving obstacles can be easily integrated with further constraints like different terrains traversability, path narrowing or fuel economy. Many classical techniques are present in literature for shortest path search with stationary obstacles, a good survey can be found in Latombe (1991). The space-time representation of obstacles moving in a 2D environment, assuming that one axis will be the time variable will be a solid structure. In Fig. 4 there is the representation of a circular obstacle moving in the*x-y* plane along *x* direction with an abrupt

As shown in Kimmel et al. (1998) the search for the time-optimal path within the aforementioned 3D time-space can be reduced to the search in a reduced set of 2D regions forming a Multivalued Distance Map (MDM). A Distance Map from a point **p** is a function that assigns to each point **q** in the considered domain a value corresponding to the "minimal length" (geodesic). The simplest approach uses a contour with a constant speed propagating from a small circle around **p** and the distance assigned to each crossed point **q** is proportional to the propagation time. If obstacles are moving the optimal path will be the minimal geodesic only when it does not cross any moving obstacles in the scene. For all the other cases a

**4. Motion planning for traveling robots in complex environments**

*<sup>R</sup>* The time evolution

#### **3.2 Liquid snowflakes surface evolution**

Assuming that the ice is radiated with an external electromagnetic radiation normal to its surface ( a focused beam with a gaussian shape) the heat will follow an exponential decay penetration, accordingly to the following heat propagation equation:

$$I\left(r,z\right) = I\_0 e^{-\left(r/r\_0\right)^2 - az} \tag{14}$$

where *I*<sup>0</sup> is the central beam intensity, *rb* is its standard deviation while *r* is the in-plane coordinate for the ice surface. Furthermore *z* is the coordinate orthogonal to the ice surface and *α* is a penetration coefficient. The formation of the liquid snowflake is then in the in-plane direction. A level-set function can the be defined where *φ* is the normalized temperature that follows heat PDE:

$$\begin{cases} \begin{array}{c} \frac{\partial \phi}{\partial t} = \nabla^2 \phi + e^{-\left(r/\beta\right)^2} \\ c\_{pr} \frac{\partial \phi}{\partial t} = k\_r \nabla^2 \phi + \alpha\_r e^{-\left(r/\beta\right)^2} \end{array} \text{for the liquid} \tag{15}$$

where *cpr* <sup>=</sup> *cps cpl* is the relative measure of heat capacities, *kr* <sup>=</sup> *ks kl* is the relative thermal conductivity and *α<sup>r</sup>* = *<sup>α</sup><sup>s</sup> <sup>α</sup><sup>l</sup>* is the relative absorbtion coefficient.

The liquid/solid interface is then associated to the zero level-set in a sort of rescaled Celsius scale. Experimental observations show that liquid snow flakes often begin as circular discs which then grow outwards in a radially symmetric manner. However, after a certain amount of time, the interface becomes unstable and small, sinusoidal perturbations with a well defined wave-number appear. In Fig. 2 it is possible to see the growing of a liquid snowflake starting from a circle and then evolving with equally spaced rippling.

Fig. 2. Liquid snow flakes formation starting from a circle and the evolving, accordingly to linear instabilities with sinusoidal perturbations

The eq. 15, accordingly to the problem symmetry can be written in polar coordinates for the quasi-steadily condition:

$$\begin{cases} \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial \phi}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 \phi}{\partial \theta^2} + e^{-\left(r/\beta\right)^2} = 0 & \text{for the liquid} \\\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial \phi}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 \phi}{\partial \theta^2} + \frac{a\_r}{k\_r}e^{-\left(r/\beta\right)^2} = 0 & \text{for the solid} \end{cases} \tag{16}$$

Even small linear perturbations must satisfy the equations above and an integer number of oscillations (first order approximation) will constitute the zero level-set. Assuming a starting 6 Will-be-set-by-IN-TECH

Assuming that the ice is radiated with an external electromagnetic radiation normal to its surface ( a focused beam with a gaussian shape) the heat will follow an exponential decay

2

2

2

*for the liquid*

<sup>−</sup>*α<sup>z</sup>* (14)

*for the solid* (15)

*kl* is the relative thermal

*I* (*r*, *z*) = *I*0*e*−(*r*/*rb*)

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>∇</sup>2*<sup>φ</sup>* <sup>+</sup> *<sup>e</sup>*−(*r*/*β*)

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> *kr*∇2*<sup>φ</sup>* <sup>+</sup> *<sup>α</sup>re*−(*r*/*β*)

*cpl* is the relative measure of heat capacities, *kr* <sup>=</sup> *ks*

*<sup>α</sup><sup>l</sup>* is the relative absorbtion coefficient.

The liquid/solid interface is then associated to the zero level-set in a sort of rescaled Celsius scale. Experimental observations show that liquid snow flakes often begin as circular discs which then grow outwards in a radially symmetric manner. However, after a certain amount of time, the interface becomes unstable and small, sinusoidal perturbations with a well defined wave-number appear. In Fig. 2 it is possible to see the growing of a liquid snowflake starting

Fig. 2. Liquid snow flakes formation starting from a circle and the evolving, accordingly to

The eq. 15, accordingly to the problem symmetry can be written in polar coordinates for the

*∂θ*<sup>2</sup> <sup>+</sup> *<sup>e</sup>*−(*r*/*β*)

Even small linear perturbations must satisfy the equations above and an integer number of oscillations (first order approximation) will constitute the zero level-set. Assuming a starting

2

= 0 *for the liquid*

<sup>=</sup> <sup>0</sup> *for the solid* (16)

where *I*<sup>0</sup> is the central beam intensity, *rb* is its standard deviation while *r* is the in-plane coordinate for the ice surface. Furthermore *z* is the coordinate orthogonal to the ice surface and *α* is a penetration coefficient. The formation of the liquid snowflake is then in the in-plane direction. A level-set function can the be defined where *φ* is the normalized temperature that

penetration, accordingly to the following heat propagation equation:

� *∂φ*

*cpr ∂φ*

from a circle and then evolving with equally spaced rippling.

linear instabilities with sinusoidal perturbations

⎧ ⎨ ⎩ 1 *r ∂ ∂r* � *r ∂φ ∂r* � + <sup>1</sup> *r*2 *∂*2*φ*

1 *r ∂ ∂r* � *r ∂φ ∂r* � + <sup>1</sup> *r*2 *∂*2*φ ∂θ*<sup>2</sup> <sup>+</sup> *<sup>α</sup><sup>r</sup> kr e*−(*r*/*β*) 2

**3.2 Liquid snowflakes surface evolution**

follows heat PDE:

where *cpr* <sup>=</sup> *cps*

conductivity and *α<sup>r</sup>* = *<sup>α</sup><sup>s</sup>*

quasi-steadily condition:

contour equal to a circle of radius *<sup>R</sup>* the curvature will be equal to *<sup>κ</sup>* <sup>=</sup> <sup>−</sup> <sup>1</sup> *<sup>R</sup>* The time evolution of the level-set accordingly to the aforemention equations is represented in Fig. 3 where a very good agreement of the zero level-set with real data is reached. For further details about analytical integration of eq. 16 refer to (Hennessy (2010)).

Fig. 3. on the top raw: different snapshots of the zero level set evolution varying the relative absorbition *αr*; instead on the bottom line different level set values (different temperatures) are presented for the larger curve of the raw above

#### **4. Motion planning for traveling robots in complex environments**

Fluid dynamics and Level-set formulation can also be fruitfully adopted in Motion planning Algorithms. In particular the flexibility of the presented technique for time optimal motion planning with moving obstacles can be easily integrated with further constraints like different terrains traversability, path narrowing or fuel economy. Many classical techniques are present in literature for shortest path search with stationary obstacles, a good survey can be found in Latombe (1991). The space-time representation of obstacles moving in a 2D environment, assuming that one axis will be the time variable will be a solid structure. In Fig. 4 there is the representation of a circular obstacle moving in the*x-y* plane along *x* direction with an abrupt slow down.

As shown in Kimmel et al. (1998) the search for the time-optimal path within the aforementioned 3D time-space can be reduced to the search in a reduced set of 2D regions forming a Multivalued Distance Map (MDM). A Distance Map from a point **p** is a function that assigns to each point **q** in the considered domain a value corresponding to the "minimal length" (geodesic). The simplest approach uses a contour with a constant speed propagating from a small circle around **p** and the distance assigned to each crossed point **q** is proportional to the propagation time. If obstacles are moving the optimal path will be the minimal geodesic only when it does not cross any moving obstacles in the scene. For all the other cases a

*∂*P (*α* (*t*))

Fluid Dynamics Without Fluids 597

Fig. 5. The time evolution of variable *α*, and the definition of the moving front velocity *VN*

**<sup>N</sup>** <sup>=</sup> (−*p*, <sup>−</sup>*q*, 1)

Since the minimal geodesics are perpendicular to the equal distance contours on the surface, we can track the optimal trajectory starting from the destination point and, considering the equal distance contour *α* crossing it we have to move backward thanks to the orthogonality

*<sup>x</sup>* + (1 + *p*2) *n*<sup>2</sup>

*<sup>∂</sup><sup>y</sup>* , accordingly to (Kimmel et al. (1995)) we can write:

<sup>1</sup> <sup>+</sup> *<sup>p</sup>*<sup>2</sup> <sup>+</sup> *<sup>q</sup>*<sup>2</sup> (21)

<sup>1</sup> <sup>+</sup> *<sup>p</sup>*<sup>2</sup> <sup>+</sup> *<sup>q</sup>*<sup>2</sup> (22)

*η* = − (**N** × **T**) (23)

*<sup>l</sup>*<sup>2</sup> (<sup>1</sup> <sup>+</sup> *<sup>p</sup>*2) <sup>+</sup> *<sup>r</sup>*<sup>2</sup> (<sup>1</sup> <sup>+</sup> *<sup>q</sup>*2) <sup>−</sup> <sup>2</sup>*pqrl* (24)

<sup>−</sup> *pqr*, *pr* <sup>+</sup> *ql*

*<sup>∂</sup><sup>x</sup>* , *<sup>l</sup>* <sup>=</sup> *<sup>∂</sup>MDMp*

*<sup>∂</sup><sup>y</sup>* and

(25)

*<sup>y</sup>* − 2*pqnxny*

Calling *p* = *<sup>∂</sup>t*(*x*,*y*)

using a normalization function

surface S, the best path *β* so that: *∂β*

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>η</sup>* <sup>=</sup> *<sup>ψ</sup>* (*p*, *<sup>q</sup>*,*r*, *<sup>l</sup>*)

and:

*<sup>∂</sup><sup>x</sup>* and *<sup>q</sup>* <sup>=</sup> *<sup>∂</sup>t*(*x*,*y*)

*VN* =

between the optimal path and the equal distance contour:

(1 + *q*2) *n*<sup>2</sup>

Using an arc length parametrization for the curve *<sup>α</sup>* we can write *<sup>r</sup>* = *<sup>∂</sup>MDMp*

 *r* 1 + *q*<sup>2</sup> 

*<sup>ψ</sup>* (*p*, *<sup>q</sup>*,*r*, *<sup>l</sup>*) <sup>=</sup> *<sup>p</sup>*<sup>2</sup> <sup>+</sup> *<sup>q</sup>*<sup>2</sup> <sup>+</sup> <sup>1</sup>

and the backtracking rule will then consist in tracking, from destination point **q** over the

− *pql*, *l*

 1 + *p*<sup>2</sup> 

In Fig. 5 there is a graphical representation of the aforementioned quantities.

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>*VN***<sup>n</sup>** (20)

Fig. 4. a space-time representation for a circular obstacle moving in the *x-y* plane along the x direction that slow down its speed

more sophisticated approach is required: MDM will keep track of multiple crossing of the propagating contour when moving obstacles are present: in particular *MDMp* (**q**, **1**) will indicate the first time at which the point **q** is reached from **p** while *MDMp* (**q**, **2**) will indicate the second time at which, accordingly to the moving objects the propagating contour can reach **q**, and so on. We will call N (**q**) the number of times that the propagating front will cross that point. Once the *MDM* is built for all points in the domain, the best path is obtained using a back-track procedure.

Consider now, given the source point **p**, a continuous surface S in the space-time domain, centered in **p** indicating the wavefront contour as a function of time: S = {(*x*, *y*, *t*(*x*, *y*))}, where a contour at a specific geodesic distance ¯*d*, *α* ¯*d* is given by:

$$\mathfrak{a}\begin{pmatrix}\vec{d}\end{pmatrix} = \begin{Bmatrix}\mathbf{q}|d\_{\mathcal{S}}\begin{pmatrix}\mathbf{p},\mathbf{q}\end{pmatrix} = \vec{d}\end{Bmatrix} \tag{17}$$

where *dg* (**p**, **q**) is the minimal geodesic distance along the surface S. Accordingly to Kimmel et al. (1995), due to the continuity and smoothness of S we can define a curvilinear abscissa *u* for *α u*, ¯*d* .

If we assume that the motion velocity is constant, *V* without loss of generality we can also assume that *V* = 1 then the geodesic distance *d* = *dg* (**p**, **q**) coincides with time *t* and *α* (*u*, *d*) = *α* (*u*, *t*) the equal distance contour evolution rule is:

$$\frac{\partial \boldsymbol{\alpha}\left(\boldsymbol{u},t\right)}{\partial t} = \mathbf{N} \times \mathbf{T} \tag{18}$$

where **T** is the unit vector tangent to *α* while **N** is the unit vector normal to the surface S. To avoid point discontinuities the starting contour *α* (*u*, 0) is a small circle around the point **p** and the projection of equation 18 on the *x-y* plane will be the velocity vector orthogonal to the actual front, *VN*:

$$V\_N = \mathcal{P}\left(\mathbf{N} \times \mathbf{T}\right) \cdot \mathbf{n} \tag{19}$$

where <sup>P</sup>() is the projector on the x-y plane and **<sup>n</sup>** <sup>=</sup> *nx*, *ny* is the unit normal to the planar curve in the x-y plane.

The evolution of the planar contour will then be:

8 Will-be-set-by-IN-TECH

Fig. 4. a space-time representation for a circular obstacle moving in the *x-y* plane along the x

more sophisticated approach is required: MDM will keep track of multiple crossing of the propagating contour when moving obstacles are present: in particular *MDMp* (**q**, **1**) will indicate the first time at which the point **q** is reached from **p** while *MDMp* (**q**, **2**) will indicate the second time at which, accordingly to the moving objects the propagating contour can reach **q**, and so on. We will call N (**q**) the number of times that the propagating front will cross that point. Once the *MDM* is built for all points in the domain, the best path is obtained using a

Consider now, given the source point **p**, a continuous surface S in the space-time domain, centered in **p** indicating the wavefront contour as a function of time: S = {(*x*, *y*, *t*(*x*, *y*))},

where *dg* (**p**, **q**) is the minimal geodesic distance along the surface S. Accordingly to Kimmel et al. (1995), due to the continuity and smoothness of S we can define a curvilinear

If we assume that the motion velocity is constant, *V* without loss of generality we can also assume that *V* = 1 then the geodesic distance *d* = *dg* (**p**, **q**) coincides with time *t* and *α* (*u*, *d*) =

*∂α* (*u*,*t*)

where **T** is the unit vector tangent to *α* while **N** is the unit vector normal to the surface S. To avoid point discontinuities the starting contour *α* (*u*, 0) is a small circle around the point **p** and the projection of equation 18 on the *x-y* plane will be the velocity vector orthogonal to the

 ¯*d* 

**<sup>q</sup>**|*dg* (**p**, **<sup>q</sup>**) <sup>=</sup> ¯*<sup>d</sup>*

is given by:

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> **<sup>N</sup>** <sup>×</sup> **<sup>T</sup>** (18)

*VN* = P (**N** × **T**) · **n** (19)

*nx*, *ny* 

(17)

is the unit normal to the planar

direction that slow down its speed

where a contour at a specific geodesic distance ¯*d*, *α*

*α* (*u*, *t*) the equal distance contour evolution rule is:

where <sup>P</sup>() is the projector on the x-y plane and **<sup>n</sup>** <sup>=</sup>

The evolution of the planar contour will then be:

*α* ¯*d* =

back-track procedure.

abscissa *u* for *α*

actual front, *VN*:

curve in the x-y plane.

 *u*, ¯*d* .

$$\frac{\partial \mathcal{P}\left(\mathbf{a}\left(t\right)\right)}{\partial t} = -V\_{\rm N}\mathbf{n} \tag{20}$$

In Fig. 5 there is a graphical representation of the aforementioned quantities.

Fig. 5. The time evolution of variable *α*, and the definition of the moving front velocity *VN* Calling *p* = *<sup>∂</sup>t*(*x*,*y*) *<sup>∂</sup><sup>x</sup>* and *<sup>q</sup>* <sup>=</sup> *<sup>∂</sup>t*(*x*,*y*) *<sup>∂</sup><sup>y</sup>* , accordingly to (Kimmel et al. (1995)) we can write:

$$\mathbf{N} = \frac{(-p\_\prime - q\_\prime 1)}{\sqrt{1 + p^2 + q^2}} \tag{21}$$

and:

$$V\_N = \sqrt{\frac{(1+q^2)\ n\_x^2 + (1+p^2)\ n\_y^2 - 2pqn\_xn\_y}{1+p^2+q^2}}\tag{22}$$

Since the minimal geodesics are perpendicular to the equal distance contours on the surface, we can track the optimal trajectory starting from the destination point and, considering the equal distance contour *α* crossing it we have to move backward thanks to the orthogonality between the optimal path and the equal distance contour:

$$\eta = -\left(\mathbf{N} \times \mathbf{T}\right) \tag{23}$$

Using an arc length parametrization for the curve *<sup>α</sup>* we can write *<sup>r</sup>* = *<sup>∂</sup>MDMp <sup>∂</sup><sup>x</sup>* , *<sup>l</sup>* <sup>=</sup> *<sup>∂</sup>MDMp <sup>∂</sup><sup>y</sup>* and using a normalization function

$$\Psi\left(p,q,r,l\right) = \frac{p^2 + q^2 + 1}{\sqrt{l^2\left(1+p^2\right) + r^2\left(1+q^2\right) - 2pqrl}}\tag{24}$$

and the backtracking rule will then consist in tracking, from destination point **q** over the surface S, the best path *β* so that:

$$\frac{\partial \beta}{\partial t} = \eta = \psi \left( p, q, r, l \right) \left[ r \left( 1 + q^2 \right) - p q l, l \left( 1 + p^2 \right) - p q r, p r + q l \right] \tag{25}$$

If there are moving obstacles along the path the *MDMp* becomes multivalued and it could be possible that some points change their index along the optimal path, to account for this all distance values must be evaluated at a candidate point.

#### **4.1 The Level-set approach**

Moving obstacles introduce new boundary constraints and dealing with them could be very complex following the previously describer approach, a simpler way could be using the level-set approach where the projection of *α* on the x-y plane is the zero level set *φ* = 0 and *φ* will be negative in the interior and positive in the exterior of the zero level set.

$$\phi\left(\mathcal{P}\left(\mathfrak{a}\left(t\right)\right),t\right) = 0\tag{26}$$

Fig. 7. (a) the terrain definition, (b) a squared object moves across the scene. (c) *α*(*u*, *t*) for different values of *t* representation (d)Back tracking from the final point toward the starting

Fluid Dynamics Without Fluids 599

In Fig. 7 there is the representation of level-set evolution, accordingly to moving objects, in an

Terrain traversability can be easily integrated in the presented approach since it is possible to define a map *F* (*x*, *y*) describing the terrain traversability according to continuous values from 1 (optimal traversability) to 0 (impassable) and the optimal path can be obtained just changing

Further velocity dependencies, e.g. from the vehicle weight reduction due to fuel combustion, can be accounted introducing a time dependance in *VN*. More complex conditions can be

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>*<sup>F</sup>* (*x*, *<sup>y</sup>*) *VN***<sup>n</sup>** (29)

*∂*P (*α* (*t*))

The steps for the best path definition can then be summarized as: 1. Evolve the PDE (eq.28) till the destination point is reached.

"egg-box" terrain where a square is moving across the scene.

2. back-track along the optimal path till the origin.

one along the estimated geodesic path

eq. 18 into

A graphical representation of this process can be seen in Fig. 6.

The next step requires the definition of the evolution law for *φ*, that, accordingly to the chain rule can be written as:

$$\nabla \phi \left( \mathcal{P} \left( \mathfrak{a} \left( t \right) \right), t \right) \frac{\partial \mathcal{P} \left( \mathfrak{a} \left( t \right) \right)}{\partial t} + \frac{\partial \phi \left( \mathcal{P} \left( \mathfrak{a} \left( t \right) \right), t \right)}{\partial t} = 0 \tag{27}$$

Furthermore, accordingly to this level-set formulation, the in-plane normal vector can be written as **n** = <sup>∇</sup>*<sup>φ</sup>* �∇*φ*� . This relation together with the Lagrangian evolution equation:

$$\frac{\partial \phi\left(\mathcal{P}\left(a\left(t\right)\right), t\right)}{\partial t} = -V\_N \left\|\nabla \phi\right\| = \sqrt{\frac{\left(1 + q^2\right)\phi\_x^2 + \left(1 + p^2\right)\phi\_y^2 - 2pq\phi\_x\phi\_y}{1 + p^2 + q^2}}\tag{28}$$

This is the Eulerial formulation for curve evolution.

Fig. 6. Implicit representation of the projection of the *α* function through the level-set.

10 Will-be-set-by-IN-TECH

If there are moving obstacles along the path the *MDMp* becomes multivalued and it could be possible that some points change their index along the optimal path, to account for this all

Moving obstacles introduce new boundary constraints and dealing with them could be very complex following the previously describer approach, a simpler way could be using the level-set approach where the projection of *α* on the x-y plane is the zero level set *φ* = 0 and *φ*

The next step requires the definition of the evolution law for *φ*, that, accordingly to the chain

Furthermore, accordingly to this level-set formulation, the in-plane normal vector can be

Fig. 6. Implicit representation of the projection of the *α* function through the level-set.

�∇*φ*� . This relation together with the Lagrangian evolution equation:

(1 + *q*2) *φ*<sup>2</sup>

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> *∂φ* (<sup>P</sup> (*<sup>α</sup>* (*t*)), *<sup>t</sup>*)

*φ* (P (*α* (*t*)), *t*) = 0 (26)

*<sup>x</sup>* + (1 + *p*2) *φ*<sup>2</sup>

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>0</sup> (27)

*<sup>y</sup>* − 2*pqφxφ<sup>y</sup>* <sup>1</sup> <sup>+</sup> *<sup>p</sup>*<sup>2</sup> <sup>+</sup> *<sup>q</sup>*<sup>2</sup> (28)

will be negative in the interior and positive in the exterior of the zero level set.

distance values must be evaluated at a candidate point.

A graphical representation of this process can be seen in Fig. 6.

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>*VN* �∇*φ*� <sup>=</sup>

This is the Eulerial formulation for curve evolution.

<sup>∇</sup>*<sup>φ</sup>* (<sup>P</sup> (*<sup>α</sup>* (*t*)), *<sup>t</sup>*) *<sup>∂</sup>*<sup>P</sup> (*<sup>α</sup>* (*t*))

**4.1 The Level-set approach**

rule can be written as:

written as **n** = <sup>∇</sup>*<sup>φ</sup>*

*∂φ* (P (*α* (*t*)), *t*)

Fig. 7. (a) the terrain definition, (b) a squared object moves across the scene. (c) *α*(*u*, *t*) for different values of *t* representation (d)Back tracking from the final point toward the starting one along the estimated geodesic path

The steps for the best path definition can then be summarized as:


In Fig. 7 there is the representation of level-set evolution, accordingly to moving objects, in an "egg-box" terrain where a square is moving across the scene.

Terrain traversability can be easily integrated in the presented approach since it is possible to define a map *F* (*x*, *y*) describing the terrain traversability according to continuous values from 1 (optimal traversability) to 0 (impassable) and the optimal path can be obtained just changing eq. 18 into

$$\frac{\partial \mathcal{P}\left(\mathfrak{a}\left(t\right)\right)}{\partial t} = -F\left(\mathfrak{x}, \mathcal{y}\right)V\_{\mathcal{N}}\mathbf{n} \tag{29}$$

Further velocity dependencies, e.g. from the vehicle weight reduction due to fuel combustion, can be accounted introducing a time dependance in *VN*. More complex conditions can be

where Γ represents the surface, **x** is the set of 3D points, *d* is the euclidean distance between the point **x** and the closest point on the surface Γ and 1 ≤ *p* ≤ ∞ is the order of the functional distance: for *p* = 2 it gives a root mean square error while for *p* = 1 we just focus on the point with the maximum distance from Γ. looking for the surface that minimizes the functional (30) is similar to enveloping the data set with a membrane with certain elastic properties and have this membrane to evolve in time until it comes to rest. This time evolution paradigm can be

Fluid Dynamics Without Fluids 601

*dp*−<sup>1</sup> (**x**)

∇*d* (**x**) · **n** +

where **n** points accordingly to the evolving point direction, *κ* indicates the mean curvature, ∇*d* (**x**) · **n** represents the surface attraction term and *d* (**x**) *κ* is the surface tension. An important consideration is that function *d* (**x**) will automatically impose that the final surface is more flexible where the 3D points dataset is denser and more rigid in low density regions: this will allow smoother regions where we have less samples and accurate fitting in highly sampled regions. Fluids, with their properties like surface tension and viscosity represent an optimal modelization for the warping surface providing us with a flexible tool for this goal. Also in this case we present a Level-set approach, in particular we will adopt an implicit function *φ* over the spatial 3D domain so that *φ* (**x**, *t*) = 0 represent the warping surface Γ (*t*) at time *t*, for the initialization we use *φ* (**x**, *t*)|*t* = 0 = *D* where *D* is the signed distance of **x** from Γ (*t*) (positive values indicate external points while negative ones represent internal points). One possible approach to make surface Γ to converge towards the cloud of points is, accordingly to Zhao et al. (2001), to use the fluids convection defining a velocity field **v** (**x**)

�

1 *p d* (**x**) *κ* �

∇*d* (**x**) · **n** +

1 *p d* (**x**) *κ* �

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> **<sup>v</sup>** (<sup>Γ</sup> (*t*)) (33)

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>∇</sup>*<sup>d</sup>* (**x**) · ∇*<sup>φ</sup>* (34)

**n** (31)

(32)

described as:

accordingly to:

*∂*Γ *<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup> ⎡ ⎣ �

Γ

*d<sup>p</sup>* (**x**) *ds*

*dp*−<sup>1</sup> (**x**)

Euler-Lagrange equation can the be used to find its minimum:

⎤ ⎦ 1 *<sup>p</sup>* −1

�

*∂*Γ (*t*)

accordingly to eq. 8, the level-set evolution in this case can be defined as:

*∂φ*

In particular we adopted, as a steering potential field, the distance function *d* (**x**) so that the velocity field **v** is oriented towards the minima of the potential field: **v** = −∇*d* (**x**). So,

A major computational drawback comes from the fact that the *φ* function at each time iteration becomes more and more compressed presenting a growing slope across the zero level-set which could introduce numerical imprecisions. To solve this problem some solutions are present in literature, for example a level-set re-initialization in a region around the zero level: this is particularly useful for Fast Marching methods which, as described in section 3 for the Narrow Band approach Sethian (1999). Other approaches, instead, suggest a different implementation for the Hamilton Jacobi PDE like Gomes & Faugeras (1999). In particular in

easily added just acting on *VN*, e.g. if the Robot has a specific shape the velocity can reflect the ability to walk through a narrow path or to turn at a specific point. In some cases the *F* (*x*, *y*) correcting function is called viscosity for it similarity to fluid property that can prevent or slow down motion in particular circumstances.

#### **5. Surface reconstruction from cloud of points**

The temporal evolution of a volumetric implicit function can also be fruitfully applied to reconstruct 3D surfaces form large set of unorganized sample points. The evolving front can be thought as the surface that separates two different fluids obeying specific fluid dynamics laws. One remarkable feature of this approach is its ability to model complex topologies using a set of intuitive tools derived from fluid physics: Global and local surface descriptors are used allowing the parallelization of the algorithm on different object parts working with different resolutions and accuracies. The problem of building surfaces from unorganized sets of 3D points has recently gained a great deal of attention. In fact, in addition to being an interesting problem of topology extraction from geometric information, its applications are becoming more and more numerous. For example, the acquisition of large sets of 3D points is becoming easier and more affordable using, for example, 3D-scanners (*Registration and fusion of intensity and range data for 3d modelling of real world scenes* (2003)). The presented approach is particularly important in applications where objects are better described by their external surface rather than by simple unorganized data (clouds of points, data slices, etc.). For example, in medical applications based on Tomography scans or NMRs it is often necessary to visualize some specific tissues such as the external surface of an organ starting from the acquired 3D points. This can be achieved by selecting the points that belong to a specific class (organ boundary, tissue, etc.) and then generating the surface from their interpolation. In most cases the definition of this surface is an ill-posed problem as there is no unique way to connect points of a dataset into a surface, therefore it is often necessary to introduce constraints for globally or locally controlling the surface behavior. As a matter of fact, the resulting surface often turns out to exhibit a complex topology due to noise in the acquired data or ambiguities in the case of non-convex objects Hoppe (1994). In order to overcome such problems, surface wrapping algorithms need to incorporate specific constraints on the quality of the data fitting (surface closeness to the acquired points), on the maximum surface curvature and roughness, on the number of resulting triangles, etc. The existing surface reconstruction methods can be classified into two broad categories: the former describes the surface as an implicit function while the latter describes it in an explicit form. Explicit (boundary) representations describe the surface in terms of point connections, and traditional approaches are based on Delaunay triangulation, Voronoi diagrams (Amenta et al. (1998)) or NURBS (Piegel & Tiller (1996)).

Here we will show how Level-set technique together with Volume of Fluid end Fluid Dynamics equations can be fruitfully applied to obtain a wrapping surface with intuitive parameters tuning. To evaluate quality of the final reconstruction an important parameter that can also be used as stopping condition for the surface evolution is an energy functional defined as:

$$E\left(\Gamma\right) = \left[\int d^p\left(\mathbf{x}\right)ds\right]^{1/p}, 1 \leqslant p \leqslant \infty \tag{30}$$

12 Will-be-set-by-IN-TECH

easily added just acting on *VN*, e.g. if the Robot has a specific shape the velocity can reflect the ability to walk through a narrow path or to turn at a specific point. In some cases the *F* (*x*, *y*) correcting function is called viscosity for it similarity to fluid property that can prevent or

The temporal evolution of a volumetric implicit function can also be fruitfully applied to reconstruct 3D surfaces form large set of unorganized sample points. The evolving front can be thought as the surface that separates two different fluids obeying specific fluid dynamics laws. One remarkable feature of this approach is its ability to model complex topologies using a set of intuitive tools derived from fluid physics: Global and local surface descriptors are used allowing the parallelization of the algorithm on different object parts working with different resolutions and accuracies. The problem of building surfaces from unorganized sets of 3D points has recently gained a great deal of attention. In fact, in addition to being an interesting problem of topology extraction from geometric information, its applications are becoming more and more numerous. For example, the acquisition of large sets of 3D points is becoming easier and more affordable using, for example, 3D-scanners (*Registration and fusion of intensity and range data for 3d modelling of real world scenes* (2003)). The presented approach is particularly important in applications where objects are better described by their external surface rather than by simple unorganized data (clouds of points, data slices, etc.). For example, in medical applications based on Tomography scans or NMRs it is often necessary to visualize some specific tissues such as the external surface of an organ starting from the acquired 3D points. This can be achieved by selecting the points that belong to a specific class (organ boundary, tissue, etc.) and then generating the surface from their interpolation. In most cases the definition of this surface is an ill-posed problem as there is no unique way to connect points of a dataset into a surface, therefore it is often necessary to introduce constraints for globally or locally controlling the surface behavior. As a matter of fact, the resulting surface often turns out to exhibit a complex topology due to noise in the acquired data or ambiguities in the case of non-convex objects Hoppe (1994). In order to overcome such problems, surface wrapping algorithms need to incorporate specific constraints on the quality of the data fitting (surface closeness to the acquired points), on the maximum surface curvature and roughness, on the number of resulting triangles, etc. The existing surface reconstruction methods can be classified into two broad categories: the former describes the surface as an implicit function while the latter describes it in an explicit form. Explicit (boundary) representations describe the surface in terms of point connections, and traditional approaches are based on Delaunay triangulation, Voronoi

slow down motion in particular circumstances.

**5. Surface reconstruction from cloud of points**

diagrams (Amenta et al. (1998)) or NURBS (Piegel & Tiller (1996)).

*E* (Γ) =

⎡ ⎣ �

Γ

defined as:

Here we will show how Level-set technique together with Volume of Fluid end Fluid Dynamics equations can be fruitfully applied to obtain a wrapping surface with intuitive parameters tuning. To evaluate quality of the final reconstruction an important parameter that can also be used as stopping condition for the surface evolution is an energy functional

> ⎤ ⎦

1/*p*

, 1 *p* ∞ (30)

*d<sup>p</sup>* (**x**) *ds*

where Γ represents the surface, **x** is the set of 3D points, *d* is the euclidean distance between the point **x** and the closest point on the surface Γ and 1 ≤ *p* ≤ ∞ is the order of the functional distance: for *p* = 2 it gives a root mean square error while for *p* = 1 we just focus on the point with the maximum distance from Γ. looking for the surface that minimizes the functional (30) is similar to enveloping the data set with a membrane with certain elastic properties and have this membrane to evolve in time until it comes to rest. This time evolution paradigm can be described as:

$$\frac{\partial \Gamma}{\partial t} = -\left[ \int\_{\Gamma} d^p \left( \mathbf{x} \right) ds \right]^{\frac{1}{p} - 1} d^{p - 1} \left( \mathbf{x} \right) \left[ \nabla d \left( \mathbf{x} \right) \cdot \mathbf{n} + \frac{1}{p} d \left( \mathbf{x} \right) \kappa \right] \mathbf{n} \tag{31}$$

Euler-Lagrange equation can the be used to find its minimum:

$$d^{p-1}\left(\mathbf{x}\right)\left[\nabla d\left(\mathbf{x}\right)\cdot\mathbf{n} + \frac{1}{p}d\left(\mathbf{x}\right)\kappa\right] \tag{32}$$

where **n** points accordingly to the evolving point direction, *κ* indicates the mean curvature, ∇*d* (**x**) · **n** represents the surface attraction term and *d* (**x**) *κ* is the surface tension. An important consideration is that function *d* (**x**) will automatically impose that the final surface is more flexible where the 3D points dataset is denser and more rigid in low density regions: this will allow smoother regions where we have less samples and accurate fitting in highly sampled regions. Fluids, with their properties like surface tension and viscosity represent an optimal modelization for the warping surface providing us with a flexible tool for this goal. Also in this case we present a Level-set approach, in particular we will adopt an implicit function *φ* over the spatial 3D domain so that *φ* (**x**, *t*) = 0 represent the warping surface Γ (*t*) at time *t*, for the initialization we use *φ* (**x**, *t*)|*t* = 0 = *D* where *D* is the signed distance of **x** from Γ (*t*) (positive values indicate external points while negative ones represent internal points). One possible approach to make surface Γ to converge towards the cloud of points is, accordingly to Zhao et al. (2001), to use the fluids convection defining a velocity field **v** (**x**) accordingly to:

$$\frac{\partial \Gamma \left( t \right)}{\partial t} = \mathbf{v} \left( \Gamma \left( t \right) \right) \tag{33}$$

In particular we adopted, as a steering potential field, the distance function *d* (**x**) so that the velocity field **v** is oriented towards the minima of the potential field: **v** = −∇*d* (**x**). So, accordingly to eq. 8, the level-set evolution in this case can be defined as:

$$\frac{\partial \phi}{\partial t} = \nabla d\left(\mathbf{x}\right) \cdot \nabla \phi \tag{34}$$

A major computational drawback comes from the fact that the *φ* function at each time iteration becomes more and more compressed presenting a growing slope across the zero level-set which could introduce numerical imprecisions. To solve this problem some solutions are present in literature, for example a level-set re-initialization in a region around the zero level: this is particularly useful for Fast Marching methods which, as described in section 3 for the Narrow Band approach Sethian (1999). Other approaches, instead, suggest a different implementation for the Hamilton Jacobi PDE like Gomes & Faugeras (1999). In particular in

Fig. 8. a 2D non-convex example, arrows indicate the steering velocity field toward the cloud of points and different level-sets are also represented. Zones A and B are zoomed in Fig.9

Fluid Dynamics Without Fluids 603

Fig. 9. zoom of zone A and B of Fig. 8 in particular different behavior of the steering field prevents external fluid from entering the convex part in B while favor its penetration in A

Further aspects from fluid dynamics can be fruitfully applied, e.g. the diffusive parameter *ξ* in eq. 35 is a delicate parameter as it is strictly connected to fluid viscosity. High values of *ξ* result in a smooth but often inaccurate surface while low values result in a surface that strictly honors the data but is sensitive to acquisition noise. As a matter of fact for laser scanner acquisitions the portions of the object with high samples density ( usually due to multiple acquisitions from different viewpoints) requires high accuracy and/or present complex topologies, while smoother or less critical areas usually exhibit a lower sample density. In order to account for that, we can establish an inverse dependency between the

> *<sup>ξ</sup>* <sup>=</sup> *<sup>α</sup>* 1 + *ρ<sup>p</sup>*

(36)

(nonconvex region)

diffusive constant *ξ* and the local point density:

Marcon et al. (2008) a quite different interpretation for the level-set function is presented that aims to conciliate Level-set upgrade, narrow band approach, discretization issues and Volume of Fluid (VoF) method: in particular, the VoF method (see Noh & Woodward (1976)) has been formulated in a variety of forms and re-introduced under different names such as the "cell method" and the "partial fractions method". The key idea behind this technique is to define a fixed computational grid and assign to each grid cell a value that describes the relative proportions of two materials contained in that cell. Our particular modeling metaphor is based on two immiscible fluids of opposite densities *ρ* = 1 and *ρ* = −1, which identify the surfaceŠs outside and inside, respectively. More specifically, cells completely filled with outer fluid take on the value +1 while cells that are only filled with inner fluid will take on the value -1. When filled with a mix of both fluids, the cell is assigned an intermediate value in the ] − 1, 1[ range, which measures the relative density of the mix (i.e. a value of zero indicates the same amount, half and half of the two fluids inside the cell). As a first step we define a bounding box that completely encloses the point cloud. We assume that this box is completely filled with inner fluid and that the whole space outside of the box is filled with outer fluid. We define the rules of evolution so that both fluids are set in motion by the presence of data points, which act as attractors. This fluid migration takes place in compliance with laws of conservation. According to such rules, the front starts propagating inward (towards the cloud of points) and stops only when the inner fluid is confined "inside" the cloud of points and the interface between the two fluids wraps the point-cloud. At each iteration step we will replace the value of points on the domain boundary B with 1 (this is equivalent to place outer fluid sources on the domain boundary), the cells with values greater than 1 will be clipped with value '1' while cells with values lower that -1 will be clipped to -1: this will be equivalent to remove excess fluids from cells without any need to introduce further fluid sources or drains. Two main contribution to fluid motion are considered: convective and diffusive: the first one is important for steering the zero level-set toward the cloud of points while the second has an important role in smoothing evolving front. Accordingly to the Gauss theorem when no sources or drains are present the fluid density will follow eq. (35):

$$\frac{d\rho}{dt} = \nabla \cdot (\mathbf{v}\rho) + \nabla \cdot (\xi \nabla \rho) \tag{35}$$

where *ρ* indicates the fluid density, **v***ρ* is the convective term indicating the amount of fluid per volume unit that is transported by the velocity field while *ξ*∇*ρ* represents the diffusive part where *ξ* is the diffusive term (further details in (Gueyffier et al. (1999))). In this approach a point will be assigned a velocity that is proportional to the vector **v** that joins that point to the closest data point. This way data points act as attractors for both internal and external fluids and, therefore, for the evolving front: each sample point represents a central vector field. Even if this attraction field has no physical equivalent, it is physically consistent thanks to its irrotational characteristic (see Marcon et al. (2008) for further details). This fact guarantees algorithmic convergence. In Figures 8 and 9 it is possible to see the direction and intensity of the steering field toward the cloud of points and different level-set curves; in particular in the magnification of zones A and B is possible to see how the steering field automatically stops external liquid from entering the object (all the internal arrows point outside) while, in the convex part, the steering field force the outer fluid to enter the convex region.

In Fig. 10 we represent the final state for the wrapping of a two dimensional cloud of points.

14 Will-be-set-by-IN-TECH

Marcon et al. (2008) a quite different interpretation for the level-set function is presented that aims to conciliate Level-set upgrade, narrow band approach, discretization issues and Volume of Fluid (VoF) method: in particular, the VoF method (see Noh & Woodward (1976)) has been formulated in a variety of forms and re-introduced under different names such as the "cell method" and the "partial fractions method". The key idea behind this technique is to define a fixed computational grid and assign to each grid cell a value that describes the relative proportions of two materials contained in that cell. Our particular modeling metaphor is based on two immiscible fluids of opposite densities *ρ* = 1 and *ρ* = −1, which identify the surfaceŠs outside and inside, respectively. More specifically, cells completely filled with outer fluid take on the value +1 while cells that are only filled with inner fluid will take on the value -1. When filled with a mix of both fluids, the cell is assigned an intermediate value in the ] − 1, 1[ range, which measures the relative density of the mix (i.e. a value of zero indicates the same amount, half and half of the two fluids inside the cell). As a first step we define a bounding box that completely encloses the point cloud. We assume that this box is completely filled with inner fluid and that the whole space outside of the box is filled with outer fluid. We define the rules of evolution so that both fluids are set in motion by the presence of data points, which act as attractors. This fluid migration takes place in compliance with laws of conservation. According to such rules, the front starts propagating inward (towards the cloud of points) and stops only when the inner fluid is confined "inside" the cloud of points and the interface between the two fluids wraps the point-cloud. At each iteration step we will replace the value of points on the domain boundary B with 1 (this is equivalent to place outer fluid sources on the domain boundary), the cells with values greater than 1 will be clipped with value '1' while cells with values lower that -1 will be clipped to -1: this will be equivalent to remove excess fluids from cells without any need to introduce further fluid sources or drains. Two main contribution to fluid motion are considered: convective and diffusive: the first one is important for steering the zero level-set toward the cloud of points while the second has an important role in smoothing evolving front. Accordingly to the Gauss theorem when no

sources or drains are present the fluid density will follow eq. (35):

*∂ρ*

convex part, the steering field force the outer fluid to enter the convex region.

In Fig. 10 we represent the final state for the wrapping of a two dimensional cloud of points.

where *ρ* indicates the fluid density, **v***ρ* is the convective term indicating the amount of fluid per volume unit that is transported by the velocity field while *ξ*∇*ρ* represents the diffusive part where *ξ* is the diffusive term (further details in (Gueyffier et al. (1999))). In this approach a point will be assigned a velocity that is proportional to the vector **v** that joins that point to the closest data point. This way data points act as attractors for both internal and external fluids and, therefore, for the evolving front: each sample point represents a central vector field. Even if this attraction field has no physical equivalent, it is physically consistent thanks to its irrotational characteristic (see Marcon et al. (2008) for further details). This fact guarantees algorithmic convergence. In Figures 8 and 9 it is possible to see the direction and intensity of the steering field toward the cloud of points and different level-set curves; in particular in the magnification of zones A and B is possible to see how the steering field automatically stops external liquid from entering the object (all the internal arrows point outside) while, in the

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> ∇ · (**v***ρ*) <sup>+</sup> ∇ · (*ξ*∇*ρ*) (35)

Fig. 8. a 2D non-convex example, arrows indicate the steering velocity field toward the cloud of points and different level-sets are also represented. Zones A and B are zoomed in Fig.9

Fig. 9. zoom of zone A and B of Fig. 8 in particular different behavior of the steering field prevents external fluid from entering the convex part in B while favor its penetration in A (nonconvex region)

Further aspects from fluid dynamics can be fruitfully applied, e.g. the diffusive parameter *ξ* in eq. 35 is a delicate parameter as it is strictly connected to fluid viscosity. High values of *ξ* result in a smooth but often inaccurate surface while low values result in a surface that strictly honors the data but is sensitive to acquisition noise. As a matter of fact for laser scanner acquisitions the portions of the object with high samples density ( usually due to multiple acquisitions from different viewpoints) requires high accuracy and/or present complex topologies, while smoother or less critical areas usually exhibit a lower sample density. In order to account for that, we can establish an inverse dependency between the diffusive constant *ξ* and the local point density:

$$\tilde{\xi} = \frac{\mathfrak{a}}{1 + \rho\_p} \tag{36}$$

Fig. 11. Evolution of the zero level-set toward a 3D cloud of points acquired by a laser scanner. The level-set starting from the boundary of the considered volume (top left) evolves

Fig. 12. a set of points presenting a high curvature region within a concavity

In some particular cases, where high curvature is present in nonconvex regions with a small entrance, external fluid could be inhibited from entering (e.g. consider the cloud of points of

Fluid Dynamics Without Fluids 605

For cases like these a further help could come from vorticity. In particular, since our steering field is based on the distance from the closest point it will be a conservative irrotational vector field and when the outer fluid will approach the concavity mouth it will be pushed to its borders without the opportunity to get inside (we underline that this condition is different from the one in Fig. 8 where the concavity of the region A does not present high curvature and external fluid is correctly steered inside the region). A possible solution consists in introducing a rotational component in the steering field, this could be obtained adding periodically a small random displacement for each point around its original position: it will result in a vorticity

towards the final cloud of points.

Fig. 12).

**5.1 Turbolence and convex regions**

Fig. 10. the wrapping of a set of 2D points (top right), the zero level-set at the end of the evolution represent the joining line, inner points have positive *ψ* while external ones have negative values

where *α* is a parameter that indicates the estimated noise in the point location (the higher the value the smoother the surface). The density *ρp* at a given point is defined as the number of samples that fall into a cube of volume *V*/*N* around that point, where *V* is the volume of the point cloud (approximated by the volume of the smallest parallelepiped that contains the point cloud) and *N* is the total number of samples. This means that if all the samples would be uniformly distributed in the considered volume (like in a cubic crystal) in each cubic cell of volume *V*/*N* there would be just one sample. In our case, since we have to define the *ρp* variable for each voxel, we define it as the sum of all the points falling inside a cube of volume *V*/*N* centered at the voxel center. In surface regions with a low point density the diffusive constant *ξ* corresponds to a smooth surface, while in high-density zones *ξ* becomes very low, to guarantee that the surface will strictly honor the sample points. The level-set evolution equation is then the following:

$$\begin{cases} \frac{\partial \phi(\mathbf{x},t)}{\partial t} = \nabla \cdot \left( \mathbf{v} \phi \left( \mathbf{x},t \right) \right) + \xi \nabla^2 \phi \left( \mathbf{x},t \right), \ t \ge 0\\ \phi \left( \mathbf{x},t \right) = +1, & t \ge 0, \mathbf{x} \in \mathcal{B} \\ \phi \left( \mathbf{x},0 \right) = -1 & \mathbf{x} \notin \mathcal{B} \end{cases} \tag{37}$$

where the second raw indicates that the boundary B of the considered domain is always kept equal to +1 assuming that those points are sources of external fluid. The third raw indicates that at the starting point all the points inside the domain are filled with the internal fluid. In Fig. 11 the level-set evolves toward the cloud of points which is a laser scanner acquisition of a plaster wolf.

16 Will-be-set-by-IN-TECH

Fig. 10. the wrapping of a set of 2D points (top right), the zero level-set at the end of the evolution represent the joining line, inner points have positive *ψ* while external ones have

to guarantee that the surface will strictly honor the sample points.

The level-set evolution equation is then the following:

*∂φ*(**x**,*t*)

⎧ ⎪⎨

⎪⎩

where *α* is a parameter that indicates the estimated noise in the point location (the higher the value the smoother the surface). The density *ρp* at a given point is defined as the number of samples that fall into a cube of volume *V*/*N* around that point, where *V* is the volume of the point cloud (approximated by the volume of the smallest parallelepiped that contains the point cloud) and *N* is the total number of samples. This means that if all the samples would be uniformly distributed in the considered volume (like in a cubic crystal) in each cubic cell of volume *V*/*N* there would be just one sample. In our case, since we have to define the *ρp* variable for each voxel, we define it as the sum of all the points falling inside a cube of volume *V*/*N* centered at the voxel center. In surface regions with a low point density the diffusive constant *ξ* corresponds to a smooth surface, while in high-density zones *ξ* becomes very low,

> *<sup>∂</sup><sup>t</sup>* <sup>=</sup> ∇ · (**v***<sup>φ</sup>* (**x**, *<sup>t</sup>*)) <sup>+</sup> *<sup>ξ</sup>*∇2*<sup>φ</sup>* (**x**, *<sup>t</sup>*), *<sup>t</sup>* <sup>0</sup> *φ* (**x**, *t*) = +1, *t* 0, **x** ∈ B *φ* (**x**, 0) = −1 **x** ∈ B /

where the second raw indicates that the boundary B of the considered domain is always kept equal to +1 assuming that those points are sources of external fluid. The third raw indicates that at the starting point all the points inside the domain are filled with the internal fluid. In Fig. 11 the level-set evolves toward the cloud of points which is a laser scanner acquisition of

(37)

negative values

a plaster wolf.

Fig. 11. Evolution of the zero level-set toward a 3D cloud of points acquired by a laser scanner. The level-set starting from the boundary of the considered volume (top left) evolves towards the final cloud of points.

#### **5.1 Turbolence and convex regions**

In some particular cases, where high curvature is present in nonconvex regions with a small entrance, external fluid could be inhibited from entering (e.g. consider the cloud of points of Fig. 12).

Fig. 12. a set of points presenting a high curvature region within a concavity

For cases like these a further help could come from vorticity. In particular, since our steering field is based on the distance from the closest point it will be a conservative irrotational vector field and when the outer fluid will approach the concavity mouth it will be pushed to its borders without the opportunity to get inside (we underline that this condition is different from the one in Fig. 8 where the concavity of the region A does not present high curvature and external fluid is correctly steered inside the region). A possible solution consists in introducing a rotational component in the steering field, this could be obtained adding periodically a small random displacement for each point around its original position: it will result in a vorticity

and

respect to *φ* will give:

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> *δ�* (*φ*)

 −

*<sup>u</sup>*<sup>+</sup> <sup>−</sup> *<sup>u</sup>*<sup>0</sup>

Fig. 13. update procedure in the narrow band points

 <sup>2</sup> <sup>−</sup> *<sup>μ</sup>* ∇*u*<sup>+</sup> <sup>2</sup> +

*∂φ*

*<sup>u</sup>*<sup>−</sup> <sup>−</sup> *<sup>u</sup>*<sup>0</sup> <sup>=</sup> *<sup>μ</sup>*∇2*u*<sup>−</sup> (42)

*<sup>∂</sup>***<sup>n</sup>** = 0 (i.e.

(43)

 <sup>∇</sup>*<sup>φ</sup>* |∇*φ*|

which correspond to the diffusive part inside the borders defined by Γ, since *<sup>∂</sup>u*+,<sup>−</sup>

no diffusion will take place across the boundary keeping edges sharp). While, deriving with

Fluid Dynamics Without Fluids 607

Accordingly to a steepest descendent approach. The *δ�* indicates that the delta of Dirac, in a discrete environment, must be calculated in a narrow band and, to prevent discontinuities *u*<sup>+</sup> has to be extended, in the narrow band, for *φ* < 0 values and also *u*<sup>−</sup> has to be extended in the narrow band for *φ* > 0. A possible approach to these extension could be the "Ghost fluid Method" (Fedkiw et al. (1999)) where the interface is treated as a moving boundary. Solving the governing equations (Terashima & Tryggvason (2009)) with the extension of the discontinuous variables across the fluids interface, (for real CFD applications this is usually entropy), will reduce smearing of discontinuities providing smooth convergence to the final image without oscillations. In Fig. 13 is represented how the contour between the two image parts is steered accordingly to the Ghost Fluid Method in the narrow band while in Fig. 14 it is possible to see the evolution of the contour across the noisy image: diffusive terms smooth uniform regions reducing noise while the contour evolution smoothly warp the objects.

*u*<sup>−</sup> − *u*<sup>0</sup>

 <sup>2</sup> + *<sup>μ</sup>* ∇*u*<sup>−</sup> <sup>2</sup> <sup>+</sup> *<sup>ν</sup>*<sup>∇</sup>

contribution: *ζ* = ∇ × **v** (**x**) the resulting effect is that inner fluid in convex regions is pushed inside and outside with the same probability while in nonconvex regions inner fluid that went inside stays there and progressively fills the region. Particular care must be placed in introducing the amount of random points displacement since convergence of the algorithm is no more guaranteed by the conservative vector field.

#### **6. Image processing and viscous noise removal**

Fluid dynamics applied together with level-set method can also provide a powerful tool for noise removal in complex images (i.e. when noise is not simply removable using intensity thresholds of histogram analysis), a Total Variation denoising approach is presented in (Osher & Fedkiw (2002)) where the regularization algorithm is based on this PDE:

$$L\mu = \lambda R\mu,\tag{38}$$

where *u* is the image represented as a 2D function (*u* (*x*, *y*)) that indicates intensity value for each point; *R* is the regularization operator and *L* is the time-space operator. A very common energy functional *E*(*u*) whose minimization provide the denoised image is the Mumford-Shah multiscale segmentation model:

$$\min\_{\boldsymbol{u}} \mathcal{E}\left(\boldsymbol{u}\right) = \frac{1}{2} \int\_{\Omega} (\boldsymbol{u} - \boldsymbol{u}\_{0})^{2} d\boldsymbol{x} + \mu \int\_{\Omega \backslash \Gamma} |\nabla \boldsymbol{u}|^{2} d\boldsymbol{x} + \nu \mathcal{H}^{2}(\Gamma) \tag{39}$$

where *u*<sup>0</sup> is the given noisy image, Ω is the image function domain and *μ* is a weighting coefficient for the average smoothness of regions divided by the contour Γ. *ν* is the weight coefficient for the total contour length expressed by the Haussdorf measure H.

The minimization of *E*(*u*) will provide an piecewise-smooth approximation of the initial *u*<sup>0</sup> image, Γ has the role of approximating the edges in the image *u*<sup>0</sup> and *u* will be smooth only in Ω\Γ; Accordingly to the level-set formulation the implicit function will be *φ* and in particular *φ*(**x**) = 0 represents the contour Γ separating region of positive values of *φ* from the region of negative values. In particular *u*+(**x**) will be the intensity function representing values of the region where *φ*(**x**) > 0 while *u*−(**x**) will represent the intensity of points in the *φ*(**x**) < 0 region. Eq. 39 then become:

$$\begin{split} \min\_{\boldsymbol{\mu}} & E\left(\boldsymbol{u}^{+}, \boldsymbol{u}^{-}, \boldsymbol{\phi}\right) = \int\_{\Omega} \left(\boldsymbol{u}^{+} - \boldsymbol{u}\_{0}\right)^{2} H\left(\boldsymbol{\phi}\right) d\mathbf{x} + \int\_{\Omega} \left(\boldsymbol{u}^{-} - \boldsymbol{u}\_{0}\right)^{2} \left(1 - H\left(\boldsymbol{\phi}\right)\right) d\mathbf{x} + \\ & + \mu \int\_{\Omega} \left|\nabla \boldsymbol{u}^{+}\right|^{2} H\left(\boldsymbol{\phi}\right) d\mathbf{x} + \mu \int\_{\Omega} \left|\nabla \boldsymbol{u}^{-}\right|^{2} \left(1 - H\left(\boldsymbol{\phi}\right)\right) d\mathbf{x} + \nu \int\_{\Omega} \left|\nabla H\left(\boldsymbol{\phi}\right)\right| \end{split} \tag{40}$$

Where *H*() represents the Heaviside step function and the last term, Ω |∇*H* (*φ*)| is the integral in the sense of geometric measure, i.e. for a 2D case, is the length of contour Γ.

The minimization can be done iteratively on *u*<sup>+</sup> , *u*<sup>−</sup> and on *φ*, in particular, minimizing with respect to *u*<sup>+</sup> and *u*−, i.e. deriving eq. 40, will give:

$$
\mu^{+} - \mu\_{0} = \mu \nabla^{2} \mu^{+} \tag{41}
$$

and

18 Will-be-set-by-IN-TECH

contribution: *ζ* = ∇ × **v** (**x**) the resulting effect is that inner fluid in convex regions is pushed inside and outside with the same probability while in nonconvex regions inner fluid that went inside stays there and progressively fills the region. Particular care must be placed in introducing the amount of random points displacement since convergence of the algorithm is

Fluid dynamics applied together with level-set method can also provide a powerful tool for noise removal in complex images (i.e. when noise is not simply removable using intensity thresholds of histogram analysis), a Total Variation denoising approach is presented in

where *u* is the image represented as a 2D function (*u* (*x*, *y*)) that indicates intensity value for each point; *R* is the regularization operator and *L* is the time-space operator. A very common energy functional *E*(*u*) whose minimization provide the denoised image is the Mumford-Shah

> 2 *dx* + *μ*

where *u*<sup>0</sup> is the given noisy image, Ω is the image function domain and *μ* is a weighting coefficient for the average smoothness of regions divided by the contour Γ. *ν* is the weight

The minimization of *E*(*u*) will provide an piecewise-smooth approximation of the initial *u*<sup>0</sup> image, Γ has the role of approximating the edges in the image *u*<sup>0</sup> and *u* will be smooth only in Ω\Γ; Accordingly to the level-set formulation the implicit function will be *φ* and in particular *φ*(**x**) = 0 represents the contour Γ separating region of positive values of *φ* from the region of negative values. In particular *u*+(**x**) will be the intensity function representing values of the region where *φ*(**x**) > 0 while *u*−(**x**) will represent the intensity of points in the *φ*(**x**) < 0


Ω\Γ

*Lu* = *λRu*, (38)

<sup>2</sup> *dx* <sup>+</sup> *<sup>ν</sup>*H2(Γ) (39)

<sup>2</sup> (<sup>1</sup> <sup>−</sup> *<sup>H</sup>* (*φ*))*d***x**<sup>+</sup>


(40)


Ω

(Osher & Fedkiw (2002)) where the regularization algorithm is based on this PDE:

(*u* − *u*0)

coefficient for the total contour length expressed by the Haussdorf measure H.

2

*H* (*φ*) *d***x** +

The minimization can be done iteratively on *u*<sup>+</sup> , *u*<sup>−</sup> and on *φ*, in particular, minimizing with

 *u*<sup>−</sup> − *u*<sup>0</sup>

Ω

*<sup>u</sup>*<sup>+</sup> <sup>−</sup> *<sup>u</sup>*<sup>0</sup> <sup>=</sup> *<sup>μ</sup>*∇2*u*<sup>+</sup> (41)

Ω

<sup>2</sup> (<sup>1</sup> <sup>−</sup> *<sup>H</sup>* (*φ*)) *<sup>d</sup>***x**+*<sup>ν</sup>*

no more guaranteed by the conservative vector field.

**6. Image processing and viscous noise removal**

min*<sup>u</sup> <sup>E</sup>* (*u*) <sup>=</sup> <sup>1</sup>

2 

Ω

multiscale segmentation model:

region. Eq. 39 then become:

*u*+, *u*−, *φ*

 = 

respect to *u*<sup>+</sup> and *u*−, i.e. deriving eq. 40, will give:

Ω

<sup>2</sup> *<sup>H</sup>* (*φ*) *<sup>d</sup>***x**+*<sup>μ</sup>*

 *<sup>u</sup>*<sup>+</sup> <sup>−</sup> *<sup>u</sup>*<sup>0</sup>

 ∇*u*<sup>−</sup> 

in the sense of geometric measure, i.e. for a 2D case, is the length of contour Γ.

Ω

Where *H*() represents the Heaviside step function and the last term,

min*<sup>u</sup> <sup>E</sup>*

Ω

 ∇*u*<sup>+</sup> 

+*μ* 

$$
\mu^- - \mu\_0 = \mu \nabla^2 \mu^- \tag{42}
$$

which correspond to the diffusive part inside the borders defined by Γ, since *<sup>∂</sup>u*+,<sup>−</sup> *<sup>∂</sup>***<sup>n</sup>** = 0 (i.e. no diffusion will take place across the boundary keeping edges sharp). While, deriving with respect to *φ* will give:

$$\frac{\partial \phi}{\partial t} = \delta\_{\text{ef}}(\phi) \left[ - \left| u^+ - u\_0 \right|^2 - \mu \left| \nabla u^+ \right|^2 + \left| u^- - u\_0 \right|^2 + \mu \left| \nabla u^- \right|^2 + \nu \nabla \left( \frac{\nabla \phi}{|\nabla \phi|} \right) \right] \tag{43}$$

Accordingly to a steepest descendent approach. The *δ�* indicates that the delta of Dirac, in a discrete environment, must be calculated in a narrow band and, to prevent discontinuities *u*<sup>+</sup> has to be extended, in the narrow band, for *φ* < 0 values and also *u*<sup>−</sup> has to be extended in the narrow band for *φ* > 0. A possible approach to these extension could be the "Ghost fluid Method" (Fedkiw et al. (1999)) where the interface is treated as a moving boundary. Solving the governing equations (Terashima & Tryggvason (2009)) with the extension of the discontinuous variables across the fluids interface, (for real CFD applications this is usually entropy), will reduce smearing of discontinuities providing smooth convergence to the final image without oscillations. In Fig. 13 is represented how the contour between the two image parts is steered accordingly to the Ghost Fluid Method in the narrow band while in Fig. 14 it is possible to see the evolution of the contour across the noisy image: diffusive terms smooth uniform regions reducing noise while the contour evolution smoothly warp the objects.

Fig. 13. update procedure in the narrow band points

Bridson, R. (2003). *Computational aspects of dynamic surfaces*, PhD thesis, University of British

Fluid Dynamics Without Fluids 609

Dobrushin, R. L., Kotecky, R. & Shlosman, S. (1992). ` *Wulff Construction: A Global Shape from*

Fedkiw, R. P., Aslam, T., Merriman, B. & Osher., S. (1999). A nonoscillatory eulerian approach

Gomes, J. & Faugeras, O. (1999). Reconciling distance functions and level sets, tech. rep. 3666,

Gueyffier, D., Li, J., Nadim, A., Scardovelli, R. & Zaleski, S. (1999). Volume-of-fluid

Hennessy, M. G. (2010). *Liquid Snowflake Formation in Superheated Ice*, PhD thesis, University of

Hills, R. N. & Roberts, P. H. (1993). A note on the kinetic conditions at a supercooled interface,

Hoppe, H. (1994). *Surface reconstruction from unorganized points*, PhD thesis, University of

Kimmel, R., Amir, A. & Bruckstein, A. M. (1995). Finding shortest paths on surfaces using

Kimmel, R., Kyriati, N. & Bruckstein, A. M. (1998). Multivauled distance maps fro motion

Libbrecht, K. G. (2005). The physics of snow crystals, *Rep. Prog. Phys. 68*, Norman Bridge

Losasso, F., F. Gibou, F. & Fedkiw, R. (2004). Simulating water and smoke with an octree data

Marcon, M., Piccarreta, L., Sarti, A. & Tubaro, S. (2008). Fast pde approach to surface

Markstein, G. H. (1951). Experimental and theoretical studies of flame front stability, *J. Aero.*

Noh, W. & Woodward, P. (1976). A simple line interface calculation, *in* Springer-Verlag (ed.),

Osher, S. & Fedkiw, R. (2002). *Level Set Methods and Dynamic Implicit Surfaces*, Springer-Verlag. Piegel, L. & Tiller, W. (1996). *The NURBS Book*, second ed. edn, Springer-Verlag, Berlin. Pimpinelli, A. & Villain, J. (1999). *Physics of Crystal Growth*, Cambridge University Press. *Registration and fusion of intensity and range data for 3d modelling of real world scenes* (2003).

Laboratory of Physics, California Institute of Technology.

*Proceedings of 5th International Conference on Fluid Dynamics*.

structure, *ACM Transactions on Graphics* 23(3): 457–462.

level sets propagation, *IEEE Transactions on Pattern Analysis and Machine Intelligence*

planning on surfaces with moving obstacles, *IEEE Trans. on Robotics and automation*

reconstruction from large cloud of points, *Computer Vision and Image Understanding*

Vol. Fourth International Conference on 3-D Digital Imaging and Modeling, 3DIM

to interfaces in multimaterial flows (the ghost fluid method), *Journal of Computational*

interface tracking with smoothed surface stress methods for three-dimensional flows,

Davis, S. H. (2001). *Theory of solidification*, Cambridge University Press.

*Local Interaction*, American Mathematical Society.

*Technical report*, INRIA Sophia Antipolis.

*Int. Comm. Heat Mass Transfer* 20: 407–416. Hobbs, P. V. (2010). *Ice Physics*, Oxford University Press.

Latombe, J. C. (1991). *Robot Motion Planning*, Kluwer.

Columbia.

Oxford.

Washington.

17(6): 635–640.

14(3): 427–436.

112: 274–285.

*Sci.* 18: 199.

Proceedings.

*Physics* 152(2): 457–492.

*J. Comput. Phys.* 152: 423–456.

Fig. 14. zero level set evolution accordingly to the Mumford-Shah energy minimization for object segmentation and image denoising

#### **7. Conclusion**

In this chapter we have shown some significant implementations of Computational Fluid Dynamics equations in completely different contexts ranging from image denoising to optimal path research for moving robots. In recent literature it is also possible to find different applications to many other research fields: The main reason for this large diffusion and reuse of CFD equations can be mainly addressed to:


All these aspects concurred in making CFD a science that goes well beyond fluids analysis providing a self consistent, well known, spread of tools to find fruitful and unlimited applications in engineering, computer science and other scientific research fields.

#### **8. References**


20 Will-be-set-by-IN-TECH

Fig. 14. zero level set evolution accordingly to the Mumford-Shah energy minimization for

In this chapter we have shown some significant implementations of Computational Fluid Dynamics equations in completely different contexts ranging from image denoising to optimal path research for moving robots. In recent literature it is also possible to find different applications to many other research fields: The main reason for this large diffusion and reuse

• the actual deep insight and knowledge both into the theoretical and computational aspect

• The accuracy and stability of modelling equations that can grant convergence in a wide

• Versatility and adaptability of different models from multiphase, to viscous or compressible fluids that provide the user with many intuitive and easy-to-tune parameters

• The close relation to physical quantities that in many cases can automatically provide the

All these aspects concurred in making CFD a science that goes well beyond fluids analysis providing a self consistent, well known, spread of tools to find fruitful and unlimited

Adalsteinsson, D. & Sethain, J. A. (1995). A fast level set method for propagating interfaces,

Amenta, N., Bern, M. & Eppstein, D. (1998). The crust and the *β*-skeleton: combinatorial curve

Biben, T., Klaus, K. & Chaouqi, M. (2005). Phase-field approach to three-dimensional vesicle

applications in engineering, computer science and other scientific research fields.

object segmentation and image denoising

of CFD equations can be mainly addressed to:

allowing an easy adaptation to different contexts.

*Journal of Computational Physics* 118(2): 269–277.

dynamics, *Phys. Rev. E* 72(4).

reconstruction, *Graph. Models Image Process* 60(2): 125–135.

**7. Conclusion**

of fluids.

**8. References**

variety of applications.

minimum energy solution.


**28** 

*Mexico* 

H. Pérez-de-Tejada

*Institute of Geophysics, UNAM* 

**Fluid Dynamics in Space Sciences** 

Much of what it is understood in nature and that it is also inherent to our common activities can be appropriately interpreted as representing evidence of a collective response that substantiates the basis of fluid dynamics. Such is the case for the behavior for a group of particles or individuals that interact with each other as they follow common trajectories. Under most circumstances their interaction takes place across distances (mean free path) that are much smaller than the size of the region where they move and, as a whole, they exhibit properties (density, speed, temperature) that represent average local values of the conglomerate. Besides the application of this view to standard problems in physics and engineering it has been intuitively suggested that it could also be relevant to describe the motion of the solar wind which at large distances from the sun travels with its particles barely executing any collisions among them. Even though the solar wind is, in fact, a (collisionless) ionized gas it maintains a collective response as it expands through interplanetary space and interacts with the planets of the solar system. The text below describes the manner in which the fluid dynamic interpretation of the solar wind was initiated and how it has expanded to

The over one million degree temperature of the solar corona is significantly larger than the nearly six thousand degrees of the solar surface and thus provides the energy source that ultimately leads to its strong outward expansion. In fact, the large amounts of energy that are delivered to the solar corona from the solar interior can be mechanically used for producing the solar wind. The prediction and theoretical description of this phenomenon was provided by E. N. Parker 1 at the University of Chicago who over 50 years ago devised a fluid dynamic interpretation of the manner in which the ionized coronal gas (plasma) is forced to expand outward reaching speeds of the order of 300-400 km/s. A remarkable aspect of this concept is that the solar wind is predicted to be an outflowing continuum gas that expands away from the sun as its density decreases to values that by the earth´s orbit are very small (~ 10 cm-3). As a result, collisions among its particles (mostly protons and electrons) become extremely rare and the solar wind rapidly becomes a collisionless plasma with an effective mean free path for the collisions of its particles that is comparable to one astronomical unit (this is the distance between the earth and the sun). Despite this constraint observations made with various spacecraft in the interplanetary medium have confirmed the existence of the solar wind and its overall collective response when it interacts with the

examine its interaction with the planets of the solar system.

**2. The solar wind as a continuous expanding gas** 

**1. Introduction** 

