inlet_distortion_effect.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "922d2ea3-6d89-48dc-8a71-5cfe91771ddb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from IPython.display import display, clear_output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09ecb14a-6c1f-4f98-a462-308b97a640f4",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38f3653a-0595-4a81-9389-7e1d5aec02f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "P01 = 160            #kPa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd56c1a1-ed61-456a-9593-6e664fcc74f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "T01 = 400            #kelvin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58a2000a-f4b6-4906-a2d8-75736ed96f0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "P_N = 100            #kPa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38272cc2-8cfb-4554-aef5-a0beff586f78",
   "metadata": {},
   "outputs": [],
   "source": [
    "gama = 1.4\n",
    "R = 0.287            #kj/kg-k\n",
    "C_v = 0.717          #kj/kg-k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf2fc7af-8d19-4fb3-9d51-8584be238761",
   "metadata": {},
   "outputs": [],
   "source": [
    "rho01 = P01/(R*T01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2256c047-6cbc-4d81-a345-fc2e1ced44af",
   "metadata": {},
   "outputs": [],
   "source": [
    "L1 = 3       #meter\n",
    "L2 = 1       #meter\n",
    "L3 = 5       #meter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39a38def-2d53-40ac-9665-4eed132519d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "D1 = 0.2     #meter\n",
    "D3 = 0.1     #meter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34228cbd-9901-4b21-abf5-dea97e0546d6",
   "metadata": {},
   "source": [
    "# Grid Points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1865662b-f641-48a1-8b26-c46a6e40852d",
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = 0.1                                        #meter\n",
    "x1 = np.linspace(0, L1, int(L1/dx+1))\n",
    "x2 = np.linspace(L1, L1+L2, int(L2/dx+1))\n",
    "x3 = np.linspace(L1+L2, L1+L2+L3, int(L3/dx+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd2316c0-e5a5-41a9-939d-ec4021f2d92a",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = int(L1/dx+1)+int(L2/dx+1)+int(L3/dx+1)-2                     #Number of Grid Points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd0ca5f3-83e1-46fa-a10f-eeddfe190369",
   "metadata": {},
   "outputs": [],
   "source": [
    "N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3aa35936-6e2c-4c43-bcfb-c25b80935236",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array(list(x1)+list(x2[1:])+list(x3[1:]), dtype=float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aefc5720-0b0a-4f39-aae9-0d1a20ee8ac7",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.zeros(len(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3627c89-f021-4fca-98a1-b8fc6e9c05bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(0, len(x)):\n",
    "    \n",
    "    if x[i] <= L1:\n",
    "        \n",
    "        A[i] = math.pi*D1**2/4\n",
    "        \n",
    "    elif x[i] <= L1+L2:\n",
    "        \n",
    "        D2 = (x[i]-L1)/L2*(D3-D1)+D1\n",
    "        A[i] = math.pi*D2**2/4\n",
    "        \n",
    "    else:\n",
    "        \n",
    "        A[i] = math.pi*D3**2/4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb2028aa-16fb-472e-8a49-7a5fc457c231",
   "metadata": {},
   "source": [
    "# suddenly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61e2e17c-a280-4fee-8027-b1c08dd24ce1",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_excel('steady.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "463a0c64-a43c-4a2d-b212-8a5d8cbf418d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#primitive variables at t = 0:\n",
    "rho_n = np.array(data['rho'])\n",
    "v_n = np.array(data['v'])\n",
    "T_n = np.array(data['T'])\n",
    "P_n = rho_n*T_n\n",
    "\n",
    "CFL = 0.7\n",
    "\n",
    "n = np.where(x == L1+L2)[0][0]          #index for throat\n",
    "A = A/A[n]\n",
    "iteration = 1\n",
    "\n",
    "#dU1/dt, dU2/dt, dU3/dt:\n",
    "dU1_dt = np.ones(len(x))\n",
    "dU2_dt = np.ones(len(x))\n",
    "dU3_dt = np.ones(len(x))\n",
    "\n",
    "#(dU1/dt)_bar, (dU2/dt)_bar:, (dU3/dt)_bar:\n",
    "dU1_dt_bar = np.ones(len(x))\n",
    "dU2_dt_bar = np.ones(len(x))\n",
    "dU3_dt_bar = np.ones(len(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "732efdec-4bdc-4510-8ea8-3063a8a39f38",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the initial plot\n",
    "fig, ax = plt.subplots()\n",
    "mass_flux_nd = [rho_n[n]*v_n[n]]\n",
    "time = [0]\n",
    "P_nd = [P_n[n]]\n",
    "line1, = ax.plot(time, mass_flux_nd)\n",
    "ax.set_xlabel('time')\n",
    "ax.set_ylabel('Non-Dimensional Mass Flux')\n",
    "ax.set_title('Live Mass Flux')\n",
    "m_dot = 1\n",
    "t = 0\n",
    "\n",
    "while (m_dot>0.01):\n",
    "    \n",
    "    print(f'this is iteration = {iteration}')\n",
    "    dt = np.min(CFL*dx/((gama*R*1000*abs(T_n)*T01)**0.5+abs(v_n)*(gama*R*1000*T01)**0.5))\n",
    "    t += dt\n",
    "    \n",
    "    if np.isnan(dt):\n",
    "        \n",
    "        break\n",
    "    #dt = 0.000015\n",
    "    print(f'dt = {dt}')\n",
    "\n",
    "    #U1, U2, U3:\n",
    "    U1 = rho_n*A\n",
    "    U2 = rho_n*v_n*A\n",
    "    U3 = rho_n*(T_n/(gama-1)+gama/2*v_n**2)*A\n",
    "\n",
    "    #F1, F2, F3:\n",
    "    F1 = U2\n",
    "    F2 = U2**2/U1+(gama-1)/gama*(U3-gama/2*U2**2/U1)\n",
    "    F3 = gama*U2*U3/U1-gama*(gama-1)/2*U2**3/U1**2\n",
    "    J2 = (gama-1)/gama*(U3-gama/2*U2**2/U1)*math.pi/2*A*(D3-D1)/L2\n",
    "\n",
    "    for i in range(0, len(x)-1):\n",
    "    \n",
    "        dU1_dt[i] = -(F1[i+1]-F1[i])/dx\n",
    "        dU2_dt[i] = -(F2[i+1]-F2[i])/dx+J2[i]\n",
    "        dU3_dt[i] = -(F3[i+1]-F3[i])/dx\n",
    "        \n",
    "    #print(f'dU1_dt_bar = {dU1_dt[i]}')\n",
    "    \n",
    "    #U1_np1_bar, U1_np1_bar, U1_np1_bar:    np1: n+1\n",
    "    U1_np1_bar = U1+dU1_dt*dt\n",
    "    U2_np1_bar = U2+dU2_dt*dt\n",
    "    U3_np1_bar = U3+dU3_dt*dt    \n",
    "\n",
    "    #primitive variables:\n",
    "    rho_bar = U1_np1_bar/A\n",
    "    v_bar = U2_np1_bar/U1_np1_bar\n",
    "    T_bar = (gama-1)*(U3_np1_bar/U1_np1_bar-gama/2*v_bar**2)\n",
    "    P_bar = rho_bar*T_bar\n",
    "\n",
    "    #F1_bar, F2_bar, F3_bar, J2_bar:\n",
    "    F1_bar = U2_np1_bar\n",
    "    F2_bar = U2_np1_bar**2/U1_np1_bar+(gama-1)/gama*(U3_np1_bar-gama/2*U2_np1_bar**2/U1_np1_bar)\n",
    "    F3_bar = gama*U2_np1_bar*U3_np1_bar/U1_np1_bar-gama*(gama-1)/2*U2_np1_bar**3/U1_np1_bar**2\n",
    "    J2_bar = (gama-1)/gama*(U3_np1_bar-gama/2*U2_np1_bar**2/U1_np1_bar)*math.pi/2*A*(D3-D1)/L2\n",
    "\n",
    "    for i in range(1, len(x)):\n",
    "    \n",
    "        dU1_dt_bar[i] = -(F1_bar[i]-F1_bar[i-1])/dx\n",
    "        dU2_dt_bar[i] = -(F2_bar[i]-F2_bar[i-1])/dx+J2_bar[i]\n",
    "        dU3_dt_bar[i] = -(F3_bar[i]-F3_bar[i-1])/dx\n",
    "        \n",
    "    print(f'dU1_dt_bar = {dU1_dt_bar[i]}')\n",
    "    \n",
    "    #average of time derivatives:\n",
    "    dU1_dt_avg = 0.5*(dU1_dt_bar+dU1_dt)\n",
    "    dU2_dt_avg = 0.5*(dU2_dt_bar+dU2_dt)\n",
    "    dU3_dt_avg = 0.5*(dU3_dt_bar+dU3_dt)\n",
    "\n",
    "    #U1_np1, U1_np2, U1_np3\n",
    "    U1_np1 = U1+dU1_dt_avg*dt\n",
    "    U2_np1 = U2+dU2_dt_avg*dt\n",
    "    U3_np1 = U3+dU3_dt_avg*dt\n",
    "\n",
    "    #decoding primitive variables:\n",
    "    rho_np1 = U1_np1/A\n",
    "    v_np1 = U2_np1/U1_np1\n",
    "    T_np1 = (gama-1)*(U3_np1/U1_np1-gama/2*(U2_np1/U1_np1)**2)\n",
    "    P_np1 = rho_np1*T_np1\n",
    "    \n",
    "    #subsonic inlet:\n",
    "    P_np1[0] = P01+0.2*P01*math.sin(2*math.pi*40*t)\n",
    " \n",
    "    #outlet:\n",
    "    P_np1[N-1] = P_N/P01\n",
    "    rho_np1[N-1] = rho_n[N-1]\n",
    "    a = (gama*R*1000*T_np1[N-1]*T01)\n",
    "    v_np1[N-1] = v_n[N-1]-dt*a*v_n[N-1]*math.log(A[N-1]/A[N-2])\n",
    "\n",
    "    #residuals:\n",
    "    drho_dt_avg = abs((rho_np1[n]-rho_n[n])/dt)\n",
    "    dv_dt_avg = abs((v_np1[n]-v_n[n])/dt)\n",
    "    dT_dt_avg = abs((T_np1[n]-T_n[n])/dt)\n",
    "    \n",
    "    rho_n = rho_np1.copy()\n",
    "    v_n = v_np1.copy()    \n",
    "    T_n = T_np1.copy() \n",
    "    P_n = P_np1.copy()\n",
    "    \n",
    "    #delta_m_dot\n",
    "    m_dot = rho_n[n]*v_n[n]\n",
    "    print(f'nd mass flux = {m_dot}')\n",
    "    \n",
    "    # Append the current iteration number and residual values to the x and y arrays\n",
    "    mass_flux_nd.append(m_dot)\n",
    "    time.append(t)\n",
    "    P_nd.append(P_n[n])\n",
    "    \n",
    "    line1.set_data(time, mass_flux_nd)\n",
    "    ax.relim()\n",
    "    ax.autoscale_view()\n",
    "    \n",
    "    # Update and display the plot after each iteration\n",
    "    clear_output(wait=True)\n",
    "    display(fig)\n",
    "    plt.pause(0.001)\n",
    "    \n",
    "    iteration += 1\n",
    "        \n",
    "    if iteration==40000:\n",
    "        break\n",
    "# Show the final plot after all iterations are completed\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce8608c2-359d-4238-9a78-086da7f979c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time, P_nd, linestyle = \"solid\")\n",
    "plt.xlabel('time[second]')\n",
    "plt.ylabel('Non-Dimensional Pressure')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d16f1b2-fff7-4c81-a2b8-5e5988995aa9",
   "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}