-
Notifications
You must be signed in to change notification settings - Fork 115
Expand file tree
/
Copy pathrun.py
More file actions
315 lines (267 loc) · 10.8 KB
/
run.py
File metadata and controls
315 lines (267 loc) · 10.8 KB
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
#!/usr/bin/env python
## \file run.py
# \brief turbulent premixed dump combustor simulation (PSI flame)
# phi=0.5, methane-air, U=40 m/s
# \version 8.1.0 "Harrier"
#
# SU2 Project Website: https://su2code.github.io
#
# The SU2 Project is maintained by the SU2 Foundation
# (http://su2foundation.org)
#
# Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md)
#
# SU2 is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License as published by the Free Software Foundation; either
# version 2.1 of the License, or (at your option) any later version.
#
# SU2 is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public
# License along with SU2. If not, see <http://www.gnu.org/licenses/>.
import sys
import pysu2
import numpy as np
# with mpi:
# $ mpirun -n 4 python run.py
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
# without mpi:
# $ python run.py
# rank = 0
# comm = 0
# reference temperature for enthalpy computations
Tref = 298.15
# flame temperature of the methane-air mixture (phi=0.5, P=5)
Tf = 1777
# laminar burning velocity of methane-air at phi=0.5, P=5
Slu = 0.232
Le = 1.0
# unburnt temperature of the methane-air mixture (phi=0.5, P=5)
Tu = 673.0
Pu = 5.0
# unburnt density at P=5
rho_u = 2.52
# unburnt thermal conductivity of methane-air (phi=0.5, P=5)
#k_u = 0.0523
# unburnt heat capacity of methane-air (phi=0.5, P=5)
cp_u = 1350.0
# ################################################################## #
# returns laminar flamespeed as function of the enthalpy deficiency #
# dh = 0 returns S_L at Tu=673K (methane-air, phi=0.5, P=5 bar) #
# ################################################################## #
def SL(val_dh):
dh = val_dh/1000.0
dh = max(-500.0, min(dh, 200.0))
S_L = max(0.05, 0.2327 + 1.177278e-3 * dh + 2.234098e-6*dh*dh + 1.56795e-9*dh*dh*dh)
#print("dH = ", dh, ", S_L =", S_L)
return S_L
# ################################################################## #
# create a function for the initial progress variable c #
# ################################################################## #
def initC(coord):
x = coord[0]
#y = coord[1]
#z = coord[2]
# location where the flame should be
flame_x = 0.012
if (x < flame_x):
C = 0.0
else:
C = 1.0
return C
# ################################################################## #
# loop over all vertices on a boundary and impose wall temperature #
# ################################################################## #
def ApplyWallTemperature(SU2Driver, marker_ids):
"""Applies an isothermal wall boundary condition on a specific wall"""
for marker_id in marker_ids:
if marker_id < 0:
continue
allCoords = SU2Driver.Coordinates()
# the top horizontal wall is 1200K, but the vertical side wall goes from 1200K down to 700K
# to match the inlet pipe temperature (which has zero heat flux BC)
for i_vertex in range(SU2Driver.GetNumberMarkerNodes(marker_id)):
i_point = SU2Driver.GetMarkerNode(marker_id, i_vertex)
coord = allCoords.Get(i_point)
Tw = 673 + 500.0 * (coord[1]-0.0125)/0.025
SU2Driver.SetMarkerCustomTemperature(marker_id, i_vertex, Tw)
# ################################################################## #
# loop over all vertices and set the species progress variable c #
# ################################################################## #
def SetInitialSpecies(SU2Driver):
allCoords = SU2Driver.Coordinates()
iSPECIESSOLVER = SU2Driver.GetSolverIndices()['SPECIES']
for iPoint in range(SU2Driver.GetNumberNodes()):
coord = allCoords.Get(iPoint)
C = initC(coord)
# now update the initial condition for the species
SU2Driver.Solution(iSPECIESSOLVER).Set(iPoint,0,C)
# ################################################################## #
# Temperature is an algebraic function of c
# ################################################################## #
def update_temperature(SU2Driver, iPoint):
# first, get the progress variable
iSPECIESSOLVER = SU2Driver.GetSolverIndices()['SPECIES']
# Note: returns a list
C = SU2Driver.Solution(iSPECIESSOLVER)(iPoint,0)
T = Tu*(1-C) + Tf*C
iFLOWSOLVER = SU2Driver.GetSolverIndices()['INC.FLOW']
# 0 1 2 3
# pressure, velocity-x, velocity-y, enthalpy
iENTH = 3
SU2Driver.Solution(iFLOWSOLVER).Set(iPoint,iENTH, cp_u*(T-Tref))
# ################################################################## #
# Source term according to Zimont
# ################################################################## #
def zimont(SU2Driver, iPoint, nDim):
iSSTSOLVER = SU2Driver.GetSolverIndices()['SST']
tke, dissipation = SU2Driver.Solution(iSSTSOLVER)(iPoint)
iSPECIESSOLVER = SU2Driver.GetSolverIndices()['SPECIES']
# get the gradient of species_0
gradc = SU2Driver.Gradient(iSPECIESSOLVER)(iPoint,0)
primindex = SU2Driver.GetPrimitiveIndices()
iDENSITY = primindex.get("DENSITY")
iMU = primindex.get("LAMINAR_VISCOSITY")
rho = SU2Driver.Primitives()(iPoint,iDENSITY)
mu = SU2Driver.Primitives()(iPoint,iMU)
nu=mu/rho
# Turbulent Flamespeed Closure with Dinkelacker correction
up = np.sqrt((2.0/3.0) * tke )
lt = (0.09**0.75) * (tke**1.5) / dissipation
Re = up*lt/nu
# assume linear relation between progress variable and heat of combustion
iTEMPERATURE = primindex.get("TEMPERATURE")
T = SU2Driver.Primitives()(iPoint,iTEMPERATURE)
C = SU2Driver.Solution(iSPECIESSOLVER)(iPoint,0)
delta_H = C * 0.0284*50.5*1.0e6
# adiabatic temperature when no heat losses:
T_adiabatic = (delta_H/cp_u) + Tu
# heat loss
dh = cp_u*(T - T_adiabatic)
#print("T_adiabatic=",T_adiabatic, ", actual T = ",T, ", heat los = dh=",dh )
# correction to SL using simple correlation for methane-air
# no heat loss is original S_L
# heat loss gives the new S_L
Ut = SL(dh) * (1.0 + (0.46/Le) * np.power(Re,0.25) * np.power(up/SL(dh),0.3) * np.power(Pu,0.2) )
norm_gradc = 0.0
for idim in range(nDim):
norm_gradc += gradc[idim]*gradc[idim]
norm_gradc = np.sqrt(norm_gradc)
Sc = rho_u * Ut * norm_gradc
return Sc
# ################################################################## #
# Get the list of solver variable names
# ################################################################## #
def getsolvar(SU2Driver):
primindex = SU2Driver.GetPrimitiveIndices()
iFLOWSOLVER = SU2Driver.GetSolverIndices()['INC.FLOW']
nVars = SU2Driver.Solution(iFLOWSOLVER).Shape()[1]
varindex = primindex.copy()
for prim in varindex.copy():
if varindex[prim] >=nVars:
del varindex[prim]
varindex = dict(sorted(varindex.items(), key=lambda item: item[1]))
return varindex
# ################################################################## #
# Main routine
# ################################################################## #
def main():
# Initialize the primal driver of SU2, this includes solver preprocessing.
try:
driver = pysu2.CSinglezoneDriver('psi.cfg', 1, comm)
except TypeError as exception:
print('A TypeError occured in pysu2.CSinglezoneDriver : ', exception)
raise
# Get the ID of the markers where the temperature is applied.
all_marker_ids = driver.GetMarkerIndices()
marker_names = ['wall_top', 'wall_side']
marker_ids = []
for name in marker_names:
marker_ids.append(all_marker_ids[name] if name in all_marker_ids else -1)
print("marker names where the temperature boundary condition will be applied:", marker_names)
print("marker IDs where the temperature boundary condition will be applied:", marker_ids)
if rank == 0:
print("\n------------------------------ Begin Solver -----------------------------")
sys.stdout.flush()
nDim = driver.GetNumberDimensions()
# index to the flow solver
# C.FLOW
# INC.FLOW
# HEAT
# FLAMELET
# SPECIES
# SA
# SST
iFLOWSOLVER = driver.GetSolverIndices()['INC.FLOW']
iSPECIESSOLVER = driver.GetSolverIndices()['SPECIES']
iSSTSOLVER = driver.GetSolverIndices()['SST']
# all the indices and the map to the names of the primitives
primindex = driver.GetPrimitiveIndices()
nElem = driver.GetNumberElements()
nVars = driver.Solution(iFLOWSOLVER).Shape()[1]
nVarsSpecies = driver.Solution(iSPECIESSOLVER).Shape()[1]
nVarsTurb = driver.Solution(iSSTSOLVER).Shape()[1]
if rank == 0:
print("Dimensions of the problem = ",nDim)
print("index of flow solver = ",iFLOWSOLVER)
print("index of turbulence solver = ",iSSTSOLVER)
print("indices of primitives=",primindex)
print("number of primitives:",len(primindex))
print("number of elements:",nElem)
print("number of flow solver variables:",nVars)
print("indices of solver variables: ", getsolvar(driver))
print("number of species solver variables:",nVarsSpecies)
print("number of turbulence solver variables:",nVarsTurb)
sys.stdout.flush()
# ### Check if we do a restart or not. ###
with open('psi.cfg') as f:
if 'RESTART_SOL= YES' in f.read():
if rank == 0:
print("restarting from file")
else:
# We can set an initial condition by calling this function:
if rank == 0:
print("Using user defined initial condition.")
SetInitialSpecies(driver)
# super important to actually push the commands.
sys.stdout.flush()
# run N iterations
N = 10
for inner_iter in range(N):
if (rank==0):
print("python iteration ", inner_iter)
# impose wall temperature
ApplyWallTemperature(driver,marker_ids)
driver.Preprocess(inner_iter)
driver.Run()
Source = driver.UserDefinedSource(iSPECIESSOLVER)
Source_h = driver.UserDefinedSource(iFLOWSOLVER)
# set the source term, per point
for i_node in range(driver.GetNumberNodes() - driver.GetNumberHaloNodes()):
# add source term:
# default TFC of Zimont: rho*Sc = rho_u * U_t * grad(c)
S = zimont(driver,i_node, nDim)
iPV = 0
Source.Set(i_node, iPV, S)
iENTH = 3
Source_h.Set(i_node, iENTH, 1.0* 0.0284*50.5*1e6*S)
# for the update of temperature, we need to update also the halo nodes
#for i_node in range(driver.GetNumberNodes()):
# # set the temperature to T = c*Tf + (1-c)*Tu
# update_temperature(driver, i_node)
driver.Postprocess()
driver.Update()
# Monitor the solver and output solution to file if required.
#driver.Monitor(inner_iter)
# Output the solution to file
driver.Output(inner_iter)
# Finalize the solver and exit cleanly.
driver.Finalize()
if __name__ == '__main__':
main()