# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ## Generating OpenFoam 6DoF Files
#
# This Jupyter notebook 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 np
import pandas as pd
import matplotlib.pyplot as plt
import os
from 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 = 5
alphaGeoAmp = 15
plungeamplitude = 0.05
period = 3.801279221
resolution = 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"

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

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

#time
time = 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 displacement
linearDisplacement = pd.DataFrame(data=time,columns={'Time'})
linearDisplacement['X'] = 0
linearDisplacement['Y'] = plungeamplitude*np.sin(time)
linearDisplacement['Z'] = 0

#rotational displacement
rotationalDisplacement = pd.DataFrame(data=time,columns={'Time'})
rotationalDisplacement['phi'] = 0
rotationalDisplacement['psi'] = 0
rotationalDisplacement['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 file

f= 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")

StarcoderdataJupyterScriptsDedupFiltered