Calculating Volumetric Flow Rate at Internal Boundary

General
1

Hi,
I am trying to calculate volumetric flow rate at an internal boundary x=0.2. My velocity, pressure, etc. are correct, but I am getting an incorrect flow rate. Also, if I refine my mesh, the flow rate value changes.

from dolfin import *
from ufl_legacy import indices, Jacobian, Min
from mshr import *

####### Domain and mesh setup #######

# Parameters defining domain geometry:
SOLID_LEFT = 0.45
SOLID_RIGHT = 0.55
SOLID_TOP = 0.5
OMEGA_H = 0.75
OMEGA_W = 1.0
REF_MARGIN = 0.1

# Define the mshr geometrical primitives for this domain:
r_Omega = Rectangle(Point(0,0),Point(OMEGA_W,OMEGA_H))
r_Omega_s = Rectangle(Point(SOLID_LEFT,0),
                      Point(SOLID_RIGHT,SOLID_TOP))
r_Omega_int = Rectangle(Point(0,0),
                      Point(0.2,OMEGA_H))                      
# Enumerate subdomain markers
FLUID_FLAG = 0
SOLID_FLAG = 1
INT_FLAG = 5
# Zero is the default flag, and does not need to be
# explicitly set for the fluid subdomain.
r_Omega.set_subdomain(SOLID_FLAG,r_Omega_s)
r_Omega.set_subdomain(INT_FLAG,r_Omega_int)
# Parameters defining refinement level:
N = 32

# Generate mesh of Omega, which will have a fitted
# subdomain for Omega_s.
mesh = generate_mesh(r_Omega, N)


# Mesh-derived quantities:
d = mesh.geometry().dim()
n_y = FacetNormal(mesh)
I = Identity(d)

# Marking internal boundary
boundaries = MeshFunction("size_t", mesh, d-1,0)
for f in facets(mesh): # For all edges in the mesh
    p0 = Vertex(mesh, f.entities(0)[0]) # save the two ending points p0 and p1
    p1 = Vertex(mesh, f.entities(0)[1])
    if near(p0.x(0), 0.2, DOLFIN_EPS) and near(p1.x(0), 0.2, DOLFIN_EPS): # check if the points lie on the circle - if yes, put a label on this     edge
        boundaries[f] = 4

# Define surface integral
dS_int = Measure("dS", domain=mesh, subdomain_data=boundaries) 

# Calculating Flow rate
# Here v is velocity. I am sure that the velocity field is correct
Q_surf = assemble((dot((v),n_y('+'))*dS_int(4)))

I also tried to visualize ‘boundaries’. They look correct.

Screenshot (1720)
Screenshot (1720)764×467 92.7 KB

Please correct me.
Thank you!

2

Please try first with

Q_surf = assemble(v[0]*dS_int(4)))

which, on paper, should be an equivalent formula because n_y is the unit vector aligned with the x-axis.
I suspect that it will give the correct result.

3

Thank you!
Looks like it is working. But, I am not sure what’s wrong with my script.

Q_surf = assemble((dot((v),n_y('+'))*dS_int(4)))
4

My understanding is that + and - sides do not have a specific meaning (e.g., in your case, + aligned with positive x and - aligned with negative x), and the + and - orientation may change from one cell to the next. Taking your picture, it’s possible that one cell had + aligned with positive x, and the cell above aligned with negative x, ending up in subtracting the flow rate from that cell, rather than adding it.

In legacy dolfinx, the orientation of +/- was guaranteed to be consistent only when 4 is at the interface between two subdomains. If you want to use that form, make sure to divide the fluid part into two subdomains, separated by the interface 4. No such guarantee is present in dolfinx.

5

Thank you for explanation!

6

There is even more to this, as you need to add the MeshFunction for cells to this, as explained in Integrating over an interior surface - #4 by MiroK

In DOLFINx, one can reduce it to a one-sided internal integral, as explained in for instance: Add one-sided integration example as example of custom integration · Issue #158 · jorgensd/dolfinx-tutorial · GitHub or use the same notion of custom integration to consistently orient the integration entities of a dS integral.