This page was generated from wta/rotating.ipynb.

Rotating domains

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.

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw

tau = 0.0003
tend = 1
omega = 2*pi

square = MoveTo(0,0).Rectangle(1,1).Face()
circo = Circle((0.5,0.5), 0.30001).Face()
circ = Circle((0.5,0.5), 0.3).Face()
holes = Circle((0.65,0.5), 0.05).Face() + Circle((0.35,0.5), 0.05).Face()

square.edges.name="outer"
circ.edges.name="gammai"
holes.edges.name="hole"
circo.edges.name="gammao"
square.edges.name="wall"
square.edges.Min(X).name="inlet"
square.edges.Max(X).name="outlet"
outer = square-circo
outer.faces.name = "outer"

circ.faces.name = "inner"

both = Compound([outer, circ-holes])
mesh = Mesh(OCCGeometry(both, dim=2).GenerateMesh(maxh=0.05)).Curve(3)
print (mesh.GetMaterials(), mesh.GetBoundaries())
Draw (mesh);

('outer', 'inner') ('wall', 'inlet', 'outlet', 'wall', 'gammao', 'gammai', 'hole', 'hole')

Define a GridFunction deformation describing the current configuration:

fesdef = VectorH1(mesh, order=3)

deformation = GridFunction(fesdef)
defold = GridFunction(fesdef)

def MeshRotation(angle, deform=None):
    mesh.UnsetDeformation()
    if not deform: deform = GridFunction(fesdef)
    rotmat = CF( (cos(angle), -sin(angle), sin(angle), cos(angle))).Reshape( (2,2))
    center = CF( (0.5, 0.5) )
    pos = CF( (x,y) )

    deform.Set( (rotmat-Id(2))*(pos-center), definedon=mesh.Materials("inner"))
    return deform

scene = Draw(mesh)
tau1 = 1e-3
for step in range(int(tend/tau1)):
    MeshRotation(step*omega*tau1, deformation)
    mesh.SetDeformation(deformation)
    scene.Redraw()

Solve for a flow potential such that \(\frac{\partial \phi}{\partial n}\) matches the normal velocity:

fest = H1(mesh, order=3, dirichlet="inlet|outlet")

festgrad = VectorH1(mesh, order=3)
gfutgrad = GridFunction(festgrad)

ut,vt = fest.TnT()

n = specialcf.normal(2)
h = specialcf.mesh_size

gfut = GridFunction(fest)

meshVelocity = (deformation-defold) / tau

at = BilinearForm(grad(ut)*grad(vt)*dx)
ft = LinearForm(fest)
ft += -InnerProduct(meshVelocity,n)*vt*ds(definedon="hole")

contactt = ContactBoundary(mesh.Boundaries("gammai"), mesh.Boundaries("gammao"), volume=True)
contactt.AddIntegrator (3/h*(ut-ut.Other())*(vt-vt.Other()))
contactt.AddIntegrator (n*grad(ut)*(vt.Other()-vt)+n*grad(vt)*(ut.Other()-ut))


def solveWind(gfut,at,ft):
    contactt.Update (deformation, bf=at, intorder=10)

    at.Assemble()
    ft.Assemble()

    # the solution field
    gfut.Set((x), BND)
    rt = ft.vec.CreateVector()
    rt.data = ft.vec - at.mat * gfut.vec
    gfut.vec.data += at.mat.Inverse(freedofs=fest.FreeDofs(), inverse = "sparsecholesky") * rt
    gfutgrad.Set(Grad(gfut))

scene = Draw(gfut)
tau1 = 2e-3
for step in range(int(tend/tau1)):
    defold.vec.data = deformation.vec
    MeshRotation(step*omega*tau1, deformation)

    mesh.SetDeformation(deformation)
    solveWind(gfut,at,ft)
    scene.Redraw()

Solve a transport problem:

fes = L2(mesh, order=3)
u,v = fes.TnT()

feshat = FacetFESpace(mesh, order=3)
uhat, vhat = feshat.TnT()
traceop = fes.TraceOperator(feshat, average=True)

mesh.SetDeformation(MeshRotation(0))

wind = -(meshVelocity - gfutgrad)

a = BilinearForm(fes)
a += -wind*u*grad(v)*dx
uup = IfPos(wind*n, u, u.Other(bnd=0))
a += wind*n*uup*v * dx(element_boundary=True) # upwind


ahat = BilinearForm(feshat)

f = LinearForm(fes)
f.Assemble()

gfu = GridFunction(fes)
gfu.Set(exp(-10**2*((x-0.15)**2 +(y-0.5)**2)))

solveWind(gfut,at,ft)
scene = Draw(gfu, min=0, max=2, order=3, autoscale=False)

contact = ContactBoundary(mesh.Boundaries("gammao"), mesh.Boundaries("gammai"), volume=False)
nc = contact.normal

term = gfutgrad*nc * IfPos (gfutgrad*nc, uhat.Trace()*(-vhat.Trace().Other()), \
                               uhat.Trace().Other()*(vhat.Trace()))

contact.AddIntegrator (term)

t = 0
cnt = 0
deformation.vec[:] = 0
w = gfu.vec.CreateVector()
gfuhat = GridFunction(feshat)
# what = gfuhat.vec.CreateVector()

invm = fes.Mass(rho=1).Inverse()

with TaskManager():
    while t < tend:
        defold.vec.data = deformation.vec
        MeshRotation(t*omega, deformation)

        contact.Update (deformation, bf=ahat, intorder=10)
        # apply the transport operator
        mesh.SetDeformation(deformation)
        solveWind(gfut,at,ft)

        gfuhat.vec[:] = traceop * gfu.vec
        w[:] = a.Apply (gfu.vec) + traceop.T * ahat.Apply(gfuhat.vec)

        gfu.vec.data -= tau * invm * w

        if cnt%10 == 0:
            mesh.SetDeformation(deformation)
            scene.Redraw()

        t += tau
        cnt +=1