# 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.

[1]:

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

[2]:

tau = 0.0003
tend = 1
omega = 2*pi

[3]:

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:

[4]:

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

[5]:

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:

[6]:

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

ut,vt = fest.TnT()

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

gfut = GridFunction(fest)

meshVelocity = (deformation-defold) / tau

ft = LinearForm(fest)
ft += -InnerProduct(meshVelocity,n)*vt*ds(definedon="hole")

contactt = ContactBoundary(mesh.Boundaries("gammai"), mesh.Boundaries("gammao"), volume=True)

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

[7]:

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:

[8]:

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))

a = BilinearForm(fes)
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)

[9]:

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

uhat.Trace().Other()*(vhat.Trace()))


[10]:

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

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

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

[ ]: