### Modified the code of Navier stokes

375

views

0

Hi everyone..

I am try to modifie the tutorial code of navier stokes . I want to move the cylinder throuhgt the channel in the time. When I run the code I dont get just one graphic, appear a lot of graphics. Someone can help me with this. the code is this:

"""

FEniCS tutorial demo program: Incompressible Navier-Stokes equations

for flow around a cylinder using the Incremental Pressure Correction

Scheme (IPCS).

u' + u . nabla(u)) - div(sigma(u, p)) = f

div(u) = 0

"""

from __future__ import print_function

from fenics import *

from mshr import *

import numpy as np

T = 5.0 # final time

num_steps = 50 # number of time steps

dt = T / num_steps # time step size

mu = 0.001 # dynamic viscosity

rho = 1 # density

# Create mesh

channel = Rectangle(Point(0, 0), Point(2.2, 0.41))

cylinder = Circle(Point(0.2, 0.2), 0.05)

domain = channel - cylinder

mesh = generate_mesh(domain, 64)

# Define function spaces

V = VectorFunctionSpace(mesh, 'P', 2)

Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries

inflow = 'near(x[0], 0)'

outflow = 'near(x[0], 2.2)'

walls = 'near(x[1], 0) || near(x[1], 0.41)'

cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'

# Define inflow profile

inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')

# Define boundary conditions

bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)

bcu_walls = DirichletBC(V, Constant((0, 0)), walls)

bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)

bcp_outflow = DirichletBC(Q, Constant(0), outflow)

bcu = [bcu_inflow, bcu_walls, bcu_cylinder]

bcp = [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 symmetric gradient

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]

# Create XDMF files for visualization output

xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')

xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')

# Create time series (for use in reaction_system.py)

timeseries_u = TimeSeries('navier_stokes_cylinder/velocity_series')

timeseries_p = TimeSeries('navier_stokes_cylinder/pressure_series')

# Save mesh to file (for use in reaction_system.py)

File('navier_stokes_cylinder/cylinder.xml.gz') << mesh

# Create progress bar

progress = Progress('Time-stepping')

set_log_level(PROGRESS)

# Time-stepping

t = 0

for n in range(num_steps):

# Update current time

t += dt

nv = 0.2 +0.1*t

channel = Rectangle(Point(0, 0), Point(2.2, 0.41))

cylinder = Circle(Point(nv, 0.2), 0.05)

domain = channel - cylinder

mesh = generate_mesh(domain, 64)

# Define function spaces

V = VectorFunctionSpace(mesh, 'P', 2)

Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries

inflow = 'near(x[0], 0)'

outflow = 'near(x[0], 2.2)'

walls = 'near(x[1], 0) || near(x[1], 0.41)'

cylinder = 'on_boundary && x[0]>0.1 && x[0]<2.1 && x[1]>0.1 && x[1]<0.3'

# Define inflow profile

inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')

# Define boundary conditions

bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)

bcu_walls = DirichletBC(V, Constant((0, 0)), walls)

bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)

bcp_outflow = DirichletBC(Q, Constant(0), outflow)

bcu = [bcu_inflow, bcu_walls, bcu_cylinder]

bcp = [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 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]

# 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')

# Plot solution

plot(u_, title='Velocity')

plot(p_, title='Pressure')

# Save solution to file (XDMF/HDF5)

xdmffile_u.write(u_, t)

xdmffile_p.write(p_, t)

# Save nodal values to file

timeseries_u.store(u_.vector(), t)

timeseries_p.store(p_.vector(), t)

# Update previous solution

u_n.assign(u_)

p_n.assign(p_)

# Update progress bar

progress.update(t / T)

print('u max:', u_.vector().array().max())

#Hold plot

interactive()
Please format properly your code and try to show us just a minimal (non)working example that reproduces your problem, not the entire code (is hard to read it).

I am try to modifie the tutorial code of navier stokes . I want to move the cylinder throuhgt the channel in the time. When I run the code I dont get just one graphic, appear a lot of graphics. Someone can help me with this. the code is this:

"""

FEniCS tutorial demo program: Incompressible Navier-Stokes equations

for flow around a cylinder using the Incremental Pressure Correction

Scheme (IPCS).

u' + u . nabla(u)) - div(sigma(u, p)) = f

div(u) = 0

"""

from __future__ import print_function

from fenics import *

from mshr import *

import numpy as np

T = 5.0 # final time

num_steps = 50 # number of time steps

dt = T / num_steps # time step size

mu = 0.001 # dynamic viscosity

rho = 1 # density

# Create mesh

channel = Rectangle(Point(0, 0), Point(2.2, 0.41))

cylinder = Circle(Point(0.2, 0.2), 0.05)

domain = channel - cylinder

mesh = generate_mesh(domain, 64)

# Define function spaces

V = VectorFunctionSpace(mesh, 'P', 2)

Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries

inflow = 'near(x[0], 0)'

outflow = 'near(x[0], 2.2)'

walls = 'near(x[1], 0) || near(x[1], 0.41)'

cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'

# Define inflow profile

inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')

# Define boundary conditions

bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)

bcu_walls = DirichletBC(V, Constant((0, 0)), walls)

bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)

bcp_outflow = DirichletBC(Q, Constant(0), outflow)

bcu = [bcu_inflow, bcu_walls, bcu_cylinder]

bcp = [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 symmetric gradient

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]

# Create XDMF files for visualization output

xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')

xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')

# Create time series (for use in reaction_system.py)

timeseries_u = TimeSeries('navier_stokes_cylinder/velocity_series')

timeseries_p = TimeSeries('navier_stokes_cylinder/pressure_series')

# Save mesh to file (for use in reaction_system.py)

File('navier_stokes_cylinder/cylinder.xml.gz') << mesh

# Create progress bar

progress = Progress('Time-stepping')

set_log_level(PROGRESS)

# Time-stepping

t = 0

for n in range(num_steps):

# Update current time

t += dt

nv = 0.2 +0.1*t

channel = Rectangle(Point(0, 0), Point(2.2, 0.41))

cylinder = Circle(Point(nv, 0.2), 0.05)

domain = channel - cylinder

mesh = generate_mesh(domain, 64)

# Define function spaces

V = VectorFunctionSpace(mesh, 'P', 2)

Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries

inflow = 'near(x[0], 0)'

outflow = 'near(x[0], 2.2)'

walls = 'near(x[1], 0) || near(x[1], 0.41)'

cylinder = 'on_boundary && x[0]>0.1 && x[0]<2.1 && x[1]>0.1 && x[1]<0.3'

# Define inflow profile

inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')

# Define boundary conditions

bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)

bcu_walls = DirichletBC(V, Constant((0, 0)), walls)

bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)

bcp_outflow = DirichletBC(Q, Constant(0), outflow)

bcu = [bcu_inflow, bcu_walls, bcu_cylinder]

bcp = [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 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]

# 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')

# Plot solution

plot(u_, title='Velocity')

plot(p_, title='Pressure')

# Save solution to file (XDMF/HDF5)

xdmffile_u.write(u_, t)

xdmffile_p.write(p_, t)

# Save nodal values to file

timeseries_u.store(u_.vector(), t)

timeseries_p.store(p_.vector(), t)

# Update previous solution

u_n.assign(u_)

p_n.assign(p_)

# Update progress bar

progress.update(t / T)

print('u max:', u_.vector().array().max())

#Hold plot

interactive()

Community: FEniCS Project

written
9 months ago by
Hernán Mella

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