ALEMO - Homework 1
François THOMAS
Table of Contents
1.0 - Background/Motivations
2.0 - Problem Statement
3.0 - Dataset Description
4.0 - Methods
4.1 - Data Parsing
import numpy as np
import pylab as plt
#IMPORTING DATA
file_5k = '/work/Files_for_Flow/Data/Retau_5k_basic_stats.txt'
file_10k = '/work/Files_for_Flow/Data/Retau_10k_basic_stats.txt'
file_20k = '/work/Files_for_Flow/Data/Retau_20k_basic_stats.txt'
file_30k = '/work/Files_for_Flow/Data/Retau_30k_basic_stats.txt'
file_40k = '/work/Files_for_Flow/Data/Retau_40k_basic_stats.txt'
#PARSING THE DATA
#FOR Re_Tau = 5K
Re_tau_5 = np.loadtxt(file_5k,delimiter=' ', usecols=(4), max_rows=1)
R_5 = np.loadtxt(file_5k,delimiter=' ', usecols=(7), skiprows=3, max_rows=1)
y_5,U_5=np.loadtxt(file_5k,delimiter=' ',usecols=(0,3),skiprows=6,unpack=True)
plt.plot(U_5,R_5-y_5,'.-',label='For Re_Tau = 5k')
#FOR Re_Tau = 10K
Re_tau_10 = np.loadtxt(file_10k,delimiter=' ', usecols=(4), max_rows=1)
R_10 = np.loadtxt(file_10k,delimiter=' ', usecols=(7), skiprows=3, max_rows=1)
y_10,U_10=np.loadtxt(file_10k,delimiter=' ',usecols=(0,3),skiprows=6,unpack=True)
plt.plot(U_10,R_10-y_10,'.-',label='For Re_Tau = 10k')
#FOR Re_Tau = 20K
Re_tau_20 = np.loadtxt(file_20k,delimiter=' ', usecols=(3), max_rows=1)
R_20 = np.loadtxt(file_20k,delimiter=' ',usecols=(7),skiprows=3,max_rows=1)
y_20 = np.loadtxt(file_20k,delimiter=' ',usecols=(0),skiprows=6)
U_20 = np.loadtxt(file_20k,delimiter=' ',usecols=(3),skiprows=6,max_rows=9)
U_20=np.append(U_20,np.loadtxt(file_20k,delimiter=' ',usecols=(2),skiprows=15))
plt.plot(U_20,R_20-y_20,'.-',label='For Re_Tau = 20k')
#FOR Re_Tau = 30K
Re_tau_30 = np.loadtxt(file_30k,delimiter=' ', usecols=(3), max_rows=1)
R_30 = np.loadtxt(file_30k,delimiter=' ',usecols=(7),skiprows=3,max_rows=1)
y_30 = np.loadtxt(file_30k,delimiter=' ',usecols=(0),skiprows=6)
U_30 = np.loadtxt(file_30k,delimiter=' ',usecols=(3),skiprows=6,max_rows=2)
U_30=np.append(U_30,np.loadtxt(file_30k,delimiter=' ',usecols=(2),skiprows=8))
plt.plot(U_30,R_30-y_30,'.-',label='For Re_Tau = 30k')
#FOR Re_Tau = 40K
Re_tau_40 = np.loadtxt(file_40k,delimiter=' ', usecols=(3), max_rows=1)
R_40 = np.loadtxt(file_40k,delimiter=' ', usecols=(7), skiprows=3, max_rows=1)
y_40,U_40=np.loadtxt(file_40k,delimiter=' ',usecols=(0,2),skiprows=6,unpack=True)
plt.plot(U_40,R_40-y_40,'.-',label='For Re_Tau = 40k')
#ADDING FIGURE LABELS
plt.title("Parsed Data")
plt.xlabel("Velocity U [m/s]")
plt.ylabel("Distance from the pipe centerline y [m]")
plt.legend(loc='upper right')
plt.grid()
Note: unfortunately I could not use a loop to parse the data of all files, as Deepnote does not allow for multiple spaces delimiters, therefore the position of the data in the files was changed from file to file as the number of spaces in between the values varies. As an alternative could not be found, the data parsing method is quite simple and repetitive.
4.2 - Preprocessing
4.3 - Data Preprocessing
4.4 - Linear Regression
from scipy import stats
#FOR 5K
a_5=stats.linregress(np.log(y_5/R_5),np.log(U_5))
print('For Re_Tau=5k: ',f'R-squared = {a_5.rvalue**2:.3f},',f'Umax = {np.exp(a_5.intercept):.3f}, ',f'n = {1/a_5.slope:.3f}')
plt.plot(np.log(y_5/R_5),np.log(U_5),'ko',label='Given Data')
plt.plot(np.log(y_5/R_5),a_5.intercept + (a_5.slope)*np.log(y_5/R_5),label='Fitted Line for $Re_{\mathcal{T}}=5k$')
#FOR 10K
a_10=stats.linregress(np.log(y_10/R_10),np.log(U_10))
print('For Re_Tau=10k:',f'R-squared = {a_10.rvalue**2:.3f},',f'Umax = {np.exp(a_10.intercept):.3f},',f'n = {1/a_10.slope:.3f}')
plt.plot(np.log(y_10/R_10),np.log(U_10),'ko')
plt.plot(np.log(y_10/R_10),a_10.intercept + (a_10.slope)*np.log(y_10/R_10),label='Fitted Line for $Re_{\mathcal{T}}=10k$')
#FOR 20K
a_20=stats.linregress(np.log(y_20/R_20),np.log(U_20))
print('For Re_Tau=20k:',f'R-squared = {a_20.rvalue**2:.3f},',f'Umax = {np.exp(a_20.intercept):.3f},',f'n = {1/a_20.slope:.3f}')
plt.plot(np.log(y_20/R_20),np.log(U_20),'ko')
plt.plot(np.log(y_20/R_20),a_20.intercept + (a_20.slope)*np.log(y_20/R_20),label='Fitted Line for $Re_{\mathcal{T}}=20k$')
#FOR 30K
a_30=stats.linregress(np.log(y_30/R_30),np.log(U_30))
print('For Re_Tau=30k:',f'R-squared = {a_30.rvalue**2:.3f},',f'Umax = {np.exp(a_30.intercept):.3f},',f'n = {1/a_30.slope:.3f}')
plt.plot(np.log(y_30/R_30),np.log(U_30),'ko')
plt.plot(np.log(y_30/R_30),a_30.intercept + (a_30.slope)*np.log(y_30/R_30),label='Fitted Line for $Re_{\mathcal{T}}=30k$')
#FOR 40K
a_40=stats.linregress(np.log(y_40/R_40),np.log(U_40))
print('For Re_Tau=40k:',f'R-squared = {a_40.rvalue**2:.3f},',f'Umax = {np.exp(a_40.intercept):.3f},',f'n = {1/a_40.slope:.3f}')
plt.plot(np.log(y_40/R_40),np.log(U_40),'ko')
plt.plot(np.log(y_40/R_40),a_40.intercept + (a_40.slope)*np.log(y_40/R_40),label='Fitted Line for $Re_{\mathcal{T}}=40k$')
plt.title('Linear Regression for various values of $Re_{\mathcal{T}}$')
plt.xlabel("$log(\\frac{y}{R})$")
plt.ylabel("$log(u(y))$")
plt.legend(loc='best',fontsize='x-small')
plt.show()
For Re_Tau=5k: R-squared = 0.941, Umax = 5.868, n = 5.156
For Re_Tau=10k: R-squared = 0.945, Umax = 10.664, n = 5.994
For Re_Tau=20k: R-squared = 0.959, Umax = 24.823, n = 6.490
For Re_Tau=30k: R-squared = 0.958, Umax = 37.627, n = 6.603
For Re_Tau=40k: R-squared = 0.980, Umax = 45.931, n = 7.451
4.5 - Interpolation
from scipy.interpolate import interp1d
Re_Tau = [Re_tau_5,Re_tau_10,Re_tau_20,Re_tau_30,Re_tau_40]
n = [1/a_5.slope,1/a_10.slope,1/a_20.slope,1/a_30.slope,1/a_40.slope]
Umax = [np.exp(a_5.intercept),np.exp(a_10.intercept),np.exp(a_20.intercept),np.exp(a_30.intercept),np.exp(a_40.intercept)]
Re_Tau_35 = 35061
f_n = interp1d(Re_Tau,n)
f_Umax = interp1d(Re_Tau,Umax)
n_35=f_n(Re_Tau_35)
Umax_35=f_Umax(Re_Tau_35)
print(f'For Re_Tau_35k, n = {n_35:.3f} and Umax = {Umax_35:.3f}')
Re_Tau = np.insert(Re_Tau,-1,Re_Tau_35)
n = np.insert(n,-1,n_35)
For Re_Tau_35k, n = 6.947 and Umax = 40.993
plt.plot(Re_Tau,n,'o',label='Coefficients $n$ for each $Re_{\mathcal{T}}$')
a_n=stats.linregress(np.log(Re_Tau),n)
plt.plot(Re_Tau,a_n.intercept + (a_n.slope)*np.log(Re_Tau),'r-', label='Fitted Curve')
print(f'R-squared = {a_n.rvalue**2:.3f}')
plt.title('Coefficients $n$ for all $Re_{\mathcal{T}}$ values')
plt.xscale('log')
plt.xlabel('$Re_{\mathcal{T}}$')
plt.ylabel('$n$')
plt.legend(loc='upper left')
plt.grid()
plt.show()
R-squared = 0.921
y = np.linspace(0,0.450,num=1000)
U_35_prdct = np.zeros(1000)
for i in range(1000):
U_35_prdct[i] = Umax_35*(y[i]/0.45)**(1/n_35)
plt.plot(U_35_prdct,0.45-y,'r-')
plt.title("Projected Velocity Curve for $Re_{\mathcal{T}}=35,061$")
plt.xlabel("Velocity U [m/s]")
plt.ylabel("Distance from the pipe centerline y [m]")
plt.grid()
4.6 - Validation
First let's compare the experimental and predicted velocity profiles on the same plot (Figure 4.6.1):
#Experimental data for Re_Tau=35k
file_35k = '/work/Files_for_Flow/Data/Retau_35k_basic_stats.txt'
Re_tau_35 = np.loadtxt(file_35k,delimiter=' ', usecols=(3), max_rows=1)
R_35 = np.loadtxt(file_35k,delimiter=' ', usecols=(7), skiprows=3, max_rows=1)
y_35,U_35=np.loadtxt(file_35k,delimiter=' ',usecols=(0,2),skiprows=6,unpack=True)
plt.plot(U_35,R_35-y_35,'b.-',label='Experimental Data')
#Predicted data
plt.plot(U_35_prdct,R_35-y,'r-',label='Predicted data')
plt.title("Experimental and Predicted Velocity Curves for $Re_{\mathcal{T}}=35,061$")
plt.xlabel("Velocity U [m/s]")
plt.ylabel("Distance from the pipe centerline y [m]")
plt.legend(loc='lower left')
plt.grid()
U_35_prdct2 = np.zeros(len(y_35))
for i in range(len(y_35)):
U_35_prdct2[i] = Umax_35*(y_35[i]/R_35)**(1/n_35)
plt.plot(U_35_prdct2,R_35-y_35,'r.-',label='Predicted Data')
plt.plot(U_35,R_35-y_35,'b.-',label='Experimental Data')
plt.fill_betweenx(R_35-y_35,U_35_prdct2,U_35,color='yellow',label='Error')
plt.title("Comparison between Experimental and Predicted curves \nfor $Re_{\mathcal{T}}=35,061$")
plt.xlabel("Velocity U [m/s]")
plt.ylabel("Distance from the pipe centerline y [m]")
plt.grid()
plt.legend(loc='lower left')
#Average Error
err1,err2 = 0,0
for i in range(len(y_35)):
err1 += abs(U_35[i] - U_35_prdct2[i])
err2 += 100*(abs(U_35[i] - U_35_prdct2[i]))/U_35_prdct2[i]
err_avg1 = err1 / len(y_35)
err_avg2 = err2 / len(y_35)
print(f'The average absolute error is {err_avg1:.3f} m/s, and the average relative error is {err_avg2:.3f} %.')
The average absolute error is 1.295 m/s, and the average relative error is 5.550 %.
5.0 - Results
6.0 - Conclusion
7.0 - Bibliography
8.0 - Appendix
8.1 - Figures
8.2 - Full Code
##############################################
DATA PARSING
##############################################
import numpy as np
import pylab as plt
#IMPORTING DATA
file_5k = '/work/Files_for_Flow/Data/Retau_5k_basic_stats.txt'
file_10k = '/work/Files_for_Flow/Data/Retau_10k_basic_stats.txt'
file_20k = '/work/Files_for_Flow/Data/Retau_20k_basic_stats.txt'
file_30k = '/work/Files_for_Flow/Data/Retau_30k_basic_stats.txt'
file_40k = '/work/Files_for_Flow/Data/Retau_40k_basic_stats.txt'
#PARSING THE DATA
#FOR Re_Tau = 5K
Re_tau_5 = np.loadtxt(file_5k,delimiter=' ', usecols=(4), max_rows=1)
R_5 = np.loadtxt(file_5k,delimiter=' ', usecols=(7), skiprows=3, max_rows=1)
y_5,U_5=np.loadtxt(file_5k,delimiter=' ',usecols=(0,3),skiprows=6,unpack=True)
plt.plot(U_5,R_5-y_5,'.-',label='For Re_Tau = 5k')
#FOR Re_Tau = 10K
Re_tau_10 = np.loadtxt(file_10k,delimiter=' ', usecols=(4), max_rows=1)
R_10 = np.loadtxt(file_10k,delimiter=' ', usecols=(7), skiprows=3, max_rows=1)
y_10,U_10=np.loadtxt(file_10k,delimiter=' ',usecols=(0,3),skiprows=6,unpack=True)
plt.plot(U_10,R_10-y_10,'.-',label='For Re_Tau = 10k')
#FOR Re_Tau = 20K
Re_tau_20 = np.loadtxt(file_20k,delimiter=' ', usecols=(3), max_rows=1)
R_20 = np.loadtxt(file_20k,delimiter=' ',usecols=(7),skiprows=3,max_rows=1)
y_20 = np.loadtxt(file_20k,delimiter=' ',usecols=(0),skiprows=6)
U_20 = np.loadtxt(file_20k,delimiter=' ',usecols=(3),skiprows=6,max_rows=9)
U_20=np.append(U_20,np.loadtxt(file_20k,delimiter=' ',usecols=(2),skiprows=15))
plt.plot(U_20,R_20-y_20,'.-',label='For Re_Tau = 20k')
#FOR Re_Tau = 30K
Re_tau_30 = np.loadtxt(file_30k,delimiter=' ', usecols=(3), max_rows=1)
R_30 = np.loadtxt(file_30k,delimiter=' ',usecols=(7),skiprows=3,max_rows=1)
y_30 = np.loadtxt(file_30k,delimiter=' ',usecols=(0),skiprows=6)
U_30 = np.loadtxt(file_30k,delimiter=' ',usecols=(3),skiprows=6,max_rows=2)
U_30=np.append(U_30,np.loadtxt(file_30k,delimiter=' ',usecols=(2),skiprows=8))
plt.plot(U_30,R_30-y_30,'.-',label='For Re_Tau = 30k')
#FOR Re_Tau = 40K
Re_tau_40 = np.loadtxt(file_40k,delimiter=' ', usecols=(3), max_rows=1)
R_40 = np.loadtxt(file_40k,delimiter=' ', usecols=(7), skiprows=3, max_rows=1)
y_40,U_40=np.loadtxt(file_40k,delimiter=' ',usecols=(0,2),skiprows=6,unpack=True)
plt.plot(U_40,R_40-y_40,'.-',label='For Re_Tau = 40k')
#ADDING FIGURE LABELS
plt.title("Parsed Data")
plt.xlabel("Velocity U [m/s]")
plt.ylabel("Distance from the pipe centerline y [m]")
plt.legend(loc='upper right')
plt.grid()
##############################################
LINEAR REGRESSION
##############################################
from scipy import stats
#FOR 5K
a_5=stats.linregress(np.log(y_5/R_5),np.log(U_5))
print('For Re_Tau=5k: ',f'R-squared = {a_5.rvalue**2:.3f},',f'Umax = {np.exp(a_5.intercept):.3f}, ',f'n = {1/a_5.slope:.3f}')
plt.plot(np.log(y_5/R_5),np.log(U_5),'ko',label='Given Data')
plt.plot(np.log(y_5/R_5),a_5.intercept + (a_5.slope)*np.log(y_5/R_5),label='Fitted Line for $Re_{\mathcal{T}}=5k$')
#FOR 10K
a_10=stats.linregress(np.log(y_10/R_10),np.log(U_10))
print('For Re_Tau=10k:',f'R-squared = {a_10.rvalue**2:.3f},',f'Umax = {np.exp(a_10.intercept):.3f},',f'n = {1/a_10.slope:.3f}')
plt.plot(np.log(y_10/R_10),np.log(U_10),'ko')
plt.plot(np.log(y_10/R_10),a_10.intercept + (a_10.slope)*np.log(y_10/R_10),label='Fitted Line for $Re_{\mathcal{T}}=10k$')
#FOR 20K
a_20=stats.linregress(np.log(y_20/R_20),np.log(U_20))
print('For Re_Tau=20k:',f'R-squared = {a_20.rvalue**2:.3f},',f'Umax = {np.exp(a_20.intercept):.3f},',f'n = {1/a_20.slope:.3f}')
plt.plot(np.log(y_20/R_20),np.log(U_20),'ko')
plt.plot(np.log(y_20/R_20),a_20.intercept + (a_20.slope)*np.log(y_20/R_20),label='Fitted Line for $Re_{\mathcal{T}}=20k$')
#FOR 30K
a_30=stats.linregress(np.log(y_30/R_30),np.log(U_30))
print('For Re_Tau=30k:',f'R-squared = {a_30.rvalue**2:.3f},',f'Umax = {np.exp(a_30.intercept):.3f},',f'n = {1/a_30.slope:.3f}')
plt.plot(np.log(y_30/R_30),np.log(U_30),'ko')
plt.plot(np.log(y_30/R_30),a_30.intercept + (a_30.slope)*np.log(y_30/R_30),label='Fitted Line for $Re_{\mathcal{T}}=30k$')
#FOR 40K
a_40=stats.linregress(np.log(y_40/R_40),np.log(U_40))
print('For Re_Tau=40k:',f'R-squared = {a_40.rvalue**2:.3f},',f'Umax = {np.exp(a_40.intercept):.3f},',f'n = {1/a_40.slope:.3f}')
plt.plot(np.log(y_40/R_40),np.log(U_40),'ko')
plt.plot(np.log(y_40/R_40),a_40.intercept + (a_40.slope)*np.log(y_40/R_40),label='Fitted Line for $Re_{\mathcal{T}}=40k$')
plt.title('Linear Regression for various values of $Re_{\mathcal{T}}$')
plt.xlabel("$log(\\frac{y}{R})$")
plt.ylabel("$log(u(y))$")
plt.legend(loc='best',fontsize='x-small')
plt.show()
##############################################
INTERPOLATION
##############################################
from scipy.interpolate import interp1d
Re_Tau = [Re_tau_5,Re_tau_10,Re_tau_20,Re_tau_30,Re_tau_40]
n = [1/a_5.slope,1/a_10.slope,1/a_20.slope,1/a_30.slope,1/a_40.slope]
Umax = [np.exp(a_5.intercept),np.exp(a_10.intercept),np.exp(a_20.intercept),np.exp(a_30.intercept),np.exp(a_40.intercept)]
Re_Tau_35 = 35061
f_n = interp1d(Re_Tau,n)
f_Umax = interp1d(Re_Tau,Umax)
n_35=f_n(Re_Tau_35)
Umax_35=f_Umax(Re_Tau_35)
print(f'For Re_Tau_35k, n = {n_35:.3f} and Umax = {Umax_35:.3f}')
Re_Tau = np.insert(Re_Tau,-1,Re_Tau_35)
n = np.insert(n,-1,n_35)
plt.plot(Re_Tau,n,'o',label='Coefficients $n$ for each $Re_{\mathcal{T}}$')
a_n=stats.linregress(np.log(Re_Tau),n)
plt.plot(Re_Tau,a_n.intercept + (a_n.slope)*np.log(Re_Tau),'r-', label='Fitted Curve')
print(f'R-squared = {a_n.rvalue**2:.3f}')
plt.title('Coefficients $n$ for all $Re_{\mathcal{T}}$ values')
plt.xscale('log')
plt.xlabel('$Re_{\mathcal{T}}$')
plt.ylabel('$n$')
plt.legend(loc='upper left')
plt.grid()
plt.show()
##############################################
PREDICTION
##############################################
y = np.linspace(0,0.450,num=1000)
U_35_prdct = np.zeros(1000)
for i in range(1000):
U_35_prdct[i] = Umax_35*(y[i]/0.45)**(1/n_35)
plt.plot(U_35_prdct,0.45-y,'r-')
plt.title("Projected Velocity Curve for $Re_{\mathcal{T}}=35,061$")
plt.xlabel("Velocity U [m/s]")
plt.ylabel("Distance from the pipe centerline y [m]")
plt.grid()
##############################################
VALIDATION
##############################################
#Experimental data for Re_Tau=35k
file_35k = '/work/Files_for_Flow/Data/Retau_35k_basic_stats.txt'
Re_tau_35 = np.loadtxt(file_35k,delimiter=' ', usecols=(3), max_rows=1)
R_35 = np.loadtxt(file_35k,delimiter=' ', usecols=(7), skiprows=3, max_rows=1)
y_35,U_35=np.loadtxt(file_35k,delimiter=' ',usecols=(0,2),skiprows=6,unpack=True)
plt.plot(U_35,R_35-y_35,'b.-',label='Experimental Data')
#Predicted data
plt.plot(U_35_prdct,R_35-y,'r-',label='Predicted data')
plt.title("Experimental and Predicted Velocity Curves for $Re_{\mathcal{T}}=35,061$")
plt.xlabel("Velocity U [m/s]")
plt.ylabel("Distance from the pipe centerline y [m]")
plt.legend(loc='lower left')
plt.grid()
U_35_prdct2 = np.zeros(len(y_35))
for i in range(len(y_35)):
U_35_prdct2[i] = Umax_35*(y_35[i]/R_35)**(1/n_35)
plt.plot(U_35_prdct2,R_35-y_35,'r.-',label='Predicted Data')
plt.plot(U_35,R_35-y_35,'b.-',label='Experimental Data')
plt.fill_betweenx(R_35-y_35,U_35_prdct2,U_35,color='yellow',label='Error')
plt.title("Comparison between Experimental and Predicted curves \nfor $Re_{\mathcal{T}}=35,061$")
plt.xlabel("Velocity U [m/s]")
plt.ylabel("Distance from the pipe centerline y [m]")
plt.grid()
plt.legend(loc='lower left')
#Average Error
err1,err2 = 0,0
for i in range(len(y_35)):
err1 += abs(U_35[i] - U_35_prdct2[i])
err2 += 100*(abs(U_35[i] - U_35_prdct2[i]))/U_35_prdct2[i]
err_avg1 = err1 / len(y_35)
err_avg2 = err2 / len(y_35)
print(f'The average absolute error is {err_avg1:.3f} m/s, and the average relative error is {err_avg2:.3f} %.')