|
|
|
@ -0,0 +1,40 @@ |
|
|
|
from fenics import * |
|
|
|
from mshr import * |
|
|
|
import numpy as np |
|
|
|
|
|
|
|
finalt = 5 #Final time |
|
|
|
steps = 5000 #number of steps |
|
|
|
dt = finalt/steps #time steps |
|
|
|
|
|
|
|
# model parameters |
|
|
|
mu = 0.001 #dynamic viscosity |
|
|
|
rho = 1 #density |
|
|
|
|
|
|
|
#creating 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) |
|
|
|
|
|
|
|
#Defining spaces |
|
|
|
V = VectorFunctionSpace(mesh,'P', 2) |
|
|
|
Q = FunctionSpace(mesh, 'P', 1) |
|
|
|
|
|
|
|
#defining 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' |
|
|
|
|
|
|
|
#defining inflow profile |
|
|
|
inflowProfile('4.0*1.5*x[1]*(0.41-x[1])/pow(0.41,2)', '0') |
|
|
|
|
|
|
|
#defining boundary conditions: |
|
|
|
|
|
|
|
bcuInflow = DirichletBC(V,Expression(inflowProfile, degree=2), inflow) |
|
|
|
bcuWalls = DirichletBC(V, Constant((0, 0)), walls) |
|
|
|
bcuCylinder = DirichletBC(V, Constant((0,0)), cylinder) |
|
|
|
bcpOutflow = DirichletBC(Q, Constant(0), outflow) |
|
|
|
bcu = [bcuInflow, bcuWalls, bcuCylinder] |
|
|
|
bcp = [bcpOutflow] |
|
|
|
|
