-
Notifications
You must be signed in to change notification settings - Fork 4
/
FVConvDiff2D.py
321 lines (254 loc) · 11.8 KB
/
FVConvDiff2D.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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
import numpy as np
from scipy.sparse import spdiags
def preprocess(n, L, Pe, problem, fvscheme):
"""
Pre-process input for simulation of steady convection-diffusion heat problem
on two-dimensional square domain using the Finite Volume Method.
Pre-processing returns the convection-diffusion system matrix A(1:n^2,1:n^2)
and the source array s(1:n,1:n).
Implementation does not include source terms and is limited to the
following two sets of convective velocity fields and boundary conditions:
1) uniform flow [u,v] = [1,0], homogeneous Neumann BC. at north and south
walls (dTn/dy=dTs/dy=0), homogeneous Dirichlet BC at west wall (Tw=0),
inhomogeneous Dirichlet BC at east wall (Te=1).
2) corner flow [u,v] = [-x,y], homogeneous Neumann BC. at north wall
(dTn/dy=0), inhomogeneous Dirichlet BC at west wall (Tw=0), inhomogeneous
Dirichlet BC at east wall (Te=1), inhomogeneous Dirichlet BC
at south wall (Ts=x)
Syntax:
-------
A, s = preprocess(n, L, Pe, problem, fvscheme)
Input:
------
n : Number of cells along x,y-axis
L : Size of square in x,y-direction
Pe : Global Peclet number
problem : Problem : 1 or 2, selects case of convective field and BCs
fvscheme : Finite volume scheme: 'cds-cds' or 'uds-cds'
Output:
-------
A : Convection-diffusion system matrix, A(1:n^2,1:n^2)
s : Source array with BC contributions, s(1:n,1:n)
"""
# Coordinate arrays
dx = L/n # cell size in x,y-direction
xf = np.arange(0., L+dx, dx) # cell face coordinate vector along x,y-axis
xc = np.arange(dx/2., L, dx) # cell center coordinates along x-axis
# Generate convective velocity field at cell faces (!)
if problem == 1: # problem 1 - uniform flow
UFw = np.ones((n, n)) # western face velocity, u = 1
UFe = np.ones((n, n)) # eastern face velocity, u = 1
VFs = np.zeros((n, n)) # southern face velocity, v = 0
VFn = np.zeros((n, n)) # northern face velocity, v = 0
elif problem == 2: # problem 2 - corner flow
UFw = np.tile(-xf[0:n], (n, 1)) # western face velocity, u = -x
UFe = np.tile(-xf[1:n+1], (n, 1)) # eastern face velocity, u = -x
VFs = np.tile(xf[0:n, np.newaxis], (1, n)) # southern face velocity, v = y
VFn = np.tile(xf[1:n+1, np.newaxis], (1, n)) # northern face velocity, v = y
else: # problem unknown
print('problem not implemented')
# Generate convective face flux matrices
Fw = dx*Pe*UFw
Fe = dx*Pe*UFe
Fs = dx*Pe*VFs
Fn = dx*Pe*VFn
# Generate diffusive flux matrix
D = -np.ones((n, n)) # -dx/dy
# Generate coefficient arrays
if fvscheme.lower() == 'cds-cds': # CDS-CDS FV-scheme applied
aW = D-Fw/2.
aE = D+Fe/2.
aS = D-Fs/2.
aN = D+Fn/2.
aP = -(aW+aE+aS+aN)+Fe-Fw+Fn-Fs
elif fvscheme.lower() == 'uds-cds': # UDS-CDS FV-scheme applied
aW = D+np.minimum(0, -Fw)
aE = D+np.minimum(0, Fe)
aS = D+np.minimum(0, -Fs)
aN = D+np.minimum(0, Fn)
aP = -(aW+aE+aS+aN)+Fe-Fw+Fn-Fs
else: # fvscheme unknown
print('fvscheme not implemented')
s = np.zeros((n, n)) # initialize source array
# Impose BCs using ghost point approach
if problem == 1: # specified Tw=0, Te=1, dT/dy|s = 0, dT/dy|n = 0
Tw = np.zeros(n)
Te = np.ones(n)
Ds = np.zeros(n)
Dn = np.zeros(n)
# Correct w domain boundary points (j,i) = [:, 0], 1st order dirichlet:
aP[:, 0] = aP[:, 0]-aW[:, 0]
s[:, 0] = s[:, 0] - 2 * Tw * aW[:, 0]
aW[:, 0] = 0
# Correct e domain boundary points (j,i) = [:, n^2-1], 1st order dirichlet:
aP[:, -1] = aP[:, -1] - aE[:, -1]
s[:, -1] = s[:, -1] - 2 * Te * aE[:, -1]
aE[:, -1] = 0
# Correct s domain boundary points (j,i) = [0, :], 2nd order NeumaNn:
aP[0, :] = aP[0, :] + aS[0, :]
s[0, :] = s[0, :] + aS[0, :] * dx * Ds
aS[0, :] = 0 # Neumann: 2nd order FD
# Correct n domain boundary points (j,i) = [n^2, :], 2nd order NeumaNn:
aP[-1, :] = aP[-1, :] + aN[-1, :]
s[-1, :] = s[-1, :] - aN[-1, :] * dx * Dn
aN[-1, :] = 0
elif problem == 2: # specified Tw=0, Te=1, Ts = x, dT/dy|n = 0
Te = np.ones(n)
Tw = np.zeros(n)
Ts = xc
Dn = np.zeros(n)
# Correct w domain boundary points (j,i) = [:, 0], 1st order dirichlet:
aP[:, 0] = aP[:, 0] - aW[:, 0]
s[:, 0] = s[:, 0] - 2 * Tw * aW[:, 0]
aW[:, 0] = 0
# Correct e domain boundary points (j,i) = [:, -1], 1st order dirichlet:
aP[:, -1] = aP[:, -1] - aE[:, -1]
s[:, -1] = s[:, -1] - 2 * Te * aE[:, -1]
aE[:, -1] = 0
# Correct s domain boundary points (j,i) = [0, :], 1st order dirichlet:
aP[0, :] = aP[0, :]-aS[0, :]
s[0, :] = s[0, :] - 2 * Ts * aS[0, :]
aS[0, :] = 0
# Correct n domain boundary points (j,i) = [n^2, :], 2nd order NeumaNn:
aP[-1, :] = aP[-1, :] + aN[-1, :]
s[-1, :] = s[-1, :] - aN[-1, :] * dx * Dn
aN[-1, :] = 0
else:
print('BCs not implemented')
# Assemble sparse 5-diagonal system matrix A (of size n^2*n^2)
offsets = np.array([-n, -1, 0, 1, n])
data = np.hstack([
np.vstack([aW.reshape(n**2, 1, order='F')[n:n**2], np.zeros((n, 1))]),
np.vstack([aS.reshape(n**2, 1, order='F')[1:n**2], np.zeros((1, 1))]),
aP.reshape(n**2, 1, order='F')[0:n**2],
np.vstack([np.zeros((1, 1)), aN.reshape(n**2, 1, order='F')[0:n**2-1]]),
np.vstack([np.zeros((n, 1)), aE.reshape(n**2, 1, order='F')[0:n**2-n]])
])
A = spdiags(data.T, offsets, n**2, n**2, format="csc")
return A, s
def postprocess(T, n, L, Pe, problem, fvscheme):
"""
Post-process results from simulation of steady convection-diffusion heat
problem on two-dimensional square domain using the Finite Volume Method.
Post-processing returns the temperature field TT(1:n+2,1:n+2) extrapolated
to the walls, the global heat conservation GHC and the net convective and
diffusive fluxes through the walls CF,DF.
Implementation does not include source terms and is limited to the
following two sets of convective velocity fields and boundary conditions:
1) uniform flow [u,v] = [1,0], homogeneous Neumann BC. at north and south
walls (dTn/dy=dTs/dy=0), homogeneous Dirichlet BC at west wall (Tw=0),
inhomogeneous Dirichlet BC at east wall (Te=1).
2) corner flow [u,v] = [-x,y], homogeneous Neumann BC. at north wall
(dTn/dy=0), inhomogeneous Dirichlet BC at west wall (Tw=0), inhomogeneous
Dirichlet BC at east wall (Te=1), inhomogeneous Dirichlet BC
at south wall (Ts=x)
Syntax:
-------
TT, GHC, CF, DF = postprocess(T, n, L, Pe, problem, fvscheme)
Input:
------
T : Temperature at cell nodes, T(1:n,1:n)
n : Number of cells along x,y-axis
L : Size of square in x,y-direction
Pe : Global Peclet number
problem : Problem : 1 or 2, selects case of convective field and BCs
fvscheme : Finite volume scheme: 'cds-cds' or 'uds-cds'
Output:
-------
TT : Temperature field extrapolated to walls, TT(1:n+2,1:n+2)
GHC : Global heat conservation, scalar (computed from wall fluxes)
CF,DF : Conv. and diff. net fluxes through walls, CF=[CFw,CFe,CFs,CFn]
"""
# Coordinate arrays
dx = L/n # cell size in x,y-direction
xf = np.arange(0., L+dx, dx) # cell face coordinate vector along x,y-axis
xc = np.arange(dx/2., L, dx) # cell center coordinates along x-axis
# Generate convective velocity field at cell faces (!)
if problem == 1: # problem 1 - uniform flow
UFw = np.ones((n, n)) # western face velocity, u = 1
UFe = np.ones((n, n)) # eastern face velocity, u = 1
VFs = np.zeros((n, n)) # southern face velocity, v = 0
VFn = np.zeros((n, n)) # northern face velocity, v = 0
elif problem == 2: # problem 2 - corner flow
UFw = np.tile(-xf[0:n], (n, 1)) # western face velocity, u = -x
UFe = np.tile(-xf[1:n+1], (n, 1)) # eastern face velocity, u = -x
VFs = np.tile(xf[0:n, np.newaxis], (1, n)) # southern face velocity, v = y
VFn = np.tile(xf[1:n+1, np.newaxis], (1, n)) # northern face velocity, v = y
else: # problem unknown
print('problem not implemented')
# Generate convective face flux matrices
Fw = dx*Pe*UFw
Fe = dx*Pe*UFe
Fs = dx*Pe*VFs
Fn = dx*Pe*VFn
# Extrapolate temperature field to walls
TT = np.zeros((n+2, n+2)) # init extended temp field
TT[1:n+1, 1:n+1] = T
if problem == 1:
Tw = np.zeros(n)
Te = np.ones(n)
Ds = np.zeros(n)
Dn = np.zeros(n)
if fvscheme.lower() == 'cds-cds': # CDS-CDS FV-scheme applied
# Correct w domain boundary points (j,i) = (:,1), 1st order dirichlet:
TT[1:n+1, 0] = Tw
# Correct e domain boundary points (j,i) = (:,end), 1st order dirichlet:
TT[1:n+1, -1] = Te
# Correct s domain boundary points (j,i) = (1,:), 2nd order Neumann:
TT[0, 1:n+1] = T[0, :] - 0.5 * dx * Ds
# Correct n domain boundary points (j,i) = (end,:), 2nd order Neumann:
TT[-1, 1:n+1] = T[-1, :] + 0.5 * dx * Dn
elif fvscheme.lower() == 'uds-cds':
# Correct w domain boundary points (j,i) = (:,1), 1st order dirichlet:
TT[1:n+1, 0] = 2.0 * Tw - T[:, 0]
# Correct e domain boundary points (j,i) = (:,end), 1st order dirichlet:
TT[1:n+1, -1] = T[:, -1]
# Correct s domain boundary points (j,i) = (1,end), 2nd order Neumann:
TT[0, 1:n+1] = T[0, :] - dx * Ds
# Correct n domain boundary points (j,i) = (end,:), 2nd order Neumann:
TT[-1, 1:n+1] = T[-1, :]
elif problem == 2:
Tw = np.zeros(n)
Te = np.ones(n)
Ts = xc
Dn = np.zeros(n)
if fvscheme.lower() == 'cds-cds':
# Correct w domain boundary points (j,i) = (:,1), 1st order dirichlet:
TT[1:n+1, 0] = Tw
# Correct e domain boundary points (j,i) = (:,end), 1st order dirichlet:
TT[1:n+1, -1] = Te
# Correct s domain boundary points (j,i) = (1,:), 1st order dirichlet:
TT[0, 1:n+1] = Ts
# Correct n domain boundary points (j,i) = (end,:), 2nd order Neumann:
TT[-1, 1:n+1] = T[-1, :] + 0.5 * dx * Dn
elif fvscheme.lower() == 'uds-cds':
# Correct w domain boundary points (j,i) = (:,1), 1st order dirichlet:
TT[1:n+1, 0] = T[:, 0]
# Correct e domain boundary points (j,i) = (:,end), 1st order dirichlet:
TT[1:n+1, -1] = 2 * Te - T[:, -1]
# Correct s domain boundary points (j,i) = (1,:), 1st order dirichlet:
TT[0, 1:n+1] = 2 * Ts - T[0, :]
# Correct n domain boundary points (j,i) = (end,:), 2nd order Neumann:
TT[-1, 1:n+1] = T[-1, :]
# Compute temperature gradients at walls
if problem == 1:
DTw = 2. / dx * (T[:, 0] - Tw)
DTe = 2. / dx * (Te - T[:, -1])
DTs = np.ones(n) * Ds
DTn = np.ones(n) * Dn
elif problem == 2:
DTw = 2. / dx * (T[:, 0] - Tw)
DTe = 2. / dx * (Te - T[:, -1])
DTs = 2. / dx * (T[0, :] - Ts)
DTn = np.ones(n) * Dn
# Convective and diffusive wall fluxes
DF = dx * np.vstack([DTw, DTe, DTs, DTn])
CFw = Fw[:, 0] * TT[1:n+1, 0]
CFe = Fe[:, -1] * TT[1:n+1, -1]
CFs = Fs[0, :] * TT[0, 1:n+1]
CFn = Fn[-1, :] * TT[-1, 1:n+1]
CF = np.vstack([CFw, CFe, CFs, CFn])
# Global heat conservation
GHC = np.sum((CF[1, :] - DF[1, :]) - (CF[0, :] - DF[0, :])) \
+ np.sum((CF[3, :] - DF[3, :]) - (CF[2, :] - DF[2, :]))
return TT, GHC, CF, DF