examples.Simulation_from_scratch
Link to the video (which does not work on GitHUB)
In this section we are generating a CFD simulation of a relatively complex tank from scratch. That does not mean we are generating the full pipeline ourselves, but we are instead using the presets provided by AQUAgpusph.
The simulation consists on a cylindrical tank with ellipsoidal caps. On the following figures it can be shown the sketch of a quarter of the tank, as well as a 3D view of the tank, both carried out with the 3D parametric free software, FreeCAD:
We are carrying this this task on 3 phases:
- We are creating the Initial Condition (IC).
- We are setting up the input XML files.
- We are running and post-processing the results.
The IC is the definition of the particles at the physical time,
In any case, we are first producing a 3D model of the geometry, that we are meshing and exporting. Then we are loading the exported meshes with a script and producing the particles definition files required by AQUAgpusph.
This contrast with many of the examples, where the geometries are simple enough to produce the IC directly by a Python script. It is worthy to incidentally remarks that, although the better the quality of the mesh the better results will be obtained, specially at the beginning of the simulation, herein we are demonstrating that the mesh is actually uncritical. In other words, the meshing process is a rather simple process that does not need the same level of detail required on other methods, like FEM or VOF.
This is not meant as a FreeCAD tutorial. If you want to learn more about FreeCAD, you can visit its documentation web page, which is actually really nice. For instance, you can start visiting the Basic Sketcher Tutorial.
We are indeed starting by creating a sketch on the X-Z plane. We then create a point around the coordinates (1000 mm, 0 mm, 0 mm). Do not worry if you do not nail the coordinates because now we are adding 2 constrains with respect to the origin:
- Constrain horizontal
- Constrain horizontal distance
When it asks about the disrtance, set 1000 mm.
Now add the 2 points at the right (1275.14 mm, 0 mm, 0 mm) and at the top (1000 mm, 0 mm, 159.6 mm) of the last one. Of course, constrain horizontal and constrain horizontal distance are used for the point at right, while constrain vertical and constrain vertical distance are used on the point at top. Once both points are added, you can select the tool "Create arc of ellipse", selecting the 3 points on the following order:
- 1000 mm, 0 mm, 0 mm
- 1275.14 mm, 0 mm, 0 mm
- 1000 mm, 0 mm, 159.6 mm
- 1275.14 mm, 0 mm, 0 mm
Now we need to add a last point on (0 mm, 0 mm, 159.6 mm). That point can use just vertical and horizontal constrains to keep in place.
Once that last point is added, you can create the horizontal line selecting the points:
- 0 mm, 0 mm, 159.6 mm
- 1000 mm, 0 mm, 159.6 mm
Now we are creating a helper vertical line between the following points:
- 0 mm, 0 mm, 0 mm
- 0 mm, 0 mm, 159.6 mm
So we can select the horizontal line and the ellipse arc, and use the Symmetry tool, selecting then the helper vertical line. After that we can remove the vertical line. The sketch should look like the following image:
The sketch is done. Feel free to press the Esc key, or the Close button on
the Tasks panel.
Now we are using the Draft workbench to convert the sketch on the wire that we are revolving to produce the tank. Let's start using the clone tool to copy the sketch object, so we can keep it safe just in case. Now we use the Downgrade tool on the new object (named "Sketch 2D"), so 4 edges are generated. Now we join them as a wire with the Upgrade tool.
Finally, we are using the Part workbench to create the final geometry. Just select the wire and use the Revolve tool, setting the direction (1, 0, 0), and leaving the box "Create Solid" unchecked.
Now select the revolved object and apply Part/Convert to solid.
Finally, hide all objects except the generated tank, and export it by using File/Export.... Select the *STEP with colors (*.step, .stp) format, and save it as tank.step.
To mesh the geometry we are using GMSH, that can be installed with a simple pip command:
pip install gmsh
Again, this is not meant to be a GMSH tutorial. We are instead focusing on the steps to be done to carry out the meshing task for this specific case.
Just launch GMSH and Open the exported step file. It is strongly recommended now to move the mouse cursor over the left point of the geometry, so we can know the dimensions of the tank:
In this example, the geometry was exported on millimeters, and thus the coordinates of the left point are (-1275.14, 0, 0). This situation can be tackled in many ways. Here in we are producing the mesh in millimeters and correcting that later on the Python script.
We are producing 2 different meshes, one for the volume and one for the surface. This is useful because, first it grants the possibility to use different resolutions on both parts, and second it evetually allows to generate different meshes for different boundaries.
Let's start with the surface. First we are setting the maximum length of the edges of the triangles. To this end click on Tools/Options. A window pops up, where we select Mesh on the left panel, and then General tab. Then we set a maximum element size of 20 (i.e. 20 mm = 0.02 m). We can close that window and proceed to the mesh generation, by selecting on the left panel Modules/Mesh/2D. After a couple of seconds you probably have a nice mesh for your boundary:
So use File/Save Mesh to export it. GMSH will generate a tank.msh file, please rename it as tank-Shell.msh.
Now we are repeating the same operation to produce the volume mesh. First remove the mesh by selecting on the left panel Modules/Geometry/Reload script. We are leaving the same 20 mm resolution, so just generate the volume mesh by clicking on Modules/Mesh/3D. Use again File/Save Mesh to export it as tank.msh.
Now you should have 2 mesh files: tank-Shell.msh and tank.msh.
To convert the mesh into particles we are using the following python script that we are calling Create.py:
import numpy as np
import meshio
# Input data
# ==========
TO_METERS = 1.e-3
# The level of the fluid in global coordinates
Z_LEVEL = 0.05
# A point inside the mesh, to check the mesh normals
IN_POINT = [0, 0, 0]
# Some physical/SPH parameters
g = 9.81
cs = 50.0
refd = 998.0
# Read the surface mesh file and create the boundary particles
# ============================================================
print("Parsing the shell...")
mesh = meshio.read("tank-Shell.msh")
verts = mesh.points
output = open("Surface.dat", "w")
output.write("# r.x, r.y, r.z, r.w")
output.write(", normal.x, normal.y, normal.z, normal.w")
output.write(", tangent.x, tangent.y, tangent.z, tangent.w")
output.write(", u.x, u.y, u.z, u.w")
output.write(", dudt.x, dudt.y, dudt.z, dudt.w")
output.write(", rho, drhodt, m, imove\n")
n_surf =0
for cell in mesh.cells:
def triangle(cell, verts):
a, b, c = [verts[i] * TO_METERS for i in cell]
r = np.mean([a, b, c], axis=0)
t = b - a
n = np.cross(b - a, c - a)
s = 0.5 * np.linalg.norm(n)
t /= np.linalg.norm(t)
n /= 2.0 * s
if np.dot(r - IN_POINT, n) < 0.0:
n = -n
return r, n, t, s
if cell.type != 'triangle':
if cell.type not in ['line', 'vertex', ]:
print(f"WARNING: Found a cell of type {cell.type}")
continue
for elem in cell.data:
r, n, t, s = triangle(elem, verts)
dens = refd
imove = -3
string = ("{} {} {} 0.0, " * 5 + "{}, {}, {}, {}\n").format(
r[0], r[1], r[2],
n[0], n[1], n[2],
t[0], t[1], t[2],
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
dens,
0.0,
s,
imove)
output.write(string)
n_surf += 1
output.close()
print(f"Done! {n_surf} boundary elements written")
print("Parsing the volume...")
mesh = meshio.read("tank.msh")
verts = mesh.points
output = open("Volume.dat", "w")
output.write("# r.x, r.y, r.z, r.w")
output.write(", normal.x, normal.y, normal.z, normal.w")
output.write(", tangent.x, tangent.y, tangent.z, tangent.w")
output.write(", u.x, u.y, u.z, u.w")
output.write(", dudt.x, dudt.y, dudt.z, dudt.w")
output.write(", rho, drhodt, m, imove\n")
n_vol =0
for cell in mesh.cells:
def tetra(cell, verts):
a, b, c, d = [verts[i] * TO_METERS for i in cell]
r = np.mean([a, b, c, d], axis=0)
matrix = np.array([a - d, b - d, c - d])
vol = np.abs(np.linalg.det(matrix)) / 6.0
return r, vol
if cell.type != 'tetra':
if cell.type not in ['triangle', 'line', 'vertex', ]:
print(f"WARNING: Found a cell of type {cell.type}")
continue
for elem in cell.data:
r, vol = tetra(elem, verts)
if r[-1] > Z_LEVEL:
# Discard the particles above the free surface
continue
p = refd * g * (Z_LEVEL - r[-1])
dens = refd + p / (cs * cs)
mass = refd * vol
imove = 1
string = ("{} {} {} 0.0, " * 5 + "{}, {}, {}, {}\n").format(
r[0], r[1], r[2],
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
0.0, 0.0, 0.0,
dens,
0.0,
mass,
imove)
output.write(string)
n_vol += 1
output.close()
print(f"Done! {n_vol} fluid particles written")
We are not going into the detail of the script. Just tke into account that it depends on meshio and numpy. Both can be installed with pip.
There are several inputs to this script:
-
TO_METERS: During the meshing process we checked that the geometry is defined on millimeters, so we must convert it to meters, multiplying by 1e-3. -
Z_LEVEL: The fluid level within the tank -
IN_POINT: A point that belongs to the inner part of the tank. On AQUAgpusph it is very important that the boundary element normals point outwards the fluid. -
g: The gravity acceleration vector length. -
cs: The spped of sound. On SPH this value is typically not the actual physical value, but a smaller one. That allows to accelerate the computation. -
refd: The fluid density of reference. In this case water, i.e.$\rho_0 = 998 ,, \mathrm{kg}/\mathrm{m}^3$ .
So feel free to execute the script typing python Create.py.
An output like the following will be printed:
Parsing the shell...
Done! 15118 boundary elements written
Parsing the volume...
Done! 77320 fluid particles written
Take note of the number of boundary elements and fluid particles written. We are using those numbers, as well as the script inputs and the meshing resolution --0.02 m--, when creating the input XML files.
As aforementioned, we are using the presets provided by AQUAgpusph to setup the pipeline, that we are then editing with our custom XML files.
To this end we start creating a Main.xml file with the following content:
<?xml version="1.0" ?>
<sphInput>
<Include file="/YOUR_AQUAGPUSPH_PATH/resources/Presets/basic.xml" />
<Include file="/YOUR_AQUAGPUSPH_PATH/resources/Presets/basic/time_scheme/midpoint.xml" />
<Include file="/YOUR_AQUAGPUSPH_PATH/resources/Presets/basic/time_scheme/midpoint/autostop.xml" />
<Include file="/YOUR_AQUAGPUSPH_PATH/resources/Presets/basic/time_scheme/midpoint/residues_profile.xml" />
<Include file="/YOUR_AQUAGPUSPH_PATH/resources/Presets/basic/domain.xml" />
<Include file="/YOUR_AQUAGPUSPH_PATH/resources/Presets/cfd.xml" />
<Include file="/YOUR_AQUAGPUSPH_PATH/resources/Presets/cfd/deltaSPH.xml" />
<Include file="/YOUR_AQUAGPUSPH_PATH/resources/Presets/cfd/deltaSPH-full.xml" />
<Include file="/YOUR_AQUAGPUSPH_PATH/resources/Presets/cfd/BIe.xml" />
<Include file="/YOUR_AQUAGPUSPH_PATH/resources/Presets/basic/timing.report.xml" />
<Include file="/YOUR_AQUAGPUSPH_PATH/resources/Presets/basic/performance.report.xml" />
<Include file="Settings.xml" />
<Include file="SPH.xml" />
<Include file="Time.xml" />
<Include file="Fluids.xml" />
<Include file="Sensors.xml" />
<Include file="Motion.xml" />
</sphInput>
where YOUR_AQUAGPUSPH_PATH shall be replaced by the path where you have your
compiled AQUAgpusph.
As can be appreciated, we are just including other XML files on Main.xml. Along this line there are 2 different parts, a first part where we are loading the needed AQUAgpusph presets, and a second part we are loading our actual input files.
We refer to the presets documentation to know more about the loaded presets. As a summary, we are setting up a Computational Fluid Dynamics (CFD) simulation, with the Euler Implicit Mid-Point (IMP) time integrator and the associated conservative wall boundary condition. That combination confeers unconditional stability to SPH. See:
Cercos-Pita, J. L., Merino-Alonso, P. E., Calderon-Sanchez, J., & Duque, D. (2023). The role of time integration in energy conservation in Smoothed Particle Hydrodynamics fluid dynamics simulations. European Journal of Mechanics-B/Fluids, 97, 78-92.
Cercos-Pita, J. L., Duque, D., Merino-Alonso, P. E., & Calderon-Sanchez, J. (2024). Boundary conditions for SPH through energy conservation. Computers & Fluids, 285, 106454.
We are anyway adding
So let's create our actual input files.
A rather standard file, you can create it with the following content:
<?xml version="1.0" ?>
<sphInput>
<Settings>
<Device platform="0" device="0" type="ALL" />
</Settings>
</sphInput>
In this file is where we are actually configuring the main parameters of the
simulation.
To this end we need to collect some data, like the
volume mesh resolution --dr=0.02--, and the
script inputs --g=9.81, cs=50--.
<?xml version="1.0" ?>
<sphInput>
<Variables>
<Variable name="g" type="vec" value="0.0, 0.0, -9.81, 0.0" />
<Variable name="dr" type="float" value="0.02" />
<Variable name="hfac" type="float" value="4.0" />
<Variable name="h" type="float" value="hfac * dr" />
<Variable name="cs" type="float" value="50" />
<Variable name="courant" type="float" value="0.5" />
<Variable name="domain_min" type="vec" value="-2.0, -2.0, -2.0, 0.0" />
<Variable name="domain_max" type="vec" value="2.0, 2.0, 2.0, 0.0" />
<Variable name="iter_midpoint_max" type="unsigned int" length="1" value="5"/>
<Variable name="relax_midpoint" type="float" length="1" value="0"/>
<Variable name="Residual_midpoint_max" type="float" length="1" value="1.e-2"/>
</Variables>
</sphInput>
Please, keep attention to these parameters, it is very important that they
match the used values before.
For instance, domain_min and domain_max should be large enought to grant
that the tank never touchs the computational domain boundaries.
The last 3 parameters are related to the IMP. In general you can set similar values and tweak afterwards assessing if it properly converged, considering to this end the midpoint.out file generated by AQUAgpusph.
This is quite dependant on the motion. Since we are modeling a 4 seconds period motion, we can consider the following XML file:
<?xml version="1.0" ?>
<sphInput>
<Timing>
<Option name="End" type="Time" value="4.0" />
<Option name="Output" type="FPS" value="30" />
</Timing>
</sphInput>
Here we define both the tank boundary and the volume inside.
To generate this file we shall collect some info from the
particles generation script, like the
density --refd=998-- and the output number of boundary elements --n=15118--
and fluid particles --n=77320--:
<?xml version="1.0" ?>
<sphInput>
<ParticlesSet n="15118">
<Scalar name="refd" value="998" />
<Scalar name="visc_dyn" value="0.000894" />
<Scalar name="delta" value="0.1" />
<Load format="FastASCII" file="Surface.dat" fields="r, normal, tangent, u, dudt, rho, drhodt, m, imove" />
<Save format="VTK" file="output_boundary" fields="r, normal, tangent, u, dudt, rho, drhodt, m, p, imove" />
</ParticlesSet>
<ParticlesSet n="77320">
<Scalar name="refd" value="998" />
<Scalar name="visc_dyn" value="0.000894" />
<Scalar name="delta" value="0.1" />
<Load format="FastASCII" file="Volume.dat" fields="r, normal, tangent, u, dudt, rho, drhodt, m, imove" />
<Save format="VTK" file="output_fluid" fields="r, normal, tangent, u, dudt, rho, drhodt, m, p, imove, shepard" />
</ParticlesSet>
</sphInput>
We have also added some other information, like the dynamic viscosity and the
It is quite important to remember the order on which the particles sets are added. Right now we have the tank surface with the set index 0, and the fluid particles with the set index 1.
On a similar fashion, we can add a third particles set --with the set index 2--, which contains 2 sernsors:
<?xml version="1.0" ?>
<sphInput>
<ParticlesSet n="2">
<Scalar name="refd" value="998.0" />
<Scalar name="visc_dyn" value="0.06700523924875554" />
<Scalar name="delta" value="0.1" />
<Load format="FastASCII" file="Sensors.dat" fields="r, normal, u, dudt, rho, drhodt, m, imove" />
<Save format="VTK" file="sensors" fields="r, normal, u, dudt, rho, drhodt, m, p, imove" />
</ParticlesSet>
<Reports>
<Report type="particles" name="Sensors" fields="p,shepard" path="sensors.out" set="2" ipf="0" fps="300"/>
</Reports>
</sphInput>
2 facts might capture our attention now.
First, we have a new <Reports>, were we are asking AQUAgpusph that export the
information collected at the sensors at a 300 Frames Per Second (FPS) pace
(please, notice that the report is aiming the particles set with the index 2).
Second, we do not have a Sensors.dat file, so we need to create it with
the following content:
1.27 0 0 0, 0 0 0 0, 0 0 0 0, 0 0 0 0, 998.0, 0, 0, 0
-1.27 0 0 0, 0 0 0 0, 0 0 0 0, 0 0 0 0, 998.0, 0, 0, 0
The sensors are just 2 particles with null mass and imove=0 flag.
We have placed them on the left and right borders of the tank.
This file is more fun! At first glance we just need to move the tank wall. Unfortunately it is a bit more complex, because we actually need to move the tank wall --with the set index 0-- and the sensors --with the set index 2--.
Let's see the Motion.xml file and we discuss how the task is done:
<?xml version="1.0" ?>
<sphInput>
<Include file="/YOUR_AQUAGPUSPH_PATH/resources/Presets/cfd/motion.xml" prefix="iset0_" />
<Include file="/YOUR_AQUAGPUSPH_PATH/resources/Presets/cfd/motion.xml" prefix="iset2_" />
<Variables>
<Variable name="motion_r" type="vec" value="0.0, 0.0, 0.0, 0.0" />
</Variables>
<Tools>
<Tool action="replace" name="iset0_cfd motion data" type="python" path="Motion.py"/>
<Tool action="insert" before="iset0_cfd motion data" type="set_scalar" name="motion_iset0" in="motion_iset" value="0"/>
<Tool action="replace" name="iset2_cfd motion data" type="python" path="Motion.py"/>
<Tool action="insert" before="iset2_cfd motion data" name="motion_iset1" type="set_scalar" in="motion_iset" value="2"/>
<Tool action="insert" after="iset0_cfd motion data,iset2_cfd motion data" name="Motion report0" type="report_screen" fields="t,motion_a,motion_dadt"/>
</Tools>
</sphInput>
Again, remember to replace YOUR_AQUAGPUSPH_PATH by the path where you have
your compiled AQUAgpusph.
As it can be appreciated we are loading the motion preset provided by
AQUAgpusph twice, with different prefixes.
Indeed we are using the prefix iset0_ to control the motion of the tank wall,
and the prefix iset1_ to control the motion of the sensors.
The first thing we do is telling AQUAgpusph to consider our own Python script
Motion.py for both motions (notice that we are replacing those tools, i.e.
action="replace").
Then we are inserting a tool before each Motion.py execution to let
AQUAgpusph which particles set is affected by the motion (notice that we are
inserting those tools, i.e. action="insert" before="...").
So it is time to give AQUAgpusph a Motion.py file:
import numpy as np
import aquagpusph as aqua
A = 5.0
T = 4.0
def main():
omega = 2.0 * np.pi / T
angle = np.radians(A)
t = aqua.get("t")
a = np.zeros(4, dtype=np.float32)
a[1] = angle * np.sin(omega * t)
dadt = np.zeros(4, dtype=np.float32)
dadt[1] = angle * omega * np.cos(omega * t)
aqua.set("motion_a", a)
aqua.set("motion_dadt", dadt)
return True
We are indeed setting up a sinusoidal motion of 5 degrees amplitude and 4
seconds period, around the y axis.
To this end we are asking the physical time from AQUAgpusph
--aqua.get("t")--, and setting the EulerXYZ angles and angular velocities
afterwards --aqua.set("motion_a", a) and aqua.set("motion_dadt", dadt)--.
Running AQUAgpusph is as simple as typing:
/YOUR_AQUAGPUSPH_PATH/bin/AQUAgpusph -d 3 -i Main.xml
Please, notice that to stop AQUAgpusph and run it again from the beginning, you should clean all the output files.
At the time of postprocessing we want to produce 2 outputs: A plot of the pressure at the sensors and an animation of the simulation.
To produce the plot we can use the following Python script:
import numpy as np
from scipy.signal import savgol_filter
import matplotlib.pyplot as plt
def readFile(filepath):
f = open(filepath, "r")
lines = f.readlines()
f.close()
data = []
for l in lines[1:-1]: # Skip the last line, which may be unready
l = l.strip()
while l.find(' ') != -1:
l = l.replace(' ', ' ')
fields = l.split(' ')
try:
data.append(map(float, fields))
except:
continue
return [list(d) for d in zip(*data)]
fig = plt.figure()
ax = fig.add_subplot(111)
data = readFile('sensors.out')
t = data[0]
p1 = [d * s * 2 for d, s in zip(data[1], data[2])]
p2 = [d * s * 2 for d, s in zip(data[3], data[4])]
ax.plot(t,
p1,
label=r'$x=1.27$',
color="blue",
linewidth=1.0)
ax.plot(t,
p2,
label=r'$x=-1.27$',
color="red",
linewidth=1.0)
ax.grid()
ax.legend(loc='best')
ax.set_xlim(0, 4)
ax.set_ylim(-500, 4000)
ax.set_autoscale_on(False)
ax.set_xlabel(r"$t \, [\mathrm{s}]$")
ax.set_ylabel(r"$p \, [\mathrm{Pa}]$")
plt.tight_layout()
plt.show()
It is worth to highlight the operation to extract the pressure,
p1 = [d * s * 2 for d, s in zip(data[1], data[2])].
Effectively, on AQUAgpusph the sensors value is renormalized with the Shepard
renormalization factor.
However, the consistent way to get the pressure along the wall, when
conservative boundaries are considered, is to use a renormalization factor of
2.
The script should produce a plot similar to the following one:
After the initial shock --caused by the sudden change from null angular velocity to the maximum one-- the pressure field quickly converges and a smooth signal is retrieved.
To produce the animation we are using Paraview. Once again, this is not meant as a tutorial of Paraview, but a summary of the basic steps to get a first visualization of the boundary elements and the fluid particles.
Just lunch your Paraview and open the file output_boundary.vtu.series. Then on the Properties panel --at left-- click on Apply. You are probably seeing nothing, so on the top panel click on Edit Color Map, and select a color that constrasts with the background, e.g. chose black if your backgorund is white.
Now open output_fluid.vtu.series, and click on Apply again. This time, on the top panel replace Solid Color by p, and if the color bar is not shown on the 3D view, click on Toogle Color Legend Visibility.
Now you can click and drag on the 3D view to get a more exciting view angle.
Then you can tweak the color scale by clicking on Edit Color Map. On the panel poping up at right, click on the folder with a heart: Choose preset. Select a preset that you like.
Press the play button and enjoy! You can also save your animation as a video!
Paraview is a quite powerful tool, so feel free to play and check the official tutorials.
On this tutorial we have traversed the process of generating a CFD simulation of a relatively complex tank filled with water.
To this end we have used the help of a 3D CAD and a meshing tool. The meshing process has been though quite simple. We just needed to tweak the mesh resolution, as the maximum length of the triangles and tetrahedra edges.
We also visited the process to generate the AQUAgpusph XML input, as well as the auxiliar files.
We have finally ran the simulation and postprocessed the results. Along this line, on the AQUAgpusph examples you can find alternative solutions to the plotters that are able to run on real time.