Hello! I am currently trying to run a simulation of flow around a glacier. I have an image in paraview for both cases for reference below. My issue is that there is no flow going around the glacier and I am confused as to why that is when it worked fine for the cylinder. I have my code below. I would appreciate the help!
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import fenics as fe
import numpy as np
# Constant Variables
T = 86400 # final time
dt = 1000 # time step size
num_steps = int(T/dt) # number of time steps
mu = 0.001 # dynamic viscosity
rho = 1 # density
# Functions
def epsilon(u): # strain rate
return sym(nabla_grad(u))
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(len(u))
# Boundary Conditions Classes
class Noslip(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
# Sub domain for inflow (right)
class Inflow(SubDomain):
def inside(self, x, on_boundary):
return x[0] > 1.0 - DOLFIN_EPS and on_boundary
# Sub domain for outflow (left)
class Outflow(SubDomain):
def inside(self, x, on_boundary):
return x[0] < DOLFIN_EPS and on_boundary
# Create Mesh
mesh = Mesh('larsen.xml')
plot(mesh)
plt.show()
# Function Space
V = VectorFunctionSpace(mesh, 'CG', 2)
Q = FunctionSpace(mesh, 'CG', 1)
# Boundary Conditions
inflow_profile = ('4.0*1500*x[1]*(735000 - x[1]) / pow(735000, 2)', '0')
bc = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
bc.set_all(0)
Inflow().mark(bc, 1)
Outflow().mark(bc, 2)
Noslip().mark(bc, 3)
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), Inflow())
bcp_outflow = DirichletBC(Q, Constant(0), Outflow())
bcu_noslip = DirichletBC(V, Constant((0, 0)), Noslip())
bcu = [bcu_inflow, bcu_noslip]
bcp = [bcp_outflow]
ds = fe.ds(subdomain_data=bc)
# Trial/Test Functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
# Time-Step Function Solutions
u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
# Variational Form
U = 0.5*(u_n + u)
n = FacetNormal(mesh)
I = Identity(mesh.geometry().dim())
f = Constant((0, 0))
k = Constant(dt)
rho = Constant(rho)
mu = Constant(mu)
nu = Constant(0.001)
# 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)
# 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
# Step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx
# Assemble Meshes
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
# Apply Boundary Conditions to Meshes
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]
# Create XDMF files for Visualization Output
xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')
xdmffile_u.parameters["flush_output"] = True
xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')
xdmffile_p.parameters["flush_output"] = True
# Create progress bar
progress = Progress('Time-stepping', num_steps)
# Solve Problem w/ Time-stepping
data = []
t = 0
for j 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, 'bicgstab', 'hypre_amg')
# Step 2: Pressure correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')
# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(), b3, 'cg', 'sor')
# Save Solution to File (XDMF/HDF5)
xdmffile_u.write(u_, t)
xdmffile_p.write(p_, t)
# Update previous solution
u_n.assign(u_)
p_n.assign(p_)
# Update progress bar
set_log_level(LogLevel.PROGRESS)
progress += 1
set_log_level(LogLevel.ERROR)
