Please how do I impose a matrix/vector as Neumann boundary condition in fenics? The examples I see are all expressions but I want to pass a vector/matrix as the boundary condition.
Thank you.
What do you mean by supplying this question. Could you please provide a mathematical formulation of this boundary condition, and your implementation for a built in mesh up to the point of where you are struggeling
The mathematical formulation
The implementation
# import dolfin
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
import sys
import pandas as pd
Time = 1.0 # final time
num_steps = 4 # number of time steps
dt = Time / num_steps # time step size
mu = 1 # kinematic viscosity
rho = 1 # density
# Create mesh and define function spaces
mesh = UnitSquareMesh(2, 2)
coordinates = mesh.coordinates()
plt.show()
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
nbc_function = Function(Q)
# Define boundaries
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 1)'
walls = 'near(x[1], 0) || near(x[1], 1)'
# Define boundary conditions for pressure
# bcu_noslip = DirichletBC(V, Constant((0, 0)), walls)
bcp_inflow = DirichletBC(Q, Constant(0.1), inflow)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
# bcu = [bcu_noslip]
bcp = [bcp_inflow, bcp_outflow]
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
# Define expressions used in variational forms
U = 0.5 * (u_n + u)
n = FacetNormal(mesh)
# f = Constant((0, 0))
k = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)
F = PETScVector()
# Define strain-rate tensor
def epsilon(u):
# return sym(nabla_grad(u))
return 0.5 * (nabla_grad(u) + nabla_grad(u).T)
# Define stress tensor
def sigma(u, p):
return 2 * mu * epsilon(u) - p * Identity(len(u))
# Define variational problem for step 1
# Set a random seed for reproducibility (optional)
np.random.seed(0)
# Create a 64x30 matrix with random components between -1 and 1
matrix = np.random.uniform(-1, 1, size=(64, 30))
# Defining the Neumann boundary condition
nbc = dot(grad(u_n) * n, v) * ds
F1 = rho * dot((u - u_) / k, v) * dx + \
rho * dot(dot(u_, nabla_grad(u_n)), v) * dx \
+ inner(sigma(U, p_n), epsilon(v)) * dx \
+ dot(p_n * n, v) * ds - dot(mu * nabla_grad(U) * n, v) * ds \
- nbc
a1 = lhs(F1)
L1 = rhs(F1)
# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q)) * dx
L2 = dot(nabla_grad(p_n), nabla_grad(q)) * dx - (1 / k) * div(u_) * q * dx
# Define variational problem for step 3
a3 = dot(u, v) * dx
L3 = dot(u_, v) * dx - k * dot(nabla_grad(p_ - p_n), v) * dx
# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
# Apply boundary conditions to matrices
# [bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]
# Initialize lists to store data
time_values = []
u_values = []
p_values = []
x_coordinates = []
y_coordinates = []
nbc_values = []
# Time-stepping
t = 0
for n in range(num_steps):
# Update current time
t += dt
# Step 1: Tentative velocity step
b1 = assemble(L1)
solve(A1, u_.vector(), b1)
# Step 2: Pressure correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2)
# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(), b3)
# Store data
time_values.append(t)
u_values.append(u_.vector().get_local()) # Store the velocity vector values
p_values.append(p_.vector().get_local()) # Store the pressure vector values
# Calculate the scalar value of nbc
nbc_scalar = assemble(nbc)
nbc_values.append(nbc_scalar.get_local()) # Store nbc scalar value
# Calculate x and y coordinates based on mesh nodes
mesh_coordinates = mesh.coordinates()
x_coordinates.append(mesh_coordinates[:, 0])
y_coordinates.append(mesh_coordinates[:, 1])
plot(u_)
plt.show()
So when I impose the neumann boundary condition as seen here it works
but when I change it to the vector it fails.
The reason I’m doing this is to be able to vary the values at the boundary to observe the flow profile for different values at the bound which is necessary for the main thing I am looking.
I hope I am a bit clear now. Thank you.
First of all, in that equation u seems to be the unknown on both sides of (1d). The only difference is the usage of n and \hat n. Could you explain what the differences between the u and the left hand side and right hand side is, as well as what the difference between the ns are?
Thank you very much for the response.
From what I know, I stand to be corrected though; the normal derivative of a vector u can be expressed as the dot product of the gradient of the vector u and the unit normal vector as expressed:
\dfrac{\partial u}{\partial n}=\nabla u\cdot \hat{n}
Though mostly this condition is specified by a function independent of the unknown variable, I don’t want to specify it a prior. But rather get it evaluated based on the behavior of the flow at the boundary during the simulation.
And so after the simulation, I would now like to provide/specify boundary values just to assess the behavior. And that is when the vector comes in.
Thank you.
This boundary condition does not make sense to me. As you are solving a fluid flow problem, the most natural condition to have on such a boundary (outflow) is either
du/dn = pn or du/dn-pn=0. See for instance chapter 3.2.1 of https://www.duo.uio.no/bitstream/handle/10852/51949/Vegard_Vinje_Thesis.pdf?sequence=1&isAllowed=y#page21
Thank you for the response.
With the response, you mentioned outflow but my intention is not imposing it on the outflow but the walls of the pipe.
Also the same boundary condition that I’m worried about can be found in the paper you recommended, specifically chapter 3.2.1 page 16:
\dfrac{\partial v}{\partial n}=n\cdot \nabla v
So please can I get more intuition or explanation to that boundary condition imposed.
Again thank you your time.
It is not a boundary condition, but a definition of the dv/dn operator.
If you have walls with a prescribed slip or noslip condition, there is no reason to add an additional boundary condition on that wall.
Thank you.
But please per what I am looking at, I cannot tell the conditions on the boundary a priori. So using that condition in the course of the simulation, the boundary condition will be evaluated.
Therefore I am not also imposing/prescribing slip or noslip boundary condition.
So please what should be the expression/function for g if
\dfrac{\partial u}{\partial n}= g
and
g\ne 0
Thank you.
If du/dn =g where g is non-zero you definitely dont have a wall, as this condition does not enforce u=0 on boundary or dot(u,n)=0 on boundary.
Applying du/dn=g Where g is known is simple with integration by parts, similar to how you would apply a Neumann condition for a Poisson problem.
Ok
In fact I appreciate your quick responses.
But please in my case the challenge is that, I don’t know g.
So please any idea on what g should be if g is not assumed to be zero (0).
Thank you.
Your problem is not well posed. You cannot have a non-defined boundary condition. It has to be related to something. The equation you claim is your boundary condition is just a way of defining the partial derivative in normal direction in a dense way. Please describe what physical problem you are trying to solve.
Thank you.
Please what I am looking at is detection of leakages in water distribution pipeline/network. And I am applying the Navier-Stokes equation. And per what I know when the neumann boundary condition \frac{\partial u}{\partial n}=0 it means that there is no flow across the boundary (i.e. the impermeability assumption). And so if the neumnn boundary condition (\frac{\partial u}{\partial n}\ne0) it implies there is a flow across the boundary.
Now how I want to detect my leakage (flow across the boundary) is based on the neumann condition (normal derivative). So I want an expression/function for the neumann condition so that depending on the value of the expression in the course of the simulation, where the value of the neumann condition is not zero, it will mean that there is a flow across the boundary.
Thank you very much once gain.
This is not correct. As I showed you in the thesis above, du/dn=0 means that water can flow out in a specific way.
Walls in Navier stokes is usually done by either using u=0 on the boundary or dot(u,n)=0 (no flow through the wall, but fluid can move along the wall.
Assuming that you have actual walls in your pipe, They should have any of these conditions. What you should then rather consider is what causes leakage. I guess it is related to the pressure at the wall (maybe the pressure gradient)?
Yes pressure has influence on the leak. And so some (if not most) techniques in literature base on the pressure difference (gradient) in detecting leakages. However changes in pressure may not be significant or obvious when the leakage is very small (what is termed background leak). And thus difference in pressure fails at some leakages. That is why I am of the view that except the upstream and the downstream, the water will flow along the walls (boundary) but not across (i.e. no penetration of pipe walls) . So that whenever the flow is not at outflow boundary but the walls and it is not along the walls but across, then it tells that there is something unusual happen (probably a leakage or an intrusion). I don’t know if I am making some sense.
In my opinion, you need to formulate the leakage as part of your boundary condition to include in the variational form.
You could for instance model the wall as porous, and use a Darcy type equation for the solid.
See for instance Coupled system, Navier-Stokes and porous medium - #2 by VegardVinje
Thank you
Will go through the recommended materials and get back to you.
I appreciate 
.
Please sorry for bringing you back to this:
I have read the materials you recommended. Really helpful, many thanks.
But please the challenge I see with the coupling is that, with my understanding, the pipe isn’t porous/permeable until a leakage occurs. So in such a case can I assume that every other portion of the boundary (pipe wall) has zero permeability except where/when there is a leakage? So that the leak area can be seen as porous so that the Darcy Law can be applied in modeling the pipe wall and the Navier-Stokes Equation for fluid in the domain, just as was done by DRØSDAL and as you suggested.
Many thanks.
That is what I would suggest. However, I am not an expert in this field, so I would strongly advice you to talk to your colleagues, supervisor or anyone that knows your domain to make sure that it can be thought of as a sensible formulation. Hopefully, there are other people on the forum that could also chime in with their opinion.
Alright I have heard you.
Thank you very much.