#!/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 . import sys import pysu2 import numpy as np # with mpi: from mpi4py import MPI comm = MPI.COMM_WORLD rank = comm.Get_rank() # without mpi: # comm = 0 Tref = 298.15 # flame temperature of the methane-air mixture (phi=0.5, P=5) Tf = 1777 # unburnt temperature of the methane-air mixture (phi=0.5, P=5) Tu = 673.0 Pu = 5.0 phi = 0.5 # 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 # P = rho*R*T # 5 = 2.55 * R * 673 # R = 0.0029 # ################################################################## # # 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 and set the species progress variable c # # ################################################################## # def SetInitialSpecies(SU2Driver): allCoords = SU2Driver.Coordinates() iSPECIESSOLVER = SU2Driver.GetSolverIndices()['SPECIES'] for iPoint in range(SU2Driver.GetNumberNodes() - SU2Driver.GetNumberHaloNodes()): 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 #prim_indices = SU2Driver.GetPrimitiveIndices() #iTemp = prim_indices['TEMPERATURE'] #SU2Driver.Primitives().Set(iPoint,iTemp, T) iFLOWSOLVER = SU2Driver.GetSolverIndices()['INC.FLOW'] 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") # laminar burning velocity of methane-air at phi=0.5, P=5 Slu = 0.232 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 Le = 1.0 Ut = Slu * (1.0 + (0.46/Le) * np.power(Re,0.25) * np.power(up/Slu,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 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("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 for inner_iter in range(2000): if (rank==0): print("python iteration ", inner_iter) driver.Preprocess(inner_iter) driver.Run() Source = driver.UserDefinedSource(iSPECIESSOLVER) # 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) Source.Set(i_node,0,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()