diff --git a/examples/Fluid_Structure_Interaction/Fluid_Structure_Interaction.ipynb b/examples/Fluid_Structure_Interaction/Fluid_Structure_Interaction.ipynb new file mode 100644 index 00000000..46f5e9e9 --- /dev/null +++ b/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 +}