{ "cells": [ { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "# MTF053 - Computer Assignment 1 (CA1) - Numerical Analysis of Fully Develooped Channel Flow\n", "\n", "In this exercise, you will study a fully developed channel flow (flow between two parallel plates) numerically. You will start with a laminar flow as that problem can be solved analytically and thus it is possible to make a comparison and get a feeling for the accuracy of the numerical method. In the second part of the assignment, you will analyze a turbulent flow numerically and compare your results to provided measured data. The numerical part will be done using `Python` codes provided below. You will need to make some minor modifications to some of the codes. In the header of each code segment, there is a description of what modifications that should be done if any. \n", "\n", "The assignemnt is divided into tasks some of which includes running numerical simulations. You can run the python codes directly in this `Jupyter NoteBook`. For some of the tasks, you are expected to write down answers to questions. This should be done in the dedicated sections in this notebook titled *\"answer to task x.y\"*.\n", "\n", "More detailed instructions, theory and background are given in the assignment instruction document [MTF053_CA1.pdf](https://courses.g3dflow.com/fluidmech/docs/MTF053_CA1.pdf).\n", "\n", "## Jupyter notebook instructions\n", "\n", "### Edit a notebook cell\n", "\n", "Double click on the cell to enter edit mode. \n", "\n", "### Update a cell\n", "\n", "With the cell marked, klick on the play symbol in at the top of the notebook\n", "\n", "### Add a new cell\n", "\n", "Click on the cell above or below which you want a new cell and click on one of the symbols at the top right corner of the marked cell. Set cell type (Markdown or Code) for the new cell by clicking the cell and select *Code* or *Markdown* in the dropdown list in the top menu of the notebook. You'll find a good documentation of Markdown [here](https://www.markdownguide.org/basic-syntax/).\n", "\n", "### Add an image in a Markdown cell\n", "\n", "If you like to include an image in one of the Markdown text cells (for example if you want to include a photo of a derivation made on paper), do as follows:\n", "\n", "1. upload your image using the *upload files* function\n", "\n", "2. include the image in the notebook by writing: \\!\\[title\\]\\(image.png\\)\n", "\n", "3. typeset the cell by klicking the play icon at the top of the notebook to see the result\n", "\n", "You can also include an image available online, add Markdown code as follows: \\!\\[title\\]\\(url\\)\n", "\n", "If it does not work to include a local image in a Markdown cell, add a new code cell and add the following code\n", "\n", "`Image(filename='filename.jpg')`\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "# Initialization\n", "---\n", "Make sure to execute the code section below before runnig any of the other code sections that follows as this will load the `python` libraries needed." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "#********************************************************\n", "#\n", "# NO NEED FOR MODIFICATIONS IN THIS NOTEBOOK SECTION\n", "#\n", "#********************************************************\n", "\n", "# IMPORT THE PYTHON LIBRARIES NEEDED \n", "\n", "import numpy as np\n", "from scipy.sparse import diags\n", "from scipy.sparse.linalg import spsolve\n", "import matplotlib.pyplot as plt\n", "from IPython.display import Image\n", "\n", "# DEFINITION OF GLOBAL VARIABLES\n", "\n", "verbose = 0 # verbose solver output activated if verbose=1 \n", "FontSize = 20\n", "\n", "# geometry\n", "\n", "h = 0.06 # channel height\n", "\n", "# fluid properties\n", "\n", "mu = 1.8e-5 # fluid viscosity (dynamic viscosity)\n", "rho = 1.2 # fluid density\n", "\n", "# log-law constants\n", "\n", "kappa = 0.41 # von Kármán constant\n", "B = 5.0 # integration constant in the log-law\n", "\n", "\n", "# DEFINITION OF HELP FUNCTIONS\n", "\n", "#========================================================\n", "# mesh generation for laminar channel flow\n", "#========================================================\n", "def generate_mesh_for_laminar_flow_simulation(N,h):\n", " \n", " # allocate y vector\n", " \n", " y = np.zeros(N+2)\n", " \n", " # calculate cell size (dy) assuming equidistant grid\n", "\n", " dy=(h/2.0)/N \n", " \n", " # add wall-normal coordinates for all cells (cell 1 to N) \n", "\n", " y[slice(1,N+1)] = np.linspace(dy/2.0,0.5*h-(dy/2.0),N)\n", "\n", " # add centerline coordinate\n", "\n", " y[-1]=h/2. \n", " \n", " return y\n", "\n", "#========================================================\n", "# get an estimate of the numerical error\n", "# for a laminar channel flow simulation \n", "# (diff from analytic)\n", "#========================================================\n", "def laminar_flow_calculate_error(U,Ua):\n", " \n", " if Ua[-2] == 0.0 :\n", " return -1\n", " \n", " # calcualte relative error (deviation in last cell)\n", " \n", " return np.abs((U[-2]-Ua[-2])/Ua[-2])\n", "\n", "#========================================================\n", "# calculate max velocity (laminar channel flow)\n", "#========================================================\n", "def laminar_flow_get_umax(h,dpdx,mu):\n", " \n", " # calculate the analytic maximum velocity (channel center)\n", " \n", " return -0.125*h**2*dpdx/mu\n", "\n", "#========================================================\n", "# calculate y+ and u+ from y, U, tau_w, rho, and mu\n", "#========================================================\n", "def calc_yplus_uplus(rho,mu,tau_w,y,U):\n", " nu=mu/rho\n", " ustar=np.sqrt(tau_w/rho)\n", " yplus=y*ustar/nu\n", " uplus=U/ustar\n", " return yplus,uplus\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 1. Laminar Channel Flow" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Task 1.1\n", "---\n", "\n", "Show that the Navier-Stokes equations\n", "\n", "$$\\dfrac{\\partial u}{\\partial x}+\\dfrac{\\partial v}{\\partial y}=0$$\n", "\n", "$$\\rho\\left(u\\dfrac{\\partial u}{\\partial x}+v\\dfrac{\\partial u}{\\partial y}\\right)=-\\dfrac{\\partial p}{\\partial x}+\\rho g_x+\\dfrac{\\partial}{\\partial x}\\left(\\mu\\dfrac{\\partial u}{\\partial x}\\right)+\\dfrac{\\partial}{\\partial y}\\left(\\mu\\dfrac{\\partial u}{\\partial y}\\right)$$\n", "\n", "$$\\rho\\left(u\\dfrac{\\partial v}{\\partial x}+v\\dfrac{\\partial v}{\\partial y}\\right)=-\\dfrac{\\partial p}{\\partial y}+\\rho g_y+\\dfrac{\\partial}{\\partial x}\\left(\\mu\\dfrac{\\partial v}{\\partial x}\\right)+\\dfrac{\\partial}{\\partial y}\\left(\\mu\\dfrac{\\partial v}{\\partial y}\\right)$$\n", "\n", "reduces to\n", "\n", "$$\\dfrac{\\partial}{\\partial y}\\left(\\mu\\dfrac{\\partial u}{\\partial y}\\right) = \\dfrac{\\partial p}{\\partial x}$$\n", "\n", "for fully developed laminar channel flow in two dimensions. Examining the equation above, it is obvious that for our specified application, it is possible to reduce the number of terms in the Navier-Stokes equations significantly, which means that you will remove terms in the equations as part of the derivations. Please note, however, that all such modifications of the equations must be justified.\n", "\n", "\n", "**Hint:** the lecture notes for chapter 4 ([MTF053_C04.pdf](https://courses.g3dflow.com/fluidmech/notes/MTF053_C04.pdf)) may be of use here.\n", "\n", "### Answer to task 1.1\n", "\n", "**write down your answer here**\n", "\n", "**Note!** derivations made uisng pen and paper or a tablet can be included as photos/images. Upload the photo or image and include as described at the top of this document.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Task 1.2\n", "---\n", "\n", "1. Solve the simplified equation from task 1.1 analytically for the above-listed boundary conditions and flow properties\n", "2. Calculate the Reynolds number\n", "\n", "Remember that we are investigating the flow in a two-dimensional channel, i.e., it is **not** a pipe flow. The correct hydraulic diameter for such a *non-circular pipe flow* can be found in [MTF053_Formulas-Tables-and-Graphs.pdf](https://courses.g3dflow.com/fluidmech/docs/MTF053_Formulas-Tables-and-Graphs.pdf). Also, the Reynolds number should be calculated using the average velocity, which can be obtained by integrating over the channel\n", "\n", "$$V=\\dfrac{1}{h}\\int_0^hu(y)dy$$\n", "\n", "\n", "### Answer to task 1.2\n", "\n", "**write down your answers here**\n", "\n", "**Note!** derivations made uisng pen and paper or a tablet can be included as photos/images. Upload the photo or image and include as described at the top of this document." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Task 1.3\n", "---\n", "\n", "Write the equation system that we get by evaluating \n", "\n", "$$u_{i+1}-2u_{i}+u_{i-1}=\\dfrac{\\partial p}{\\partial x}\\dfrac{\\Delta y^2}{\\mu}$$\n", "\n", "in all of the $N$ cells on matrix form, i.e., write the equation above as $A\\mathbf{u}=\\mathbf{b}$, where $\\mathbf{u}$ is a vector with $N$ elements. What does the matrix $A$ and the vector $\\mathbf{b}$ look like?\n", "\n", "To solve this rather simple problem, all we have to do now is to invert the matrix $A$ and then the velocity values are calculated as $\\mathbf{u}=A^{-1}\\mathbf{b}$\n", "\n", "### Answers to task 1.3\n", "\n", "**write down your answers here**\n", "\n", "**Note!** derivations made uisng pen and paper or a tablet can be included as photos/images. Upload the photo or image and include as described at the top of this document.\n" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "---\n", "## Task 1.4\n", "---\n", "\n", "Below you'll find a `Python` code that solves the equation system derived in Task 1.3 using the `Python` functions `scipy.sparse.diags` and `scipy.sparse.linalg.spsolve`.\n", "\n", "Before running the solver, you must make some minor changes in the codes:\n", "\n", "1. in the function `laminar_channel_flow_solver`, specify the boundary entries in the coefficient matrix $A$\n", "2. in the function `get_analytic_laminar_velocity_profile`, implement calculation of analytical solution for laminar channel flow\n", "3. in the main code for task 1.4, specify number of cells\n", "\n", "With these changes made it is now time to run the code. To execute the code in a code cell, click on the code cell to make it active and then click on the play button at the top of the notebook. Make sure to execute the code cells in sequence. The functions `laminar_channel_flow_solver` and `get_analytic_laminar_velocity_profile` must be executed before running the task 1.4 main code.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Laminar flow solver\n", "Below you find a template for a laminar flow solver written in `python`. Before running the solver you will have to make a few minor changes (indicated in the code)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "\n", "# FUNCTION: laminar_channel_flow_solver(y,mu,rho,dpdx)\n", "\n", "#********************************************************\n", "#\n", "# CODE MODIFICATION NEEDED: \n", "#\n", "# 1. IMPLEMENT BOUNDARY CONDITIONS\n", "#\n", "#********************************************************\n", "\n", "#========================================================\n", "# laminar channel flow solver\n", "# numerical calculation of a laminar boundary layer \n", "# velocity profile using a finite-volume approach\n", "#========================================================\n", "#\n", "# 0 1 2 N-1 N N+1\n", "# ------------------- - - - -------------------\n", "# | | | | | | | |\n", "# | o | o | o | | o | o | o |\n", "# | G1 | | | | | | G2 |\n", "# ------------------- - - - -------------------\n", "# ^ ^ \n", "# \t y=0. y=h/2 \n", "# u=0. du/dy=0\n", "#\n", "# The computational domain contains N cells ranging \n", "# from cell 1 to cell N. G1 and G2 are ghost cells \n", "# (cells outside of the computational domain used for\n", "# boundary condition implementation).\n", "#========================================================\n", "\n", "def laminar_channel_flow_solver(y,mu,rho,dpdx):\n", " \n", " N = len(y) - 2\n", " \n", " dy=2.0*(y[1]-y[0])\n", " \n", " # allocate U vector\n", " \n", " U = np.zeros(N+2)\n", "\n", " # generate matrix diagonal vector and off-diagonal vectors\n", "\n", " diagonal = -2.0*np.ones(N)\n", " off_diagonal = np.ones(N-1)\n", "\n", " # set boundary conditions (update diagonal elements)\n", "\n", " diagonal[ 0] = 0 # <= ADD LEFT-END BOUNDARY CONDITION HERE\n", " diagonal[-1] = 0 # <= ADD RIGHT-END BOUNDARY CONDITION HERE\n", " \n", " if diagonal[0] == 0 and diagonal[-1] == 0 :\n", " return U\n", "\n", " # generate sparse matrix (A) (note! the sparse library functions \n", " # are used in order to make the code more efficient)\n", "\n", " A = diags([diagonal,off_diagonal,off_diagonal],[0, -1, 1],format=\"csr\")\n", "\n", " # generate right-hand-side vector (b)\n", "\n", " b = (dpdx*dy**2/mu)*np.ones(N)\n", "\n", " # solve linear equation system AU=b => U=inv(A)b\n", "\n", " U[slice(1,N+1)] = spsolve(A,b)\n", "\n", " # update boundary velocity\n", "\n", " U[-1] = U[-2]\n", " \n", " return U\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Running the laminar flow solver\n", "Here you define variables necessary to run the laminar flow solver. Make sure to run the code above first since the function defined there will be used in the code below." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "\n", "# FUNCTION: get_analytic_laminar_velocity_profile(y,h,dpdx,mu)\n", "\n", "#********************************************************\n", "#\n", "# CODE MODIFICATION NEEDED: \n", "#\n", "# 1. CALCULATE ANALYTIC SOLUTION (LAMINAR CHANNEL FLOW)\n", "#\n", "#********************************************************\n", "\n", "#========================================================\n", "# function that calculates the analytic solution for \n", "# laminar channel flow\n", "#========================================================\n", "def get_analytic_laminar_velocity_profile(y,h,dpdx,mu):\n", " \n", " # calculate the velocity for all y-coordinates (analytic solution)\n", " \n", " Ua = y*0.0 # <= CALCULATE ANALYTIC SOLUTION HERE U(y)\n", " \n", " return Ua\n" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "\n", "# TASK 1.4 main code\n", "\n", "#********************************************************\n", "#\n", "# CODE MODIFICATION NEEDED: \n", "#\n", "# 1. SPECIFY CELL COUNT\n", "#\n", "#********************************************************\n", "\n", "#========================================================\n", "# running the laminar channel flow solver\n", "#========================================================\n", "\n", "dpdx = -0.01 # pressure gradient\n", "N = 0 # number of cells <= SET CELL COUNT HERE\n", "\n", "if N > 0 :\n", " \n", " # generate mesh\n", "\n", " y_lam = generate_mesh_for_laminar_flow_simulation(N,h)\n", "\n", " # run laminar channel flow solver \n", "\n", " U_lam = laminar_channel_flow_solver(y_lam,mu,rho,dpdx)\n", "\n", " # calculate analytic solution (laminar flow)\n", "\n", " U_analytic = get_analytic_laminar_velocity_profile(y_lam,h,dpdx,mu) \n", "\n", " # plot laminar flow velocity profile\n", "\n", " fig = plt.figure(num=1, figsize=(15, 10), dpi=80, facecolor='w', edgecolor='k')\n", " ax = fig.add_subplot(111)\n", " ax.plot(U_lam,y_lam,linewidth=2,label='numeric' )\n", " ax.plot(U_analytic,y_lam,linewidth=2,label='analytic')\n", " ax.set_xlabel('axial velocity (u [m/s])',fontsize=FontSize)\n", " ax.set_ylabel('wall-normal coordinate (y [m])',fontsize=FontSize)\n", " ax.set_title('Velocity profile - laminar channel flow',fontsize=FontSize)\n", " ax.legend(fontsize=FontSize)\n", " ax.set_xlim(0,round(laminar_flow_get_umax(h,dpdx,mu)*100)/100)\n", " ax.set_ylim(0,0.5*h)\n", " for label in (ax.get_xticklabels() + ax.get_yticklabels()):\n", " label.set_fontsize(FontSize)\n", " ax.grid()\n", " \n", " plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "---\n", "## Task 1.5\n", "---\n", "\n", "For this task, you'll find a `Python` code that runs the laminar flow solver set up in the previous task and calculates the error ($e$) according to\n", "\n", "$$e=\\left|\\dfrac{u_{numeric}-u_{analytic}}{u_{analytic}}\\right|$$\n", "\n", "for 10 different values of $N$ (number of iterations) spanning from $N=10$ to $N=1000$. Note that the `Python` function `numpy.logspace` is used to distribute the cell counts. After running the laminar flow solver with 10 different cell counts and calculating the error for each case, a `loglog` plot showing the error as a function of cell count is generated.\n", "\n", "1. How does the error decay with increased cell count?\n", "\n", " **write down your answer here**\n", " \n", "2. When calculating the error, why is $u_{analytic}$ not calculated at the center of the duct?\n", "\n", " **write down your answer here**\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "\n", "# TASK 1.5 main code\n", "\n", "#********************************************************\n", "#\n", "# NO CODE MODIFICATIONS NEEDED IN THIS NOTEBOOK SECTION\n", "#\n", "#********************************************************\n", "\n", "# case-specific properties\n", "\n", "dpdx = -0.01 # pressure gradient\n", "\n", "# generate a vector with cell counts spanning from 10 to 1000\n", "\n", "NN = np.int_(np.logspace(1,3,10))\n", "\n", "# allocate vectors\n", "\n", "err = []\n", "y_lam = []\n", "U_lam = []\n", "U_analytic = []\n", "\n", "# run the laminar flow solver for each cell count N in the array NN\n", "\n", "for N in NN:\n", " \n", " # generate mesh\n", "\n", " y_lam = generate_mesh_for_laminar_flow_simulation(N,h)\n", "\n", " # run laminar channel flow solver \n", "\n", " U_lam = laminar_channel_flow_solver(y_lam,mu,rho,dpdx)\n", "\n", " # get analytic solution (laminar flow)\n", "\n", " U_analytic = get_analytic_laminar_velocity_profile(y_lam,h,dpdx,mu)\n", "\n", " # calculate error\n", " \n", " error = laminar_flow_calculate_error(U_lam,U_analytic)\n", "\n", " # append error to the error array \n", " \n", " err.append(error)\n", " \n", " \n", "# plot velocity results\n", "\n", "if err[0] > 0.0 :\n", "\n", " fig = plt.figure(num=1, figsize=(15, 10), dpi=80, facecolor='w', edgecolor='k')\n", "\n", " ax = fig.add_subplot(121)\n", " ax.plot(U_lam ,y_lam,linewidth=2,label='numeric' )\n", " ax.plot(U_analytic,y_lam,linewidth=2,label='analytic')\n", " ax.set_xlabel('axial velocity (u [m/s])',fontsize=FontSize)\n", " ax.set_ylabel('wall-normal coordinate (y [m])',fontsize=FontSize)\n", " ax.set_title('Velocity profile - laminar channel flow',fontsize=FontSize)\n", " ax.legend(fontsize=FontSize)\n", " ax.set_xlim(0,round(laminar_flow_get_umax(h,dpdx,mu)*100)/100)\n", " ax.set_ylim(0,0.5*h)\n", " for label in (ax.get_xticklabels() + ax.get_yticklabels()):\n", " label.set_fontsize(FontSize)\n", " ax.grid()\n", "\n", " ax = fig.add_subplot(122)\n", " ax.loglog(NN,err,'-o',linewidth=2)\n", " ax.set_xlabel('number of cells',fontsize=FontSize)\n", " ax.set_ylabel('relative error [%]',fontsize=FontSize)\n", " ax.set_title('Relative error as a function of cell count',fontsize=FontSize)\n", " ax.set_xlim(NN[0],NN[-1])\n", " ax.set_ylim(err[-1],err[0])\n", " for label in (ax.get_xticklabels() + ax.get_yticklabels()):\n", " label.set_fontsize(FontSize)\n", " ax.grid(which='both')\n", " \n", " plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Task 2.1\n", "---\n", "\n", "Replace $u$, $v$, and $p$ in the Navier-Stokes equations\n", "\n", "$$\\dfrac{\\partial u}{\\partial x}+\\dfrac{\\partial v}{\\partial y}=0$$\n", "\n", "$$\\rho\\left(u\\dfrac{\\partial u}{\\partial x}+v\\dfrac{\\partial u}{\\partial y}\\right)=-\\dfrac{\\partial p}{\\partial x}+\\rho g_x+\\dfrac{\\partial}{\\partial x}\\left(\\mu\\dfrac{\\partial u}{\\partial x}\\right)+\\dfrac{\\partial}{\\partial y}\\left(\\mu\\dfrac{\\partial u}{\\partial y}\\right)$$\n", "\n", "$$\\rho\\left(u\\dfrac{\\partial v}{\\partial x}+v\\dfrac{\\partial v}{\\partial y}\\right)=-\\dfrac{\\partial p}{\\partial y}+\\rho g_y+\\dfrac{\\partial}{\\partial x}\\left(\\mu\\dfrac{\\partial v}{\\partial x}\\right)+\\dfrac{\\partial}{\\partial y}\\left(\\mu\\dfrac{\\partial v}{\\partial y}\\right)$$\n", "\n", "with\n", "\n", "$$u=\\overline{u}+u^\\prime,\\ v=\\overline{v}+v^\\prime,\\ p=\\overline{p}+p^\\prime$$\n", "\n", "and derive the the Reynolds-Averaged Navier-Stokes (RANS) equations. Use the fact that the time average of a fluctuation is zero. Then, use the same arguments as in Task 1.1 to show that the equations can be simplified to the following relation for our specific problem\n", "\n", "\n", "$$0=-\\dfrac{\\partial \\overline{p}}{\\partial x}+\\dfrac{\\partial}{\\partial y}\\left(\\mu\\dfrac{\\partial \\overline{u}}{\\partial y}-\\rho\\overline{u^\\prime v^\\prime}\\right)$$\n", "\n", "\n", "**Hint:** make use of the derivation available in the lecture notes for Chapter 6 ([MTF053_C06.pdf](https://courses.g3dflow.com/fluidmech/notes/MTF053_C06.pdf)).\n", "\n", "### Answer to task 2.1\n", "\n", "**present your derivations here**\n", "\n", "**Note!** derivations made uisng pen and paper or a tablet can be included as photos/images. Upload the photo or image and include as described at the top of this document.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Task 2.2\n", "---\n", "\n", "In the task 2.2 main code below, set the number of cells to $N=1000$. Run the turbulent flow solver for 50, 100, 200, and 300 iterations and monitor the development of the solution.\n", "\n", "1. How does the number of iterations affect the velocity profile? \n", "\n", " **write your answer here**\n", "\n", "2. Suggest a measure of convergence for the numerical simulation and estimate the number of iterations needed to reach convergence.\n", "\n", " **write your answer here**\n", "\n", "3. When converged, the velocity profile has significantly lower velocities than the parabolic profile that was used as a starting guess. What is the reason for that? Suggest an alternative starting guess.\n", "\n", " **write your answer here**\n", "\n", "4. Calculate the Reynolds number for the turbulent flow.\n", "\n", " **write your answer here**\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Turbulent flow solver\n", "Below you find a turbulent flow solver for fully developed channel flow written in `python`." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "\n", "# FUNCTION: turbulent_flow_solver(y,U,Niter,mu,rho,dpdx,u0,verbose)\n", "\n", "#********************************************************\n", "#\n", "# NO CODE MODIFICATIONS NEEDED IN THIS NOTEBOOK SECTION\n", "#\n", "#********************************************************\n", "\n", "#========================================================\n", "# numerical calculation of a turbulent boundary layer \n", "# velocity profile using a finite-volume approach\n", "#========================================================\n", "#\n", "# 0 1 2 N-1 N N+1\n", "# ------------------- - - - -------------------\n", "# | | | | | | | |\n", "# | o | o | o | | o | o | o |\n", "# | G1 | | | | | | G2 |\n", "# ------------------- - - - -------------------\n", "# ^ ^ \n", "# y=y0 y=h/2 \n", "# u=u0 du/dy=0\n", "#\n", "# The computational domain contains N cells ranging \n", "# from cell 1 to cell N. G1 and G2 are ghost cells \n", "# (cells outside of the computational domain used for\n", "# boundary condition implementation).\n", "#========================================================\n", "\n", "def turbulent_flow_solver(y,U,Niter,mu,rho,dpdx,u0,verbose):\n", " \n", " N = len(y)-2\n", " \n", " relax = 0.5 # under relaxation\n", "\n", " # cell centers (excluding ghost cells): u_(1) to y_(N)\n", " y0=y[slice(1,N+1)]\n", "\n", " # right-shifted cell centers: y_(2) to y_(N+1) \n", " # (right ghost cell (G2) included) \n", " yp = y[slice(2,N+2)]\n", "\n", " # left-shifted cell centers: y_(0) to y_(N-1)\n", " # (left ghost cell (G1) included) \n", " ym = y[slice(0,N)]\n", " \n", " # cell center velocity vector: U_(1) to U_(N)\n", " U0=U[slice(1,N+1)]\n", "\n", " # right-shifted velocity vector: U_(2) to U_(N+1)\n", " # (right ghost cell (G2) included) \n", " Up=U[slice(2,N+2)]\n", "\n", " # left-shifted velocity vector: U_(0) to U_(N-1)\n", " # (left ghost cell (G1) included)\n", " Um=U[slice(0,N)] \n", " \n", " # note! U0, Up, and Um are not copies of U but references to U, which \n", " # means that the values in U0, Up, and Um will be updated when U changes\n", "\n", " # generate right-hand-side vector (b)\n", " b0 = dpdx*dy**2\n", " b = b0*np.ones(N)\n", "\n", " # implement boundary conditions\n", " U[ 0]=u0\n", " U[-1]=U[-2]\n", " \n", " # update velocity vector iteratively until converged \n", " # (increase the number of iterations if not converged) \n", " \n", " for i in range(Niter):\n", "\n", " # store the velocity vector from the previous iteration\n", " Uold = U.copy()\n", "\n", " # generate matrix diagonal and off-diagonal elements\n", "\n", " off_diagonal_p = (rho*kappa**2* 0.25*(yp+y0)**2*(Up-U0)/(yp-y0)+mu)\n", " diagonal = -(rho*kappa**2*(0.25*(yp+y0)**2*(Up-U0)/(yp-y0)+\\\n", " 0.25*(y0+ym)**2*(U0-Um)/(y0-ym))+2*mu)\n", " off_diagonal_m = (rho*kappa**2* 0.25*(y0+ym)**2*(U0-Um)/(y0-ym)+mu)\n", "\n", " # assemble a sparse matrix (A)\n", " A = diags([diagonal,off_diagonal_m[slice(1,N)],off_diagonal_p[slice(0,N-1)]],[0, -1, 1],format=\"csr\")\n", "\n", " # update right-hand-side vector (boundary cell values)\n", " b[ 0] = b0 - off_diagonal_m[ 0]*u0\n", " b[-1] = b0 - off_diagonal_p[-1]*U0[-1]\n", "\n", " # calculate a new velocity vector solving the linear system AU=b\n", " U[slice(1,N+1)] = relax*spsolve(A,b)+(1.-relax)*Uold[slice(1,N+1)]\n", "\n", " # set the velocity in the right ghost cell (G2) to \n", " # U_(N), which ensures that du/dy=0 at the center \n", " U[-1] = U[-2]\n", " \n", " if verbose :\n", " # calculate max difference from previous solution\n", " dUmax = np.max(np.abs(U-Uold))\n", " print('iteration %4d of %4d, max change: %e' %(i+1,Niter,dUmax)) \n", " \n", " return U\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Running the turbulent flow solver\n", "Here you define variables necessary to run the turbulent flow solver. Make sure to run the code above first since the function defined there will be used in the code below." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "\n", "# TASK 2.2 main code\n", "\n", "#********************************************************\n", "#\n", "# CODE MODIFICATION NEEDED: \n", "#\n", "# 1. SPECIFY CELL COUNT\n", "# 2. RUN THE TURBULENT FLOW SOLVER FOR THE NUMBER OF \n", "# ITERATIONS SPECIFIED IN THE TASK (50, 100, 200, 300)\n", "# THE 50 ITERATIONS CASE IS GIVEN AS AN EXAMPLE\n", "# 3. ADD THE RESULTS IN THE PLOT AT THE END OF THE SCRIPT\n", "#\n", "#********************************************************\n", "\n", "# case-specific properties\n", "\n", "dpdx = -8.0 # pressure gradient\n", "\n", "# solver settings\n", "\n", "N = 0 # number of cells <= SPECIFY CELL COUNT HERE\n", "\n", "if N > 0 :\n", "\n", " # we will assume that we are in the log-region and thus the \n", " # lower boundary will coincide with the lower end of the \n", " # log region (here assumed to be at y+=30). The velocity \n", " # at y+=30 is calculated using the log law.\n", "\n", " yplus0 = 30\n", " ustar_turb = np.sqrt(-dpdx*h/(2.0*rho))\n", " y0 = yplus0*mu/(rho*ustar_turb)\n", " u0 = ustar_turb*np.log(yplus0)/kappa+B*ustar_turb\n", "\n", " # calculate cell size (dy) assuming equidistant grid\n", "\n", " dy = (0.5*h-y0)/(N+0.5)\n", "\n", " # generate cell center coordinate vectors\n", "\n", " # cell centers (including ghost cells): y_(0) to y_(N+1)\n", " y_turb = np.linspace(y0,0.5*h+0.5*dy,N+2)\n", "\n", " # initialize velocity vector (U) \n", " # here a laminar solution is used as the initial condition\n", "\n", " U_turb = get_analytic_laminar_velocity_profile(y_turb,h,dpdx,mu)\n", "\n", " #********************************************************\n", " #\n", " # RUN THE TURBULENT FLOW SOLVER WITH THE SPECIFIED INPUT\n", " #\n", " #********************************************************\n", "\n", " # 50 ITERATIONS\n", "\n", " U_turb = turbulent_flow_solver(y_turb,U_turb,50,mu,rho,dpdx,u0,verbose)\n", " U_turb_50 = U_turb.copy()\n", "\n", " # <= RUN THE TURBULENT FLOW SOLVER FOR 100, 200 AND 300 ITERATIONS\n", "\n", " # NOTE! YOU SINCE YOU WILL SEND THE PREVIOUS SOLUTION TO THE SOLVER \n", " # AS START SOLUTION, YOU SHOULD NOT SPECIFY THE TOTAL NUMBER OF ITERATIONS\n", " # BUT RATHER THE NUMBER OF ITERATIONS TO BE ADDED. SO IF YOU WANT TO RUN \n", " # 100 ITERATIONS STARTING FROM U_turb THE PREVIOUS SOLUTION YOU SHOULD \n", " # SPECIFY 50 ITERATIONS\n", "\n", " # set last y-coordinate to channel center\n", " y_turb[-1] = h/2. \n", "\n", " # plot velocity profiles\n", "\n", " fig = plt.figure(num=1, figsize=(15, 10), dpi=80, facecolor='w', edgecolor='k')\n", " ax = fig.add_subplot(111)\n", " ax.plot(U_turb_50,y_turb,linewidth=2,label='50 iterations')\n", " # <= ADD PLOTS FOR 100, 200 & 300 ITERATIONS HERE\n", " ax.set_xlabel('axial velocity (u [m/s])',fontsize=FontSize)\n", " ax.set_ylabel('wall-normal coordinate (y [m])',fontsize=FontSize)\n", " ax.set_title('Velocity profile - turbulent channel flow',fontsize=FontSize)\n", " ax.legend(fontsize=FontSize)\n", " for label in (ax.get_xticklabels() + ax.get_yticklabels()):\n", " label.set_fontsize(FontSize)\n", " ax.grid()\n", " \n", " plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Task 2.3\n", "---\n", "\n", "In the provided `Python` function `turbulent_flow_solve`, there is a variable called `relax` used as follows:\n", "\n", "`U = relax*spsolve(A,b)+(1.-relax)*Uold`\n", "\n", "where `U` is a vector containing the velocity values in all cells, `A` is a coefficient matrix, `b` is the right-hand side vector, `spsolve` is a `Python` function solving linear equation systems (in our case `AU=b`), and `Uold` is a vector containing the velocity values from the previous iteration. What does the `relax` variable do\n", "\n", "### Answer to task 2.3\n", "\n", "**write down your answer here**\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Task 2.4\n", "---\n", "\n", "$u^+$ is a non-dimensional velocity and $y^+$ is a non-dimensional wall-normal coordinate, respectively. The definitions of $u^+$, $y^+$, and the shear velocity $u^\\ast$ are repeated here for convenience\n", "\n", "$$y^+=\\dfrac{yu^\\ast}{\\nu}$$\n", "\n", "$$u^+=\\dfrac{\\bar{u}}{u^\\ast}$$\n", "\n", "$$u^\\ast\\equiv\\sqrt{\\dfrac{\\tau_w}{\\rho}}$$\n", "\n", "1. Why is it common to present the boundary layer velocity profile in terms of $y^+$ and $u^+$ rather than $y$ and $\\bar{u}$?\n", "\n", " **write your answer here**\n", "\n", "2. What range of $y^+$ values corresponds to the viscous sublayer?\n", "\n", " **write your answer here**\n", "\n", "3. What range of $y^+$ values defines the log-region?\n", "\n", " **write your answer here**\n", " \n", "4. How does the turbulence viscosity $\\mu_t$ compare to the fluid viscosity $\\mu$ in the viscous sublayer and in the log-region, respectively?\n", "\n", " **write your answer here**\n", "\n", "5. How does the total shear stress vary with distance from the wall in the viscous sublayer and in the log-region, respectively?\n", "\n", " **write your answer here**\n", "\n", "6. How are the velocity profiles characterized mathematically in the viscous sublayer and in the log-region, respectively?\n", "\n", " **write your answer here**\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Task 2.5\n", "---\n", "\n", "In the task 2.5 main code below, you have been provided with average velocity and corresponding wall-normal coordinates obtained from measurements in a turbulent boundary layer over a **flat plate** in a wind tunnel. \n", "\n", "\n", "1. Obtain an estimate of the shear velocity ($u^\\ast$) from the provided data. Follow the procedure described in [MTF053_CA1.pdf](https://courses.g3dflow.com/fluidmech/docs/MTF053_CA1.pdf) and make changes to the `Python` code below where indicated.\n", "\n", " What $u^\\ast$ value do you get?\n", " \n", " **write your answer here**\n", " \n", "2. With the estimate of $u^\\ast$ from the experimental data, you can now convert the provided $\\overline{u}$ and $y$ values to $u^+$ and $y^+$ values. Compare these values with the corresponding data from your numerical simulation by plotting $u^+$ as a function of $y^+$ for the experiments and the numerical simulation in the same plot (use `semilogx` instead of `plot`). \n", "\n", " Compare the two profiles and comment on any differences and try to draw some conclusions about the accuracy of the numerical solution and/or the measurements based on your observations.\n", " \n", " **write your answers here**\n", "\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "\n", "# TASK 2.5 main code\n", "\n", "#********************************************************\n", "#\n", "# CODE MODIFICATIONS NEEDED: \n", "# \n", "# 1. SPECIFY y+ DEFINING THE LOWER AND UPPER END \n", "# OF THE LOG-REGION\n", "# 2. IMPLEMENT A NEWTON-RAPHSON SOLVER ACCORDING TO \n", "# THE DESCRIPTION ABOVE\n", "#\n", "#********************************************************\n", "\n", "#========================================================\n", "# provided measured boundary-layer data\n", "#========================================================\n", "\n", "y_exp = np.array([0.00015, 0.000231, 0.0003, 0.000481, 0.0006, 0.000706, 0.000931, 0.001181, 0.001406, 0.001656, 0.001881, 0.002131, 0.002356, 0.002581, 0.002831, 0.003056, 0.003306, 0.003775, 0.00415, 0.004525, 0.0049, 0.005275, 0.00605, 0.0068, 0.0072, 0.0076, 0.0083, 0.0095, 0.01057, 0.012, 0.014, 0.016, 0.019, 0.022, 0.025, 0.03])\n", "U_exp = np.array([3.051, 3.531, 3.917, 4.579, 4.879, 5.086, 5.385, 5.637, 5.805, 5.928, 6.05, 6.139, 6.225, 6.32, 6.424, 6.463, 6.531, 6.657, 6.753, 6.828, 6.905, 6.985, 7.148, 7.265, 7.357, 7.417, 7.52, 7.693, 7.829, 7.984, 8.233, 8.449, 8.683, 8.897, 9.034, 9.121])\n", "\n", "#========================================================\n", "# calculate y+ and u+ for the turbulent-channel-flow \n", "# simulation velocity profile\n", "#========================================================\n", "\n", "y_turb_defined = True\n", "\n", "try: y_turb\n", "except NameError: \n", " y_turb_defined = False\n", "\n", "if y_turb_defined :\n", " yplus_turb = y_turb*ustar_turb*rho/mu\n", " uplus_turb = U_turb/ustar_turb \n", "\n", "#========================================================\n", "# calculate y+ and u+ from the provided experimental data\n", "#========================================================\n", "\n", "#--------------------------------------------------------\n", "# step 1: estimate u* \n", "#--------------------------------------------------------\n", "\n", "yplus_low = 0 # <= SPECIFY THE y+ VALUE DEFINING THE LOWER END OF THE LOG-REGION HERE\n", "yplus_high = 0 # <= SPECIFY THE y+ VALUE DEFINING THE UPPER END OF THE LOG-REGION HERE\n", "\n", "ustar_exp = 1.0\n", "\n", "if yplus_high > yplus_low :\n", "\n", " # go through all locations in the boundary layer until a \n", " # location within the log-region is found\n", "\n", " for j in range(len(U_exp)):\n", "\n", " # calculate rough estimates of tau_w, u*, y+, and u+ \n", " # make an estimate of by calculating the wall shear stress as tau_w=mu*dy/du\n", "\n", " tau_w=mu*U_exp[j]/y_exp[j]\n", "\n", " # calculate y+ and u+ from tau_w\n", "\n", " yplus,uplus=calc_yplus_uplus(rho,mu,tau_w,y_exp[j],U_exp[j])\n", "\n", " # check if the rough y+ estimate is within the log-region\n", "\n", " if( yplus > yplus_low/3 and yplus < yplus_high ): \n", "\n", " dtau_w=10.*tau_w\n", "\n", " while( abs(dtau_w) > 0.0001*tau_w ):\n", "\n", " # do a Newton-Raphson iteration\n", "\n", " # <= IMPLEMENT NEWTON-RAPHSON ITERATION HERE (UPDATE dtau_w)\n", "\n", " # update tau_w\n", "\n", " tau_w=tau_w+dtau_w\n", "\n", " # update y+, and u+ \n", "\n", " yplus,uplus=calc_yplus_uplus(rho,mu,tau_w,y_exp[j],U_exp[j])\n", "\n", " # again, verify that the estimate of y+ is within the log region\n", " # if the calculated y+ is within the log region stop, if not move \n", " # one step further out from the wall and try again\n", "\n", " if( yplus > yplus_low and yplus < yplus_high ):\n", " ustar_exp = np.sqrt(tau_w/rho)\n", " print('estimate of u*: %f m/s (found at j=%d, y+=%f)' %(ustar_exp,j,yplus) )\n", " break\n", "\n", "#--------------------------------------------------------\n", "# step 2: calculate y+ and u+ using the estimate of u*\n", "#--------------------------------------------------------\n", "\n", "yplus_exp = (rho*ustar_exp/mu)*y_exp\n", "uplus_exp = U_exp/ustar_exp\n", "\n", "#========================================================\n", "# plot velocity profiles\n", "#========================================================\n", "\n", "if ustar_exp != 1.0 :\n", "\n", " fig = plt.figure(num=1, figsize=(15, 10), dpi=80, facecolor='w', edgecolor='k')\n", " ax = fig.add_subplot(111)\n", "\n", " ax.semilogx(yplus_turb,uplus_turb,linewidth=2,label='numeric')\n", " ax.semilogx(yplus_exp,uplus_exp,linewidth=2,label='measured')\n", " ax.set_xlabel(r'$y^+$',fontsize=FontSize)\n", " ax.set_ylabel(r'$u^+$',fontsize=FontSize)\n", " ax.set_title('Velocity profile (non-dimensional) - turbulent channel flow',fontsize=FontSize)\n", " ax.legend(fontsize=FontSize)\n", " for label in (ax.get_xticklabels() + ax.get_yticklabels()):\n", " label.set_fontsize(FontSize)\n", " ax.grid(which='both')\n", " \n", " plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.6" } }, "nbformat": 4, "nbformat_minor": 4 }