D = 2.5e-2 # m
nu = 15.e-6 #m^2/s
rho = 1.2 #kg/m^3
mu = nu/rho 

R = D/2

Re = 500.

Ub = Re * nu / D
print("Bulk velocity= %2.2f m/s" %Ub)
flowrate = Ub * np.pi * R**2
print("Flow rate = %1.4e m^3/s" %flowrate)

import numpy as np

n = 30
r = np.linspace(0,R,n)

U = 2 * Ub * (1 - np.power(r,2)/R**2)
import matplotlib.pyplot as plt

plt.plot(r,U,linewidth = 2)
plt.xlabel(r"$r$ (m)", fontdict = fontlabel)
plt.ylabel(r"$U(r)$ (m/s)", fontdict = fontlabel)
plt.xlim(0,R)
plt.show()

"""
Explanation: <h1>Laminar Pipe Flow</h1>

In this first tutorial, you will simulate a laminar pipe flow using the cad file <FONT FACE="courier" style="color:blue">CAD>PipeD0.025L2.5VOF.iges</FONT> which is the <b>Volume of Fluid</b> of the CAD below. The pipe inner diameter is 2.5cm and 250cm long.
<img src="CAD/Pipe.png">
The Junior course of Fluid Mechanics has taught you that:
<ul>
<li> The entrance length for laminar flow is given by
$$
\frac{L_e}{D}=0.05Re=0.05\frac{U_{b}D}{\nu}
$$
</li>
<li> The pressure gradient is constant in the fully developed region ($x > L_e$)</li>
<li> The velocity profile in the fully developed region is governed by the reduced streamwise momentum equation 
$$
0=-\frac{dP}{dx}+\frac{\mu}{r}\frac{d}{dr}r\frac{dU}{dr}
$$
</li>
<li>The solution of this equation is
$$
U(r) = \frac{R^2}{4\mu}\left(-\frac{dP}{dx}\right)\left(1-\frac{r^2}{R^2}\right)=2U_b\left(1-\frac{r^2}{R^2}\right)
$$
</li>
<li>Hence, 
$$
U_b=\frac{R^2}{8\mu}\left(-\frac{dP}{dx}\right)
$$
<li> Finally, the wall shear stress is proportional to the pressure drop:
$$
-\frac{dP}{dx}=\frac{2\tau_w}{R}
$$
</ul>

<p class='alert alert-danger'>
Using <a href="https://www.simscale.com">https://www.simscale.com</a>, Simulate a laminar flow. Verify and Validate (V & V) your simulation.
</p>
End of explanation
"""


tauw = - mu * (U[n-1] - U[n-2])/(r[n-1] - r[n-2])
print("wall shear stress = %1.4f Pa" %tauw)

mdpdx = 8 * mu *Ub / R**2
tauw_e = R * mdpdx / 2
print("Exact wall shear stress = %1.4f Pa" %tauw_e)

import csv, sys
filename = 'Data/LaminarPipe/Coarse/centerline.csv'
with open(filename, newline='') as f:
    reader = csv.reader(f)
    try:
        for row in reader:
            print(row)
    except csv.Error as e:
        sys.exit('file {}, line {}: {}'.format(filename, reader.line_num, e))



data = np.genfromtxt(filename,delimiter=',',skip_header=1)
P = data[:,0]
U = data[:,1]
x = data[:,5]
plt.plot(x,P,linewidth = 2)
plt.xlabel(r"$x$ (m)", fontdict = fontlabel)
plt.ylabel(r"$P$ (Pa)", fontdict = fontlabel)
plt.show()
plt.plot(x,U,linewidth = 2)
plt.xlabel(r"$x$ (m)", fontdict = fontlabel)
plt.ylabel(r"$U_c$ (m/s)", fontdict = fontlabel)
plt.show()

filename = 'Data/LaminarPipe/Coarse/cross-section-x=0.025.csv'
data = np.genfromtxt(filename,delimiter=',',skip_header=1)
U_0_025 = data[:,1]
x_0_025 = data[:,5]
filename = 'Data/LaminarPipe/Coarse/cross-section-x=0.25.csv'
data = np.genfromtxt(filename,delimiter=',',skip_header=1)
U_0_25 = data[:,1]
x_0_25 = data[:,5]
filename = 'Data/LaminarPipe/Coarse/cross-section-x=1.25.csv'
data = np.genfromtxt(filename,delimiter=',',skip_header=1)
U_1_25 = data[:,1]
x_1_25 = data[:,5]
filename = 'Data/LaminarPipe/Coarse/cross-section-x=2.25.csv'
data = np.genfromtxt(filename,delimiter=',',skip_header=1)
U_2_25 = data[:,1]
x_2_25 = data[:,5]
plt.plot(x_0_025,U_0_025,linewidth = 2,label = r"$x=D$")
plt.plot(x_0_25,U_0_25,linewidth = 2,label = r"$x=10D$")
plt.plot(x_1_25,U_1_25,linewidth = 2,label = r"$x=50D$")
plt.plot(x_2_25,U_2_25,linestyle = 'dashed' ,linewidth = 2,label = r"$x=100D$")
plt.xlabel(r"$x$ (m)", fontdict = fontlabel)
plt.ylabel(r"$U$ (m/s)", fontdict = fontlabel)
plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=2, shadow=True, fancybox=True)
plt.show()

"""
Explanation: The wall shear stress is defined:
$$
\tau_w = \mu \frac{dU}{dr}
$$
End of explanation
"""


