Application (Preccinsta Combustion Chamber)

Post-processing of a CFD simulation with an unstructured multi-element dataset.

This tutorial will help you build a post-processing on the Preccinsta test case. Please download the associated data at https://www.cerfacs.fr/antares/downloads/application2_tutorial_data.tgz

untar the archive in the working directory

copy the app2_combustion.ipynb in the directory application2

Preliminary steps

import os, numpy, and antares packages

create a directory named OUTPUT

import os

import numpy as np

import antares

if not os.path.isdir('OUTPUT'):
    os.makedirs('OUTPUT')

Part I: Read Input Data

The first step is to create a base from the data of a CFD simulation of the Preccinsta combustion chamber. This data structure will be used in the following post-processing.

The first step is to create a Base named base by reading the grid from files MESH/mesh.mesh.h5 (hdf_avbp format) as shared data.

reader = antares.Reader('hdf_avbp')
reader['filename'] = os.path.join('MESH', 'mesh.mesh.h5')
reader['shared'] = True
base = reader.read()

Now add to your existing base the solution stored in the file SOLUT/sol_ite0001000.h5 (hdf_avbp format).

reader = antares.Reader('hdf_avbp')
reader['base'] = base
reader['filename'] = os.path.join('SOLUT', 'sol_ite0001000.h5')
reader.read()

Part II: Extract Skin Data

This part details the steps to extract the skin of the combustion chamber.

To extract the walls of the Preccinsta combustion chamber, create a Base skin_base from the Family Patches which contains the boundaries of your simulation.

print(base.families)

skin_base = base[base.families['Patches']]

To remove the atmospheric boundaries and only keep the combustion chamber wall, use the threshold treatment (see documentation to know how it works) to create a Base thres_skin_base that contains the cells matching the following conditions:

x < 0.11
-0.045 < y < 0.045
-0.045 < z < 0.

Printing the grid_points attribute of the resulting Base object shows that the size of the extracted dataset is much smaller.

clip = antares.Treatment('threshold')
clip['base'] = skin_base
clip['variables'] = ['x', 'y', 'z']
clip['threshold'] = [(None, 0.11), (-0.045, 0.045), (-0.045, 0.)]
clip_skin_base = clip.execute()

Finally, to visualize the extracted skin data, write the base thres_skin_base in VTK binary format (compatible with Paraview software) in the directory OUTPUT/SKIN/.

writer = antares.Writer('bin_vtk')
writer['base'] = clip_skin_base
writer['filename'] = os.path.join('OUTPUT', 'SKIN', 'Patches_<zone>')
writer.dump()

Part III: Create Cut

This part will teach you how to extract the tangential velocity on an transversal geometrical cut.

We want to visualize the tangential velocity on a geometrical cut. First, to reduce the size of the dataset, remove the atmospheric domain by keeping only the cells matching the following conditions:

x < 0.11
-0.045 < y < 0.045
-0.045 < z < 0.045

The base obtained after the threshold will be named thres_base.

clip = antares.Treatment('threshold')
clip['base'] = base
clip['variables'] = ['x', 'y', 'z']
clip['threshold'] = [(None, 0.11), (-0.045, 0.045), (-0.045, 0.045)]
clip_base = clip.execute()

print(clip_base.grid_points)

Use the cut treatment (see documentation for more information) to perform a geometrical plane cut in your domain at a constant z value of 0. The base obtained after cut will be named cut_base.

plane_cut = antares.Treatment('acut')
plane_cut['base'] = clip_base
plane_cut['type'] = 'plane'
plane_cut['origin'] = [0., 0., 0.]
plane_cut['normal'] = [0., 0., 1.]
cut_base = plane_cut.execute()

Use the function compute() to compute the tangential velocity (the velocity minus the axial contribution) on your base. To finally visualize the cut, write the base cut_base in VTK binary format (compatible with Paraview software) in directory OUTPUT/ZCUT/.

cut_base.set_computer_model('internal')
cut_base.compute('tangential_velocity=(velocity_Y**2+velocity_Z**2)** 0.5')

writer = antares.Writer('hdf_antares')
writer['base'] = cut_base
writer['filename'] = os.path.join('OUTPUT', 'zcut')
writer.dump()

Part IV: Plot Data Over Line

The objective of this part is to plot the Mach number distribution along the y-axis at the swirler outlet.

We want now to vizualize Mach number distribution along the line at x = 0.005 and z = 0.0. So first, use the cut treatment on the cut_base (from part III) to perform a cut at x = 0.005. Then use the merge treatment (see documentation to learn how it works) to merge and the unwrapline treatment to reorganize the line points. The base obtained will be named line_base.

plane_cut = antares.Treatment('cut')
plane_cut['base'] = cut_base
plane_cut['type'] = 'plane'
plane_cut['origin'] = [0.005, 0., 0.]
plane_cut['normal'] = [1., 0., 0.]
line_base = plane_cut.execute()

merge = antares.Treatment('merge')
merge['base'] = line_base
merge['duplicates_detection'] = True
# merge['tolerance_decimals'] = 12
line_base1 = merge.execute()

unw = antares.Treatment('unwrapline')
unw['base'] = line_base1
unw['sorting_variable'] = 'y'
line_base = unw.execute()

To visualize data stored in the base line_base use matplotlib. Plot the Mach number distribution versus y.

import matplotlib.pyplot as plt
plt.plot(line_base[0][0]['y'], line_base[0][0]['Mach'],lw=2.5)
plt.grid()
plt.xlabel('{}'.format('y'), fontsize=18)
plt.ylabel('{}'.format('Mach'), fontsize=18)
plt.autoscale(axis='x')
plt.autoscale(axis='y')
plt.show()
plt.close()

Write the base line_base in directory OUTPUT/LINE/ in the column format.

writer = antares.Writer('column')
writer['base'] = line_base
writer['filename'] = os.path.join('OUTPUT', 'line.dat')
writer.dump()

Part V: Compute Average Quantities

Finally, it can be interesting to compute the average thickening and efficiency over the flame front. In fact it allows to evaluate the quality of your simulation relatively to grid refinement and flame front resolution.

As we want to compute the average thickening and efficiency over the flame front, the first step is to extract this flame front data from the 3D dataset. It can be assimilated to the isosurface of temperature 1500 Kelvin. To do so, use the isosurface treatment to create the base iso_base from clip_base (from part III).

iso = antares.Treatment('isosurface')
iso['base'] = clip_base
iso['variable'] = 'temperature'
iso['value'] = 1500.
iso_base = iso.execute()

writer = antares.Writer('hdf_antares')
writer['base'] = iso_base
writer['filename'] = os.path.join('OUTPUT', 'isoT')
writer.dump()

To compute surface averages, we need to know the cell surfaces. see Base method compute_cell_normal().

Now that the variable ‘surface’ is available in the base, use it to compute the product of Efficiency and Thickening and surface (needed for surface averaging)

Finally use the numpy function sum() to integrate data and compute the average efficiency and thickening.

t = antares.Treatment('integration')
t['base'] = iso_base
integral_base = t.execute()

effic = integral_base[0][0]['Efficiency'][0]
thick = integral_base[0][0]['Thickening'][0]
surface = integral_base[0][0]['volume'][0]

effic /= surface
thick /= surface

print('\n >>>  Isosurface average :')
print('        - Efficiency : {:.2f}'.format(effic))
print('        - Thickening : {:.2f}'.format(thick))
# >>>  Isosurface average :
#        - Efficiency : 3.55
#        - Thickening : 11.62