3D Stokes flow, meshed in gmsh

mesh
1

pressure
pressure1014×412 13.2 KB

Hi everyone,
I got stuck trying to simulate Stokes flow against a cylinder inside a rectangular domain in 3D.
The geometry is designed in FreeCAD and meshed in GMSH (3D mesh), then i’ve converted the file with success through ubuntu in .xml; i obtained three files, geometry.xml, geometry_physical_region.xml and geometry_facet_region.xml. Simulation works but not at all levels (see attached images), it seems that doesn’t consider any element volume inside rectangle, despite the bcs are well setted and recognized from Fenics.

# Create mesh
mesh = Mesh("geom3.xml")
markers = fn.MeshFunction("size_t", mesh, 'geom3_physical_region.xml')
sub_domains = fn.MeshFunction('size_t', mesh, 'geom3_facet_region.xml')
dx = fn.Measure('dx', domain=mesh, subdomain_data=markers)

# Define function space
P2 = VectorElement('P', tetrahedron, 2)
P1 = FiniteElement('P', tetrahedron, 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)

# Define boundary conditions
bcu_walls = DirichletBC(W.sub(0), Constant((0, 0, 0)), sub_domains, 16)
bcp_inflow = DirichletBC(W.sub(1), Constant(50), sub_domains, 17)
bcp_outflow = DirichletBC(W.sub(1), Constant(0), sub_domains, 16)
bcu_cylinder = DirichletBC(W.sub(0), Constant((0, 0, 0)), sub_domains, 1)

bcu = [bcp_inflow, bcu_cylinder, bcu_walls, bcp_outflow]

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
a = inner(grad(u), grad(v))*dx - p*div(v)*dx + div(u)*q*dx
f = Constant((0, 0, 0))
L = dot(f, v)*dx

# Compute solution
w = Function(W)
solve(a == L, w, bcu, solver_parameters={'linear_solver':'mumps'})

# Split the mixed solution 
(u, p) = w.split(True)

# Save for PARAVIEW
vel = File('S_velocity.pvd')
vel << u
pre = File('S_pressure.pvd')
pre << p
![pressure|690x280](upload://v22oTJSnVdEnwfIwiAkLap6oleX.png)
![velocity|690x254](upload://vZklXqVwmjh0mMyvoPBRtlXIFZN.png)
2

Without a minimal reproducable example it is tricky to help you.

When you say it does not work at all levels, what do you imply?

You have only attached one figure,

Note that it is common to solve this problem in a symmetric fashion, see for instance: https://fenicsproject.org/olddocs/dolfin/2019.1.0/python/demos/stokes-iterative/demo_stokes-iterative.py.html#strong-formulation

3

Ok sorry, I’ll try to better explain what I mean.
I’m implying that it seems as if the flow is only captured at the initial instant and does not evolve, starting from the inlet with some initial condition on pressure or inflow velocity, and ending at the end of the domain with different condition in all the volume between this two surface nothing happen in term of flux.
Then the obstacle (cylinder) will not be in contact with the fluid flow field. How can I solve ?
The attached image represents the velocity field starting from the inlet (the side near the obstacle)

velocity
velocity885×327 17.6 KB

4

Start by solving the problem with a built in mesh (BoxMesh), to make sure your code works as it should . You can Also compare your implementation with the code in the link i supplied above.

5

I’ve did what you write, and I get the confirm that the code works well. What I want to reproduce it’s this situation but in a 3D domain (see attached image). The code used can be the same for both simulation in terms of solution of PDEs and it works well, I got the confirm that also with the code linked by you I obtain the same results. Is it a non correct use of ParaView ? should I apply some particular filter to the image ?
myplot

6

That is fine. As it is clearly something to do with how you generate your mesh. Thus, you need to supply how you generated the mesh and converted to a dolfin compatible format.

Without this information, it is impossible to tell what has gone wrong, as one cannot reproduce your results.

7

I’ve generated the mesh using GMSH, the physical group are attached in the images and also the mesh.
For the conversion I’ve follow the command dolfin-convert in ubuntu. I attached three photo which show the mesh and the geometry properties with a clipped part of mesh.

group
group768×258 12.3 KB

mesh
mesh1257×300 30.3 KB

Maybe, I can send you the file if you’re agree.
Many thanks for the patience

8

You should post your geo file here. Note that your mesh is incredibly coarse, you only have one non-Zero degree of freedom in the interior of your channel. Flow is not going to be able to pass around an object with such a coarse grid. Please adjust the mesh size

9

Sure, but the system does not allow me to insert .geo files, such that .msh or .xml
I can send all the files by mail if it is not a problem, I got an istitutional adress.

10

You can include it in plaintext using 3x` encapsulation such as

Add Gmsh geo file here

See for instance:

12

I’ve created a drive folder, where I’ve posted codes for simulation of Navier-Stokes and Stokes flow, on the previously geometry. Also all the files relative to the geometry and the mesh are presend.
Everyone can acces to this folder by clicking the following link:
https://drive.google.com/drive/folders/1_pOWo9-23KmrvvI6uviZc6Hwmsn-pCQT?usp=sharing

Any suggestions and improvements are welcome.

13

First of all; you are meshing the interior of the obstacle. If you are not solving a Fluid-structure interaction problem (FSI), you should not mesh the interior of the obstacle.
See for instance: Using the GMSH Python API to generate complex meshes | Jørgen S. Dokken
on how to use boolean operators with gmsh.

Secondly. you are using the deprecated xml format to convert your mesh to something to read with dolfin. I would suggest using meshio as follows (pip3 install --no-binary=h5py h5py meshio --user):

import meshio

# Convert mesh
mesh_from_file = meshio.read("geom_mini.msh")


def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
                           "name_to_read": [cell_data]})
    return out_mesh


triangle_mesh = create_mesh(mesh_from_file, "triangle")
meshio.write("facet_mesh.xdmf", triangle_mesh)

tetra_mesh = create_mesh(mesh_from_file, "tetra")
meshio.write("mesh.xdmf", tetra_mesh)

# Read mesh
mesh = Mesh()
mvc_v = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    infile.read(mvc_v, "name_to_read")
markers = cpp.mesh.MeshFunctionSizet(mesh, mvc_v)

mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("facet_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
sub_domains = cpp.mesh.MeshFunctionSizet(mesh, mvc)

Finally, you are not imposing any boundary conditions on the side walls. Thus you allow flow to flow to the sides of your geometry.

15

I’ve installed by terminal the line you post, it works fine. I then install h5py from python interpreter and then i get the following error:

3
31740×573 68.5 KB

16

See for instance: Problem with trying to get dolfin, meshio and GMSH python interface to coexist - #4 by BillS

Note that you can use the deprecated mesh format for now.
The actual issue with your code is the lack of boundary conditions on the side walls.

17

Dear dokken, I’ve try to implement all your suggestions but my initial problem still continue to exist.
I attach both the pressure and velocity field, as you can see the flow does not across the fluid domain, it remains on the inlet boundary.

pressure
pressure856×487 70.7 KB

velocity
velocity838×507 54.9 KB

I upload the code in this folder: FEniCS - Google Drive
Many thanks and i hope you can solve this problem.