So given this formulation, you have that the peak velocity should be in the center, thus
\frac{p}{x}=\beta, u_x =\frac{\beta}{8}H^2
Why would you want to have an inlet velocity of 1 if your channel is only 0.001 long? This will likely cause issues with the CFL condition, which is C=\frac{u\Delta t}{\Delta x}=\frac{\beta \Delta t H N}{8}<C_{max}, Here I have used that \Delta x=\frac{H}{N} where N are the number of elements in each direction.
If u_x=1, C=\frac{\Delta t N}{H}, thus you need a exteremly small time step to ensure stability (Say you have 20 elements in each direction, then \delta t\simeq \frac{H}{N}=10^{-4})
Apologies for not saying that: the desired velocity is 0.05, so I will need to change the length for that. In the tutorial, the maximum velocity is 1 and the length is 1. You mention inlet velocity, but I thought that is the maximum velocity, not the inlet, which is reached later in the simulation and further down the channel. And still you can see that for a smaller height, \beta is getting very large…
You are right, inlet/maximum velocity got mixed in my previous post.
If u=0.05, then \beta = \frac{8\cdot 0.05}{H^2}=400000. Which is a relatively large number, but I am not sure why you would expect a smaller number? You have an maximum velocity that is 50 times as large as the length of the channel (Compared to the scaled version of the problem where maximum velocity and length of the channel has ratio 1).
β here is the pressure difference, following the book. The length of the channel does not appear at all in the calculation of β, which is another thing I don’t understand. So is there any point increasing the length?
The problem in the book starts as a scaled version, but then it has dimensions, so I do not really understand if I should consider it scaled or not.
By using that β, the maximum velocity alternates (at the beginning) between 200 and 0.0002, which shows that something is wrong.
The length does go directly in to the formulation, but the mesh size does (as I showed above with the CFL condition). Might I ask what time step (\Delta t) you are using when you are trying to solve the problem with that given \beta? I tried to explain how \beta, u_x and H^2 are related in the post above.
As the tutorial goes through how to do the argumentation for the scaled problem, you can follow the same steps to do an analysis of the unscaled problem. I would suggest you sit down and write out the mathematical steps to make sure you understand them.
I analysed the steps of the code very carefully and the only thing I do not understand is the jump between scaled and unscaled version. In the book, the whole discussion is about the scaled problem and suddenly, with a choice of every parameter as 1, Re is gone and the code is written for the unscaled problem. This is very confusing and I cannot tell if it consistent or not with anything.
Thank you for the explanation about β, good to confirm we both get the same formula.
Maybe I need to explain the problem I am trying to solve more detailed:
I want solve the flow in a channel where the height will be way smaller than the length (so not a unit square). I am used to non-dimensionalise the equations, but I have no problem using dimensions therefore my parameters are μ=0.001, ρ=1000, Η=10^{-6}, L can vary but it will definitely be an order of magnitude larger than H and I want a maximum velocity of 0.05 (all these values in m, s, kg, N,…).
By trying to just implement these parameters in the code given from the tutorials after I change the rectangle mesh and add a subdomain (which has no effect on the flow at this point so please ignore it), I get all these problems mention above.
So my main questions are:
- Should I carry on with dimensions or formulate the actual non-dim equations? Because in this tutorial, the system is dimensional and just use everything as 1 to make in non-dim(?).
- Now that I have the correct pressure difference to achieve this velocity, should I change the length to take care of CLF condition?
To answer your question, Dt=0.02 (and ρ=μ=1 for this discussion but that needs to change).
Apologies for the long message and sorry if I come across as short, difficult via texts 
Thank you for your time again! 
1.I would continue with the dimensional problem. But have you tried to follow the steps to scale the equation? I think giving that a go would help you quite a lot. If you find this particular problem, tricky, I would suggest starting with the scaling section 2.4.1 and ensure that you understand the steps.
2. Look at the 2 dimensional CFL condition in: Courant–Friedrichs–Lewy condition - Wikipedia
To make sure your scheme is stable (for an explicit time stepping scheme), you need that u\Delta t/\Delta x \leq 1 where u is the maximum velocity, \Delta x the mesh size in x direction, \Delta t the time step. The size of the element is currently dependent on the height (as I assume you do not have elements that are very stretched in the x-direction.) This means that \Delta x = H/N where N are the number of elements in a NxN unit square mesh. Inserting this into the 1D CFL condition gives:
\frac{u \Delta t N}{H}\leq 1. Thus if you use a small H, you have to keep u\Delta t N of the same order as H to ensure stability. In your example above, you had H=10^{-3}, u=0.05, \Delta t=0.02, which gives \frac{0.05\cdot 0.02 \cdot N}{10^{-3}}\leq 1, N \leq 1, which is a very coarse discretization of our domain (1x1 unit square). Therefore you should reduce the \Delta t, such that you can use a mesh with a finer discretization than 1x1.
If you would have used H=10^{-6}, you would get N\leq \frac{1}{1000} which is not a valid discretization.
1a. Maybe I am confused with the term “scaled”. What exactly do you mean? By scaling, I assume non-dimensionalise, which I understand for the current problem, therefore I cannot see why we cannot write the system with the Reynolds number term.
1b. So do you mean just scale the lengths?
2. To create the mesh (as I said in a different post), I use the generate_mesh function as
# Define geometries
domain = Rectangle(Point(0,0),Point(L,H))
#circle = Circle(Point(-0.5,0.5),r)
circles = [[r,x0, y0]]; circles = np.array(circles)
for circle in circles:
circleD = Circle(Point(circle[1:]), circle[0], 60)
# Set subdomains
domain.set_subdomain(1,circleD)
# Create mesh
mesh=generate_mesh(domain,60)
So N=60, but does that mean that I have 60 points across H and L? As I said in the post above, when L is larger than H, there are less and less points between 0 and H, and I can end up with points only on the boundaries. Also, I can’t tell if the mesh is uniform or not. I definitely not have a unit square since L will be larger than H.
So using CFL for N=60, u=0.05 and H=0.001 I get a \Delta=3\times 10^{-4} which I think is very small? And I tried that too (just have increasing the number of time steps) with no luck, meaning the maximum velocity (printed at every time step) is not stable (alternates between 0.2 and 9e-5).
\rho and \mu has been set to one in this example to illustrate the process. If you want to use \rho \neq 1, \mu\neq 1, you need to include the Reynolds number in the analysis that is performed to determine \beta. Similarly, if u\neq 1, H\neq 1, it should be kept in the following analysis.
As mshr lacks documentation (so it is unclear what resolution 60 means) and is no longer maintained, I suggest using some other software (as indicated in the other post) to generate a mesh where you can control the mesh resolution.
You claim that this is very small, but as I’ve said multiple times, this is related to the CFL condition, which you cannot escape.
It would also be useful if you actually look at the flow solution at each time step, does it look physically correct?
1a. I thought that the tutorial presented in the book, although it jumps from scaled to non scaled, has the code for the non-scaled problem, therefore \rho and \mu can change without including the Reynolds number. So given the code from the book, since I want to use different parameters and nothing as 1, should I reformulate the code to include Reynolds number or just solve the system with the the parameters with dimensions? That changes a lot the whole discussion which was based on that code…
2. That is understandable. By saying very small, comes from experience, where I have modelled micrometre droplets, however all my equations and models were non-dimensional so maybe I judge from that.
2. Given that the maximum velocity jumps between these two extreme values, I think there is a stability issue there, which I can’t explain given that I fulfil the CFL condition. I am trying to print the flow field like
p=plot(u_) plt.colorbar(p) plt.show()
but maybe of the L and H difference, I can’t see anything and I don’t know to fix that since I am not familiar with that plotting method. By decreasing the length, while I get the same jumping, the flow seems okay (sorry for the ugly figure 
)
The scaled equation (on its strong form) only has one parameter, namely the Reynolds number. The code has been written on non-scaled form, including \mu, \rho and the actual physical dimensions of the problem, but set such that Re=1.
If you change \rho and fixing all other parameters, you are implicitly changing the Reynolds number. You do not need to reformulate the equations.
Instead of using matplotlib to visualize the velocity field, I would suggest saving it to either a pvd file using File or and xdmf-file using XDMFFile, so that you can look at the evolution of the properties, instead of snap-shots. Then you will be able to see how the flow varies in time, and see if it converges (How does it vary when it goes from one solution to the other).
As the snapshots you present does not have a time stamp attached to them, it is impossible to say if it jumps between two solutions or not.
I print the solution at every time step. If it helps, I can print them again with a time-stamp if you want. These two figures are consecutive time steps and you can see the change in the colour bar. Since I satisfy the CFL, why there is this instability? The flow seems okay to me.
So what happens at the next time steps? 0.0019 etc? The flow has changed direction between the two time steps? Do you have a variable inlet condition?
It keeps swapping between these values (next time step is not 0.0019). The code is the same as the book and I have not changed the initial condition but just the value Dp.
Could you create a minimal code that reproduces this? Are you able to reproduce it using a RectangleMesh (built in mesh from dolfin)?
import dolfin
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
import matplotlib
from matplotlib import rc # use LaTeX text
# Plotting parameters
matplotlib.rcParams['contour.negative_linestyle'] = 'solid'
matplotlib.rcParams['text.usetex'] = True
rc('text', usetex=True)
rc('font', family = 'serif')
# Problem parameters
T = 10.0 # final time
num_steps = int(3.0e4) # number of time steps
dt = T / num_steps # time step size
mu = 1 # dynamic viscosity [N s/m2]
rho = 1 # density [kg/m3]
L = 1e-2 # channel length [m]
H = 1e-3 # channel height [m]
r = 0.25e-3 # cell radius [m]
x0,y0 = 0.5e-3,0.5e-3 # cell centre [m]
Dp = 4e+3 # pressure difference
# Define geometries
domain = Rectangle(Point(0,0),Point(L,H))
circles = [[r,x0, y0]]; circles = np.array(circles)
for circle in circles:
circleD = Circle(Point(circle[1:]), circle[0],60)
# Set subdomains
domain.set_subdomain(1,circleD)
# Create mesh
mesh=generate_mesh(domain,60)
#mesh=RectangleMesh(Point(0,0),Point(L,H),16,16)
coor = mesh.coordinates()
xcoor, ycoor = coor[:,0], coor[:,1]
# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
# Define boundaries
inflow = 'near(x[0], 0)'
class Outflow(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],float(L))
outflow = Outflow()
class Walls(SubDomain):
def inside(self, x, on_boundary):
return near(x[1],0) or near(x[1],float(H))
walls = Walls()
# Define boundary conditions
bcu_noslip = DirichletBC(V, Constant((0, 0)), walls)
bcp_inflow = DirichletBC(Q, Constant(Dp), 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)
# Define strain-rate tensor
def epsilon(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(len(u))
# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx + \
rho*dot(dot(u_n, 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 \
- dot(f, v)*dx
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]
# Time-stepping
t = 0
xc_new,yc_new = np.zeros(len(xc)),np.zeros(len(yc))
plot_time=0
for n in range(num_steps):
# Update current time
t += dt
# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
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)
# Plot solution
print('t = %.5f' % t)
print('max u:', u_.vector().get_local().max())
print('max x:', np.max(xc))
print('L = %.3f' % L )
# Update previous solution
u_n.assign(u_)
p_n.assign(p_)
plt.figure()
p=plot(u_)
plt.colorbar(p)
plt.title(r'Time =$%.5f$' % t)
# plt.show()
plot_time=plot_time+1
filename="PNG/domain_%03d.png"% plot_time
plt.savefig(filename,dpi=400)
plt.close()
This is for the generate_mesh case.
For the RectangleMesh remove subdomain command domain.set_subdomain(1,circleD) and uncomment mesh=RectangleMesh(Point(0,0),Point(L,H),16,16).
Update: apologies, remove unnecessary lines as mentioned.
Please reduce the code to only the lines required to reproduce the issue. That means removing everything that is not required, for instance the computation of “cell points”, all the plotting and marking of a circle and tracking of boundaries.
I would like to help you, but I do not want to spend time removing unused lines of code and code that does not to relate to the issue.
1 Like
See updated comment above 
The issue is that you are initializing a time dependent flow with a huge pressure gradient, (with initial state as no flow). If you let the simulation run for more iterations, it will eventually stabilize and converge to the known solution.
You can get similar stability issues by starting with a non-moving fluid, and initializing a large inlet velocity. This can be handled by either:
- Gradually increase the pressure gradient to obtain a steady flow
- Solve a Stokes-problem to get an initial guess for the velocity profile in the fluid (then solve it as a coupled problem and split the solution into a velocity and pressure variable).
I will ignore for now the \Delta p difference, as I will need to think about the flow rate in my channel.
Right now, I have the same instability issue (the alternate values) when I use a very low Dp. I try with the following values:
\mu=\rho=1, L=10^{-1},10^{-2}.10^{-3}, H=10^{-3},\Delta p=10^{-3},10^{-4},10^{-5}, down to 10^{-12} and \Delta t=0.02,0.001.
For all these values (different combinations), same instability. When I use \rho=1000 and \mu=0.001, it is absolutely fine. The velocity decreases very fast and reaches its steady state quite fast.
I don’t think it has to do with Re, since the velocity is very small. Why is this case now with such small values of \Delta p compared to before?