OpenFOAM postprocessing: a few shortcuts

Run-time postprocessing

For pump performance calculation (and in many other cases as well), one usually needs data about pressure rise (head), flow, and power. If cavitation is an issue, maximum and minimum values for pressure come in handy, also maximum and minimum velocity can tell about convergence or numerical issues. All of these values are best calculated at runtime so it’s best to put the settings in controlDict.functionObjects.

// impeller forces for torque/power/viscous losses/axial and radial forces/mechanical losses
forces_impeller
{
    type    forces;
    libs    ("libforces.so");
    patches (impeller_walls);

    rho     rhoInf;
    rhoInf  985; // depends on temperature, also check viscosity in constant/transportProperties

    CofR (0 0 0);
}

// pressure difference between inlet and outlet for head calculation
delta_p
{
    type    fieldValueDelta;
    libs    ("libfieldFunctionObjects.so");
    operation   subtract;

    region1
    {
        type        surfaceFieldValue;
        operation   average;
        fields      (p);
        writeFields no;
        regionType  patch;
        name        outlet;
    }

    region2
    {
        type        surfaceFieldValue;
        operation   average;
        fields      (p);
        writeFields no;
        regionType  patch;
        name        inlet;
    }
}

// flow (is specified in 0/U but here it's added for an additional check)
Q
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");

    writeFields     false;
    surfaceFormat   foam;

    regionType      patch;
    name            outlet;

    operation       areaIntegrate;

    fields
    (
        U
    );
}

// maximum values: maximum velocity should be around impeller circumferential velocity
// difference between minimum pressure and inlet pressure is NPSH
maxValues
{
    type        fieldMinMax;
    libs        ("libfieldFunctionObjects.so");

    writeToFile no;
    log         yes;
    location    no;
    mode        magnitude;
    fields      (U p);
}

Values are calculated after every time step and results are saved in postProcessing/ directory.

Post-postprocessing

There is also a handful of other quantities of interest like yPlus for mesh quality check and Q and vorticity for turbulence/losses visualisation. The problem with these quantities are that the whole field gets saved each timestep. If your writeInterval in controlDict is not 1 OpenFOAM makes a mess with timestep directories. Now every timestep has a directory, some only have yPlus and Q fields and others have all calculated quantities. This makes working with Paraview difficult so it’s better to postprocess these fields separately.

mpirun -np 4 pimpleFoam -parallel -postProcess -latestTime -func yPlus
mpirun -np 4 pimpleFoam -parallel -postProcess -latestTime -func Q
mpirun -np 4 pimpleFoam -parallel -postProcess -latestTime -func vorticity

The postProcessing/ directory

If for any reason the simulation is interrupted and then restarted from latestTime (check your controlDict), new files are added to postProcessing directory. Every added file gets an additional timestep timestamp added to its name. With many restarts the files multiply and make gathering data difficult.

Gathering data

This is a script for reading postprocessing data from all created files mentioned above. It gathers the values in a single Python dictionary and pickles it to a file. It can then be read in Paraview with programmable filter/source for visualization or read from another Python script for plotting and so on.

from subprocess import Popen, PIPE
import os

import re

import pickle

# case times
times = []

# data to be written
data = {
    # [time]: {
    #   delta_p: [number],
    #   Q: [number]
    #   forces_impeller: [F_pressure, F_viscous, M_pressure, M_viscous]
    # }   
}

for t in times:
    data[t] = {}

def read_pp_data(pp_object, parse_function):
    for froot, _, fnames in os.walk(os.path.join('postProcessing', pp_object)):
        for fname in fnames:
            filename = os.path.join(froot, fname)
            if not os.path.isfile(filename):
                continue

            with open(filename, 'r') as f:
                print filename
                
                # skip lines with comment
                for line in f:
                    if line[0] == '#':
                        continue

                    # timestep
                    time = line.split()[0]

                    # add to an existing entry or create a new one
                    if time not in times:
                        times.append(time)
                        data[time] = {}

                    # add data to collection
                    data[time][pp_object] = parse_function(line)

# parsing files: each value has different format
split_re = re.compile('[ \t\n\r()]+') # split: whitespace, parenthesis

def parse_delta_p(line): # delta_p function object
    values = line.split()

    return float(values[1])


def parse_Q(line): # surfaceFieldValue areaIntegrate (flow)
    return float(split_re.split(line)[2])

def parse_forces_impeller(line): # forces functionObject
    # forces format:
    # 0                   1 2 3    4 5 6   7 8 9            10 11 12 13 14 15 16 17 18
    # Time         forces(pressure viscous porous)	moments(pressure viscous porous)
    forces = split_re.split(line)

    # interested in: pressure and viscous forces and moments
    indexes = (1, 2, 3, 4, 5, 6, 10, 11, 12, 13, 14, 15)
    return [float(forces[i]) for i in indexes]
    

read_pp_data('delta_p', parse_delta_p)
read_pp_data('Q', parse_Q)
read_pp_data('forces_impeller', parse_forces_impeller)

# change text numbers to floats
for k in data.keys():
    data[float(k)] = data[k]
    del data[k]

# save the data to be read in paraView (?)
with open('postProcessing.pickle', 'wb') as f:
    pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)

Plotting collected data

Loading a pickled dictionary is a piece of cake. I usually cut the first few results away since timesteps usually don’t converge and therefore values are unphysical.

import pickle
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats

import collect_pp # import the above script to avoid running it manually

# start only after transients decay
t_start = 0.001 # here interesting data starts

# pickled data format:
# 0.000444: {
#   'Q': -0.000222930392,
#   'delta_p': 200.0, 
#   'forces_impeller': [-70.6219451, -0.307142376, 0.116242584, 0.0130158996, -0.0167682449, 0.00424048528, -0.320480923, -0.00105645189, -0.00185268198, -0.0417415834, -4.51203576e-05, -0.00011419315]
# },

with open('postProcessing.pickle', 'rb') as f:
    data = pickle.load(f)

# times
times = data.keys()
times.sort()

n_start = np.searchsorted(np.array(times), t_start)
times = times[n_start:]

p_H = [] # head
p_Q = [] # flow
p_Mp = [] # pressure torque
p_Mv = [] # viscous torque
p_Fx = [] # axial force

for t in times:
    try:
        p_H.append(data[t]['delta_p']) 
        p_Q.append(-data[t]['Q']) 
        Mp.append(-data[t]['forces_impeller'][6]) 
        Mv.append(-data[t]['forces_impeller'][9]) 
        Fx.append(data[t]['forces_impeller'][0]) 
        
        # time
        p_times.append(t)
    except KeyError:
        continue

Filtering stuff

Ocassionally a timestep doesn’t converge and so calculated forces/moments/pressures are way too high or too low. That ruins the day for graphing and these values should be filtered. Actually, it would be better to find out why calculation doesn’t converge in the first place but I haven’t had the time yet for a proper ‘debugging’ session. So I use a little trim utility from scipy:

# filter out the outliers; weird unexplained results that occur randomly
# first, determine which are they. 
def filter_data(values):
    # use trimmed values for average or it doesn't 'mean' anything
    avg = np.average(scipy.stats.trim_mean(values, 0.1))
    
    return np.array(np.where(np.abs(values - avg) > abs(0.5*avg)))[0]

i_filter = np.concatenate((
        filter_data(p_H),
        filter_data(p_Q),
        filter_data(p_Mp),
        filter_data(p_Mv),
        filter_data(p_Fx)))

i_filter = list(set(i_filter))

# delete all data with given indexes
p_times = np.delete(p_times, i_filter)
p_H = np.delete(p_H, i_filter)
p_Q = np.delete(p_Q, i_filter)
p_Mp = np.delete(p_Mp, i_filter)
p_Mv = np.delete(p_Mv, i_filter)
p_Fx = np.delete(p_Fx, i_filter)

Now you’re ready to do whatever with the data. Like plotting:

# plot results
fig, ax1 = plt.subplots()

ax1.plot(p_times, p_H, 'b-')
ax1.set_xlabel('Time (s)')

ax1.set_ylabel('H [dm]', color='b')
ax1.tick_params('y', colors='b')

ax2 = ax1.twinx()

ax2.plot(p_times, p_Mp, 'r.')
ax2.set_ylabel('M [Nm]', color='r')
ax2.tick_params('y', colors='r')

Enjoy!

One Comment

  1. Pingback: OpenFOAM and Rotating Machinery: Relative Velocity | Damogran Labs

Leave a Reply