{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# Rotating domains" ] }, { "cell_type": "markdown", "id": "1", "metadata": {}, "source": [ "We model configurations with rotating sub-domains, like electric motors, or wind mills. We generate independent meshes for the components, and glue them together using Nitsche's method." ] }, { "cell_type": "code", "execution_count": null, "id": "2", "metadata": {}, "outputs": [], "source": [ "from netgen.occ import *\n", "from ngsolve import *\n", "from ngsolve.webgui import Draw" ] }, { "cell_type": "code", "execution_count": null, "id": "3", "metadata": {}, "outputs": [], "source": [ "tau = 0.0003\n", "tend = 1\n", "omega = 2*pi" ] }, { "cell_type": "code", "execution_count": null, "id": "4", "metadata": {}, "outputs": [], "source": [ "square = MoveTo(0,0).Rectangle(1,1).Face()\n", "circo = Circle((0.5,0.5), 0.30001).Face()\n", "circ = Circle((0.5,0.5), 0.3).Face()\n", "holes = Circle((0.65,0.5), 0.05).Face() + Circle((0.35,0.5), 0.05).Face()\n", "\n", "square.edges.name=\"outer\"\n", "circ.edges.name=\"gammai\"\n", "holes.edges.name=\"hole\"\n", "circo.edges.name=\"gammao\"\n", "square.edges.name=\"wall\"\n", "square.edges.Min(X).name=\"inlet\"\n", "square.edges.Max(X).name=\"outlet\"\n", "outer = square-circo\n", "outer.faces.name = \"outer\"\n", "\n", "circ.faces.name = \"inner\"\n", "\n", "both = Compound([outer, circ-holes])\n", "mesh = Mesh(OCCGeometry(both, dim=2).GenerateMesh(maxh=0.05)).Curve(3)\n", "print (mesh.GetMaterials(), mesh.GetBoundaries())\n", "Draw (mesh);" ] }, { "cell_type": "markdown", "id": "5", "metadata": {}, "source": [ "Define a GridFunction deformation describing the current configuration:" ] }, { "cell_type": "code", "execution_count": null, "id": "6", "metadata": {}, "outputs": [], "source": [ "fesdef = VectorH1(mesh, order=3)\n", "\n", "deformation = GridFunction(fesdef)\n", "defold = GridFunction(fesdef)\n", "\n", "def MeshRotation(angle, deform=None):\n", " mesh.UnsetDeformation()\n", " if not deform: deform = GridFunction(fesdef)\n", " rotmat = CF( (cos(angle), -sin(angle), sin(angle), cos(angle))).Reshape( (2,2))\n", " center = CF( (0.5, 0.5) )\n", " pos = CF( (x,y) )\n", "\n", " deform.Set( (rotmat-Id(2))*(pos-center), definedon=mesh.Materials(\"inner\"))\n", " return deform" ] }, { "cell_type": "code", "execution_count": null, "id": "7", "metadata": {}, "outputs": [], "source": [ "scene = Draw(mesh)\n", "tau1 = 1e-3\n", "for step in range(int(tend/tau1)):\n", " MeshRotation(step*omega*tau1, deformation)\n", " mesh.SetDeformation(deformation)\n", " scene.Redraw()" ] }, { "cell_type": "markdown", "id": "8", "metadata": {}, "source": [ "Solve for a flow potential such that $\\frac{\\partial \\phi}{\\partial n}$ matches the normal velocity: " ] }, { "cell_type": "code", "execution_count": null, "id": "9", "metadata": {}, "outputs": [], "source": [ "fest = H1(mesh, order=3, dirichlet=\"inlet|outlet\")\n", "\n", "festgrad = VectorH1(mesh, order=3)\n", "gfutgrad = GridFunction(festgrad)\n", "\n", "ut,vt = fest.TnT()\n", "\n", "n = specialcf.normal(2)\n", "h = specialcf.mesh_size\n", "\n", "gfut = GridFunction(fest)\n", " \n", "meshVelocity = (deformation-defold) / tau\n", "\n", "at = BilinearForm(grad(ut)*grad(vt)*dx)\n", "ft = LinearForm(fest)\n", "ft += -InnerProduct(meshVelocity,n)*vt*ds(definedon=\"hole\")\n", "\n", "contactt = ContactBoundary(mesh.Boundaries(\"gammai\"), mesh.Boundaries(\"gammao\"), volume=True)\n", "contactt.AddIntegrator (3/h*(ut-ut.Other())*(vt-vt.Other()))\n", "contactt.AddIntegrator (n*grad(ut)*(vt.Other()-vt)+n*grad(vt)*(ut.Other()-ut))\n", "\n", "\n", "def solveWind(gfut,at,ft):\n", " contactt.Update (deformation, bf=at, intorder=10)\n", "\n", " at.Assemble()\n", " ft.Assemble()\n", "\n", " # the solution field \n", " gfut.Set((x), BND)\n", " rt = ft.vec.CreateVector()\n", " rt.data = ft.vec - at.mat * gfut.vec\n", " gfut.vec.data += at.mat.Inverse(freedofs=fest.FreeDofs(), inverse = \"sparsecholesky\") * rt\n", " gfutgrad.Set(Grad(gfut))" ] }, { "cell_type": "code", "execution_count": null, "id": "10", "metadata": {}, "outputs": [], "source": [ "scene = Draw(gfut)\n", "tau1 = 2e-3\n", "for step in range(int(tend/tau1)):\n", " defold.vec.data = deformation.vec\n", " MeshRotation(step*omega*tau1, deformation)\n", "\n", " mesh.SetDeformation(deformation)\n", " solveWind(gfut,at,ft)\n", " scene.Redraw()" ] }, { "cell_type": "markdown", "id": "11", "metadata": {}, "source": [ "Solve a transport problem:" ] }, { "cell_type": "code", "execution_count": null, "id": "12", "metadata": {}, "outputs": [], "source": [ "fes = L2(mesh, order=3) \n", "u,v = fes.TnT()\n", "\n", "feshat = FacetFESpace(mesh, order=3)\n", "uhat, vhat = feshat.TnT()\n", "traceop = fes.TraceOperator(feshat, average=True)\n", "\n", "mesh.SetDeformation(MeshRotation(0))\n", "\n", "wind = -(meshVelocity - gfutgrad)\n", "\n", "a = BilinearForm(fes) \n", "a += -wind*u*grad(v)*dx\n", "uup = IfPos(wind*n, u, u.Other(bnd=0))\n", "a += wind*n*uup*v * dx(element_boundary=True) # upwind\n", "\n", "\n", "ahat = BilinearForm(feshat)\n", "\n", "f = LinearForm(fes)\n", "f.Assemble()\n", "\n", "gfu = GridFunction(fes)\n", "gfu.Set(exp(-10**2*((x-0.15)**2 +(y-0.5)**2)))\n", "\n", "solveWind(gfut,at,ft)\n", "scene = Draw(gfu, min=0, max=2, order=3, autoscale=False)" ] }, { "cell_type": "code", "execution_count": null, "id": "13", "metadata": {}, "outputs": [], "source": [ "contact = ContactBoundary(mesh.Boundaries(\"gammao\"), mesh.Boundaries(\"gammai\"), volume=False)\n", "nc = contact.normal\n", "\n", "term = gfutgrad*nc * IfPos (gfutgrad*nc, uhat.Trace()*(-vhat.Trace().Other()), \\\n", " uhat.Trace().Other()*(vhat.Trace()))\n", "\n", "contact.AddIntegrator (term)" ] }, { "cell_type": "code", "execution_count": null, "id": "14", "metadata": {}, "outputs": [], "source": [ "t = 0\n", "cnt = 0\n", "deformation.vec[:] = 0\n", "w = gfu.vec.CreateVector()\n", "gfuhat = GridFunction(feshat)\n", "# what = gfuhat.vec.CreateVector()\n", "\n", "invm = fes.Mass(rho=1).Inverse()\n", "\n", "with TaskManager():\n", " while t < tend:\n", " defold.vec.data = deformation.vec\n", " MeshRotation(t*omega, deformation)\n", " \n", " contact.Update (deformation, bf=ahat, intorder=10)\n", " # apply the transport operator\n", " mesh.SetDeformation(deformation)\n", " solveWind(gfut,at,ft)\n", " \n", " gfuhat.vec[:] = traceop * gfu.vec\n", " w[:] = a.Apply (gfu.vec) + traceop.T * ahat.Apply(gfuhat.vec)\n", " \n", " gfu.vec.data -= tau * invm * w\n", " \n", " if cnt%10 == 0:\n", " mesh.SetDeformation(deformation)\n", " scene.Redraw() \n", "\n", " t += tau\n", " cnt +=1 \n" ] }, { "cell_type": "code", "execution_count": null, "id": "15", "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.13.4" } }, "nbformat": 4, "nbformat_minor": 5 }