import numpy as np import matplotlib.pyplot as plt from scipy import interpolate from scipy import stats ################################################## #----------IMPORT DATA FROM TXT FILE-------------# ################################################## Re=[5,10,20,30,40] r_2=[] n=[] U_max=[] Re_exact=[] a=[] b=[] y=[] u=[] for i in range(0,5): file1 = f'/work/Files_for_Flow/Data/Retau_{Re[i]}k_basic_stats.txt' if Re[i]<=10: Re_tau = np.loadtxt(file1,delimiter = ' ', usecols=(4),max_rows=1) y_val,u_val = np.loadtxt(file1, delimiter = ' ',usecols=(0,3),skiprows=6, unpack = True) elif Re[i]==20: Re_tau = np.loadtxt(file1,delimiter = ' ', usecols=(3),max_rows=1) y_val = np.loadtxt(file1, delimiter = ' ',usecols=(0),skiprows=6, unpack = True) u_val1 = np.loadtxt(file1, delimiter = ' ',usecols=(3),skiprows=6, unpack = True, max_rows=9) u_val2 =np.loadtxt(file1, delimiter = ' ',usecols=(2),skiprows=15, unpack = True) u_val = np.concatenate((u_val1, u_val2)) elif Re[i]==30: Re_tau = np.loadtxt(file1,delimiter = ' ', usecols=(3),max_rows=1) y_val = np.loadtxt(file1, delimiter = ' ',usecols=(0),skiprows=6, unpack = True) u_val1 = np.loadtxt(file1, delimiter = ' ',usecols=(3),skiprows=6, unpack = True, max_rows=2) u_val2=np.loadtxt(file1, delimiter = ' ',usecols=(2),skiprows=8, unpack = True) u_val = np.concatenate((u_val1, u_val2)) else: Re_tau = np.loadtxt(file1,delimiter = ' ', usecols=(3),max_rows=1) y_val, u_val = np.loadtxt(file1, delimiter = ' ',usecols=(0,2),skiprows=6, unpack = True) R = np.loadtxt(file1,delimiter = ' ', skiprows=3, usecols=(7),max_rows=1) y.append(y_val) u.append(u_val)

######################################################################## #----------LINEAR REGRESSION OF POWER LAW VELOCITY PROFILE-------------# ########################################################################

import numpy as np import matplotlib.pyplot as plt from scipy import interpolate from scipy import stats ################################################## #----------IMPORT DATA FROM TXT FILE-------------# ################################################## Re=[5,10,20,30,40] r_2=[] n=[] U_max=[] Re_exact=[] a=[] b=[] y=[] u=[] for i in range(0,5): file1 = f'/work/Files_for_Flow/Data/Retau_{Re[i]}k_basic_stats.txt' if Re[i]<=10: Re_tau = np.loadtxt(file1,delimiter = ' ', usecols=(4),max_rows=1) y_val,u_val = np.loadtxt(file1, delimiter = ' ',usecols=(0,3),skiprows=6, unpack = True) elif Re[i]==20: Re_tau = np.loadtxt(file1,delimiter = ' ', usecols=(3),max_rows=1) y_val = np.loadtxt(file1, delimiter = ' ',usecols=(0),skiprows=6, unpack = True) u_val1 = np.loadtxt(file1, delimiter = ' ',usecols=(3),skiprows=6, unpack = True, max_rows=9) u_val2 =np.loadtxt(file1, delimiter = ' ',usecols=(2),skiprows=15, unpack = True) u_val = np.concatenate((u_val1, u_val2)) elif Re[i]==30: Re_tau = np.loadtxt(file1,delimiter = ' ', usecols=(3),max_rows=1) y_val = np.loadtxt(file1, delimiter = ' ',usecols=(0),skiprows=6, unpack = True) u_val1 = np.loadtxt(file1, delimiter = ' ',usecols=(3),skiprows=6, unpack = True, max_rows=2) u_val2=np.loadtxt(file1, delimiter = ' ',usecols=(2),skiprows=8, unpack = True) u_val = np.concatenate((u_val1, u_val2)) else: Re_tau = np.loadtxt(file1,delimiter = ' ', usecols=(3),max_rows=1) y_val, u_val = np.loadtxt(file1, delimiter = ' ',usecols=(0,2),skiprows=6, unpack = True) R = np.loadtxt(file1,delimiter = ' ', skiprows=3, usecols=(7),max_rows=1) y.append(y_val) u.append(u_val) ######################################################################## #----------LINEAR REGRESSION OF POWER LAW VELOCITY PROFILE-------------# ######################################################################## reg = stats.linregress(np.log(y[i]/R), np.log(u[i])) a.append(reg.slope) b.append(reg.intercept) r_2.append(reg.rvalue) n.append(1/a[i]) U_max.append(np.exp(b[i])) Re_exact.append(Re_tau) i=0 #Choose for which Re we want to plot x_reg=np.linspace(min(np.log(y[i]/R)),max(np.log(y[i]/R)),len(y[i])) ############################################################ #----------INTERPOLATION OF n(Re) and Umax(Re)-------------# ############################################################ reg_n = stats.linregress(np.log(Re_exact),n) a_n=(reg_n.slope) b_n=(reg_n.intercept) x_n=np.linspace(min(np.log(Re_exact)),max(np.log(Re_exact)),30) f_U_max=interpolate.interp1d(Re_exact, U_max, kind='quadratic') x_U_max=np.linspace(min(Re_exact),max(Re_exact),len(Re_exact)) ############################################### #----------VALIDATION WITH RE=35k-------------# ############################################### n_35k=a_n*np.log(35061)+b_n U_max_35k=f_U_max(35061) file_35k = '/work/Files_for_Flow/Data/Retau_35k_basic_stats.txt' y_35k_exact,u_35k_exact = np.loadtxt(file_35k, delimiter = ' ',usecols=(0,2), skiprows=6, unpack = True) u_35k = U_max_35k*(y_35k_exact/R)**(1/n_35k) Err_prct=(sum(abs(u_35k-u_35k_exact))/sum((u_35k_exact)))*100 print("Error is:",round(Err_prct,1),"%") print("n(Re=35k)=",n_35k) print("Umax(Re=35k)=",U_max_35k,'m/s') ################################################## #----------PLOT LINEAR REGRESSION----------------# ################################################## plt.figure(f"Linear regression in log-log scale for Re={Re[i]}k") plt.title(f"Linear regression in log-log scale for Re={Re[i]}k") plt.xlabel('ln(y/R)[m]')#Set a label to horizontal axis plt.ylabel('ln(u)[m/s]')#Set a label to vertical axis plt.grid() plt.scatter(np.log(y[i]/R),np.log(u[i]),label=f"Points for Re={Re[i]}k",s=10,c="r"); plt.plot(x_reg,a[i]*x_reg+b[i],label="y=ax+b",linewidth=2); plt.legend() ######################################################################################## #----------COMPARAISON BTW VELOCITY PROFILE EXPERIMENTAL AND THEORETICAL--------------# ######################################################################################## plt.figure() plt.title(f"Comparaison of theoretical and experimental velocity profile for Re={Re[i]}k") plt.xlabel('u/Umax[-]')#Set a label to horizontal axis plt.ylabel('(1-y)/R[-]')#Set a label to vertical axis plt.grid() x=np.linspace(0,1,len(u[i])) plt.scatter(u[i]/max(u[i]),1-y[i]/R,label=f"Experimental velocity profile for Re={Re}k", s=10,c="r"); plt.plot(x,1-x**n[i],label=f"Theoretical velocity profile for Re={Re[i]}k",linewidth=2) plt.legend() ##################################################### #--------EVOLUTION OF n IN FUNCTION OF Re-----------# ##################################################### plt.figure("$n$ in function of $Re$") plt.title("$n$ in function of $Re$") plt.xlabel('$Re$[-]')#Set a label to horizontal axis plt.ylabel('$n$[-]')#Set a label to vertical axis plt.grid() plt.scatter(Re_exact,n, label="$n$ in function of $Re$",s=10,c="r") plt.semilogx(np.exp(x_n),a_n*x_n+b_n, label="Linear regression of $n(ln(Re))$", color="orange",linewidth=1); plt.text(12000,5.7, f"With $Re=${Re}k",style = 'italic') plt.legend() ######################################################## #--------EVOLUTION OF Umax IN FUNCTION OF Re-----------# ######################################################## plt.figure("$U_{max}$ in function of $Re$") plt.title("$U_{max}$ in function of $Re$") plt.xlabel('$Re$[-]')#Set a label to horizontal axis plt.ylabel('$U_{max}$ [m/s]')#Set a label to vertical axis plt.grid() plt.scatter(Re_exact,U_max, label="$U_{max}$ in function of $Re$",s=10,c="r") plt.plot(x_U_max,f_U_max(x_U_max), label="Interpolation of $U_{max}(Re)$", color="orange",linewidth=1); plt.text(25000,22, f"With $Re=${Re}k",style = 'italic') plt.legend() ################################################################################################## #----------COMPARAISON BTW VELOCITY PROFILE EXPERIMENTAL AND THEORETICAL FOR Re=35k--------------# ################################################################################################## plt.figure("Comparaison of theoretical and experimental velocity profile for Re=35k") plt.title("Comparaison of theoretical and experimental velocity profile for Re=35k") plt.xlabel('u/Umax[-]')#Set a label to horizontal axis plt.ylabel('(1-y)/R[-]')#Set a label to vertical axis plt.grid() x=np.linspace(0,1,len(u_35k_exact)) plt.scatter(u_35k_exact/max(u_35k_exact),1-y_35k_exact/R, label="Experimental velocity profile for Re=35k",s=10,c="r"); plt.plot(x,1-x**n_35k,label="Theoretical velocity profile for Re=35k",linewidth=2) plt.legend()

import numpy as np import matplotlib.pyplot as plt from scipy import interpolate from scipy import stats ################################################## #----------IMPORT DATA FROM TXT FILE-------------# ################################################## Re=[5,10,20,30,40] r_2_vect=[] n_vect=[] U_max=[] Re_exact=[] a=[] b=[] a_vect=[] b_vect=[] y=[] u=[] for i in range(0,5): file1 = f'/work/Files_for_Flow/Data/Retau_{Re[i]}k_basic_stats.txt' if Re[i]<=10: Re_tau = np.loadtxt(file1,delimiter = ' ', usecols=(4),max_rows=1) y_val,u_val = np.loadtxt(file1, delimiter = ' ',usecols=(0,3),skiprows=6, unpack = True) elif Re[i]==20: Re_tau = np.loadtxt(file1,delimiter = ' ', usecols=(3),max_rows=1) y_val = np.loadtxt(file1, delimiter = ' ',usecols=(0),skiprows=6, unpack = True) u_val1 = np.loadtxt(file1, delimiter = ' ',usecols=(3),skiprows=6, unpack = True, max_rows=9) u_val2 =np.loadtxt(file1, delimiter = ' ',usecols=(2),skiprows=15, unpack = True) u_val = np.concatenate((u_val1, u_val2)) elif Re[i]==30: Re_tau = np.loadtxt(file1,delimiter = ' ', usecols=(3),max_rows=1) y_val = np.loadtxt(file1, delimiter = ' ',usecols=(0),skiprows=6, unpack = True) u_val1 = np.loadtxt(file1, delimiter = ' ',usecols=(3),skiprows=6, unpack = True, max_rows=2) u_val2=np.loadtxt(file1, delimiter = ' ',usecols=(2),skiprows=8, unpack = True) u_val = np.concatenate((u_val1, u_val2)) else: Re_tau = np.loadtxt(file1,delimiter = ' ', usecols=(3),max_rows=1) y_val, u_val = np.loadtxt(file1, delimiter = ' ',usecols=(0,2),skiprows=6, unpack = True) R = np.loadtxt(file1,delimiter = ' ', skiprows=3, usecols=(7),max_rows=1) y.append(y_val) u.append(u_val) ######################################################################## #----------LINEAR REGRESSION OF POWER LAW VELOCITY PROFILE-------------# ######################################################################## n=1 diff=[2,1] r_2=0 while diff[-2]>diff[-1] and n<20: reg = stats.linregress(y[i],u[i]**n) a.append(reg.slope) b.append(reg.intercept) r_2_vect.append(reg.rvalue) n+=0.1 diff.append(1-r_2_vect[-1]) a_vect.append(a[-1]) b_vect.append(b[-1]) n_vect.append(round(n, 2)-0.1) U_max.append((a_vect[i]*R)**(1/n_vect[i])) Re_exact.append(Re_tau) i=0 #Choose for which Re we want to plot x_reg=np.linspace(0,R,len(y[i])) ############################################################ #----------INTERPOLATION OF n(Re) AND Umax(Re)-------------# ############################################################ reg_n = stats.linregress(np.log(Re_exact),n_vect) a_n=(reg_n.slope) b_n=(reg_n.intercept) x_n=np.linspace(min(np.log(Re_exact)),max(np.log(Re_exact)),len(Re_exact)) f_U_max=interpolate.interp1d(Re_exact, U_max, kind='quadratic') x_U_max=np.linspace(min(Re_exact),max(Re_exact),len(Re_exact)) ############################################### #----------VALIDATION WITH RE=35k-------------# ############################################### n_35k=a_n*np.log(35061)+b_n U_max_35k=f_U_max(35061) file_35k = '/work/Files_for_Flow/Data/Retau_35k_basic_stats.txt' y_35k_exact,u_35k_exact = np.loadtxt(file_35k, delimiter = ' ',usecols=(0,2), skiprows=6, unpack = True) u_35k = U_max_35k*(y_35k_exact/R)**(1/n_35k) Err_prct=(sum(abs(u_35k-u_35k_exact))/sum((u_35k_exact)))*100 print("Error is:",round(Err_prct,1),"%") print("n(Re=35k)=",n_35k) print("Umax(Re=35k)=",U_max_35k,'m/s') ################################################## #----------PLOT LINEAR REGRESSION----------------# ################################################## plt.figure(f"Linear regression for Re={Re[i]}k") plt.title(f"Linear regression for Re={Re[i]}k") plt.xlabel('y [m]')#Set a label to horizontal axis plt.ylabel('u^n [m/s]')#Set a label to vertical axis plt.grid() plt.scatter(y[i],u[i]**n_vect[i],label=f"Points for Re={Re[i]}k", s=10,c="r"); plt.plot(x_reg,a_vect[i]*x_reg+b_vect[i],label="y=ax+b",linewidth=2); plt.legend() ######################################################################################## #----------COMPARAISON BTW VELOCITY PROFILE EXPERIMENTAL AND THEORETICAL--------------# ######################################################################################## plt.figure() plt.title(f"Comparaison of theoretical and experimental velocity profile for Re={Re[i]}k") plt.xlabel('u/Umax[-]')#Set a label to horizontal axis plt.ylabel('(1-y)/R[-]')#Set a label to vertical axis plt.grid() x=np.linspace(0,1,len(u[i])) plt.scatter(u[i]/max(u[i]),1-y[i]/R,label=f"Experimental velocity profile for Re={Re}k", s=10,c="r"); plt.plot(x,1-x**n_vect[i],label=f"Theoretical velocity profile for Re={Re[i]}k",linewidth=2) plt.legend() ##################################################### #--------EVOLUTION OF n IN FUNCTION OF Re-----------# ##################################################### plt.figure("$n$ in function of $Re$") plt.title("$n$ in function of $Re$") plt.xlabel('$Re$[-]')#Set a label to horizontal axis plt.ylabel('$n$[-]')#Set a label to vertical axis plt.grid() plt.scatter(Re_exact,n_vect, label="$n$ in function of $Re$",s=10,c="r") plt.semilogx(np.exp(x_n),a_n*x_n+b_n, label="Linear regression of $n(ln(Re))$", color="orange",linewidth=1); plt.text(12000,5.7, f"With $Re=${Re}k",style = 'italic') plt.legend() ######################################################## #--------EVOLUTION OF Umax IN FUNCTION OF Re-----------# ######################################################## plt.figure("$U_{max}$ in function of $Re$") plt.title("$U_{max}$ in function of $Re$") plt.xlabel('$Re$[-]')#Set a label to horizontal axis plt.ylabel('$U_{max}$ [m/s]')#Set a label to vertical axis plt.grid() plt.scatter(Re_exact,U_max, label="$U_{max}$ in function of $Re$", s=10,c="r") plt.plot(x_U_max,f_U_max(x_U_max), label="Interpolation of $U_{max}(Re)$", color="orange",linewidth=1); plt.text(25000,22, f"With $Re=${Re}k",style = 'italic') plt.legend() ################################################################################################## #----------COMPARAISON BTW VELOCITY PROFILE EXPERIMENTAL AND THEORETICAL FOR Re=35k--------------# ################################################################################################## plt.figure() plt.title("Comparaison of theoretical and experimental velocity profile for Re=35k") plt.xlabel('u/Umax[-]')#Set a label to horizontal axis plt.ylabel('(1-y)/R[-]')#Set a label to vertical axis plt.grid() x=np.linspace(0,1,len(u_35k_exact)) plt.scatter(u_35k_exact/max(u_35k_exact),1-y_35k_exact/R, label="Experimental velocity profile for Re=35k",s=10,c="r"); plt.plot(x,1-x**n_35k,label="Theoretical velocity profile for Re=35k",linewidth=2) plt.legend()