Fluid Structure Interaction Example

examples/Fluid_Structure_Interaction/Fluid_Structure_Interaction.ipynb

@@ -0,0 +1,345 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "under-sierra",
+   "metadata": {},
+   "source": [
+    "# Fluid Structure Interaction with SpatialPy"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "absolute-stamp",
+   "metadata": {},
+   "source": [
+    "## Definition of the model"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "caroline-collective",
+   "metadata": {},
+   "source": [
+    "### Imports and definitions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "promotional-profile",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "import sys\n",
+    "sys.path.insert(1, \"../..\")\n",
+    "import math\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "import spatialpy"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "matched-batch",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class Density(spatialpy.BoundaryCondition):\n",
+    "    def __init__(self, rhoF, rhoB):\n",
+    "        self.rhoF = rhoF\n",
+    "        self.rhoB = rhoB\n",
+    "        \n",
+    "    def expression(self):\n",
+    "        bcstr = \"if(me->type == 1 || me->type == 3){\"\n",
+    "        bcstr += f\"me->rho = {self.rhoF};\"\n",
+    "        bcstr += \"}\"\n",
+    "        bcstr += \"if(me->type == 2){\"\n",
+    "        bcstr += f\"me->rho = {self.rhoB};\"\n",
+    "        bcstr += \"}\"\n",
+    "        return bcstr"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "biblical-parallel",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class Teleport(spatialpy.BoundaryCondition):\n",
+    "    def __init__(self, nu, v0):\n",
+    "        self.nu = nu\n",
+    "        self.v0 = v0\n",
+    "        \n",
+    "    def expression(self):\n",
+    "        bcstr = \"if(me->x[0] > system->xhi || me->x[0] < system->xlo){\"\n",
+    "        bcstr += \"me->x[0] = system->xlo + 1e-6;\"\n",
+    "        bcstr += f\"me->v[0] = {self.v0};\"\n",
+    "        bcstr += f\"me->nu = {self.nu};\"\n",
+    "        bcstr += \"}\"\n",
+    "        bcstr += \"if(me->x[1] > system->yhi){\"\n",
+    "        bcstr += \"me->v[1]= 0.0;\"\n",
+    "        bcstr += \"me->x[1]= 99e-6;\"\n",
+    "        bcstr += \"}\"\n",
+    "        bcstr += \"if(me->x[1] < system->ylo){\"\n",
+    "        bcstr += \"me->v[1]= 0.0;\"\n",
+    "        bcstr += \"me->x[1]= 1e-6;\"\n",
+    "        bcstr += \"}\"\n",
+    "        bcstr += \"me->x[2] = 0;\"\n",
+    "        return bcstr"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "direct-tutorial",
+   "metadata": {},
+   "source": [
+    "### Model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "emerging-government",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class FluidStructureInteraction(spatialpy.Model):\n",
+    "    FLUID = 1\n",
+    "    WALLS = 2\n",
+    "    BEAM = 3\n",
+    "    \n",
+    "    def __init__(self, model_name=\"Fluid Structure Interaction\", Ra=1e4):\n",
+    "        spatialpy.Model.__init__(self, model_name)\n",
+    "        \n",
+    "        # System Constants\n",
+    "        Lxint = 300e-6\n",
+    "        Lyint = 100e-6\n",
+    "        Lbzint = -50e-6\n",
+    "        Nwall = 3\n",
+    "        rhof0 = 1000              # Density of fluid\n",
+    "        rhob0 = 7850              # Density of beam\n",
+    "        nu = 1e-3                 # Inlet viscosity\n",
+    "        v0 = 0.0333               # Inlet velocity\n",
+    "        Pratio = 0.33             # Poisson's Ratio\n",
+    "        E = 2e5                   # Young's modulus\n",
+    "        G = E/(2*(1+Pratio))      # Shear Modulus\n",
+    "        K = E/(3*(1-2*Pratio))    # Bulk Modulus\n",
+    "        cb0 = math.sqrt(K/Pratio) # Artificial speed of sound of beam\n",
+    "        cf0 = 10*v0               # Artificial speed of sound of fluid\n",
+    "        \n",
+    "        # Discretization\n",
+    "        Npx = 60\n",
+    "        \n",
+    "        # Compute Domain Bounds\n",
+    "        deltaf = Lyint/Npx\n",
+    "        deltab = 0.6*deltaf\n",
+    "        Wthick = Nwall*deltaf\n",
+    "        Lx = Lxint-Lbzint\n",
+    "        Ly = Lyint+2*Wthick\n",
+    "        yminint = 0.0\n",
+    "        ymaxint = Lyint\n",
+    "        ymin = yminb = -Wthick\n",
+    "        ymax = Lyint+Wthick\n",
+    "        xminb = 100e-6\n",
+    "        xmaxb = 105e-6\n",
+    "        ymaxb = 50e-6\n",
+    "        cy = Lyint/2\n",
+    "        \n",
+    "        # Compute Volume and Mass per Particle\n",
+    "        vol = Lx*Ly\n",
+    "        Wvol = 2*Wthick*Lx\n",
+    "        Fvol = (xmaxb-xminb)*(ymaxb-yminb)\n",
+    "        Bvol = vol-Wvol-Fvol\n",
+    "        mPFP = (Fvol*rhof0)/(Npx*((0-(-50e-6))/deltaf))\n",
+    "        mPBP = (Bvol*rhob0)/(((xmaxb-xminb)/deltab)*((ymaxb-ymin+deltab)/deltab))\n",
+    "        \n",
+    "        # Domain\n",
+    "        eps = 1e-12\n",
+    "        h = 3*deltaf\n",
+    "        \n",
+    "        domain = spatialpy.Domain(\n",
+    "            0, xlim=(Lbzint, Lxint), ylim=(ymin, ymax), zlim=(0,0), gravity=[0, -1, 0]\n",
+    "        )\n",
+    "        # Fluid and Wall Regions\n",
+    "        for x in np.arange(Lbzint, Lxint, deltaf):\n",
+    "            # Fluid Region\n",
+    "            for y in np.arange(yminint, ymaxint, deltaf):\n",
+    "                if x < 0:\n",
+    "                    domain.add_point(\n",
+    "                        [x, y, 0], type=self.FLUID, vol=Fvol, mass=mPFP, nu=nu, fixed=False\n",
+    "                    )\n",
+    "            # Top Wall Region\n",
+    "            for y in np.arange(ymin, yminint, deltaf):\n",
+    "                domain.add_point(\n",
+    "                    [x, y, 0], type=self.WALLS, vol=Wvol, mass=mPFP, nu=nu, fixed=True\n",
+    "                )\n",
+    "            # Bottom Wall Region\n",
+    "            for y in np.arange(ymaxint, ymax, deltaf):\n",
+    "                domain.add_point(\n",
+    "                    [x, y, 0], type=self.WALLS, vol=Wvol, mass=mPFP, nu=nu, fixed=True\n",
+    "                )\n",
+    "        # Beam Region\n",
+    "        for x in np.arange(xminb, xmaxb, deltab):\n",
+    "            for y in np.arange(ymin, ymaxb, deltab):\n",
+    "                domain.add_point(\n",
+    "                    [x, y, 0], type=self.BEAM, vol=Bvol, mass=mPBP, nu=nu, fixed=False\n",
+    "                )\n",
+    "        domain.dimensions = 2\n",
+    "        self.add_domain(domain)\n",
+    "        \n",
+    "        # Static Domain\n",
+    "        self.staticDomain = False\n",
+    "                \n",
+    "        # Boundary Conditions\n",
+    "        self.add_boundary_condition(Density(rhof0, rhob0))\n",
+    "        self.add_boundary_condition(Teleport(nu, v0))\n",
+    "                \n",
+    "        # Timespan\n",
+    "        dt = 1e-8\n",
+    "        nt = 1000000000\n",
+    "        freq_results = 10000\n",
+    "        self.timespan(np.arange(0, nt*dt+dt, freq_results*dt), timestep_size=dt)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "progressive-secretary",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model = FluidStructureInteraction()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "liable-swaziland",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model.domain.plot_types()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "american-simulation",
+   "metadata": {},
+   "source": [
+    "## Running model and processing the results"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "technological-gallery",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from spatialpy import Solver\n",
+    "solver = Solver(model=model, debug_level=0)\n",
+    "%time solver.compile()\n",
+    "print(solver.build_dir)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "modular-roller",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%time results = solver.run()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "lesbian-timothy",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "results.plot_property(\"type\",t_ndx=0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "homeless-cooler",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "results.plot_property(\"type\",t_ndx=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fundamental-rhythm",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#results.plot_property(\"type\", animated=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "certified-condition",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "points, properties = results.read_step(2)\n",
+    "for i, point in enumerate(points):\n",
+    "    if properties['type'][i] == 1:\n",
+    "        print(point)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "superior-conditioning",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(model.domain.xlim)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "disturbed-proposition",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(model.domain.ylim)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}