slide_set_09_cylinder_with_circulation.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Slide Set 9 - Cylinder flow\n",
    "\n",
    "This script plots streamlines and potential lines for the flow past a cylinder. It results from the superposition of a dipole and a parallel flow.\n",
    "\n",
    "## Preparation\n",
    "\n",
    "To prepare the plot, several helper functions have to be defined. This is done in the following cell.\n",
    "First, the two libraries numpy and matplotlib.pyplot have to be imported in order to user their functionality. Note, these libraries are helpful in many scenarios, when data should be plotted. Next, helper classes that provide the complex velocity for the source-/sink-flow and the parallel flow are implemented. Finally, a function providing the meshgrid for the calculations is implemented."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Define important colors\n",
    "TUMOrange = \"#E37222\"\n",
    "TUMBlue = \"#3070b3\"\n",
    "TUMGreen = \"#A2AD00\"\n",
    "\n",
    "# Define helper classes to get the complex velocity of elementary flows\n",
    "class ParallelFlow:\n",
    "\n",
    "    def __init__(self, velocity_x: float, velocity_y: float) -> None:\n",
    "        \"\"\"Initializer for the parallel flow class. Stores relevant quantities as class members.\n",
    "\n",
    "        Args:\n",
    "            velocity_x (float): Cartesion component of the velocity pointing in x-direction.\n",
    "            velocity_y (float): Cartesion component of the velocity pointing in y-direction.\n",
    "        \"\"\"\n",
    "        self.velocity_x = velocity_x\n",
    "        self.velocity_y = velocity_y\n",
    "\n",
    "    def get_complex_velocity(self, z: np.array) -> np.array:\n",
    "        \"\"\"Return the complex velocity of the parallel flow.\n",
    "\n",
    "        Args:\n",
    "            z (np.array): Complex numpy array containing the meshgrid of the complex variable z.\n",
    "\n",
    "        Returns:\n",
    "            np.array: The complex velocity as a complex numpy array.\n",
    "        \"\"\"\n",
    "        return self.velocity_x - 1.0j * self.velocity_y\n",
    "\n",
    "\n",
    "class DipoleFlow:\n",
    "\n",
    "    def __init__(self, dipole_moment, translation_x = 0.0, translation_y = 0.0) -> None:\n",
    "        \"\"\"Initializer for the dipole class. Stores relevant quantities as class members.\n",
    "\n",
    "        Args:\n",
    "            dipole_moment (float): The dipole moment.\n",
    "            translation_x (float, optional): The translation of the singularity in x direction. Defaults to 0.0.\n",
    "            translation_y (float, optional): The translation of the singularity in y direction. Defaults to 0.0.\n",
    "        \"\"\"\n",
    "        self.dipole_moment = dipole_moment\n",
    "        self.translation_x = translation_x\n",
    "        self.translation_y = translation_y\n",
    "        self.factor = self.dipole_moment / np.pi\n",
    "\n",
    "    def get_complex_velocity(self, z):\n",
    "        \"\"\"Returns the complex velocity of the source-/sink-flow.\n",
    "\n",
    "        Args:\n",
    "            z (np.array): Complex numpy array containing the meshgrid of the complex variable z.\n",
    "\n",
    "        Returns:\n",
    "            np.array: The complex velocity as a complex numpy array.\n",
    "        \"\"\"\n",
    "        z_translated = z - self.translation_x - 1.0j * self.translation_y\n",
    "        return - self.factor / (z_translated*z_translated)\n",
    "    \n",
    "def get_grid_quantities(x_min: float, x_max: float, y_min: float, y_max: float, resolution_x: int, resolution_y: int):\n",
    "    \"\"\"Returns the meshgrid used for the calculations.\n",
    "\n",
    "    Args:\n",
    "        x_min (float): The minimum x value.\n",
    "        x_max (float): The maximum x value.\n",
    "        y_min (float): The minimum y value.\n",
    "        y_max (float): The maximum y value.\n",
    "        resolution_x (int): The number of points to discretize in x direction.\n",
    "        resolution_y (int): The number of points to discretize in y direction.\n",
    "\n",
    "    Returns:\n",
    "        List[np.array]: The numpy array containing the meshgrid.\n",
    "    \"\"\"\n",
    "    x_linspace = np.linspace(x_min, x_max, resolution_x)\n",
    "    y_linspace = np.linspace(y_min, y_max, resolution_y)\n",
    "    X_meshgrid, Y_meshgrid = np.meshgrid(x_linspace, y_linspace)\n",
    "    Z_meshgrid = X_meshgrid + 1.0j * Y_meshgrid\n",
    "    return X_meshgrid, Y_meshgrid, Z_meshgrid\n",
    "\n",
    "class PotentialVortexFlow:\n",
    "\n",
    "    def __init__(self, gamma: float, translation_x: float = 0.0, translation_y: float = 0.0) -> None:\n",
    "        \"\"\"Initializer for the vortex flow class. Stores relevant quantities as class members.\n",
    "\n",
    "        Args:\n",
    "            gamma (float): The strength of the sourcce flow. It is positive for source flow and negative for sink flow.\n",
    "            translation_x (float, optional): The translation of the singularity in x direction. Defaults to 0.0.\n",
    "            translation_y (float, optional): The translation of the singularity in y direction. Defaults to 0.0.\n",
    "        \"\"\"\n",
    "        self.translation_x = translation_x\n",
    "        self.translation_y = translation_y\n",
    "        self.gamma = gamma\n",
    "        self.factor = self.gamma / (2.0 * np.pi)\n",
    "\n",
    "    def get_complex_velocity(self, z: np.array) -> np.array:\n",
    "        \"\"\"Returns the complex velocity of the vortex flow.\n",
    "\n",
    "        Args:\n",
    "            z (np.array): Complex numpy array containing the meshgrid of the complex variable z.\n",
    "\n",
    "        Returns:\n",
    "            np.array: The complex velocity as a complex numpy array.\n",
    "        \"\"\"\n",
    "        z = z - self.translation_x - 1.0j * self.translation_y\n",
    "        return - 1.0j * self.factor / z"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definition of domain\n",
    "\n",
    "Here, relevant information to describe the domain is provided. It is used to generate the meshgrid for the calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_min = -4.0\n",
    "x_max = 4.0\n",
    "y_min = -4.0\n",
    "y_max = 4.0\n",
    "resolution_per_unit_length = 50\n",
    "resolution_x = ( x_max - x_min ) * resolution_per_unit_length\n",
    "resolution_y = ( y_max - y_min ) * resolution_per_unit_length\n",
    "\n",
    "X_meshgrid, Y_meshgrid, Z_meshgrid = get_grid_quantities(x_min, x_max, y_min, y_max, int(resolution_x), int(resolution_y))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculation of the complex velocity\n",
    "\n",
    "Initializes the relevant source, sink and parallel flow objects to calculate the complex velocity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_infty = 2.0\n",
    "radius = 0.5\n",
    "theta_stagnation_point = 0.0 * np.pi\n",
    "\n",
    "dipole_moment = u_infty * np.pi * radius**2\n",
    "\n",
    "gamma = np.imag(np.exp(1.0j * theta_stagnation_point))\n",
    "circulation = 4.0 * np.pi * u_infty * radius * gamma\n",
    "\n",
    "dipole = DipoleFlow(dipole_moment)\n",
    "inflow = ParallelFlow(u_infty, 0.0)\n",
    "potential_vortex = PotentialVortexFlow(circulation)\n",
    "\n",
    "complex_velocity = inflow.get_complex_velocity(Z_meshgrid) + dipole.get_complex_velocity(Z_meshgrid) + potential_vortex.get_complex_velocity(Z_meshgrid)\n",
    "u = np.real(complex_velocity)\n",
    "v = - np.imag(complex_velocity)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting of the results\n",
    "\n",
    "Plot the streamlines and potential lines.\n",
    "Streamlines are plotted in TUMBlue.\n",
    "Potential lines are plotted in TUMOrange.\n",
    "The stagnation streamline is plotted in red."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "x_linspace = np.linspace(x_min, x_max, 25)\n",
    "y_linspace = np.linspace(x_min, x_max, 25)\n",
    "\n",
    "# Plot the streamlines.\n",
    "\n",
    "# The streamline seeds describe all the points where streamlines start. Thus, only streamline the go through these points are drawn.\n",
    "# This may heavily influence the look of the plot!\n",
    "# In this case only points to the left (*xmin refers to the left, replace with e.g. 0.0 to also obtain streamlines inside the oval)\n",
    "# of the oval are choosen. Thus, no streamlines inside the oval show.\n",
    "# In this case, they are equally distributed along the y direction.\n",
    "streamline_seeds = np.stack([np.ones_like(y_linspace)*x_min, y_linspace], axis=-1)\n",
    "\n",
    "ax.streamplot(\n",
    "    X_meshgrid, Y_meshgrid, u, v,\n",
    "    density=[1.0, 1.0],\n",
    "    minlength=1.0, maxlength=20.0,\n",
    "    broken_streamlines=False,\n",
    "    start_points=streamline_seeds,\n",
    "    color=TUMBlue,\n",
    ")\n",
    "\n",
    "# Plot the potential lines.\n",
    "\n",
    "# The potential line seed describe all points where potential lines start. Here, similar remarks as above apply.\n",
    "potential_line_seeds_1 = np.stack([x_linspace, np.ones_like(x_linspace)* 0.9 * y_min], axis=-1)\n",
    "potential_line_seeds_2 = np.stack([x_linspace, np.ones_like(x_linspace)* 0.9 * y_max], axis=-1)\n",
    "potential_line_seeds = np.concatenate([potential_line_seeds_1, potential_line_seeds_2], axis=0)\n",
    "ax.streamplot(\n",
    "    X_meshgrid, Y_meshgrid, -v, u,\n",
    "    density=[0.2, 0.2],\n",
    "    minlength=0.0, maxlength=20.0,\n",
    "    broken_streamlines=False,\n",
    "    start_points=potential_line_seeds,\n",
    "    arrowstyle = \"-\",\n",
    "    color=TUMOrange,\n",
    ")\n",
    "\n",
    "# ax.plot(\n",
    "#     [x_min, x_max], [0.0, 0.0],\n",
    "#     color=\"red\",\n",
    "#     linewidth=3.0,\n",
    "# )\n",
    "# Plot circle\n",
    "phi_linspace = np.linspace(0.0, 2.0*np.pi, 100)\n",
    "ax.plot(\n",
    "    radius * np.cos(phi_linspace), radius * np.sin(phi_linspace),\n",
    "    color=\"red\",\n",
    "    linewidth=3.0,\n",
    ")\n",
    "phi_array = np.array([theta_stagnation_point, np.pi-theta_stagnation_point])\n",
    "ax.scatter(\n",
    "    radius * np.cos(phi_array), radius * np.sin(phi_array),\n",
    "    color=TUMGreen,\n",
    "    linewidth=3.0,\n",
    ")\n",
    "\n",
    "ax.set_xlim([0.5 * x_min, 0.5 * x_max])\n",
    "ax.set_ylim([0.5 * y_min, 0.5 * y_max])\n",
    "ax.set_aspect('equal', 'box')\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting streamlines and potential lines without prescribing any seed, and by that not imposing prior knowledge, is also possible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# Plot the streamlines.\n",
    "ax.streamplot(\n",
    "    X_meshgrid, Y_meshgrid, u, v,\n",
    "    density=[1.0, 1.0],\n",
    "    color=TUMBlue,\n",
    ")\n",
    "\n",
    "# Plot the potential lines.\n",
    "ax.streamplot(\n",
    "    X_meshgrid, Y_meshgrid, -v, u,\n",
    "    density=[0.2, 0.2],\n",
    "    arrowstyle = \"-\",\n",
    "    color=TUMOrange,\n",
    ")\n",
    "\n",
    "ax.set_aspect('equal', 'box')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pressure distribution along the cylinder's surface.\n",
    "\n",
    "Here, we calculate and plot the pressure coefficient along the surface of the cylinder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution_around_cylinder = 200\n",
    "\n",
    "phi_linspace = np.linspace(0.0, np.pi, resolution_around_cylinder)\n",
    "x_circle = radius * np.cos(phi_linspace)\n",
    "y_circle = radius * np.sin(phi_linspace)\n",
    "z_circle = x_circle + 1.0j * y_circle\n",
    "complex_velocity_surface = inflow.get_complex_velocity(z_circle) + dipole.get_complex_velocity(z_circle) + potential_vortex.get_complex_velocity(z_circle)\n",
    "q_squared_surface = np.absolute(complex_velocity_surface)**2\n",
    "pressure_coefficient_surface_upper = 1.0 - q_squared_surface / u_infty**2\n",
    "\n",
    "phi_linspace = np.linspace(np.pi, 2.0 * np.pi, resolution_around_cylinder)\n",
    "x_circle = radius * np.cos(phi_linspace)\n",
    "y_circle = radius * np.sin(phi_linspace)\n",
    "z_circle = x_circle + 1.0j * y_circle\n",
    "complex_velocity_surface = inflow.get_complex_velocity(z_circle) + dipole.get_complex_velocity(z_circle) + potential_vortex.get_complex_velocity(z_circle)\n",
    "q_squared_surface = np.absolute(complex_velocity_surface)**2\n",
    "pressure_coefficient_surface_lower = 1.0 - q_squared_surface / u_infty**2\n",
    "\n",
    "plt.close()\n",
    "fig, ax = plt.subplots()\n",
    "# ax.set_aspect('equal', 'box')\n",
    "\n",
    "ax.scatter(np.real(z_circle), -pressure_coefficient_surface_upper, c=TUMGreen)\n",
    "ax.scatter(np.real(z_circle), -pressure_coefficient_surface_lower, c=TUMOrange)\n",
    "\n",
    "x_s1 = np.real(radius * np.exp(1.0j * theta_stagnation_point))\n",
    "x_s2 = np.real(radius * np.exp(1.0j * (np.pi - theta_stagnation_point)))\n",
    "\n",
    "ax.scatter([x_s1, x_s2], [-1.0, -1.0], c=TUMBlue, marker=\"*\", s=200)\n",
    "\n",
    "ax.legend([\"Upper half of the cylinder\", \"Lower half of the cylinder\", \"Stagnation points\"])\n",
    "ax.set_xlabel(\"x coordinate of point on the cylinder\")\n",
    "ax.set_ylabel(\"- C_p value\")\n",
    "ax.set_ylim(ymin=-1.2)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "plot",
   "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.10.14"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}