# Data Science for Science & Engineering: Creating a CFD 6 Degree of Freedom Kinematics OpenFoam File

In this short article I will show you how to generate a simple 6 degree of freedom kinematic file for OpenFoam. The data file we will generate specifies the magnitude and direction of motion for a series of time points.

It is a combination of a linear displacement vector and a rotation vector about the specified center of gravity (CofG) of the object, which is defined on the dynamicMeshDict file.

This data file is later linearly interpolated by OpenFoam when running your CFD simulation, so your time points do not need to match exactly the timesteps you have set in OpenFoam. However, more time points will smooth your dynamics, increasing the accuracy of your simulation, specially if, as seen in the case below, you have a sinusoidal or other non-linear motion.

The format of the data file we want to create for OpenFoam is shown below:

`4 //number of data points in the file //Position formatting is not important. File is based on the character sequence only. //Vectors are not relative. Each vector is total displacement and total rotation.(//(time_point ( (linear displacement vector) (rotation vector roll-yaw-pitch) ) )//(seconds ( (following unit system, usually meters) (degrees) ) )(0 ( (0.25 0.50 1.0) (0.220 0.30 0.40) ) )(0.25 ( (0.50 1.0 2.0) (0.60 0.60 0.60) ) )(0.75 ( (0.75 5.0 1.0) (1.2 2.4 5.0) ) )(10.0 ( (0.1 6.0 1.0) (5.0 3.0 5.5) ) ))`

Each time point first specifies the reference time, and the following bracket groups contains the motion, first linear and then angular displacement in their own lists.

We will begin by importing the libraries we need: Numpy, Pandas, MatPlotLib to visualize our output, OS to import paths, and IPython to clean up our display.

`import numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport osfrom IPython.display import clear_output`

We will now add our parameters for the kinematics. In this case we are making a sinusoidal pitch-plunge movement, for 5 periods. The geometric angle of attack (our pitch angle) will be 15 degrees, and our plunge amplitude will be 0.05m. It is a 2D motion, so most of our degrees of freedom are actually just going to be 0.

`numberofperiods = 5alphaGeoAmp = 15plungeamplitude = 0.05period = 3.801279221resolution = 1000 #this is our time resolution, the number of data      points our file will have`

We are going to be saving this file as “pitchplunge.dat”

`#outputfileoutputfilename = os.getcwd()outputfilename += '/pitchplunge.dat'`

We generate our time vector by using the numpy function linspace, with our set resolution:

`#timetime = np.linspace(0,period*numberofperiods,resolution)`

And now it’s just a question of creating our X,Y and Z vectors for our linear and rotational displacement. We will store them in a Pandas DataFrame (our motion matrix).

`#linear displacementlinearDisplacement = pd.DataFrame(data=time,columns={'Time'})linearDisplacement['X'] = 0linearDisplacement['Y'] = plungeamplitude*np.sin(time)linearDisplacement['Z'] = 0#rotational displacementrotationalDisplacement = pd.DataFrame(data=time,columns={'Time'})rotationalDisplacement['phi'] = 0rotationalDisplacement['psi'] = 0rotationalDisplacement['theta'] = -alphaGeoAmp * np.sin(time)`

Finally we can save our motion file, adding our different brackets and required formatting as per OpenFoam’s requirements. We will use the f.write file writer, and build each line with a for loop.

`#saving motion filef= open(outputfilename,"w+")f.write(str(resolution) + '\n(\n')for i in range(0,resolution,1):        line = '('    line += str(linearDisplacement['Time'].iloc[i])    line += '(('    line += str(linearDisplacement['X'].iloc[i])    line += ' '    line += str(linearDisplacement['Y'].iloc[i])    line += ' '    line += str(linearDisplacement['Z'].iloc[i])    line += ')('    line += str(rotationalDisplacement['phi'].iloc[i])    line += ' '    line += str(rotationalDisplacement['psi'].iloc[i])    line += ' '    line += str(rotationalDisplacement['theta'].iloc[i])    line += ')))\n'        f.write(line)        if i%1000==0: #this is a tiny add-on I like to use as a test and map progress        print('Completing save ',round(i/resolution*100,2), '%')        clear_output(wait=True)            f.write(')')`

We will now plot our linear displacement to double check our motion:

`plt.plot(linearDisplacement['Time']/period, linearDisplacement['Y'])plt.xlim(0,numberofperiods)plt.ylabel("Linear Displacement /m")plt.xlabel("t/T")plt.title("Linear Displacement against non-dimensional time")`

And now our rotational displacement:

`plt.plot(rotationalDisplacement['Time']/period,rotationalDisplacement['theta'])plt.xlim(0,numberofperiods)plt.ylabel("Angular displacement /rad")plt.xlabel("t/T")plt.title("Angular displacement against non-dimensional time")`

Hope you enjoyed! You can find a video for this here, and a Jupyter Notebook to follow along here.

Scout for the IX Hispana studying vortices beyond Hadrian’s Wall.

## More from Pedro Hernández Gelado

Scout for the IX Hispana studying vortices beyond Hadrian’s Wall.