### How to integrate normal flux on boundary

109

views

0

The following implements a natural convection model applied to a heat-driven cavity benchmark in FEniCS 2017.2.0, and solves until an approximately steady state:

```
import fenics
""" Solution space """
mesh = fenics.UnitSquareMesh(20, 20)
P1, P2 = fenics.FiniteElement('P', mesh.ufl_cell(), 1), fenics.VectorElement('P', mesh.ufl_cell(), 2)
mixed_element = fenics.MixedElement([P1, P2, P1])
W = fenics.FunctionSpace(mesh, mixed_element)
psi_p, psi_u, psi_theta = fenics.TestFunctions(W)
w = fenics.Function(W)
p, u, theta = fenics.split(w)
""" Benchmark parameters """
mu, Pr, Ra = fenics.Constant(1.), fenics.Constant(0.71), fenics.Constant(1.e6)
ghat = fenics.Constant((0., -1.))
theta_h, theta_c = fenics.Constant(0.5), fenics.Constant(-0.5)
""" Initial values """
w_n = fenics.interpolate(
fenics.Expression(
("0.", "0.", "0.", "theta_h + x[0]*(theta_c - theta_h)"),
theta_h = theta_h, theta_c = theta_c, element = mixed_element),
W)
p_n, u_n, theta_n = fenics.split(w_n)
""" Finite difference time discretization: Backward Euler"""
Delta_t = fenics.Constant(0.001)
u_t, theta_t = (u - u_n)/Delta_t, (theta - theta_n)/Delta_t
""" Variational form of the governing equations """
inner, dot, grad, div, sym = \
fenics.inner, fenics.dot, fenics.grad, fenics.div, fenics.sym
mass = -psi_p*div(u)
f_B = Ra/Pr*theta*ghat
momentum = dot(psi_u, u_t + dot(grad(u), u) + f_B) - div(psi_u)*p \
+ 2.*mu*inner(sym(grad(psi_u)), sym(grad(u)))
energy = psi_theta*theta_t + dot(grad(psi_theta), 1./Pr*grad(theta) - theta*u)
F = (mass + momentum + energy)*fenics.dx
""" Penalty method stabilization """
gamma = fenics.Constant(1.e-7)
F += -psi_p*gamma*p*fenics.dx
""" Boundary conditions """
class HotWall(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[0], 0.)
class ColdWall(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[0], 1.)
class Walls(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary
W_u, W_theta = W.sub(1), W.sub(2)
boundary_conditions = [
fenics.DirichletBC(W_u, (0., 0.), Walls()),
fenics.DirichletBC(W_theta, theta_h, HotWall()),
fenics.DirichletBC(W_theta, theta_c, ColdWall())]
""" Solve until approximate steady state """
w.assign(w_n)
fenics.solve(F == 0., w, boundary_conditions)
for timestep in range(4):
w_n.assign(w)
Delta_t.assign(2.*Delta_t.values()[0])
fenics.solve(F == 0., w, boundary_conditions)
```

For reference, the temperature solution look like this:

`fenics.plot(theta)`

Now I wish to compute something like the heat flux on the cold wall, the difficult parts of which are captured in the following:

$f=\int_{\Gamma_c}^{ }\partial_n\theta ds=\int_{\Gamma_c}\left(\mathbf{n}\cdot\nabla\right)\theta ds=\int_{\Gamma_c}\nabla\theta\cdot\mathbf{n}ds$`ƒ` =∫_{Γc}∂_{n}`θ``d``s`=∫_{Γc}(** n**·∇)

`θ`

`d`

`s`=∫

_{Γc}∇

`θ`·

`n``d`

`s`

**How does one compute this integral?**

With minor modification from the FEniCS book, I tried

```
boundary_parts = fenics.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
cold_wall_marker = 2
ColdWall().mark(boundary_parts, cold_wall_marker)
nhat = fenics.FacetNormal(mesh)
integrate, ds, dot, grad = fenics.assemble, fenics.ds, fenics.dot, fenics.grad
cold_wall_heat_flux = integrate(dot(grad(theta), nhat)*ds(cold_wall_marker))
print(cold_wall_heat_flux)
```

which runs but prints "0.0", while I expect a non-zero value.

Thanks!

Community: FEniCS Project

written
7 weeks ago by
Miguel

### 2 Answers

2

Probably a bug in the user code along the lines of following

which prints

Thanks for the demonstration. Between your answer and looking at

which prints

-7.507112749003061

7.507112749003009

which has the qualities I expected, e.g. approximately opposite flux for the hot and cold walls.

```
from dolfin import *
mesh = UnitSquareMesh(3, 3)
ff = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
print(assemble(1*ds(subdomain_id=0, domain=mesh)))
print(assemble(1*ds(subdomain_id=0, domain=mesh, subdomain_data=ff)))
```

which prints

```
0.0
4.000000000000001
```

`help(fenics.ds)`

again, I came up with the following:```
mesh_function = fenics.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
cold_wall_id = 2
hot_wall_id = 3
ColdWall().mark(mesh_function, cold_wall_id)
HotWall().mark(mesh_function, hot_wall_id)
nhat = fenics.FacetNormal(mesh)
integrate, dot, grad = fenics.assemble, fenics.dot, fenics.grad
ds = fenics.ds(
domain = mesh, subdomain_data = mesh_function, subdomain_id = cold_wall_id)
cold_wall_heat_flux = integrate(dot(grad(theta), nhat)*ds)
print(cold_wall_heat_flux)
ds = fenics.ds(
domain = mesh, subdomain_data = mesh_function, subdomain_id = hot_wall_id)
hot_wall_heat_flux = integrate(dot(grad(theta), nhat)*ds)
print(hot_wall_heat_flux)
```

which prints

-7.507112749003061

7.507112749003009

which has the qualities I expected, e.g. approximately opposite flux for the hot and cold walls.

written
7 weeks ago by
Alexander G. Zimmerman

0

Hi,

I had a similar problem (with electric potential/current).

Besides the definition of the measure ds which Jan Blechta poitned out, there is an other problem.

The problem is (in my understanding) if you use grad on the solution you differentiate the the basis functions (going from linear to constant and therefore discontinous). This gives you some interpolation errors. You will this if you solve a poisson equation with two different Dirichlet BCs and a zero source term, and than integrate the in and ouflow. They will differ (for my problem about 20%).

I solved this by imposing the Dirichlet BCs through Lagrange Multipliers, the Lagrange Multipliers than have the meaning of the flow through that boundary. This worked for me and had far less performance impact than using higher degree elements (which also gives you better results for the flow interpolation).

Greetings

Moritz
Thanks for the warning. In the foreseeable future I will only be using this for regression testing, so thankfully I am okay with the interpolation errors in this case.

I had a similar problem (with electric potential/current).

Besides the definition of the measure ds which Jan Blechta poitned out, there is an other problem.

The problem is (in my understanding) if you use grad on the solution you differentiate the the basis functions (going from linear to constant and therefore discontinous). This gives you some interpolation errors. You will this if you solve a poisson equation with two different Dirichlet BCs and a zero source term, and than integrate the in and ouflow. They will differ (for my problem about 20%).

I solved this by imposing the Dirichlet BCs through Lagrange Multipliers, the Lagrange Multipliers than have the meaning of the flow through that boundary. This worked for me and had far less performance impact than using higher degree elements (which also gives you better results for the flow interpolation).

Greetings

Moritz

written
7 weeks ago by
Alexander G. Zimmerman

Please login to add an answer/comment or follow this question.