-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathBoussenesq.py
139 lines (111 loc) · 3.77 KB
/
Boussenesq.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
from dolfin import *
import mshr
# Define domain
center = Point(0.2, 0.2)
radius = 0.05
L = 2.2
W = 0.41
geometry = mshr.Rectangle(Point(0.0, 0.0), Point(L, W)) \
- mshr.Circle(center, radius, 10)
# Build mesh
mesh = mshr.generate_mesh(geometry, 50)
# Construct facet markers
bndry = FacetFunction("size_t", mesh)
for f in facets(mesh):
mp = f.midpoint()
if near(mp[0], 0.0): # inflow
bndry[f] = 1
elif near(mp[0], L): # outflow
bndry[f] = 2
elif near(mp[1], 0.0) or near(mp[1], W): # walls
bndry[f] = 3
elif mp.distance(center) <= radius: # cylinder
bndry[f] = 5
# Build function spaces (Taylor-Hood)
V = VectorFunctionSpace(mesh, "CG", 2)
P = FunctionSpace(mesh, "CG", 1)
E = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([V, P, E])
# No-slip boundary condition for velocity on walls and cylinder - boundary id 3
noslip = Constant((0, 0))
bcv_walls = DirichletBC(W.sub(0), noslip, bndry, 3)
vc= Expression(("-0.5*t*cos(atan2(x[0]-0.2,x[1]-0.2))","0.5*t*sin(atan2(x[0]-0.2,x[1]-0.2))"),t=0)
bcv_cylinder = DirichletBC(W.sub(0), vc, bndry, 5)
bce_cylinder = DirichletBC(W.sub(2), Constant(1.0), bndry, 5)
# Inflow boundary condition for velocity - boundary id 1
v_in = Expression(("1.5 * 4.0 * x[1] * (0.41 - x[1]) / ( 0.41 * 0.41 )", "0.0"))
bcv_in = DirichletBC(W.sub(0), v_in, bndry, 1)
# Collect boundary conditions
bcs = [bcv_cylinder, bcv_walls, bcv_in, bce_cylinder]
# Facet normal, identity tensor and boundary measure
n = FacetNormal(mesh)
I = Identity(mesh.geometry().dim())
ds = Measure("ds", subdomain_data=bndry)
nu = Constant(0.001)
dt = 0.1
t_end = 10
theta=0.5 # Crank-Nicholson timestepping
k=0.01
# Define unknown and test function(s)
(v_, p_, e_) = TestFunctions(W)
# current unknown time step
w = Function(W)
(v, p, e) = split(w)
# previous known time step
w0 = Function(W)
(v0, p0, e0) = split(w0)
def a(v,u) :
D = sym(grad(v))
return (inner(grad(v)*v, u) + inner(2*nu*D, grad(u)))*dx
def b(q,v) :
return inner(div(v),q)*dx
def c(v,e,g) :
return ( inner(k*grad(e),grad(g)) + inner(v,grad(e))*g )*dx
# Define variational forms without time derivative in previous time
F0_eq1 = a(v0,v_) + b(p,v_)
F0_eq2 = b(p_,v)
F0_eq3 = c(v0,e0,e_)
F0 = F0_eq1 + F0_eq2 + F0_eq3
# variational form without time derivative in current time
F1_eq1 = a(v,v_) + b(p,v_)
F1_eq2 = b(p_,v)
F1_eq3 = c(v,e,e_)
F1 = F1_eq1 + F1_eq2 + F1_eq3
#combine variational forms with time derivative
#
# dw/dt + F(t) = 0 is approximated as
# (w-w0)/dt + (1-theta)*F(t0) + theta*F(t) = 0
#
F = (inner((v-v0),v_)/dt + inner((e-e0),e_)/dt)*dx + (1.0-theta)*F0 + theta*F1
# Create files for storing solution
name="a"
vfile = XDMFFile(mpi_comm_world(),"results_%s/v.xdmf" % name)
pfile = XDMFFile(mpi_comm_world(),"results_%s/p.xdmf" % name)
efile = XDMFFile(mpi_comm_world(),"results_%s/e.xdmf" % name)
vfile.parameters["flush_output"] = True
pfile.parameters["flush_output"] = True
efile.parameters["flush_output"] = True
J = derivative(F, w)
problem=NonlinearVariationalProblem(F,w,bcs,J)
solver=NonlinearVariationalSolver(problem)
# Time-stepping
t = dt
while t < t_end:
print ("t =", t)
vc.t=t
# Compute
begin("Solving ....")
solver.solve()
end()
# Extract solutions:
(v, p, e) = w.split()
v.rename("v", "velocity")
p.rename("p", "pressure")
e.rename("e", "temperature")
# Save to file
vfile << v
pfile << p
efile << e
# Move to next time step
w0.assign(w)
t += dt