Solving Stokes flow on Mesh with Two Materials

General
1

Hello,

I am trying to solve the Stokes equation for a flow in a channel. The mesh consists of two parts (solid representing the walls of the channel and fluid representing the region where I want to solve the Stokes flow)

image
image782×782 65.2 KB

The code provides zero values for velocity and pressure when I try to define the simulation to run in the fluid region by using dx(simulationRegion) , when I remove (simulationRegion) , the code below works but provides incorrect results because it does not distinguish between solid and fluid. How to correctly tell the model run only in the fluid region ?

Note: bcs_Dirichlet includes the interface mesh elements between the solid and fluid (internal wall elements)
Here is my code:

    for idx, material in enumerate(meshNGSolve.GetMaterials()):
        domainType = domainDetails[idx].getDomainType()
        if domainType == 'Fluid':
            simulationRegion = material
    V = VectorH1(meshNGSolve, order=order,
                 dirichlet=bcs_Dirichlet)
    Q = H1(meshNGSolve, order=order - 1)
    fes = V * Q
    u, p = fes.TrialFunction()
    v, q = fes.TestFunction()

    bodyForceCoeffFun = CoefficientFunction((bodyForce[0], bodyForce[1], bodyForce[2]))
    stokes = viscosity * InnerProduct(grad(u), grad(v)) + div(u) * q + div(v) * p - 1e-10 * p * q

    a = BilinearForm(fes, symmetric=True, check_unused=False)
    a += stokes * dx(simulationRegion)
    a.Assemble()

    # right hand side
    f = LinearForm(fes)
    f += -bodyForceCoeffFun * v * dx(simulationRegion)
    f.Assemble()

    gfu = GridFunction(fes)
    dicBCValues = self.getDicBCValues()

    values_list = [ (0.0, 0.0, 0.0)]

    gfu.components[0].Set(CoefficientFunction(values_list),definedon=meshNGSolve.Boundaries(bcs_Dirichlet))

    inv_stokes = a.mat.Inverse(freedofs=fes.FreeDofs(), inverse="sparsecholesky")
    gfu.vec.data = inv_stokes * f.vec
2

Hi,
it is hard to tell what the problem could be without an executable minimal example including the mesh. Have you tried out to define e.g.

cf = mesh.MaterialCF({ simulationRegion :  1 }, default=0)
Draw(cf, mesh)

and test if everything is zero or is it one on the fluid domain? What do you get when printing simulationRegion, print(simulationRegion)? When you set and draw the initial data, is it nonzero at the boundaries of the fluid domain?

You might face the problem of still inverting the matrix on the whole domain. You could either restrict the fes.FreeDofs() array only containing dofs on the fluid domain, or directly tell the spaces to live only on the fluid domain with

H1(mesh, order=..., definedon = "fluid-domain-name")

Best,
Michael

3

Hi Michael,

Thank you very much for responding. I took some time to put some effort on the problem. First, I uploaded my exported netgen mesh. Note that it was exported using:
netgenMesh.Save(‘netgenMesh_channelSimulation_twoMaterials.vol’)

netgenMesh_channelSimulation_twoMaterials.vol (203.2 KB)

After defining this I get the correct output of 1 in the fluid region and 0 in the solid region:

cf = meshNGSolve.MaterialCF({ materialFluid : 1, materialSolid : 0 })
p1 = dictOfMaxMinCoordinates[‘minXcoordinates’] + 0.00051
p2 = dictOfMaxMinCoordinates[‘minYcoordinates’] + 0.00051
p3 = dictOfMaxMinCoordinates[‘minZcoordinates’]
p4 = dictOfMaxMinCoordinates[‘minXcoordinates’]
p5 = dictOfMaxMinCoordinates[‘minYcoordinates’]
p6 = dictOfMaxMinCoordinates[‘minZcoordinates’]
cf(meshNGSolve(p1, p2, p3))
cf(meshNGSolve(p4, p5, p6))

print(simulationRegion) gives ‘Fluid’.

Also, by directly telling the spaces to live only on the fluid domain I get nan values in the grid function for velocity and pressure.
print(gfu.components[1](meshNGSolve(p1, p2, p3)))
nan
print(gfu.components[0](meshNGSolve(p1, p2, p3)))
(nan, nan, nan)
On the other hand, if I do not do that, I get zero values as expected.
print(gfu.components[1](meshNGSolve(p1, p2, p3)))
0.0
print(gfu.components[0](meshNGSolve(p1, p2, p3)))
(0.0, 0.0, 0.0)

Finally, when I try to import the .vol netgen mesh, it does not have any more the same attributes (face descriptors, material names, etc …). Even when asking for coordinates I am getting:
netgenMesh = ImportMesh(‘netgenMesh_channelSimulation_twoMaterials.vol.vol’)
meshCoordinates = netgenMesh.Coordinates()
ValueError: PyMemoryView_FromBuffer(): info->buf must not be NULL

Just incase I am also providing another exported format :

netgenMesh_channelSimulation_twoMaterials.gmsh (236.5 KB)

To be more clear am putting my test script which could also show the steps that I am following which provided the nan and zero values described above.
Pressure induced flow 3D channel_twoMaterials.py (16.5 KB)

A final question, how to restrict the fes.FreeDofs() array only containing dofs on the fluid domain ?
Thank you again for your help.

4

Hi,

the definedon-flag of the spaces was not enough.
Preparing the appropriate freedof bitarray worked:

freedofs = fes.FreeDofs() & fes.GetDofs(meshNGSolve.Materials(materialFluid))
inv = a.mat.Inverse(freedofs=freedofs, inverse="sparsecholesky")

I changed the following things (specifying Dirichlet bc, Compress the space to get rid of unused dofs, and set again the materialFluid in the integrators)

V = Compress(VectorH1(meshNGSolve, order=2, definedon = materialFluid, dirichlet='InternalBC'))
Q = Compress(H1(meshNGSolve, order=1, definedon = materialFluid, dirichlet='InternalBC'))

a += stokes * dx(materialFluid)
f += -imposed_pressure_gradient * v * dx(materialFluid)

Best,
Michael

5

Thank you Michael. That worked as you suggested. Have a great day.