diff --git a/REGGIE.md b/REGGIE.md index 589591f98..d62e19251 100644 --- a/REGGIE.md +++ b/REGGIE.md @@ -439,7 +439,7 @@ Overview of the test cases performed every week. | ** | Flow_N2_70degCone | ** | 2D axisymmetric 70 degree cone | nProcs=6 | Surface Sampling, includes CalcSurfaceImpact and adaptive wall temperature | [Link](regressioncheck/WEK_DSMC/Flow_N2_70degCone/readme.md) | | ** | fully_periodic_3D | ** | Periodic boundary conditions in all three directions | nProcs=10,20,30 | Check whether particles end up outside of the domain | [Link](regressioncheck/WEK_DSMC/fully_periodic_3D/readme.md) | | ** | Surface_Sticking_Coefficient | ** | Channel flow with a sticking coefficient model | nProcs=5 | Surface sampling | [Link](regressioncheck/WEK_DSMC/Surface_Sticking_Coefficient/readme.md) | -| 5 | Flow_N2-O2_70degCone | [BGK](regressioncheck/WEK_BGKFlow/builds.ini) | 2D axisymmetric 70 degree cone with a N2-O2 mixture | nProcs=6 | | [Link](regressioncheck/WEK_DSMC/Flow_N2-O2_70degCone/readme.md) | -| ** | Flow_N2_70degCone | ** | 2D axisymmetric 70 degree cone | nProcs=6 | | [Link](regressioncheck/WEK_DSMC/Flow_N2_70degCone/readme.md) | +| 5 | Flow_N2_70degCone | [BGK](regressioncheck/WEK_BGKFlow/builds.ini) | 2D axisymmetric 70 degree cone | nProcs=6 | | [Link](regressioncheck/WEK_DSMC/Flow_N2_70degCone/readme.md) | | ** | MultiSpec_Supersonic_Couette_Ar-He | ** | Supersonic Couette flow with an Ar-He mixture | nProcs=5 | Temperature | [Link](regressioncheck/WEK_DSMC/MultiSpec_Supersonic_Couette_Ar-He/readme.md) | +| ** | MultiSpec_Supersonic_Couette_CO2-N2 | ** | Supersonic Couette flow with a CO2-N2 mixture | nProcs=5 | Temperature | [Link](regressioncheck/WEK_DSMC/MultiSpec_Supersonic_Couette_CO2-N2/readme.md) | | 6 | Flow_N2_70degCone | [FP](regressioncheck/WEK_FPFlow/builds.ini) | 2D axisymmetric 70 degree cone | nProcs=6 | Surface Sampling, includes CalcSurfaceImpact | [Link](regressioncheck/WEK_DSMC/Flow_N2_70degCone/readme.md) | diff --git a/docs/documentation/references.bib b/docs/documentation/references.bib index 1c72dd77e..4fcca4049 100644 --- a/docs/documentation/references.bib +++ b/docs/documentation/references.bib @@ -66,6 +66,15 @@ @article{Pfeiffer2018b volume = {30}, year = {2018} } +@ARTICLE{Pfeiffer2021, +author = {Pfeiffer, Marcel and Mirza, Asim and Nizenkov, Paul}, +year = "2021", +journal = "Physics of Fluids", +volume = "33", +pages = "036106", +doi = "10.1063/5.0037915", +title = {{Multi-species modeling in the particle-based ellipsoidal statistical Bhatnagar–Gross–Krook method for monatomic gas species}}, +} @article{Jun2019, author = {Jun, Eunji and Pfeiffer, Marcel and Mieussens, Luc and Gorji, M. Hossein}, @@ -336,7 +345,7 @@ @article{Abe1994 } @phdthesis{Farbar2010, -author = {Farbar, Erin D.}, +author = {Farbar, Erin D.}, title = {Kinetic Simulation of Rarefied and Weakly Ionized Hypersonic Flow Fields.}, school = {University of Michigan, Horace H. Rackham School of Graduate Studies}, year = 2010, diff --git a/docs/documentation/userguide/features-and-models/Bhatnagar-Gross-Krook.md b/docs/documentation/userguide/features-and-models/Bhatnagar-Gross-Krook.md index 720795cfc..05833f388 100644 --- a/docs/documentation/userguide/features-and-models/Bhatnagar-Gross-Krook.md +++ b/docs/documentation/userguide/features-and-models/Bhatnagar-Gross-Krook.md @@ -2,7 +2,7 @@ # Bhatnagar-Gross-Krook Collision Operator The implementation of the BGK-based collision operator is based on the publications by {cite}`Pfeiffer2018a` and {cite}`Pfeiffer2018b`. -It allows the simulation of gas flows in the continuum and transitional regime, where the DSMC method is computationally too expensive. +It allows the simulation of gas flows in the continuum and transitional regimes, where the DSMC method is computationally too expensive. The collision integral is hereby approximated by a relaxation process: $$ \left.\frac{\partial f}{\partial t}\right|_{\mathrm{coll}} \approx \nu(f^t-f), $$ @@ -11,74 +11,59 @@ where $f^t$ is the target distribution function and $\nu$ the relaxation frequen The current implementation supports: -- 3 different methods (i.e. different target distribution functions): Ellipsoidal Statistical, Shakov, and standard BGK -- Single species, monatomic and polyatomic gases -- Multi species, monatomic and diatomic gas mixtures +- Three different BGK methods (i.e. different target distribution functions): Ellipsoidal Statistical, Shakov, and standard BGK +- Single species: monatomic, diatomic, and polyatomic gases +- Gas mixtures with an arbitrary number of monatomic, diatomic, and polyatomic species - Thermal non-equilibrium with rotational and vibrational excitation (continuous or quantized treatment) -- 2D/Axisymmetric simulations +- 2D axisymmetric simulations - Variable time step (adaption of the distribution according to the maximal relaxation factor and linear scaling) Relevant publications of the developers: - Implementation and comparison of the ESBGK, SBGK, and Unified models in PICLas for atomic species {cite}`Pfeiffer2018a` -- Extension of the modelling to diatomic species including quantized vibrational energy treatment, validation of ESBGK with the Mach -20 hypersonic flow measurements of the heat flux on a $70^\circ$ cone {cite}`Pfeiffer2018b` -- Simulation of a nozzle expansion (including the pressure chamber) with ESBGK, SBGK and coupled ESBGK-DSMC, comparison to experimental -measurements {cite}`Pfeiffer2019a`,{cite}`Pfeiffer2019b` -- Extension to polyatomic molecules, simulation of the carbon dioxide hypersonic flow around a flat-faced cylinder, comparison of -ESBGK, SBGK and DSMC regarding the shock structure and heat flux {cite}`Pfeiffer2019c` -- Implemention of Brull's multi-species modelling using Wilke's mixture rules and collision integrals for the calculation of -transport coefficients (under review) +- Extension of the modeling to diatomic species including quantized vibrational energy treatment, validation of ESBGK with the Mach 20 hypersonic flow measurements of the heat flux on a $70^\circ$ cone {cite}`Pfeiffer2018b` +- Simulation of a nozzle expansion (including the pressure chamber) with ESBGK, SBGK and coupled ESBGK-DSMC, comparison to experimental measurements {cite}`Pfeiffer2019a`,{cite}`Pfeiffer2019b` +- Extension to polyatomic molecules, simulation of the carbon dioxide hypersonic flow around a flat-faced cylinder, comparison of ESBGK, SBGK and DSMC regarding the shock structure and heat flux {cite}`Pfeiffer2019c` +- Implemention of Brull's multi-species modeling for monatomic gas mixtures using Wilke's mixture rules and collision integrals for the calculation of transport coefficients {cite}`Pfeiffer2021` +- Extension of the implementation of Brull's ESBGK multi-species model to diatomic gas mixtures using Wilke's mixture rules (under review) +- Extension of the ESBGK method to multi-species modeling of polyatomic molecules, based on the ESBGK models of Mathiaud, Mieussens, Pfeiffer, and Brull, and including internal energies with multiple vibrational degrees of freedom, using Wilke's mixture rules and collision integrals in comparison (under review) To enable the simulation with the BGK module, the respective compiler setting has to be activated: PICLAS_TIMEDISCMETHOD = BGK-Flow -A parameter file and species initialization file is required, analogous to the DSMC setup. It is recommended to utilize a previous -DSMC parameter file to ensure a complete simulation setup. To enable the simulation with the BGK methods, select the BGK method, -ES (`= 1`), Shakov (`= 2`), and standard BGK (`= 3`): +A parameter file and species initialization file is required, analogous to the DSMC setup. It is recommended to utilize a previous DSMC parameter file to ensure a complete simulation setup. To enable the simulation with the BGK methods, select the BGK method, ES (`= 1`), Shakov (`= 2`), and standard BGK (`= 3`): Particles-BGK-CollModel = 1 -The **recommended method is ESBGK**. If the simulation contains a gas mixture, a choice for the determination of the transport -coefficients is available. The first model uses Wilke's mixture rules (`= 1`) to calculate the gas mixture viscosity and thermal -conductivity. The second model utilizes collision integrals (derived for the VHS model, `= 2`) to calculate these mixture properties. -While both allow mixtures with three or more components, only the implementation of Wilke's mixing rules allows diatomic molecules. +The **recommended method is ESBGK**. If the simulation contains a gas mixture, a choice for the determination of the transport coefficients is available. The first model uses Wilke's mixture rules (`= 1`) to calculate the gas mixture viscosity and thermal conductivity. The **recommended second model utilizes collision integrals** (derived for the VHS model, `= 2`) to calculate these mixture properties. - Particles-BGK-MixtureModel = 1 + Particles-BGK-MixtureModel = 2 -The vibrational excitation can be controlled with the following flags, including the choice between continuous and quantized -vibrational energy. Quantized vibrational energy levels are currently only available for the single-species implementation. +The vibrational excitation can be controlled with the following flags, including the choice between continuous and quantized vibrational energy. Particles-BGK-DoVibRelaxation = T Particles-BGK-UseQuantVibEn = T -An octree cell refinement until the given number of particles is reached can be utilized, which corresponds to an equal refinement -in all three directions (x,y,z): +An octree cell refinement can be utilized until the given number of particles is reached, which corresponds to an equal refinement in all three directions (x,y,z): Particles-BGK-DoCellAdaptation = T Particles-BGK-MinPartsPerCell = 10 -It is recommended to utilize at least between 7 and 10 particles per (sub)cell. To enable the cell refinement above a certain number -density, the following option can be utilized +It is recommended to utilize at least between 7 and 10 particles per (sub)cell. To enable the cell refinement above a certain number density, the following option can be utilized: Particles-BGK-SplittingDens = 1E23 -A coupled BGK-DSMC simulation can be enabled, where the BGK method will be utilized if the number density $[\text{m}^{-3}]$ is -above a certain value: +A coupled BGK-DSMC simulation can be enabled, where the BGK method will be utilized if the number density $[\text{m}^{-3}]$ is above a certain value: Particles-CoupledBGKDSMC = T Particles-BGK-DSMC-SwitchDens = 1E22 -The flag `Particles-DSMC-CalcQualityFactors` controls the output of quality factors such as mean/maximal relaxation factor (mean: -average over a cell, max: maximal value within the octree), max rotational relaxation factor, which are defined as +The flag `Particles-DSMC-CalcQualityFactors` controls the output of quality factors such as mean/maximum relaxation factor (mean: average over a cell, max: maximal value within the octree), max. rotational relaxation factor, which are defined as $$ \frac{\Delta t}{\tau} < 1,$$ -where $\Delta t$ is the chosen time step and $1/\tau$ the relaxation frequency. The time step should be chosen as such that the -relaxation factors are below unity. The `BGK_DSMC_Ratio` gives the percentage of the sampled time during which the BGK model was -utilized. In a couple BGK-DSMC simulation this variable indicates the boundary between BGK and DSMC. However, a value below 1 can -occur for pure BGK simulations due to low particle numbers, when an element is skipped. +where $\Delta t$ is the chosen time step and $1/\tau$ the relaxation frequency. The time step should be chosen as such that the relaxation factors are below unity. The `BGK_DSMC_Ratio` gives the percentage of the sampled time during which the BGK model was utilized. In a coupled BGK-DSMC simulation this variable indicates the boundary between BGK and DSMC. However, a value below 1 can occur for pure BGK simulations due to low particle numbers, when an element is skipped. An option is available to utilize a moving average for the variables used in the calculation of the relaxation frequency: @@ -88,7 +73,6 @@ The purpose is to increase the sample size and reduce the noise for steady gas f Particles-BGK-MovingAverageFac = 0.01 -between zero and one must be defined with which the old $M^n$ and newly sampled moments $M$ are weighted -to define the moments for the next time step $M^{n+1}$: +between zero and one must be defined with which the old $M^n$ and newly sampled moments $M$ are weighted to define the moments for the next time step $M^{n+1}$: $$ M^{n+1}=f M+(1-f) M^n.$$ diff --git a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/70degCone2D_Set1_DSMCSurfState_000.00200000000000000_reference.h5 b/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/70degCone2D_Set1_DSMCSurfState_000.00200000000000000_reference.h5 deleted file mode 100644 index 6144f0287..000000000 Binary files a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/70degCone2D_Set1_DSMCSurfState_000.00200000000000000_reference.h5 and /dev/null differ diff --git a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/DSMC.ini b/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/DSMC.ini deleted file mode 100644 index 996695049..000000000 --- a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/DSMC.ini +++ /dev/null @@ -1,21 +0,0 @@ -! =============================================================================== ! -! Species1, N2 -! =============================================================================== ! -Part-Species1-SpeciesName = N2 -Part-Species1-InteractionID = 2 -Part-Species1-Tref = 273 -Part-Species1-dref = 4.17E-10 -Part-Species1-omega = 0.24 -Part-Species1-CharaTempVib = 3393.3 -Part-Species1-Ediss_eV = 9.79 -! =============================================================================== ! -! Species2, O2 -! =============================================================================== ! -Part-Species2-SpeciesName = O2 -Part-Species2-InteractionID = 2 -Part-Species2-Tref = 273 -Part-Species2-dref = 3.98E-10 -Part-Species2-omega = 0.24 -Part-Species2-CharaTempVib = 2272.8 -Part-Species2-Ediss_eV = 5.16 - diff --git a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/analyze.ini b/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/analyze.ini deleted file mode 100644 index 030f7f682..000000000 --- a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/analyze.ini +++ /dev/null @@ -1,7 +0,0 @@ -! hdf5 diff -h5diff_file = 70degCone2D_Set1_DSMCSurfState_000.00200000000000000.h5 -h5diff_reference_file = 70degCone2D_Set1_DSMCSurfState_000.00200000000000000_reference.h5 -h5diff_data_set = SurfaceData -h5diff_tolerance_value = 15E-2 -h5diff_tolerance_type = relative -h5diff_max_differences = 5 diff --git a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/mesh_70degCone2D_Set1_noWake_mesh.h5 b/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/mesh_70degCone2D_Set1_noWake_mesh.h5 deleted file mode 100755 index 021626708..000000000 Binary files a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/mesh_70degCone2D_Set1_noWake_mesh.h5 and /dev/null differ diff --git a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/parameter.ini b/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/parameter.ini deleted file mode 100644 index 03ab0e44e..000000000 --- a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/parameter.ini +++ /dev/null @@ -1,156 +0,0 @@ -! =============================================================================== ! -! EQUATION (linearscalaradvection) -! =============================================================================== ! -IniExactFunc = 0 - -! =============================================================================== ! -! DISCRETIZATION -! =============================================================================== ! -N = 1 ! Polynomial degree -NAnalyze = 1 ! Number of analyze points - -! =============================================================================== ! -! MESH -! =============================================================================== ! -MeshFile = mesh_70degCone2D_Set1_noWake_mesh.h5 -useCurveds = F -TrackingMethod = triatracking -! =============================================================================== ! -! OUTPUT / VISUALIZATION -! =============================================================================== ! -ProjectName = 70degCone2D_Set1 -Logging = F -printRandomSeeds = F -IterDisplayStep = 100 -! =============================================================================== ! -! CALCULATION -! =============================================================================== ! -tend = 2.0E-3 ! End time -Analyze_dt = 5.0E-4 ! Timestep of analyze outputs -CFLscale = 0.2 ! Scaling of theoretical CFL number -c0 = 299792458. -eps = 8.8541878176E-12 -mu = 12.566370614e-7 -! =============================================================================== ! -! LOAD BALANCE -! =============================================================================== ! -DoLoadBalance = T -PartWeightLoadBalance = T -DoInitialAutoRestart = T -InitialAutoRestart-PartWeightLoadBalance = T -LoadBalanceMaxSteps = 2 ! one initial and one at 5E-4 -! =============================================================================== ! -! BOUNDARIES -! =============================================================================== ! -Part-nBounds = 5 -Part-Boundary1-SourceName = IN -Part-Boundary1-Condition = open -Part-Boundary2-SourceName = OUT -Part-Boundary2-Condition = open -Part-Boundary3-SourceName = WALL -Part-Boundary3-Condition = reflective -Part-Boundary3-WallTemp = 300. -Part-Boundary3-TransACC = 1. -Part-Boundary3-MomentumACC = 1. -Part-Boundary3-VibACC = 1. -Part-Boundary3-RotACC = 1. -Part-Boundary4-SourceName = SYMAXIS -Part-Boundary4-Condition = symmetric_axis -Part-Boundary5-SourceName = ROTSYM -Part-Boundary5-Condition = symmetric -Part-FIBGMdeltas = (/0.001,0.001,0.01/) -! =============================================================================== ! -! PARTICLES -! =============================================================================== ! -Part-maxParticleNumber = 500000 -Part-nSpecies = 2 -! =============================================================================== ! -! Species1 N2 -! =============================================================================== ! -Part-Species1-ChargeIC=0 -Part-Species1-MassIC=4.65E-26 -Part-Species1-MacroParticleFactor = 1E10 - -Part-Species1-nInits = 1 -Part-Species1-Init1-SpaceIC = cell_local -Part-Species1-Init1-velocityDistribution = maxwell_lpn -Part-Species1-Init1-VeloIC = 1502.57 -Part-Species1-Init1-VeloVecIC = (/1.,0.,0./) -Part-Species1-Init1-PartDensity = 2.775E+020 -Part-Species1-Init1-MWTemperatureIC = 13.3 -Part-Species1-Init1-TempVib = 13.3 -Part-Species1-Init1-TempRot = 13.3 - -Part-Species1-nSurfaceFluxBCs = 1 -Part-Species1-Surfaceflux1-BC = 1 -Part-Species1-Surfaceflux1-velocityDistribution = maxwell_lpn -Part-Species1-Surfaceflux1-VeloIC = 1502.57 -Part-Species1-Surfaceflux1-VeloVecIC = (/1.,0.,0./) -Part-Species1-Surfaceflux1-PartDensity = 2.775E+020 -Part-Species1-Surfaceflux1-MWTemperatureIC = 13.3 -Part-Species1-Surfaceflux1-TempVib = 13.3 -Part-Species1-Surfaceflux1-TempRot = 13.3 -! =============================================================================== ! -! Species2 O2 -! =============================================================================== ! -Part-Species2-ChargeIC=0 -Part-Species2-MassIC=5.31400E-26 -Part-Species2-MacroParticleFactor = 1E10 - -Part-Species2-nInits = 1 -Part-Species2-Init1-SpaceIC = cell_local -Part-Species2-Init1-velocityDistribution = maxwell_lpn -Part-Species2-Init1-VeloIC = 1502.57 -Part-Species2-Init1-VeloVecIC = (/1.,0.,0./) -Part-Species2-Init1-PartDensity = 0.925E+020 -Part-Species2-Init1-MWTemperatureIC = 13.3 -Part-Species2-Init1-TempVib = 13.3 -Part-Species2-Init1-TempRot = 13.3 - -Part-Species2-nSurfaceFluxBCs = 1 -Part-Species2-Surfaceflux1-BC = 1 -Part-Species2-Surfaceflux1-velocityDistribution = maxwell_lpn -Part-Species2-Surfaceflux1-VeloIC = 1502.57 -Part-Species2-Surfaceflux1-VeloVecIC = (/1.,0.,0./) -Part-Species2-Surfaceflux1-PartDensity = 0.925E+020 -Part-Species2-Surfaceflux1-MWTemperatureIC = 13.3 -Part-Species2-Surfaceflux1-TempVib = 13.3 -Part-Species2-Surfaceflux1-TempRot = 13.3 -! =============================================================================== ! -! DSMC -! =============================================================================== ! -UseDSMC = T -ManualTimeStep= 2.0000E-07 -Particles-HaloEpsVelo = 12.000E+03 -Particles-NumberForDSMCOutputs = 1 -Part-TimeFracForSampling = 0.5 -Particles-DSMC-CalcSurfaceVal = T -Particles-DSMC-CalcQualityFactors = T -Particles-DSMCReservoirSim = F -Particles-DSMC-CollisMode = 2 !(1:elast coll, 2: elast + rela, 3:chem) -Part-NumberOfRandomSeeds = 2 -Particles-RandomSeed1 = 1 -Particles-RandomSeed2 = 2 -Particles-ModelForVibrationEnergy = 0 !(0:SHO, 1:TSHO) -Particles-DSMC-UseOctree = T -Particles-DSMC-UseNearestNeighbour = T -Particles-OctreePartNumNode = 40 -Particles-OctreePartNumNodeMin = 28 -Particles-MPIWeight = 1000 -! Symmetry -Particles-Symmetry-Order = 2 -Particles-Symmetry2DAxisymmetric = T -! Radial Weighting -Particles-RadialWeighting = T -Particles-RadialWeighting-PartScaleFactor = 60 -Particles-RadialWeighting-CloneMode = 2 -Particles-RadialWeighting-CloneDelay = 5 -! BGK-Flow -Particles-BGK-CollModel = 1 -Particles-BGK-MixtureModel = 1 -Particles-BGK-DoVibRelaxation = T -Particles-BGK-UseQuantVibEn = F -! BGK Refinement -Particles-BGK-DoCellAdaptation = T -Particles-BGK-MinPartsPerCell = 12 -Particles-BGK-SplittingDens = 3.8E20 diff --git a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/readme.md b/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/readme.md deleted file mode 100644 index 7d59cfd97..000000000 --- a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/readme.md +++ /dev/null @@ -1,4 +0,0 @@ -# BGK-Flow - Hypersonic flow around a 70° Cone (Axisymmetric) -* Simulation of a hypersonic N2/O2 flow around a 70° blunted cone at M = 20 -* Test case based on Allègre, J., Bisch, D., & Lengrand, J. C. (1997). Experimental Rarefied Heat Transfer at Hypersonic Conditions over 70-Degree Blunted Cone. Journal of Spacecraft and Rockets, 34(6), 724–728. https://doi.org/10.2514/2.3302 -* Comparison of the heat flux with a reference surface state file (produced with DSMC) diff --git a/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_Ar-He/CouetteFlow_DSMCState_001.000000_CollInt_ref.h5 b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_Ar-He/CouetteFlow_DSMCState_001.000000_CollInt_ref.h5 index 6091c7653..a4322bc38 100644 Binary files a/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_Ar-He/CouetteFlow_DSMCState_001.000000_CollInt_ref.h5 and b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_Ar-He/CouetteFlow_DSMCState_001.000000_CollInt_ref.h5 differ diff --git a/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_Ar-He/CouetteFlow_DSMCState_001.000000_Wilke_ref.h5 b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_Ar-He/CouetteFlow_DSMCState_001.000000_Wilke_ref.h5 index 6aa323bb9..3cc6b5f34 100644 Binary files a/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_Ar-He/CouetteFlow_DSMCState_001.000000_Wilke_ref.h5 and b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_Ar-He/CouetteFlow_DSMCState_001.000000_Wilke_ref.h5 differ diff --git a/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/CouetteFlow_DSMCState_001.000000_ref.h5 b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/CouetteFlow_DSMCState_001.000000_ref.h5 new file mode 100644 index 000000000..7a1197320 Binary files /dev/null and b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/CouetteFlow_DSMCState_001.000000_ref.h5 differ diff --git a/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/DSMC.ini b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/DSMC.ini new file mode 100644 index 000000000..e606a8d77 --- /dev/null +++ b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/DSMC.ini @@ -0,0 +1,26 @@ +! =============================================================================== ! +! Species1, CO2 +! =============================================================================== ! +Part-Species1-SpeciesName = CO2 +Part-Species1-PolyatomicMol = true +Part-Species1-InteractionID = 2 +Part-Species1-Tref = 273 +Part-Species1-dref = 5.10E-10 +Part-Species1-omega = 0.24 +Part-Species1-NumOfAtoms = 3 +Part-Species1-LinearMolec = true +Part-Species1-CharaTempVib1 = 959.66 +Part-Species1-CharaTempVib2 = 959.66 +Part-Species1-CharaTempVib3 = 1918.6 +Part-Species1-CharaTempVib4 = 3382 +Part-Species1-Ediss_eV = 5.45 +! =============================================================================== ! +! Species2, N2 +! =============================================================================== ! +Part-Species2-SpeciesName = N2 +Part-Species2-InteractionID = 2 +Part-Species2-Tref = 273 +Part-Species2-dref = 4.17E-10 +Part-Species2-omega = 0.24 +Part-Species2-CharaTempVib = 3393.3 +Part-Species2-Ediss_eV = 9.79 diff --git a/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/analyze.ini b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/analyze.ini new file mode 100644 index 000000000..550e3909b --- /dev/null +++ b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/analyze.ini @@ -0,0 +1,8 @@ +! hdf5 diff +h5diff_file = CouetteFlow_DSMCState_001.000000.h5 +h5diff_reference_file = CouetteFlow_DSMCState_001.000000_ref.h5 +h5diff_data_set = ElemData +h5diff_tolerance_value = 4 +h5diff_tolerance_type = absolute +h5diff_max_differences = 600 ! Allowed differences for 3x100 densities and 3x100 vibrational temperatures +h5diff_one_diff_per_run = T diff --git a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/command_line.ini b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/command_line.ini similarity index 76% rename from regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/command_line.ini rename to regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/command_line.ini index 11fbc4524..a2534cbf8 100644 --- a/regressioncheck/WEK_BGKFlow/Flow_N2-O2_70degCone/command_line.ini +++ b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/command_line.ini @@ -1,2 +1,2 @@ -MPI=6 +MPI=5 cmd_suffix=DSMC.ini diff --git a/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/externals.ini b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/externals.ini new file mode 100644 index 000000000..38d2ee749 --- /dev/null +++ b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/externals.ini @@ -0,0 +1,9 @@ +! --- Externals Tool Reggie +MPI = 1 , 1 ! Single execution +externalbinary = ./bin/piclas2vtk , ./bin/hopr ! Relative binary path in build directory +externaldirectory = post-vtk-DSMC-conversion , ./pre-hopr ! Directory name, where the files are located for the external tool reggie +externalruntime = post , pre ! Run after piclas is completed (post: after, pre: before) +cmd_suffix = ../CouetteFlow_DSMCState_001.000000.h5 , ! Suffix for the binary execution +cmd_pre_execute = cp\s-r\s../pre-hopr\s. , ! "\s" resembles a white space character in the command (simply using " " is not allowed) + +nocrosscombination:MPI,externalbinary,externaldirectory,externalruntime,cmd_suffix,cmd_pre_execute diff --git a/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/parameter.ini b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/parameter.ini new file mode 100644 index 000000000..aa3a184ae --- /dev/null +++ b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/parameter.ini @@ -0,0 +1,136 @@ +! =============================================================================== ! +! EQUATION (linearscalaradvection) +! =============================================================================== ! +IniExactFunc = 0 +! =============================================================================== ! +! DISCRETIZATION +! =============================================================================== ! +N = 1 ! Polynomial degree +NAnalyze = 1 ! Number of analyze points +CFLscale = 0.2 ! Scaling of theoretical CFL number +! =============================================================================== ! +! POSTI +! =============================================================================== ! +NVisu = 1 ! Number of visualization points +TimeStampLength = 10 +! =============================================================================== ! +! MESH +! =============================================================================== ! +MeshFile = pre-hopr/tunnel_mesh.h5 +useCurveds = F +! =============================================================================== ! +! OUTPUT / VISUALIZATION +! =============================================================================== ! +ProjectName = CouetteFlow +Logging = F +WriteErrorFiles = F +! =============================================================================== ! +! CALCULATION +! =============================================================================== ! +IterDisplayStep = 1000 +Part-AnalyzeStep = 1000 + +tend = 1.0 +Analyze_dt = 0.25 + +ManualTimeStep = 1E-5 + +Particles-NumberForDSMCOutputs=1 +Part-TimeFracForSampling=0.75 +! =============================================================================== ! +! LOAD BALANCE +! =============================================================================== ! +DoLoadBalance = T +PartWeightLoadBalance = T + +! Initial load balance +DoInitialAutoRestart = T +InitialAutoRestart-PartWeightLoadBalance = T +LoadBalanceMaxSteps = 2 +Load-DeviationThreshold = 1E-9 +! =============================================================================== ! +! ESBGK +! =============================================================================== ! +Particles-BGK-CollModel = 1 ! 1: ESBGK, 2: SBGK, 3: BGK +Particles-BGK-MixtureModel = 2 +! =============================================================================== ! +! BOUNDARIES +! =============================================================================== ! +Part-nBounds=6 +Part-Boundary1-SourceName=BC_periodicx+ +Part-Boundary1-Condition=periodic +Part-Boundary2-SourceName=BC_periodicx- +Part-Boundary2-Condition=periodic +Part-Boundary3-SourceName=BC_periodicy+ +Part-Boundary3-Condition=reflective +Part-Boundary3-MomentumACC=1. +Part-Boundary3-TransACC=1. +Part-Boundary3-WallTemp=273. +Part-Boundary3-WallVelo=(/350,0.,0./) +Part-Boundary4-SourceName=BC_periodicy- +Part-Boundary4-Condition=reflective +Part-Boundary4-MomentumACC=1. +Part-Boundary4-TransACC=1. +Part-Boundary4-WallTemp=273. +Part-Boundary4-WallVelo=(/-350,0.,0./) +Part-Boundary5-SourceName=BC_periodicz+ +Part-Boundary5-Condition=periodic +Part-Boundary6-SourceName=BC_periodicz- +Part-Boundary6-Condition=periodic +Part-nPeriodicVectors=2 +Part-FIBGMdeltas=(/0.005,0.005,0.005/) +! =============================================================================== ! +! PARTICLES +! =============================================================================== ! +Part-maxParticleNumber=50000 +Part-nSpecies=2 +Part-Species1-MacroParticleFactor=1E11 +Part-Species2-MacroParticleFactor=1E11 +! =============================================================================== ! +! Species1 CO2 +! =============================================================================== ! +Part-Species1-MassIC = 7.306E-26 +Part-Species1-ChargeIC=0 + +Part-Species1-nInits=1 +Part-Species1-Init1-SpaceIC=cuboid +Part-Species1-Init1-PartDensity=6.5E19 +Part-Species1-Init1-velocityDistribution=maxwell_lpn +Part-Species1-Init1-MWTemperatureIC=273. +Part-Species1-Init1-TempVib=273. +Part-Species1-Init1-TempRot=273. +Part-Species1-Init1-BasePointIC=(/0.0,-0.5,0.0/) +Part-Species1-Init1-BaseVector1IC=(/0.005,0.0,0.0/) +Part-Species1-Init1-BaseVector2IC=(/0.,1.0,0.0/) +Part-Species1-Init1-CuboidHeightIC=0.005 +Part-Species1-Init1-VeloIC=0.0 +Part-Species1-Init1-VeloVecIC=(/1.,0.,0./) +! =============================================================================== ! +! Species2 N2 +! =============================================================================== ! +Part-Species2-ChargeIC=0 +Part-Species2-MassIC=4.65E-26 ! N2 Molecular Mass + +Part-Species2-nInits=1 +Part-Species2-Init1-SpaceIC=cuboid +Part-Species2-Init1-PartDensity=6.5E19 +Part-Species2-Init1-velocityDistribution=maxwell_lpn +Part-Species2-Init1-MWTemperatureIC=273. +Part-Species2-Init1-TempVib=273. +Part-Species2-Init1-TempRot=273. +Part-Species2-Init1-BasePointIC=(/0.0,-0.5,0.0/) +Part-Species2-Init1-BaseVector1IC=(/0.005,0.0,0.0/) +Part-Species2-Init1-BaseVector2IC=(/0.,1.0,0.0/) +Part-Species2-Init1-CuboidHeightIC=0.005 +Part-Species2-Init1-VeloIC=0.0 +Part-Species2-Init1-VeloVecIC=(/1.,0.,0./) +! =============================================================================== ! +! DSMC +! =============================================================================== ! +Particles-DSMC-CalcQualityFactors=true +UseDSMC=true +Particles-DSMC-CollisMode=2 !(1:elast coll, 2: elast + rela, 3:chem) +Part-NumberOfRandomSeeds =2 +Particles-RandomSeed1= 1 +Particles-RandomSeed2= 2 +Particles-HaloEpsVelo = 9000 diff --git a/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/post-vtk-DSMC-conversion/parameter.ini b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/post-vtk-DSMC-conversion/parameter.ini new file mode 100644 index 000000000..26b9fd923 --- /dev/null +++ b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/post-vtk-DSMC-conversion/parameter.ini @@ -0,0 +1 @@ +NVisu = 1 \ No newline at end of file diff --git a/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/pre-hopr/hopr.ini b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/pre-hopr/hopr.ini new file mode 100755 index 000000000..17388af01 --- /dev/null +++ b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/pre-hopr/hopr.ini @@ -0,0 +1,39 @@ +! =============================================================================== ! +! PREPROC +! =============================================================================== ! +projectname=tunnel +mode=1 ! 1 Cartesian 2 gambit file 3 CGNS +useCurveds=F +DebugVisu=T +!=============================================================================== ! +! MESH +!=============================================================================== ! + Mode =1 ! 1 Cartesian 2 gambit file 3 CGNS + nZones =1 ! number of zones + Corner =(/0.,-0.5,0.,,0.005,-0.5,0.,,0.005,0.5,0.,,0.,0.5,0.,,0.,-0.5,0.005,,0.005,-0.5,0.005,,0.005,0.5,0.005,,0.,0.5,0.005/) + nElems =(/1,100,1/) + BCIndex =(/5,3,2,4,1,6/) ! Indices of UserDefinedBoundaries + elemtype =108 ! Elementform (108: Hexaeder) + useCurveds =F ! T if curved boundaries defined + SpaceQuandt =1. ! characteristic length of the mesh + ConformConnect=T + +!=============================================================================== ! +! BOUNDARY CONDITIONS +!=============================================================================== ! + nUserDefinedBoundaries=6 + BoundaryName=BC_periodicx+ ! Periodic (+vv1) + BoundaryType=(/1,0,0,1/) ! Periodic (+vv1) + BoundaryName=BC_periodicx- ! Periodic (-vv1) + BoundaryType=(/1,0,0,-1/) ! Periodic (-vv1) + BoundaryName=BC_periodicy+ ! Periodic (+vv2) + BoundaryType=(/4,0,0,0/) ! Periodic (+vv2) + BoundaryName=BC_periodicy- ! Periodic (-vv2) + BoundaryType=(/4,0,0,0/) ! Periodic (-vv2) + BoundaryName=BC_periodicz+ ! Periodic (+vv3) + BoundaryType=(/1,0,0,2/) ! Periodic (+vv3) + BoundaryName=BC_periodicz- ! Periodic (-vv3) + BoundaryType=(/1,0,0,-2/) ! Periodic (-vv3) + nVV=2 ! Anzahl der Verschiebungsvektoren für periodische RB (=Anzahl periodische Ränder) + VV=(/0.005,0.,0./) ! Verschiebungsvektor 1 (x-Richtung) + VV=(/0.,0.,0.005/) ! Verschiebungsvektor 3 (z-Richtung) diff --git a/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/readme.md b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/readme.md new file mode 100644 index 000000000..21d214e4d --- /dev/null +++ b/regressioncheck/WEK_BGKFlow/MultiSpec_Supersonic_Couette_CO2-N2/readme.md @@ -0,0 +1,3 @@ +# BGK Multispecies - Supersonic Couette flow +* Simulation of supersonic Couette flow: upper (y+) and lower (y-) boundaries with a wall velocity of +/-350 m/s, periodic boundary conditions in x and z +* CO2-N2 mixture diff --git a/src/particles/bgk/bgk_adaptation.f90 b/src/particles/bgk/bgk_adaptation.f90 index 798a570be..f9d39cb2c 100644 --- a/src/particles/bgk/bgk_adaptation.f90 +++ b/src/particles/bgk/bgk_adaptation.f90 @@ -51,6 +51,7 @@ SUBROUTINE BGK_octree_adapt(iElem) USE MOD_FP_CollOperator ,ONLY: FP_CollisionOperator USE MOD_BGK_Vars ,ONLY: BGKInitDone,BGK_MeanRelaxFactor,BGK_MeanRelaxFactorCounter,BGK_MaxRelaxFactor USE MOD_BGK_Vars ,ONLY: BGK_QualityFacSamp, BGK_MaxRotRelaxFactor, BGK_PrandtlNumber, BGK_ExpectedPrandtlNumber +USE MOD_BGK_Vars ,ONLY: BGK_Viscosity, BGK_ThermalConductivity USE MOD_FPFlow_Vars ,ONLY: FPInitDone, FP_PrandtlNumber, FP_QualityFacSamp USE MOD_FPFlow_Vars ,ONLY: FP_MaxRelaxFactor, FP_MaxRotRelaxFactor, FP_MeanRelaxFactor, FP_MeanRelaxFactorCounter USE MOD_part_tools ,ONLY: GetParticleWeight @@ -75,7 +76,7 @@ SUBROUTINE BGK_octree_adapt(iElem) IF(DSMC%CalcQualityFactors) THEN IF(BGKInitDone) THEN BGK_MeanRelaxFactorCounter = 0; BGK_MeanRelaxFactor = 0.; BGK_MaxRelaxFactor = 0.; BGK_MaxRotRelaxFactor = 0. - BGK_PrandtlNumber=0.; BGK_ExpectedPrandtlNumber=0. + BGK_PrandtlNumber=0.; BGK_ExpectedPrandtlNumber=0.; BGK_Viscosity=0.; BGK_ThermalConductivity=0. END IF IF(FPInitDone) THEN FP_MeanRelaxFactorCounter = 0; FP_MeanRelaxFactor = 0.; FP_MaxRelaxFactor = 0.; FP_MaxRotRelaxFactor = 0.; FP_PrandtlNumber = 0. @@ -207,6 +208,8 @@ SUBROUTINE BGK_octree_adapt(iElem) BGK_QualityFacSamp(5,iElem) = BGK_QualityFacSamp(5,iElem) + BGK_MaxRotRelaxFactor BGK_QualityFacSamp(6,iElem) = BGK_QualityFacSamp(6,iElem) + BGK_PrandtlNumber BGK_QualityFacSamp(7,iElem) = BGK_QualityFacSamp(7,iElem) + BGK_ExpectedPrandtlNumber + BGK_QualityFacSamp(8,iElem) = BGK_QualityFacSamp(8,iElem) + BGK_Viscosity + BGK_QualityFacSamp(9,iElem) = BGK_QualityFacSamp(9,iElem) + BGK_ThermalConductivity END IF IF(FPInitDone) THEN FP_QualityFacSamp(1,iElem) = FP_QualityFacSamp(1,iElem) + FP_MeanRelaxFactor @@ -440,6 +443,7 @@ SUBROUTINE BGK_quadtree_adapt(iElem) USE MOD_FP_CollOperator ,ONLY: FP_CollisionOperator USE MOD_BGK_Vars ,ONLY: BGKInitDone,BGK_MeanRelaxFactor,BGK_MeanRelaxFactorCounter,BGK_MaxRelaxFactor USE MOD_BGK_Vars ,ONLY: BGK_QualityFacSamp, BGK_MaxRotRelaxFactor, BGK_PrandtlNumber, BGK_ExpectedPrandtlNumber +USE MOD_BGK_Vars ,ONLY: BGK_Viscosity, BGK_ThermalConductivity USE MOD_FPFlow_Vars ,ONLY: FPInitDone, FP_PrandtlNumber, FP_QualityFacSamp USE MOD_FPFlow_Vars ,ONLY: FP_MaxRelaxFactor, FP_MaxRotRelaxFactor, FP_MeanRelaxFactor, FP_MeanRelaxFactorCounter USE MOD_part_tools ,ONLY: GetParticleWeight @@ -469,7 +473,7 @@ SUBROUTINE BGK_quadtree_adapt(iElem) IF(DSMC%CalcQualityFactors) THEN IF(BGKInitDone) THEN BGK_MeanRelaxFactorCounter = 0; BGK_MeanRelaxFactor = 0.; BGK_MaxRelaxFactor = 0.; BGK_MaxRotRelaxFactor = 0. - BGK_PrandtlNumber=0.; BGK_ExpectedPrandtlNumber=0. + BGK_PrandtlNumber=0.; BGK_ExpectedPrandtlNumber=0.; BGK_Viscosity=0.; BGK_ThermalConductivity=0. END IF IF(FPInitDone) THEN FP_MeanRelaxFactorCounter = 0; FP_MeanRelaxFactor = 0.; FP_MaxRelaxFactor = 0.; FP_MaxRotRelaxFactor = 0.; FP_PrandtlNumber = 0. @@ -596,6 +600,8 @@ SUBROUTINE BGK_quadtree_adapt(iElem) BGK_QualityFacSamp(5,iElem) = BGK_QualityFacSamp(5,iElem) + BGK_MaxRotRelaxFactor BGK_QualityFacSamp(6,iElem) = BGK_QualityFacSamp(6,iElem) + BGK_PrandtlNumber BGK_QualityFacSamp(7,iElem) = BGK_QualityFacSamp(7,iElem) + BGK_ExpectedPrandtlNumber + BGK_QualityFacSamp(8,iElem) = BGK_QualityFacSamp(8,iElem) + BGK_Viscosity + BGK_QualityFacSamp(9,iElem) = BGK_QualityFacSamp(9,iElem) + BGK_ThermalConductivity END IF IF(FPInitDone) THEN FP_QualityFacSamp(1,iElem) = FP_QualityFacSamp(1,iElem) + FP_MeanRelaxFactor diff --git a/src/particles/bgk/bgk_colloperator.f90 b/src/particles/bgk/bgk_colloperator.f90 index a5ff618e2..fb1a8cf74 100644 --- a/src/particles/bgk/bgk_colloperator.f90 +++ b/src/particles/bgk/bgk_colloperator.f90 @@ -29,7 +29,7 @@ MODULE MOD_BGK_CollOperator !----------------------------------------------------------------------------------------------------------------------------------- ! Private Part --------------------------------------------------------------------------------------------------------------------- ! Public Part ---------------------------------------------------------------------------------------------------------------------- -PUBLIC :: BGK_CollisionOperator, ARShakhov, CalcTEquiPoly, CalcTEqui +PUBLIC :: BGK_CollisionOperator, ARShakhov !=================================================================================================================================== CONTAINS @@ -56,7 +56,7 @@ SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume, AveragingVal USE MOD_TimeDisc_Vars ,ONLY: dt USE MOD_BGK_Vars ,ONLY: SpecBGK, BGKDoVibRelaxation, BGKMovingAverage USE MOD_BGK_Vars ,ONLY: BGK_MeanRelaxFactor, BGK_MeanRelaxFactorCounter, BGK_MaxRelaxFactor, BGK_MaxRotRelaxFactor -USE MOD_BGK_Vars ,ONLY: BGK_PrandtlNumber +USE MOD_BGK_Vars ,ONLY: BGK_PrandtlNumber, BGK_Viscosity, BGK_ThermalConductivity USE MOD_part_tools ,ONLY: GetParticleWeight #ifdef CODE_ANALYZE USE MOD_Globals ,ONLY: abort,unit_stdout,myrank @@ -74,12 +74,13 @@ SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume, AveragingVal !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES REAL :: vBulk(3), u0ij(3,3), u2, V_rel(3), dtCell -REAL :: alpha, CellTemp, dens, InnerDOF, NewEn, OldEn, Prandtl, relaxfreq, TEqui +REAL :: alpha, alphaRot(nSpecies), CellTemp, dens, InnerDOF, NewEn, OldEn, Prandtl, relaxfreq +REAL :: dynamicvis, thermalcond INTEGER, ALLOCATABLE :: iPartIndx_NodeRelax(:),iPartIndx_NodeRelaxTemp(:),iPartIndx_NodeRelaxRot(:),iPartIndx_NodeRelaxVib(:) -INTEGER :: iLoop, iPart, nRelax, iPolyatMole, nXiVibDOF -REAL, ALLOCATABLE :: Xi_vib_DOF(:), VibEnergyDOF(:,:) +INTEGER :: iLoop, iPart, nRelax, nXiVibDOF +REAL, ALLOCATABLE :: Xi_vib_DOF(:,:), VibEnergyDOF(:,:) INTEGER :: iSpec, nSpec(nSpecies), jSpec, nRotRelax, nVibRelax -REAL :: OldEnRot, NewEnRot, NewEnVib +REAL :: OldEnRot, NewEnRot(nSpecies), NewEnVib(nSpecies) REAL :: TotalMass, u2Spec(nSpecies), u2i(3), vBulkAll(3) REAL :: SpecTemp(nSpecies) #ifdef CODE_ANALYZE @@ -88,11 +89,13 @@ SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume, AveragingVal REAL,PARAMETER :: RelMomTol=1e-6 ! Relative tolerance applied to conservation of momentum before/after reaction REAL,PARAMETER :: RelEneTol=1e-12 ! Relative tolerance applied to conservation of energy before/after reaction #endif /* CODE_ANALYZE */ -REAL :: totalWeightSpec(nSpecies), totalWeight, partWeight, CellTemptmp -REAL :: EVibSpec(nSpecies), ERotSpec(nSpecies), Xi_VibSpec(nSpecies), Xi_RotSpec(nSpecies),Xi_Vib_oldSpec(nSpecies) -REAL :: TVibSpec(nSpecies), TRotSpec(nSpecies), RotExpSpec(nSpecies), VibExpSpec(nSpecies) -REAL :: collisionfreqSpec(nSpecies),rotrelaxfreqSpec(nSpecies), vibrelaxfreqSpec(nSpecies), Xi_RotTotal -INTEGER :: nVibRelaxSpec(nSpecies), nRotRelaxSpec(nSpecies) +REAL :: totalWeightSpec(nSpecies), totalWeight, totalWeight2, partWeight, CellTemptmp, MassIC_Mixture +REAL :: EVibSpec(nSpecies), Xi_VibSpec(nSpecies), Xi_VibSpecNew(nSpecies) +REAL :: ERotSpec(nSpecies), Xi_RotSpec(nSpecies), Xi_RotTotal +REAL :: CellTempRel, TEqui +REAL :: TVibSpec(nSpecies), TRotSpec(nSpecies), VibRelaxWeightSpec(nSpecies), RotRelaxWeightSpec(nSpecies) +REAL :: collisionfreqSpec(nSpecies), rotrelaxfreqSpec(nSpecies), vibrelaxfreqSpec(nSpecies) +REAL :: betaR(nSpecies), betaV(nSpecies) !=================================================================================================================================== #ifdef CODE_ANALYZE ! Momentum and energy conservation check: summing up old values @@ -104,6 +107,7 @@ SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume, AveragingVal Momentum_old(1:3) = Momentum_old(1:3) + PartState(4:6,iPart)*Species(iSpec)%MassIC*partWeight Energy_old = Energy_old + DOTPRODUCT(PartState(4:6,iPart))*0.5*Species(iSpec)%MassIC*partWeight IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN + ! Add internal energies (vibration, rotation) for molecules and molecular ions Energy_old = Energy_old + (PartStateIntEn(1,iPart) + PartStateIntEn(2,iPart))*partWeight END IF END DO @@ -112,9 +116,10 @@ SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume, AveragingVal IF(nPart.LE.2) RETURN ! 1.) Moment calculation: Summing up the relative velocities and their squares -CALL CalcMoments(nPart, iPartIndx_Node, nSpec, vBulkAll, totalWeight, totalWeightSpec, TotalMass, u2, u2Spec, u0ij, & +CALL CalcMoments(nPart, iPartIndx_Node, nSpec, vBulkAll, totalWeight, totalWeight2, totalWeightSpec, TotalMass, u2, u2Spec, u0ij, & u2i, OldEn, EVibSpec, ERotSpec, CellTemp, SpecTemp, dtCell) -IF((CellTemp.LE.0).OR.(MAXVAL(nSpec(:)).EQ.1).OR.(totalWeight.LE.0.0)) RETURN + +IF((CellTemp.LE.0.0).OR.(MAXVAL(nSpec(:)).EQ.1).OR.(totalWeight.LE.0.0)) RETURN IF(UseVarTimeStep) THEN dtCell = dt * dtCell / totalWeight @@ -126,41 +131,31 @@ SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume, AveragingVal ! totalWeight contains the weighted particle number dens = totalWeight / NodeVolume ELSE + ! MPF is the same for all species dens = totalWeight * Species(1)%MacroParticleFactor / NodeVolume END IF IF (BGKMovingAverage) THEN CALL DoAveraging(dens, u2, u0ij, u2i, CellTemp, AveragingValues) END IF -! Calculation of the rotational and vibrational degrees of freedom for molecules -nXiVibDOF=0 ! Initialize -IF (nSpecies.EQ.1) THEN - IF((SpecDSMC(1)%InterID.EQ.2).OR.(SpecDSMC(1)%InterID.EQ.20)) THEN - IF(SpecDSMC(1)%PolyatomicMol) THEN - iPolyatMole = SpecDSMC(1)%SpecToPolyArray - nXiVibDOF = PolyatomMolDSMC(iPolyatMole)%VibDOF - ALLOCATE(Xi_vib_DOF(nXiVibDOF)) - Xi_vib_DOF(:) = 0. - END IF - END IF -END IF -CALL CalcInnerDOFs(nSpec, EVibSpec, ERotSpec, totalWeightSpec, TVibSpec, TRotSpec, InnerDOF, Xi_VibSpec, Xi_Vib_oldSpec & - ,Xi_RotSpec) +CALL CalcInnerDOFs(nSpec, EVibSpec, ERotSpec, totalWeightSpec, TVibSpec, TRotSpec, InnerDOF, Xi_VibSpec, Xi_RotSpec) ! 2.) Calculation of the relaxation frequency of the distribution function towards the target distribution function -CALL CalcGasProperties(nSpec, dens, InnerDOF, totalWeightSpec, totalWeight, TotalMass, u2Spec, u0ij, u2, SpecTemp, CellTemp, & - Xi_VibSpec, Xi_RotSpec, Prandtl, relaxfreq) +CALL CalcGasProperties(nSpec, dens, InnerDOF, totalWeightSpec, totalWeight, TotalMass, u2Spec, SpecTemp, CellTemp, Xi_VibSpec, & + Xi_RotSpec, Prandtl, relaxfreq, dynamicvis, thermalcond, MassIC_Mixture) IF(DSMC%CalcQualityFactors) THEN BGK_MeanRelaxFactor = BGK_MeanRelaxFactor + relaxfreq * dtCell BGK_MeanRelaxFactorCounter = BGK_MeanRelaxFactorCounter + 1 BGK_MaxRelaxFactor = MAX(BGK_MaxRelaxFactor,relaxfreq*dtCell) BGK_PrandtlNumber = BGK_PrandtlNumber + Prandtl + BGK_Viscosity = BGK_Viscosity + dynamicvis + BGK_ThermalConductivity = BGK_ThermalConductivity + thermalcond END IF ! 3.) Treatment of molecules: determination of the rotational and vibrational relaxation frequency using the collision frequency, -! which is not the same as the relaxation frequency of distribution function, calculated above. +! which is not the same as the relaxation frequency of distribution function, calculated above IF(ANY(SpecDSMC(:)%InterID.EQ.2).OR.ANY(SpecDSMC(:)%InterID.EQ.20)) THEN collisionfreqSpec = 0.0 DO iSpec = 1, nSpecies @@ -170,79 +165,117 @@ SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume, AveragingVal ELSE CellTemptmp = CellTemp END IF - collisionfreqSpec(iSpec) = collisionfreqSpec(iSpec) + SpecBGK(iSpec)%CollFreqPreFactor(jSpec) * totalWeightSpec(iSpec)*totalWeightSpec(jSpec) & - *Dens *CellTemptmp**(-CollInf%omega(iSpec,jSpec) +0.5) /(totalWeight*totalWeight) + ! Sum up collision frequencies of species i with itself and the other species + ! S. Chapman and T.G. Cowling, "The mathematical Theory of Non-Uniform Gases", Cambridge University Press, 1970, S. 87f + ! For SpecBGK(iSpec)%CollFreqPreFactor(jSpec) see bgk_init.f90 + ! VHS according to M. Pfeiffer, "Extending the particle ellipsoidal statistical Bhatnagar-Gross-Krook method to diatomic + ! molecules including quantized vibrational energies", Phys. Fluids 30, 116103 (2018), Eq. (18) + collisionfreqSpec(iSpec) = collisionfreqSpec(iSpec) + SpecBGK(iSpec)%CollFreqPreFactor(jSpec) * totalWeightSpec(jSpec) & + * (Dens / totalWeight) *CellTemptmp**(-CollInf%omega(iSpec,jSpec) +0.5) END DO END DO + ! Calculate relaxation frequencies of rotation and vibration with relaxation properties + ! M. Pfeiffer, "Extending the particle ellipsoidal statistical Bhatnagar-Gross-Krook method to diatomic molecules including + ! quantized vibrational energies", Phys. Fluids 30, 116103 (2018) + ! N.E. Gimelshein et. al, "Vibrational relaxation rates in the direct simulation Monte Carlo method", Phys. Fluids 14, 4452 (2018) + ! relaxfreqSpec = collisionfreqSpec / collision number Z with RelaxProb = 1/Z rotrelaxfreqSpec(:) = collisionfreqSpec(:) * DSMC%RotRelaxProb vibrelaxfreqSpec(:) = collisionfreqSpec(:) * DSMC%VibRelaxProb - RotExpSpec=0.; VibExpSpec=0. - - IF(SpecDSMC(1)%PolyatomicMol) THEN - CALL CalcTEquiPoly(nPart, CellTemp, TRotSpec(1), TVibSpec(1), nXiVibDOF, Xi_vib_DOF, Xi_Vib_oldSpec(1), RotExpSpec(1), VibExpSpec(1), & - TEqui, rotrelaxfreqSpec(1), vibrelaxfreqSpec(1), dtCell) - Xi_VibSpec(1) = SUM(Xi_vib_DOF(1:PolyatomMolDSMC(iPolyatMole)%VibDOF)) - ELSE - CALL CalcTEquiMulti(nPart, nSpec, CellTemp, TRotSpec, TVibSpec, Xi_VibSpec, Xi_Vib_oldSpec, RotExpSpec, VibExpSpec, & - TEqui, rotrelaxfreqSpec, vibrelaxfreqSpec, dtCell) - END IF + IF(DSMC%CalcQualityFactors) THEN BGK_MaxRotRelaxFactor = MAX(BGK_MaxRotRelaxFactor,MAXVAL(rotrelaxfreqSpec(:))*dtCell) END IF END IF -! 4.) Determine the number of particles undergoing a relaxation (including vibration and rotation) +! 4.) Determine the relaxation temperatures and the number of particles undergoing a relaxation (including vibration + rotation) ALLOCATE(iPartIndx_NodeRelax(nPart), iPartIndx_NodeRelaxTemp(nPart)) iPartIndx_NodeRelax = 0; iPartIndx_NodeRelaxTemp = 0 ALLOCATE(iPartIndx_NodeRelaxRot(nPart),iPartIndx_NodeRelaxVib(nPart)) iPartIndx_NodeRelaxRot = 0; iPartIndx_NodeRelaxVib = 0 -CALL DetermineRelaxPart(nPart, iPartIndx_Node, relaxfreq, dtCell, RotExpSpec, VibExpSpec, nRelax, nRotRelax, nVibRelax, & - nRotRelaxSpec, nVibRelaxSpec, iPartIndx_NodeRelax, iPartIndx_NodeRelaxTemp, iPartIndx_NodeRelaxRot, & - iPartIndx_NodeRelaxVib, vBulk, OldEnRot, OldEn) -IF ((nRelax.EQ.0).AND.(nRotRelax.EQ.0).AND.(nVibRelax.EQ.0)) RETURN +! Allocate Xi_vib_DOF +IF(BGKDoVibRelaxation) THEN + IF(DSMC%NumPolyatomMolecs.GT.0) THEN + nXiVibDOF = MAXVAL(PolyatomMolDSMC(:)%VibDOF) + ALLOCATE(Xi_vib_DOF(DSMC%NumPolyatomMolecs,nXiVibDOF)) + END IF +END IF + +CALL CalcTRelax(ERotSpec, Xi_RotSpec, Xi_VibSpec, EVibSpec, totalWeightSpec, totalWeight, totalWeight2, dtCell, CellTemp, & + TRotSpec, TVibSpec, relaxfreq, rotrelaxfreqSpec, vibrelaxfreqSpec, Xi_VibSpecNew, Xi_vib_DOF, nXiVibDOF, CellTempRel, TEqui, & + betaR, betaV) + +CALL DetermineRelaxPart(nPart, iPartIndx_Node, relaxfreq, dtCell, nRelax, nRotRelax, nVibRelax, & + RotRelaxWeightSpec, VibRelaxWeightSpec, iPartIndx_NodeRelax, iPartIndx_NodeRelaxTemp, iPartIndx_NodeRelaxRot, & + iPartIndx_NodeRelaxVib, vBulk, OldEnRot, OldEn, rotrelaxfreqSpec, vibrelaxfreqSpec, betaR, betaV) +! Allocate VibEnergyDOF IF(BGKDoVibRelaxation) THEN - IF(SpecDSMC(1)%PolyatomicMol) THEN - ALLOCATE(VibEnergyDOF(nVibRelax,PolyatomMolDSMC(iPolyatMole)%VibDOF)) - END IF + IF(DSMC%NumPolyatomMolecs.GT.0) THEN + ALLOCATE(VibEnergyDOF(nVibRelax,nXiVibDOF)) + VibEnergyDOF = 0.0 + END IF +END IF + +! Return if no particles are undergoing a relaxation +IF ((nRelax.EQ.0).AND.(nRotRelax.EQ.0).AND.(nVibRelax.EQ.0)) RETURN + +! 5.) Determine the new rotational and vibrational states of molecules undergoing a relaxation +IF(ANY(SpecDSMC(:)%InterID.EQ.2).OR.ANY(SpecDSMC(:)%InterID.EQ.20)) THEN + CALL RelaxInnerEnergy(nVibRelax, nRotRelax, iPartIndx_NodeRelaxVib, iPartIndx_NodeRelaxRot, nXiVibDOF, Xi_vib_DOF, & + Xi_VibSpecNew, Xi_RotSpec, VibEnergyDOF, TEqui, NewEnVib, NewEnRot) END IF -! 5.) Determine the new rotational and vibrational state of molecules undergoing a relaxation -CALL RelaxInnerEnergy(nPart, nVibRelax, nRotRelax, iPartIndx_NodeRelaxVib, iPartIndx_NodeRelaxRot, nXiVibDOF, Xi_vib_DOF, Xi_VibSpec, & - Xi_RotSpec , TEqui, VibEnergyDOF, NewEnVib, NewEnRot) ! 6.) Sample new particle velocities from the target distribution function, depending on the chosen model -CALL SampleFromTargetDistr(nRelax, iPartIndx_NodeRelax, Prandtl, u2, u0ij, u2i, vBulkAll, CellTemp, vBulk) +CALL SampleFromTargetDistr(nRelax, iPartIndx_NodeRelax, Prandtl, u2, u0ij, u2i, vBulkAll, CellTempRel, CellTemp, vBulk, & + MassIC_Mixture) NewEn = 0. + +! Calculation of the new bulk velocity vBulk = vBulk/TotalMass + +! Loop over all relaxing particles for calculation of the new energy DO iLoop = 1, nRelax iPart = iPartIndx_NodeRelax(iLoop) iSpec = PartSpecies(iPart) partWeight = GetParticleWeight(iPart) + ! Thermal velocity of all relaxing particles V_rel(1:3) = PartState(4:6,iPart) - vBulk(1:3) + ! Sum up kinetic energies NewEn = NewEn + (V_rel(1)**(2.) + V_rel(2)**(2.) + V_rel(3)**(2.))*0.5*Species(iSpec)%MassIC*partWeight END DO +! Loop over all non-relaxing particles for calculation of the new energy DO iLoop = 1, nPart-nRelax iPart = iPartIndx_NodeRelaxTemp(iLoop) iSpec = PartSpecies(iPart) partWeight = GetParticleWeight(iPart) + ! Thermal velocity of all non-relaxing particles V_rel(1:3) = PartState(4:6,iPart) - vBulk(1:3) + ! Sum up kinetic energies NewEn = NewEn + (V_rel(1)**(2.) + V_rel(2)**(2.) + V_rel(3)**(2.))*0.5*Species(iSpec)%MassIC*partWeight END DO ! 7.) Vibrational energy of the molecules: Ensure energy conservation by scaling the new vibrational states with the factor alpha IF(ANY(SpecDSMC(:)%InterID.EQ.2).OR.ANY(SpecDSMC(:)%InterID.EQ.20)) THEN - CALL EnergyConsVib(nPart, nVibRelax, nVibRelaxSpec, iPartIndx_NodeRelaxVib, NewEnVib, OldEn, nXiVibDOF, Xi_VibSpec, VibEnergyDOF, TEqui) + CALL EnergyConsVib(nPart, totalWeight, nVibRelax, VibRelaxWeightSpec, iPartIndx_NodeRelaxVib, NewEnVib, OldEn, nXiVibDOF, & + VibEnergyDOF, Xi_VibSpecNew, TEqui) END IF +! Remaining vibrational + translational energy + old rotational energy for translation and rotation OldEn = OldEn + OldEnRot -! 8.) Determine the new particle state and ensure energy conservation by scaling the new velocities with the factor alpha. Xi_RotTotal = 0.0 DO iSpec = 1, nSpecies - Xi_RotTotal = Xi_RotTotal + Xi_RotSpec(iSpec)*nRotRelaxSpec(iSpec) + ! Sum of relaxing rotational degrees of freedom + Xi_RotTotal = Xi_RotTotal + Xi_RotSpec(iSpec)*RotRelaxWeightSpec(iSpec) END DO -alpha = SQRT(OldEn/NewEn*(3.*(nPart-1.))/(Xi_RotTotal+3.*(nPart-1.))) + +! 8.) Determine the new particle state and ensure energy conservation by scaling the new velocities with the factor alpha +! Calculation of scaling factor alpha +alpha = SQRT(OldEn/NewEn*(3.*(totalWeight-1.))/(Xi_RotTotal+3.*(totalWeight-1.))) +! Calculation of the final particle velocities with vBulkAll (average flow velocity before relaxation), scaling factor alpha, +! the particle velocity PartState(4:6,iPart) after the relaxation but before the energy conservation and vBulk (average value of +! the latter) DO iLoop = 1, nRelax iPart = iPartIndx_NodeRelax(iLoop) PartState(4:6,iPart) = vBulkAll(1:3) + alpha*(PartState(4:6,iPart)-vBulk(1:3)) @@ -263,10 +296,20 @@ SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume, AveragingVal END IF ! 9.) Rotation: Scale the new rotational state of the molecules to ensure energy conservation -IF ( (nRotRelax.GT.0)) alpha = OldEn/NewEnRot*(Xi_RotTotal/(Xi_RotTotal+3.*(nPart-1.))) +DO iSpec = 1, nSpecies + ! Calculate scaling factor alpha per species, see F. Hild, M. Pfeiffer, "Multi-species modeling in the particle-based ellipsoidal + ! statistical Bhatnagar-Gross-Krook method including internal degrees of freedom", subitted to Phys. Fluids, August 2023 + IF (NewEnRot(iSpec).GT.0.0) THEN + alphaRot(iSpec) = OldEn/NewEnRot(iSpec)*(Xi_RotSpec(iSpec)*RotRelaxWeightSpec(iSpec)/(Xi_RotTotal+3.*(totalWeight-1.))) + ELSE + alphaRot(iSpec) = 0.0 + END IF +END DO DO iLoop = 1, nRotRelax iPart = iPartIndx_NodeRelaxRot(iLoop) - PartStateIntEn( 2,iPart) = alpha*PartStateIntEn( 2,iPart) + iSpec = PartSpecies(iPart) + ! Scaling of rotational energy with factor alpha + PartStateIntEn( 2,iPart) = alphaRot(iSpec)*PartStateIntEn( 2,iPart) END DO ! CODE ANALYZE: Compare the old momentum and energy of the cell with the new, abort if relative difference is above the limits @@ -324,8 +367,9 @@ SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume, AveragingVal END SUBROUTINE BGK_CollisionOperator -SUBROUTINE CalcMoments(nPart, iPartIndx_Node, nSpec, vBulkAll, totalWeight, totalWeightSpec, TotalMass, u2, u2Spec, & - u0ij, u2i, OldEn, EVibSpec, ERotSpec, CellTemp, SpecTemp, dtCell) + +SUBROUTINE CalcMoments(nPart, iPartIndx_Node, nSpec, vBulkAll, totalWeight, totalWeight2, totalWeightSpec, TotalMass, u2, u2Spec, & + u0ij, u2i, OldEn, EVibSpec, ERotSpec, CellTemp, SpecTemp, dtCell) !=================================================================================================================================== !> Moment calculation: Summing up the relative velocities and their squares !=================================================================================================================================== @@ -335,6 +379,7 @@ SUBROUTINE CalcMoments(nPart, iPartIndx_Node, nSpec, vBulkAll, totalWeight, tota USE MOD_BGK_Vars ,ONLY: BGKDoVibRelaxation, BGKCollModel USE MOD_part_tools ,ONLY: GetParticleWeight USE MOD_Globals_Vars ,ONLY: BoltzmannConst +USE MOD_Globals ,ONLY: DOTPRODUCT ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -345,14 +390,16 @@ SUBROUTINE CalcMoments(nPart, iPartIndx_Node, nSpec, vBulkAll, totalWeight, tota INTEGER, INTENT(OUT) :: nSpec(nSpecies) REAL, INTENT(OUT) :: u2Spec(nSpecies),u0ij(3,3), OldEn, EVibSpec(nSpecies), ERotSpec(nSpecies), u2i(3), u2 REAL, INTENT(OUT) :: CellTemp, SpecTemp(nSpecies), totalWeightSpec(nSpecies) -REAL, INTENT(OUT) :: vBulkAll(3), totalWeight, TotalMass, dtCell +REAL, INTENT(OUT) :: vBulkAll(3), totalWeight, totalWeight2, TotalMass, dtCell !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES INTEGER :: iLoop, iPart, iSpec, fillMa1, fillMa2 -REAL :: V_rel(1:3), vmag2, partWeight, EnerTotal, totalWeightSpec2(nSpecies), vBulkSpec(3,nSpecies) -REAL :: tempweight, tempweight2, tempmass, vBulkTemp(3), totalWeight2, totalWeight3 +REAL :: V_rel(1:3), vmag2, EnerTotal, ThermEner, totalWeightSpec2(nSpecies), vBulkSpec(3,nSpecies) +REAL :: partWeight, tempweight, tempweight2, tempmass, vBulkTemp(3), totalWeight3 +LOGICAL :: validSpec(nSpecies) !=================================================================================================================================== totalWeightSpec = 0.0; totalWeightSpec2=0.0; vBulkAll=0.0; TotalMass=0.0; vBulkSpec=0.0; nSpec=0; dtCell=0.0 +! Loop over all simulation particles to sum up particle velocities to calculate bulk velocities DO iLoop = 1, nPart iPart = iPartIndx_Node(iLoop) partWeight = GetParticleWeight(iPart) @@ -362,31 +409,36 @@ SUBROUTINE CalcMoments(nPart, iPartIndx_Node, nSpec, vBulkAll, totalWeight, tota vBulkAll(1:3) = vBulkAll(1:3) + PartState(4:6,iPart)*Species(iSpec)%MassIC*partWeight TotalMass = TotalMass + Species(iSpec)%MassIC*partWeight vBulkSpec(1:3,iSpec) = vBulkSpec(1:3,iSpec) + PartState(4:6,iPart)*partWeight - nSpec(iSpec) = nSpec(iSpec) + 1 + nSpec(iSpec) = nSpec(iSpec) + 1 ! Count number of simulation particles per species IF(UseVarTimeStep) THEN dtCell = dtCell + PartTimeStep(iPart)*partWeight END IF END DO totalWeight = SUM(totalWeightSpec) totalWeight2 = SUM(totalWeightSpec2) -IF ((MAXVAL(nSpec(:)).EQ.1).OR.(totalWeight.LE.0.0)) RETURN + +! Calculate total bulk velocity and bulk velocities per species vBulkAll(1:3) = vBulkAll(1:3) / TotalMass DO iSpec = 1, nSpecies - IF (nSpec(iSpec).GT.0) vBulkSpec(:,iSpec) = vBulkSpec(:,iSpec) /totalWeightSpec(iSpec) + IF (nSpec(iSpec).GT.0) vBulkSpec(:,iSpec) = vBulkSpec(:,iSpec) / totalWeightSpec(iSpec) END DO totalWeight3 = 0.; u2Spec=0.0; u0ij=0.0; u2i=0.0; OldEn=0.0; EVibSpec=0.0; ERotSpec=0.0 +! Loop over all simulation particles to sum up relative velocities DO iLoop = 1, nPart iPart = iPartIndx_Node(iLoop) partWeight = GetParticleWeight(iPart) iSpec = PartSpecies(iPart) + ! Calculate thermal velocity with bulk velocity of the species V_rel(1:3)=PartState(4:6,iPart)-vBulkSpec(1:3,iSpec) vmag2 = V_rel(1)**2 + V_rel(2)**2 + V_rel(3)**2 + ! Summing up thermal velocities (squared) of all particles per species u2Spec(iSpec) = u2Spec(iSpec) + vmag2*partWeight + ! Calculate thermal velocity with bulk velocity of the gas V_rel(1:3)=PartState(4:6,iPart)-vBulkAll(1:3) vmag2 = V_rel(1)**2 + V_rel(2)**2 + V_rel(3)**2 - IF (BGKCollModel.EQ.1) THEN + IF (BGKCollModel.EQ.1) THEN ! ESBGK DO fillMa1 =1, 3 DO fillMa2 =fillMa1, 3 u0ij(fillMa1, fillMa2)= u0ij(fillMa1, fillMa2) & @@ -394,34 +446,48 @@ SUBROUTINE CalcMoments(nPart, iPartIndx_Node, nSpec, vBulkAll, totalWeight, tota END DO END DO END IF - IF (BGKCollModel.EQ.2) THEN + IF (BGKCollModel.EQ.2) THEN ! Shakhov u2i(1:3) = u2i(1:3) + V_rel(1:3)*vmag2 * partWeight*Species(iSpec)%MassIC totalWeight3 = totalWeight3 + partWeight*partWeight*partWeight END IF + + ! Sum up old energy of thermal velocities and sum up internal energies --> E_T OldEn = OldEn + 0.5*Species(iSpec)%MassIC * vmag2*partWeight IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN IF(BGKDoVibRelaxation) THEN + ! EVib without zero-point energy EVibSpec(iSpec) = EVibSpec(iSpec) + (PartStateIntEn(1,iPart) - SpecDSMC(iSpec)%EZeroPoint) * partWeight END IF ERotSpec(iSpec) = ERotSpec(iSpec) + PartStateIntEn(2,iPart) * partWeight END IF END DO -u0ij = u0ij* totalWeight / (TotalMass*(totalWeight - totalWeight2/totalWeight)) -IF (BGKCollModel.EQ.2) THEN +IF (BGKCollModel.EQ.1) THEN ! ESBGK + ! Pressure tensor + u0ij = u0ij* totalWeight / (TotalMass*(totalWeight - totalWeight2/totalWeight)) +END IF +IF (BGKCollModel.EQ.2) THEN ! Shakhov + ! Heatflux u2i = u2i*totalWeight**3/(TotalMass*(totalWeight**3-3.*totalWeight*totalWeight2+2.*totalWeight3)) END IF -IF (nSpecies.GT.1) THEN +! Calculation of cell temperature and bulk velocity (squared) +IF (nSpecies.GT.1) THEN ! mixture + validSpec = .FALSE. SpecTemp = 0.0 EnerTotal = 0.0 tempweight = 0.0; tempweight2 = 0.0; tempmass = 0.0; vBulkTemp = 0.0 DO iSpec = 1, nSpecies + ! At least two particles and non-zero squared thermal velocity needed for a valid species IF ((nSpec(iSpec).GE.2).AND.(.NOT.ALMOSTZERO(u2Spec(iSpec)))) THEN + validSpec = .TRUE. + ! Calculation of the species temperature --> translational temperatures of the different species SpecTemp(iSpec) = Species(iSpec)%MassIC * u2Spec(iSpec) & /(3.0*BoltzmannConst*(totalWeightSpec(iSpec) - totalWeightSpec2(iSpec)/totalWeightSpec(iSpec))) + ! Thermal energy EnerTotal = EnerTotal + 3./2.*BoltzmannConst*SpecTemp(iSpec) * totalWeightSpec(iSpec) - vmag2 = vBulkSpec(1,iSpec)**(2.) + vBulkSpec(2,iSpec)**(2.) + vBulkSpec(3,iSpec)**(2.) + vmag2 = DOTPRODUCT(vBulkSpec(1:3,iSpec)) + ! Add kinetic energy EnerTotal = EnerTotal + totalWeightSpec(iSpec) * Species(iSpec)%MassIC / 2. * vmag2 tempweight = tempweight + totalWeightSpec(iSpec) tempweight2 = tempweight2 + totalWeightSpec2(iSpec) @@ -429,22 +495,31 @@ SUBROUTINE CalcMoments(nPart, iPartIndx_Node, nSpec, vBulkAll, totalWeight, tota vBulkTemp(1:3) = vBulkTemp(1:3) + vBulkSpec(1:3,iSpec)*totalWeightSpec(iSpec) * Species(iSpec)%MassIC END IF END DO - - vBulkTemp(1:3) = vBulkTemp(1:3) / tempmass - vmag2 = vBulkTemp(1)*vBulkTemp(1) + vBulkTemp(2)*vBulkTemp(2) + vBulkTemp(3)*vBulkTemp(3) - EnerTotal = EnerTotal - tempmass / 2. * vmag2 - CellTemp = 2. * EnerTotal / (3.*tempweight*BoltzmannConst) - u2 = 3. * CellTemp * BoltzmannConst * (tempweight - tempweight2/tempweight) / tempmass -ELSE + IF (ANY(validSpec)) THEN + vBulkTemp(1:3) = vBulkTemp(1:3) / tempmass + ! Squared bulk velocity of the mixture + vmag2 = DOTPRODUCT(vBulkTemp(1:3)) + ! EnerTotal = kinetic energy (tempmass / 2. * vmag2) + thermal energy (3. * tempweight * BoltzmannConst * CellTemp / 2) + ThermEner = EnerTotal - tempmass / 2. * vmag2 + ! Calculation of the cell temperature from the thermal energy --> translational temperature of the mixture + CellTemp = 2. * ThermEner / (3.*tempweight*BoltzmannConst) + ! Mean squared thermal velocity c^2 of a particle, calculated with the cell temperature and the density-averaged mass + u2 = 3. * CellTemp * BoltzmannConst * (tempweight - tempweight2/tempweight) / tempmass + ELSE ! only one part per species or cloned species with u2spec = 0 because PartState(4:6) = vBulkAll + u2 = OldEn / (TotalMass*(1. - totalWeight2/totalWeight**2)) * 2. ! variance-free + CellTemp = TotalMass/totalWeight * u2 / (3.0*BoltzmannConst) + END IF +ELSE ! single species gas u2 = u2Spec(1) / (totalWeight - totalWeight2/totalWeight) CellTemp = Species(1)%MassIC * u2 / (3.0*BoltzmannConst) END IF END SUBROUTINE CalcMoments + SUBROUTINE DoAveraging(dens, u2, u0ij, u2i, CellTemp, AverageValues) !=================================================================================================================================== -!> Moment calculation: Summing up the relative velocities and their squares +!> BGK moving average !=================================================================================================================================== ! MODULES USE MOD_Particle_Vars ,ONLY: Species @@ -505,18 +580,18 @@ SUBROUTINE DoAveraging(dens, u2, u0ij, u2i, CellTemp, AverageValues) dens = AverageValues(1) u2 = AverageValues(2) END IF - CellTemp = Species(1)%MassIC * u2 / (3.0*BoltzmannConst) +CellTemp = Species(1)%MassIC * u2 / (3.0*BoltzmannConst) END SUBROUTINE DoAveraging -SUBROUTINE CalcInnerDOFs(nSpec, EVibSpec, ERotSpec, totalWeightSpec, TVibSpec, TRotSpec, InnerDOF, Xi_VibSpec, Xi_Vib_oldSpec & - ,Xi_RotSpec) + +SUBROUTINE CalcInnerDOFs(nSpec, EVibSpec, ERotSpec, totalWeightSpec, TVibSpec, TRotSpec, InnerDOF, Xi_VibSpec, Xi_RotSpec) !=================================================================================================================================== !> Determine the internal degrees of freedom and the respective temperature (rotation/vibration) for diatomic/polyatomic species !=================================================================================================================================== ! MODULES USE MOD_Particle_Vars ,ONLY: nSpecies -USE MOD_DSMC_Vars ,ONLY: SpecDSMC, PolyatomMolDSMC +USE MOD_DSMC_Vars ,ONLY: SpecDSMC USE MOD_BGK_Vars ,ONLY: BGKDoVibRelaxation USE MOD_Globals_Vars ,ONLY: BoltzmannConst USE MOD_Particle_Analyze_Tools ,ONLY: CalcTVibPoly @@ -528,41 +603,34 @@ SUBROUTINE CalcInnerDOFs(nSpec, EVibSpec, ERotSpec, totalWeightSpec, TVibSpec, T REAL, INTENT(IN) :: EVibSpec(nSpecies), ERotSpec(nSpecies), totalWeightSpec(nSpecies) !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES -REAL, INTENT(OUT) :: TVibSpec(nSpecies), TRotSpec(nSpecies), InnerDOF, Xi_VibSpec(nSpecies), Xi_Vib_oldSpec(nSpecies) +REAL, INTENT(OUT) :: TVibSpec(nSpecies), TRotSpec(nSpecies), InnerDOF, Xi_VibSpec(nSpecies) REAL, INTENT(OUT) :: Xi_RotSpec(nSpecies) !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER :: iPolyatMole, iSpec, iDOF -REAL :: exparg +INTEGER :: iSpec !=================================================================================================================================== -Xi_VibSpec=0.; InnerDOF=0.; Xi_RotSpec=0.; Xi_Vib_oldSpec=0.; TVibSpec=0.; TRotSpec=0. +Xi_VibSpec=0.; InnerDOF=0.; Xi_RotSpec=0.; TVibSpec=0.; TRotSpec=0. DO iSpec = 1, nSpecies IF (nSpec(iSpec).EQ.0) CYCLE + ! Only for molecules IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN IF(BGKDoVibRelaxation) THEN - IF(SpecDSMC(iSpec)%PolyatomicMol) THEN - iPolyatMole = SpecDSMC(iSpec)%SpecToPolyArray - TVibSpec(iSpec) = CalcTVibPoly(EVibSpec(iSpec)/totalWeightSpec(iSpec), 1) - IF (TVibSpec(iSpec).GT.0.0) THEN - DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF - exparg = PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)/TVibSpec(iSpec) - IF(CHECKEXP(exparg))THEN - Xi_VibSpec(iSpec) = Xi_VibSpec(iSpec) + 2.*exparg/(EXP(exparg) - 1.) - ELSE - Xi_VibSpec(iSpec) = 0. - END IF ! CHECKEXP(exparg) - END DO - END IF - ELSE - TVibSpec(iSpec)=EVibSpec(iSpec) / (totalWeightSpec(iSpec)*BoltzmannConst*SpecDSMC(iSpec)%CharaTVib) - IF (TVibSpec(iSpec).GT.0.0) THEN - TVibSpec(iSpec)= SpecDSMC(iSpec)%CharaTVib/LOG(1. + 1./(TVibSpec(iSpec))) - Xi_VibSpec(iSpec) = 2.* EVibSpec(iSpec) / (totalWeightSpec(iSpec)*BoltzmannConst*TVibSpec(iSpec)) - END IF + IF(SpecDSMC(iSpec)%PolyatomicMol) THEN ! polyatomic + ! Calculation of the vibrational temperature (zero-point search) for polyatomic molecules + TVibSpec(iSpec) = CalcTVibPoly(EVibSpec(iSpec)/totalWeightSpec(iSpec) + SpecDSMC(iSpec)%EZeroPoint, iSpec) + ELSE ! diatomic + ! Calculation of vibrational temperature and DOFs from Pfeiffer, Physics of Fluids 30, 116103 (2018), "Extending the + ! particle ellipsoidal statistical Bhatnagar-Gross-Krook method to diatomic molecules including quantized vibrational + ! energies" + ! TVibSpec = vibrational energy without zero-point energy + TVibSpec(iSpec) = EVibSpec(iSpec) / (totalWeightSpec(iSpec)*BoltzmannConst*SpecDSMC(iSpec)%CharaTVib) + IF (TVibSpec(iSpec).GT.0.0) TVibSpec(iSpec) = SpecDSMC(iSpec)%CharaTVib/LOG(1. + 1./(TVibSpec(iSpec))) END IF - Xi_Vib_oldSpec(iSpec) = Xi_VibSpec(iSpec) + IF (TVibSpec(iSpec).GT.0.0) Xi_VibSpec(iSpec) = 2.* EVibSpec(iSpec) / (totalWeightSpec(iSpec)*BoltzmannConst*TVibSpec(iSpec)) END IF Xi_RotSpec(iSpec) = SpecDSMC(iSpec)%Xi_Rot + ! Calculation of rotational temperature from Pfeiffer, Physics of Fluids 30, 116103 (2018), "Extending the particle ellipsoidal + ! statistical Bhatnagar-Gross-Krook method to diatomic molecules including quantized vibrational energies" TRotSpec(iSpec) = 2.*ERotSpec(iSpec)/(Xi_RotSpec(iSpec)*totalWeightSpec(iSpec)*BoltzmannConst) END IF InnerDOF = InnerDOF + Xi_RotSpec(iSpec) + Xi_VibSpec(iSpec) @@ -570,8 +638,9 @@ SUBROUTINE CalcInnerDOFs(nSpec, EVibSpec, ERotSpec, totalWeightSpec, TVibSpec, T END SUBROUTINE CalcInnerDOFs -SUBROUTINE CalcGasProperties(nSpec, dens, InnerDOF, totalWeightSpec, totalWeight, TotalMass, u2Spec, u0ij, u2, SpecTemp, CellTemp, & - Xi_VibSpec, Xi_RotSpec, Prandtl, relaxfreq) + +SUBROUTINE CalcGasProperties(nSpec, dens, InnerDOF, totalWeightSpec, totalWeight, TotalMass, u2Spec, SpecTemp, CellTemp, & + Xi_VibSpec, Xi_RotSpec, Prandtl, relaxfreq, dynamicvis, thermalcond, MassIC_Mixture) !=================================================================================================================================== !> Calculate the reference dynamic viscosity, Prandtl number and the resulting relaxation frequency of the distribution function !=================================================================================================================================== @@ -586,38 +655,45 @@ SUBROUTINE CalcGasProperties(nSpec, dens, InnerDOF, totalWeightSpec, totalWeight ! INPUT VARIABLES INTEGER, INTENT(IN) :: nSpec(nSpecies) REAL, INTENT(IN) :: totalWeightSpec(nSpecies), totalWeight, TotalMass, u2Spec(nSpecies), SpecTemp(nSpecies), CellTemp -REAL, INTENT(IN) :: u0ij(3,3), u2, Xi_VibSpec(nSpecies), Xi_RotSpec(nSpecies), dens, InnerDOF +REAL, INTENT(IN) :: Xi_VibSpec(nSpecies), Xi_RotSpec(nSpecies), dens, InnerDOF !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES -REAL, INTENT(OUT) :: Prandtl, relaxfreq +REAL, INTENT(OUT) :: Prandtl, relaxfreq, dynamicvis, thermalcond, MassIC_Mixture !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER :: iSpec, jSpec, INFO -REAL :: MolarFraction(1:nSpecies), MassFraction(1:nSpecies), MassIC_Mixture +INTEGER :: iSpec, jSpec +REAL :: MolarFraction(1:nSpecies), DOFFraction(1:nSpecies), MassFraction(1:nSpecies) REAL :: PrandtlCorrection, dynamicvisSpec(nSpecies), thermalcondSpec(nSpecies), Phi(nSpecies) -REAL :: dynamicvis, thermalcond, C_P, nu, A(3,3), W(3), Theta, CellTempSpec(nSpecies+1), Work(100) +REAL :: TotalDOFWeight, C_P, CellTempSpec(nSpecies+1) !=================================================================================================================================== -IF (nSpecies.GT.1) THEN +MassIC_Mixture = TotalMass / totalWeight +IF (nSpecies.GT.1) THEN ! gas mixture MolarFraction(1:nSpecies) = totalWeightSpec(1:nSpecies) / totalWeight - MassIC_Mixture = TotalMass / totalWeight MassFraction(1:nSpecies) = MolarFraction(1:nSpecies) * Species(1:nSpecies)%MassIC / MassIC_Mixture + DOFFraction(1:nSpecies) = totalWeightSpec(1:nSpecies) * (5.+Xi_RotSpec(1:nSpecies)+Xi_VibSpec(1:nSpecies)) + TotalDOFWeight = SUM(DOFFraction) PrandtlCorrection = 0. C_P = 0.0 DO iSpec = 1, nSpecies IF (nSpec(iSpec).EQ.0) CYCLE - PrandtlCorrection = PrandtlCorrection + MolarFraction(iSpec)*MassIC_Mixture/Species(iSpec)%MassIC - IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN - C_P = C_P + ((5. + (Xi_VibSpec(iSpec)+Xi_RotSpec(iSpec)))/2.) * BoltzmannConst / Species(iSpec)%MassIC * MassFraction(iSpec) - ELSE - C_P = C_P + (5./2.) * BoltzmannConst / Species(iSpec)%MassIC * MassFraction(iSpec) - END IF + ! Correction of Pr for calculation of relaxation frequency, see alpha - F. Hild, M. Pfeiffer, "Multi-species modeling in the + ! particle-based ellipsoidal statistical Bhatnagar-Gross-Krook method including internal degrees of freedom", subitted to Phys. + ! Fluids, August 2023 + PrandtlCorrection = PrandtlCorrection + DOFFraction(iSpec)*MassIC_Mixture/Species(iSpec)%MassIC/TotalDOFWeight + C_P = C_P + ((5. + (Xi_VibSpec(iSpec)+Xi_RotSpec(iSpec)))/2.) * BoltzmannConst / Species(iSpec)%MassIC * MassFraction(iSpec) END DO SELECT CASE(BGKMixtureModel) + ! Both cases are described in Pfeiffer et. al., Physics of Fluids 33, 036106 (2021), + ! "Multi-species modeling in the particle-based ellipsoidal statistical Bhatnagar-Gross-Krook method for monatomic gas species" + ! Extension to poyatomic mixtures according to F. Hild, M. Pfeiffer, "Multi-species modeling in the particle-based ellipsoidal + ! statistical Bhatnagar-Gross-Krook method including internal degrees of freedom", subitted to Phys. Fluids, August 2023 + CASE (1) ! Wilke's mixing rules DO iSpec = 1, nSpecies ! Dynamic viscosity per species + ! Omega = OmegaVHS + 0.5 IF ((nSpec(iSpec).GE.2).AND.(.NOT.ALMOSTZERO(u2Spec(iSpec)))) THEN ! Species temperature: Sufficient number of particles per species are available dynamicvisSpec(iSpec) = 30.*SQRT(Species(iSpec)%MassIC* BoltzmannConst*CollInf%Tref(iSpec,iSpec)/Pi) & @@ -630,14 +706,17 @@ SUBROUTINE CalcGasProperties(nSpec, dens, InnerDOF, totalWeightSpec, totalWeight *CollInf%Tref(iSpec,iSpec)**(CollInf%omega(iSpec,iSpec) + 0.5)*CellTemp**(-CollInf%omega(iSpec,iSpec) - 0.5)) END IF ! Thermal conductivity per species (Eucken's formula with a correction by Hirschfelder for the internal degrees of freedom) - IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN + IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN ! inner DOF + ! Istomin et. al., "Eucken correction in high-temperature gases with electronic excitation", J. Chem. Phys. 140, + ! 184311 (2014) thermalcondspec(iSpec) = 0.25 * (15. + 2. * (Xi_VibSpec(iSpec)+Xi_RotSpec(iSpec)) * 1.328) & * dynamicvisSpec(iSpec) * BoltzmannConst / Species(iSpec)%MassIC - ELSE + ELSE ! atoms thermalcondspec(iSpec) = 0.25 * 15. * dynamicvisSpec(iSpec) * BoltzmannConst / Species(iSpec)%MassIC END IF END DO Phi= 0.0 + ! Calculation of factor phi, depending on mass ratios and ratios of dynamic viscosities DO iSpec = 1, nSpecies DO jSpec = 1, nSpecies Phi(iSpec) = Phi(iSpec) & @@ -649,11 +728,13 @@ SUBROUTINE CalcGasProperties(nSpec, dens, InnerDOF, totalWeightSpec, totalWeight END DO dynamicvis = 0.0 thermalcond = 0.0 + ! Sum up dynamic viscosities and thermal conductivities of species DO iSpec = 1, nSpecies IF (nSpec(iSpec).EQ.0) CYCLE dynamicvis = dynamicvis + REAL(totalWeightSpec(iSpec)) * dynamicvisSpec(iSpec) / Phi(iSpec) thermalcond = thermalcond + REAL(totalWeightSpec(iSpec)) * thermalcondspec(iSpec) / Phi(iSpec) END DO + CASE(2) ! Collision integrals (VHS) DO iSpec = 1, nSpecies IF ((nSpec(iSpec).LT.2).OR.ALMOSTZERO(u2Spec(iSpec))) THEN @@ -663,27 +744,29 @@ SUBROUTINE CalcGasProperties(nSpec, dens, InnerDOF, totalWeightSpec, totalWeight END IF END DO CellTempSpec(nSpecies+1) = CellTemp - CALL CalcViscosityThermalCondColIntVHS(CellTempSpec(1:nSpecies+1), MolarFraction(1:nSpecies),dens, dynamicvis, thermalcond) + CALL CalcViscosityThermalCondColIntVHS(CellTempSpec(1:nSpecies+1), MolarFraction(1:nSpecies),dens, Xi_RotSpec, Xi_VibSpec, & + dynamicvis, thermalcond) END SELECT + ! Calculation of Prandtl number Prandtl = C_P*dynamicvis/thermalcond*PrandtlCorrection + IF(DSMC%CalcQualityFactors) BGK_ExpectedPrandtlNumber = BGK_ExpectedPrandtlNumber + Prandtl - A = u0ij - CALL DSYEV('N','U',3,A,3,W,Work,100,INFO) - Theta = u2 / 3. - nu = 1.-1./Prandtl - nu= MAX(nu,-Theta/(W(3)-Theta)) - Prandtl = 1./(1.-nu) + ! Calculation of relaxation frequency relaxfreq = Prandtl*dens*BoltzmannConst*CellTemp/dynamicvis -ELSE + +ELSE ! single species gas + ! Calculation of reference dynamic viscosity dynamicvis = 30.*SQRT(Species(1)%MassIC* BoltzmannConst*CollInf%Tref(1,1)/Pi) & / (4.*(4.- 2.*CollInf%omega(1,1)) * (6. - 2.*CollInf%omega(1,1))* CollInf%dref(1,1)**2.) - Prandtl =2.*(InnerDOF + 5.)/(2.*InnerDOF + 15.) - IF (BGKCollModel.EQ.1) THEN + ! Calculation of Prandtl number: Pr = cp*mu/K with inner DOF, atoms: Pr = 2/3 + Prandtl = 2.*(InnerDOF + 5.)/(2.*InnerDOF + 15.) + ! Calculation of relaxation frequency using the exponential ansatz of the viscosity mu + IF (BGKCollModel.EQ.1) THEN ! ESBGK: relaxfreq nu = Pr*n*kB*T/mu relaxfreq = Prandtl*dens*BoltzmannConst*CollInf%Tref(1,1)**(CollInf%omega(1,1) + 0.5) & /dynamicvis*CellTemp**(-CollInf%omega(1,1) +0.5) - ELSE + ELSE ! relaxfreq nu = n*kB*T/mu relaxfreq = dens*BoltzmannConst*CollInf%Tref(1,1)**(CollInf%omega(1,1) + 0.5) & /dynamicvis*CellTemp**(-CollInf%omega(1,1) +0.5) END IF @@ -691,9 +774,252 @@ SUBROUTINE CalcGasProperties(nSpec, dens, InnerDOF, totalWeightSpec, totalWeight END SUBROUTINE CalcGasProperties -SUBROUTINE DetermineRelaxPart(nPart, iPartIndx_Node, relaxfreq, dtCell, RotExpSpec, VibExpSpec, nRelax, nRotRelax, nVibRelax, & - nRotRelaxSpec, nVibRelaxSpec, iPartIndx_NodeRelax, iPartIndx_NodeRelaxTemp, iPartIndx_NodeRelaxRot, & - iPartIndx_NodeRelaxVib, vBulk, OldEnRot, OldEn) + +SUBROUTINE CalcTRelax(ERotSpec, Xi_RotSpec, Xi_VibSpec, EVibSpec, totalWeightSpec, totalWeight, totalWeight2, dtCell, CellTemp, & + TRotSpec, TVibSpec, relaxfreq, rotrelaxfreqSpec, vibrelaxfreqSpec, Xi_VibSpecNew, Xi_vib_DOF, nXiVibDOF, CellTempRel, TEqui, & + betaR, betaV) +!=================================================================================================================================== +!> Calculate the relaxation energies and temperatures +!=================================================================================================================================== +! MODULES +USE MOD_Particle_Vars ,ONLY: nSpecies +USE MOD_DSMC_Vars ,ONLY: PolyatomMolDSMC, SpecDSMC, DSMC +USE MOD_BGK_Vars ,ONLY: BGKDoVibRelaxation +USE MOD_Globals_Vars ,ONLY: BoltzmannConst +USE MOD_Globals ,ONLY: abort +! IMPLICIT VARIABLE HANDLING + IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +INTEGER, INTENT(IN) :: nXiVibDOF +REAL, INTENT(IN) :: TRotSpec(nSpecies), ERotSpec(nSpecies), Xi_RotSpec(nSpecies) +REAL, INTENT(IN) :: TVibSpec(nSpecies), EVibSpec(nSpecies), Xi_VibSpec(nSpecies) +REAL, INTENT(IN) :: totalWeightSpec(nSpecies), totalWeight, totalWeight2, CellTemp, dtCell +REAL, INTENT(IN) :: relaxfreq, rotrelaxfreqSpec(nSpecies), vibrelaxfreqSpec(nSpecies) +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +REAL, INTENT(OUT) :: Xi_vib_DOF(DSMC%NumPolyatomMolecs,nXiVibDOF), Xi_VibSpecNew(nSpecies), CellTempRel, TEqui +REAL, INTENT(OUT) :: betaR(nSpecies), betaV(nSpecies) +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER :: iSpec, iDOF, iPolyatMole, i, j +REAL :: RotFracSpec(nSpecies), VibFracSpec(nSpecies) +REAL :: ERotSpecMean(nSpecies), EVibSpecMean(nSpecies), ETransRelMean +REAL :: ERotTtransSpecMean(nSpecies), EVibTtransSpecMean(nSpecies) +REAL :: TEqui_Old, TEqui_Old2, TEquiNum, TEquiDenom, exparg +REAL :: eps_prec=1.0E-0 +!=================================================================================================================================== +! According to J. Mathiaud et. al., "An ES-BGK model for diatomic gases with correct relaxation rates for internal energies", +! European Journal of Mechanics - B/Fluids, 96, pp. 65-77, 2022 +! For implentation, see F. Hild, M. Pfeiffer, "Multi-species modeling in the particle-based ellipsoidal statistical Bhatnagar-Gross- +! Krook method including internal degrees of freedom", subitted to Phys. Fluids, August 2023 + +RotFracSpec=0.0; VibFracSpec=0.0; Xi_vib_DOF=0.0; Xi_VibSpecNew=0.0; betaR=1.0; betaV=1.0 +ERotSpecMean=0.0; ERotTtransSpecMean=0.0; EVibSpecMean=0.0; EVibTtransSpecMean=0.0; ETransRelMean=0.0; CellTempRel=0.0 + +DO iSpec = 1, nSpecies + IF(totalWeightSpec(iSpec).GT.0.) THEN + IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN ! molecules + ! Mean rotational energy per particle of a species + ERotSpecMean(iSpec) = ERotSpec(iSpec)/totalWeightSpec(iSpec) + ! Mean rotational energy per particle of a species for the mixture translational temperature, ERot(Ttrans) + ERotTtransSpecMean(iSpec) = CellTemp * Xi_RotSpec(iSpec) * BoltzmannConst / 2. + ! Calculate number of rotational relaxing molecules + RotFracSpec(iSpec) = totalWeightSpec(iSpec)*(rotrelaxfreqSpec(iSpec)/relaxfreq)*(1.-EXP(-relaxfreq*dtCell)) + + IF(BGKDoVibRelaxation) THEN + ! Mean vibrational energy per particle of a species + EVibSpecMean(iSpec) = EVibSpec(iSpec)/totalWeightSpec(iSpec) + IF(SpecDSMC(iSpec)%PolyatomicMol) THEN ! polyatomic + iPolyatMole = SpecDSMC(iSpec)%SpecToPolyArray + ! Loop over all vibrational DOF + DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF + ! Mean vibrational energy per particle of a species for the mixture translational temperature, EVib(Ttrans) + EVibTtransSpecMean(iSpec) = EVibTtransSpecMean(iSpec) + BoltzmannConst*PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF) & + / (EXP(PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)/CellTemp) - 1.) + END DO + ELSE ! diatomic + ! Mean vibrational energy per particle of a species for the mixture translational temperature, EVib(Ttrans) + EVibTtransSpecMean(iSpec) = BoltzmannConst * SpecDSMC(iSpec)%CharaTVib / (EXP(SpecDSMC(iSpec)%CharaTVib/CellTemp) - 1.) + END IF + ! Calculate number of vibrational relaxing molecules + VibFracSpec(iSpec) = totalWeightSpec(iSpec)*(vibrelaxfreqSpec(iSpec)/relaxfreq)*(1.-EXP(-relaxfreq*dtCell)) + END IF + + ! Mean translational energy per particle to satisfy the Landau-Teller equation + ETransRelMean = ETransRelMean + (3./2. * BoltzmannConst * CellTemp - (rotrelaxfreqSpec(iSpec)/relaxfreq) * & + (ERotTtransSpecMean(iSpec)-ERotSpecMean(iSpec))) * totalWeightSpec(iSpec)/totalWeight + IF (BGKDoVibRelaxation) THEN + ETransRelMean = ETransRelMean - (vibrelaxfreqSpec(iSpec)/relaxfreq)*(EVibTtransSpecMean(iSpec)-EVibSpecMean(iSpec)) * & + totalWeightSpec(iSpec)/totalWeight + END IF + ELSE ! atomic + ! Mean translational energy per particle to satisfy the Landau-Teller equation + ETransRelMean = ETransRelMean + 3./2. * BoltzmannConst * CellTemp * totalWeightSpec(iSpec)/totalWeight + END IF + END IF +END DO + +! Calculation of the cell temperature with ETransRelMean to satisfy the Landau-Teller equation +IF (ETransRelMean.GT.0.0) THEN + CellTempRel = 2. * ETransRelMean / (3. * BoltzmannConst) +ELSE + CellTempRel = CellTemp +END IF + +! Calculation of equilibrium temperature for relaxation and energy conservation +TEqui_Old = 0.0 +TEquiNum = 3.*(totalWeight - totalWeight2/totalWeight)*CellTemp +TEquiDenom = 3.*(totalWeight - totalWeight2/totalWeight) +! Sum up over all species +DO iSpec = 1, nSpecies + IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN + TEquiNum = TEquiNum + Xi_RotSpec(iSpec)*RotFracSpec(iSpec)*TRotSpec(iSpec) + TEquiDenom = TEquiDenom + Xi_RotSpec(iSpec)*RotFracSpec(iSpec) + IF(BGKDoVibRelaxation) THEN + TEquiNum = TEquiNum + Xi_VibSpec(iSpec)*VibFracSpec(iSpec)*TVibSpec(iSpec) + TEquiDenom = TEquiDenom + Xi_VibSpec(iSpec)*VibFracSpec(iSpec) + END IF + END IF +END DO +TEqui = TEquiNum/TEquiDenom + +i=0 +! Solving of equation system until accuracy eps_prec is reached +outerLoop: DO WHILE ( ABS( TEqui - TEqui_Old ) .GT. eps_prec ) + i = i + 1 + DO iSpec = 1, nSpecies + IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN + ! if difference small: equilibrium, no beta + IF (ABS(TRotSpec(iSpec)-TEqui).GT.1E-3) THEN + betaR(iSpec) = (TRotSpec(iSpec)-CellTemp)/(TRotSpec(iSpec)-TEqui) + IF (betaR(iSpec).LT.0.0) THEN + betaR(iSpec) = 1. + END IF + ! new calculation of number of rotational relaxing molecules with betaR + RotFracSpec(iSpec) = totalWeightSpec(iSpec)*(rotrelaxfreqSpec(iSpec)/relaxfreq)*(1.-EXP(-relaxfreq*dtCell))*betaR(iSpec) + END IF + IF(BGKDoVibRelaxation) THEN + ! if difference small: equilibrium, no beta + IF (ABS(TVibSpec(iSpec)-TEqui).GT.1E-3) THEN + betaV(iSpec) = (TVibSpec(iSpec)-CellTemp)/(TVibSpec(iSpec)-TEqui) + IF (betaV(iSpec).LT.0.0) THEN + betaV(iSpec) = 1. + END IF + ! new calculation of number of vibrational relaxing molecules + VibFracSpec(iSpec) = totalWeightSpec(iSpec)*(vibrelaxfreqSpec(iSpec)/relaxfreq)*(1.-EXP(-relaxfreq*dtCell))*betaV(iSpec) + END IF + + ! new calculation of the vibrational degrees of freedom per species + IF(SpecDSMC(iSpec)%PolyatomicMol) THEN + iPolyatMole = SpecDSMC(iSpec)%SpecToPolyArray + ! Loop over all vibrational degrees of freedom to calculate them using TEqui + DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF + exparg = PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)/TEqui + ! Check if the exponent is within the range of machine precision for calculation of vibrational degrees of freedom + IF(CHECKEXP(exparg)) THEN + IF(exparg.gt.0.) THEN ! positive overflow: exp -> inf + Xi_vib_DOF(iSpec,iDOF) = 2.*exparg/(EXP(exparg)-1.) + ELSE ! negative overflow: exp -> 0 + Xi_vib_DOF(iSpec,iDOF) = 2.*exparg/(-1.) + END IF ! exparg.gt.0. + ELSE + Xi_vib_DOF(iSpec,iDOF) = 0.0 + END IF ! CHECKEXP(exparg) + END DO + Xi_VibSpecNew(iSpec) = SUM(Xi_vib_DOF(iSpec,1:PolyatomMolDSMC(iPolyatMole)%VibDOF)) + ELSE ! diatomic + exparg = SpecDSMC(iSpec)%CharaTVib/TEqui + ! Check if the exponent is within the range of machine precision for calculation of vibrational degrees of freedom + IF(CHECKEXP(exparg)) THEN + Xi_VibSpecNew(iSpec) = 2.*SpecDSMC(iSpec)%CharaTVib/TEqui/(EXP(exparg)-1.) + ELSE + Xi_VibSpecNew(iSpec) = 0.0 + END IF ! CHECKEXP(exparg) + END IF + END IF + END IF + END DO + TEqui_Old = TEqui + TEqui_Old2 = TEqui + ! new calculation of equilibrium temperature with new RotFracSpec, new VibFracSpec, new VibDOF(TEqui) in denominator + TEquiNum = 3.*(totalWeight - totalWeight2/totalWeight)*CellTemp + TEquiDenom = 3.*(totalWeight - totalWeight2/totalWeight) + ! Sum up over all species + DO iSpec = 1, nSpecies + IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN + TEquiNum = TEquiNum + Xi_RotSpec(iSpec)*RotFracSpec(iSpec)*TRotSpec(iSpec) + TEquiDenom = TEquiDenom + Xi_RotSpec(iSpec)*RotFracSpec(iSpec) + IF(BGKDoVibRelaxation) THEN + TEquiNum = TEquiNum + Xi_VibSpec(iSpec)*VibFracSpec(iSpec)*TVibSpec(iSpec) + TEquiDenom = TEquiDenom + Xi_VibSpecNew(iSpec)*VibFracSpec(iSpec) + END IF + END IF + END DO + TEqui = TEquiNum/TEquiDenom + IF(BGKDoVibRelaxation) THEN + j=0 + ! accuracy eps_prec not reached yet + innerLoop: DO WHILE ( ABS( TEqui - TEqui_Old2 ) .GT. eps_prec ) + j=j+1 + ! mean value of old and new equilibrium temperature + TEqui = (TEqui + TEqui_Old2) * 0.5 + DO iSpec = 1, nSpecies + IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN + ! new calculation of the vibrational degrees of freedom per species + IF(SpecDSMC(iSpec)%PolyatomicMol) THEN + iPolyatMole = SpecDSMC(iSpec)%SpecToPolyArray + ! Loop over all vibrational degrees of freedom to calculate them using TEqui + DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF + exparg = PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)/TEqui + ! Check if the exponent is within the range of machine precision for calculation of vibrational degrees of freedom + IF(CHECKEXP(exparg)) THEN + IF(exparg.gt.0.) THEN ! positive overflow: exp -> inf + Xi_vib_DOF(iSpec,iDOF) = 2.*exparg/(EXP(exparg)-1.) + ELSE ! negative overflow: exp -> 0 + Xi_vib_DOF(iSpec,iDOF) = 2.*exparg/(-1.) + END IF ! exparg.gt.0. + ELSE + Xi_vib_DOF(iSpec,iDOF) = 0.0 + END IF ! CHECKEXP(exparg) + END DO + Xi_VibSpecNew(iSpec) = SUM(Xi_vib_DOF(iSpec,1:PolyatomMolDSMC(iPolyatMole)%VibDOF)) + ELSE ! diatomic + exparg = SpecDSMC(iSpec)%CharaTVib/TEqui + ! Check if the exponent is within the range of machine precision for calculation of vibrational degrees of freedom + IF(CHECKEXP(exparg)) THEN + Xi_VibSpecNew(iSpec) = 2.*SpecDSMC(iSpec)%CharaTVib/TEqui/(EXP(exparg)-1.) + ELSE + Xi_VibSpecNew(iSpec) = 0.0 + END IF ! CHECKEXP(exparg) + END IF + END IF + END DO + TEqui_Old2 = TEqui + ! new calculation of equilibrium temperature with new RotFracSpec, new VibFracSpec, new VibDOF(TEqui) in denominator + TEquiNum = 3.*(totalWeight - totalWeight2/totalWeight)*CellTemp + TEquiDenom = 3.*(totalWeight - totalWeight2/totalWeight) + ! Sum up over all species + DO iSpec = 1, nSpecies + IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN + TEquiNum = TEquiNum + Xi_RotSpec(iSpec)*RotFracSpec(iSpec)*TRotSpec(iSpec) + & + Xi_VibSpec(iSpec)*VibFracSpec(iSpec)*TVibSpec(iSpec) + TEquiDenom = TEquiDenom + Xi_RotSpec(iSpec)*RotFracSpec(iSpec) + Xi_VibSpecNew(iSpec)*VibFracSpec(iSpec) + END IF + END DO + TEqui = TEquiNum/TEquiDenom + IF (j.EQ.30) EXIT innerLoop + END DO innerLoop + END IF + IF (i.EQ.30) EXIT outerLoop +END DO outerLoop + +END SUBROUTINE CalcTRelax + + +SUBROUTINE DetermineRelaxPart(nPart, iPartIndx_Node, relaxfreq, dtCell, nRelax, nRotRelax, nVibRelax, & + RotRelaxWeightSpec, VibRelaxWeightSpec, iPartIndx_NodeRelax, iPartIndx_NodeRelaxTemp, iPartIndx_NodeRelaxRot, & + iPartIndx_NodeRelaxVib, vBulk, OldEnRot, OldEn, rotrelaxfreqSpec, vibrelaxfreqSpec, betaR, betaV) !=================================================================================================================================== !> Determine the number of particles undergoing a relaxation (including vibration and rotation) !=================================================================================================================================== @@ -707,52 +1033,71 @@ SUBROUTINE DetermineRelaxPart(nPart, iPartIndx_Node, relaxfreq, dtCell, RotExpSp !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES INTEGER, INTENT(IN) :: nPart, iPartIndx_Node(nPart) -REAL, INTENT(IN) :: relaxfreq, dtCell, RotExpSpec(nSpecies), VibExpSpec(nSpecies) +REAL, INTENT(IN) :: relaxfreq, dtCell, rotrelaxfreqSpec(nSpecies), vibrelaxfreqSpec(nSpecies) +REAL, INTENT(IN) :: betaR(nSpecies), betaV(nSpecies) !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES -INTEGER, INTENT(OUT) :: nRelax, iPartIndx_NodeRelax(nPart), iPartIndx_NodeRelaxTemp(nPart) -INTEGER, INTENT(OUT) :: iPartIndx_NodeRelaxRot(nPart), iPartIndx_NodeRelaxVib(nPart) -INTEGER, INTENT(OUT) :: nRotRelax, nVibRelax, nRotRelaxSpec(nSpecies), nVibRelaxSpec(nSpecies) -REAL, INTENT(OUT) :: vBulk(3), OldEnRot +INTEGER, INTENT(OUT) :: iPartIndx_NodeRelax(:), iPartIndx_NodeRelaxTemp(:) +INTEGER, INTENT(OUT) :: iPartIndx_NodeRelaxRot(:), iPartIndx_NodeRelaxVib(:) +INTEGER, INTENT(OUT) :: nRelax, nRotRelax, nVibRelax +REAL, INTENT(OUT) :: vBulk(3), OldEnRot, RotRelaxWeightSpec(nSpecies), VibRelaxWeightSpec(nSpecies) !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT-OUTPUT VARIABLES REAL, INTENT(INOUT) :: OldEn !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER :: iPart, nNotRelax, iSpec, iLoop -REAL :: ProbAddPartTrans, iRan, partWeight +INTEGER :: iPart, iSpec, iLoop, nNotRelax +REAL :: ProbAddPartTrans, iRan, partWeight, ProbAddPartRot(nSpecies), ProbAddPartVib(nSpecies) !=================================================================================================================================== -nVibRelaxSpec =0; nRotRelaxSpec =0; nRelax=0; nNotRelax=0; vBulk=0.0; nRotRelax=0; nVibRelax=0; OldEnRot=0.0 +VibRelaxWeightSpec=0.0; RotRelaxWeightSpec=0.0; nRelax=0; nNotRelax=0; vBulk=0.0; nRotRelax=0; nVibRelax=0; OldEnRot=0.0 + +! Calculate probability of relaxation of a particle towards the target distribution function ProbAddPartTrans = 1.-EXP(-relaxfreq*dtCell) +! Calculate probabilities of relaxation of a particle in the rotation and vibration +! See F. Hild, M. Pfeiffer, "Multi-species modeling in the particle-based ellipsoidal statistical Bhatnagar-Gross-Krook method +! including internal degrees of freedom", subitted to Phys. Fluids, August 2023 +ProbAddPartRot(:) = ProbAddPartTrans * rotrelaxfreqSpec(:)/relaxfreq*betaR(:) +ProbAddPartVib(:) = ProbAddPartTrans * vibrelaxfreqSpec(:)/relaxfreq*betaV(:) + +! Loop over all simulation particles DO iLoop = 1, nPart iPart = iPartIndx_Node(iLoop) iSpec = PartSpecies(iPart) partWeight = GetParticleWeight(iPart) CALL RANDOM_NUMBER(iRan) + ! Count particles that are undergoing a relaxation IF (ProbAddPartTrans.GT.iRan) THEN nRelax = nRelax + 1 iPartIndx_NodeRelax(nRelax) = iPart + ! Count particles that are not undergoing a relaxation ELSE nNotRelax = nNotRelax + 1 iPartIndx_NodeRelaxTemp(nNotRelax) = iPart + ! Sum up velocities of non-relaxing particles for bulk velocity vBulk(1:3) = vBulk(1:3) + PartState(4:6,iPart)*Species(iSpec)%MassIC*partWeight END IF + + ! For molecules: relaxation of inner DOF IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN - !Rotation + ! Rotation CALL RANDOM_NUMBER(iRan) - IF ((1.-RotExpSpec(iSpec)).GT.iRan) THEN + ! Count particles that are undergoing a relaxation, in total and per species + IF (ProbAddPartRot(iSpec).GT.iRan) THEN nRotRelax = nRotRelax + 1 - nRotRelaxSpec(iSpec) = nRotRelaxSpec(iSpec) + 1 + RotRelaxWeightSpec(iSpec) = RotRelaxWeightSpec(iSpec) + partWeight iPartIndx_NodeRelaxRot(nRotRelax) = iPart + ! Sum up total rotational energy OldEnRot = OldEnRot + PartStateIntEn(2,iPart) * partWeight END IF ! Vibration IF(BGKDoVibRelaxation) THEN CALL RANDOM_NUMBER(iRan) - IF ((1.-VibExpSpec(iSpec)).GT.iRan) THEN + ! Count particles that are undergoing a relaxation, in total and per species + IF (ProbAddPartVib(iSpec).GT.iRan) THEN nVibRelax = nVibRelax + 1 - nVibRelaxSpec(iSpec) = nVibRelaxSpec(iSpec) + 1 + VibRelaxWeightSpec(iSpec) = VibRelaxWeightSpec(iSpec) + partWeight iPartIndx_NodeRelaxVib(nVibRelax) = iPart + ! Sum up total vibrational energy of all relaxing particles, considering zero-point energy, and add to translational energy OldEn = OldEn + (PartStateIntEn(1,iPartIndx_NodeRelaxVib(nVibRelax)) - SpecDSMC(iSpec)%EZeroPoint) * partWeight END IF END IF @@ -761,14 +1106,15 @@ SUBROUTINE DetermineRelaxPart(nPart, iPartIndx_Node, relaxfreq, dtCell, RotExpSp END SUBROUTINE DetermineRelaxPart -SUBROUTINE RelaxInnerEnergy(nPart, nVibRelax, nRotRelax, iPartIndx_NodeRelaxVib, iPartIndx_NodeRelaxRot, nXiVibDOF, Xi_vib_DOF, Xi_VibSpec, & - Xi_RotSpec , TEqui, VibEnergyDOF, NewEnVib, NewEnRot) + +SUBROUTINE RelaxInnerEnergy(nVibRelax, nRotRelax, iPartIndx_NodeRelaxVib, iPartIndx_NodeRelaxRot, nXiVibDOF, Xi_vib_DOF, & + Xi_VibSpec, Xi_RotSpec, VibEnergyDOF, TEqui, NewEnVib, NewEnRot) !=================================================================================================================================== !> Determine the new rotational and vibrational energy of relaxing particles !=================================================================================================================================== ! MODULES USE MOD_Particle_Vars ,ONLY: PartSpecies, nSpecies -USE MOD_DSMC_Vars ,ONLY: SpecDSMC, PartStateIntEn, PolyatomMolDSMC +USE MOD_DSMC_Vars ,ONLY: SpecDSMC, PartStateIntEn, PolyatomMolDSMC, DSMC USE MOD_BGK_Vars ,ONLY: BGKDoVibRelaxation USE MOD_part_tools ,ONLY: GetParticleWeight USE MOD_Globals_Vars ,ONLY: BoltzmannConst @@ -776,10 +1122,11 @@ SUBROUTINE RelaxInnerEnergy(nPart, nVibRelax, nRotRelax, iPartIndx_NodeRelaxVib, IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES -INTEGER, INTENT(IN) :: nPart,nXiVibDOF -INTEGER, INTENT(IN) :: nVibRelax, nRotRelax, iPartIndx_NodeRelaxVib(nPart), iPartIndx_NodeRelaxRot(nPart) -REAL, INTENT(IN) :: Xi_vib_DOF(nXiVibDOF), TEqui, Xi_VibSpec(nSpecies), Xi_RotSpec(nSpecies) -REAL, INTENT(INOUT) :: NewEnVib, NewEnRot +INTEGER, INTENT(IN) :: nXiVibDOF +INTEGER, INTENT(IN) :: nVibRelax, nRotRelax, iPartIndx_NodeRelaxVib(nVibRelax), iPartIndx_NodeRelaxRot(nRotRelax) +REAL, INTENT(IN) :: Xi_vib_DOF(DSMC%NumPolyatomMolecs,nXiVibDOF), Xi_VibSpec(nSpecies), Xi_RotSpec(nSpecies) +REAL, INTENT(IN) :: TEqui +REAL, INTENT(INOUT) :: NewEnVib(nSpecies), NewEnRot(nSpecies) !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES REAL, INTENT(OUT) :: VibEnergyDOF(nVibRelax,nXiVibDOF) @@ -788,41 +1135,53 @@ SUBROUTINE RelaxInnerEnergy(nPart, nVibRelax, nRotRelax, iPartIndx_NodeRelaxVib, INTEGER :: iLoop, iPart, iDOF, iPolyatMole, iSpec REAL :: partWeight, iRan !=================================================================================================================================== -! VIB Relaxation NewEnVib = 0.0; NewEnRot=0.0 IF(BGKDoVibRelaxation) THEN + ! Loop over all particles undergoing a relaxation in the vibration DO iLoop = 1, nVibRelax iPart = iPartIndx_NodeRelaxVib(iLoop) iSpec = PartSpecies(iPart) partWeight = GetParticleWeight(iPart) + ! polyatomic, more than one vibrational DOF IF(SpecDSMC(iSpec)%PolyatomicMol) THEN - iPolyatMole = SpecDSMC(iSpec)%SpecToPolyArray - PartStateIntEn(1,iPart) = 0.0 - DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF - CALL RANDOM_NUMBER(iRan) - VibEnergyDOF(iLoop,iDOF) = - LOG(iRan)*Xi_vib_DOF(iDOF)/2.*TEqui*BoltzmannConst - PartStateIntEn(1,iPart) = PartStateIntEn(1,iPart)+VibEnergyDOF(iLoop,iDOF) - END DO + iPolyatMole = SpecDSMC(iSpec)%SpecToPolyArray + PartStateIntEn(1,iPart) = 0.0 + ! Sum up the new vibrational energy over all DOFs, see M. Pfeiffer et. al., "Extension of Particle-based BGK Models to + ! Polyatomic Species in Hypersonic Flow around a Flat-faced Cylinder", AIP Conference Proceedings 2132, 100001 (2019) + DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF + CALL RANDOM_NUMBER(iRan) + VibEnergyDOF(iLoop,iDOF) = - LOG(iRan)*Xi_vib_DOF(iPolyatMole,iDOF)/2.*TEqui*BoltzmannConst + PartStateIntEn(1,iPart) = PartStateIntEn(1,iPart)+VibEnergyDOF(iLoop,iDOF) + END DO + ! ELSE: diatomic, only one vibrational DOF, calculate new vibrational energy according to M. Pfeiffer, "Extending the particle + ! ellipsoidal statistical Bhatnagar-Gross-Krook method to diatomic molecules including quantized vibrational energies", + ! Phys. Fluids 30, 116103 (2018) ELSE CALL RANDOM_NUMBER(iRan) PartStateIntEn( 1,iPart) = -LOG(iRan)*Xi_VibSpec(iSpec)/2.*TEqui*BoltzmannConst END IF - NewEnVib = NewEnVib + PartStateIntEn(1,iPart) * partWeight + ! Sum up new vibrational energy per species + NewEnVib(iSpec) = NewEnVib(iSpec) + PartStateIntEn(1,iPart) * partWeight END DO END IF -! ROT Relaxation +! Loop over all particles undergoing a relaxation in the rotation DO iLoop = 1, nRotRelax iPart = iPartIndx_NodeRelaxRot(iLoop) iSpec = PartSpecies(iPart) partWeight = GetParticleWeight(iPart) CALL RANDOM_NUMBER(iRan) + ! Calculate new rotational energy according to M. Pfeiffer et. al., "Extension of Particle-based BGK Models to Polyatomic Species + ! in Hypersonic Flow around a Flat-faced Cylinder", AIP Conference Proceedings 2132, 100001 (2019) PartStateIntEn( 2,iPart) = -Xi_RotSpec(iSpec) / 2. * BoltzmannConst*TEqui*LOG(iRan) - NewEnRot = NewEnRot + PartStateIntEn( 2,iPart) * partWeight + ! Sum up new rotational energy per species + NewEnRot(iSpec) = NewEnRot(iSpec) + PartStateIntEn( 2,iPart) * partWeight END DO END SUBROUTINE RelaxInnerEnergy -SUBROUTINE SampleFromTargetDistr(nRelax, iPartIndx_NodeRelax, Prandtl, u2, u0ij, u2i, vBulkAll, CellTemp, vBulk) + +SUBROUTINE SampleFromTargetDistr(nRelax, iPartIndx_NodeRelax, Prandtl, u2, u0ij, u2i, vBulkAll, CellTempRel, CellTemp, vBulk, & + MassIC_Mixture) !=================================================================================================================================== !> Sample new particle velocities from the target distribution function, depending on the chosen model !=================================================================================================================================== @@ -831,12 +1190,13 @@ SUBROUTINE SampleFromTargetDistr(nRelax, iPartIndx_NodeRelax, Prandtl, u2, u0ij, USE MOD_BGK_Vars ,ONLY: BGKCollModel, ESBGKModel USE MOD_part_tools ,ONLY: GetParticleWeight USE MOD_Globals_Vars ,ONLY: BoltzmannConst +USE MOD_Globals ,ONLY: abort ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES INTEGER, INTENT(IN) :: nRelax, iPartIndx_NodeRelax(:) -REAL, INTENT(IN) :: Prandtl, u2, u0ij(3,3), u2i(3), vBulkAll(3), CellTemp +REAL, INTENT(IN) :: Prandtl, u2, u0ij(3,3), u2i(3), vBulkAll(3), CellTempRel, CellTemp, MassIC_Mixture REAL, INTENT(INOUT) :: vBulk(3) !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES @@ -845,92 +1205,115 @@ SUBROUTINE SampleFromTargetDistr(nRelax, iPartIndx_NodeRelax, Prandtl, u2, u0ij, INTEGER :: iPart, fillMa1, fillMa2, INFO, iLoop, iSpec REAL :: iRanPart(3, nRelax), A(3,3), KronDelta, SMat(3,3), W(3), Work(100), tempVelo(3), partWeight !=================================================================================================================================== +! According to M. Pfeiffer, "Particle-based fluid dynamics: Comparison of different Bhatnagar-Gross-Krook models and the direct +! simulation Monte Carlo method for hypersonic flows", Phys. Fluids 30, 106106 (2018) and F. Hild, M. Pfeiffer, "Multi-species +! modeling in the particle-based ellipsoidal statistical Bhatnagar-Gross-Krook method including internal degrees of freedom", +! subitted to Phys. Fluids, August 2023 IF (nRelax.GT.0) THEN SELECT CASE(BGKCollModel) - CASE (1) ! Ellipsoidal Statistical + CASE (1) ! Ellipsoidal Statistical BGK IF (ESBGKModel.EQ.1) THEN - !! Approximated Solution - DO fillMa1 =1, 3 - DO fillMa2 =fillMa1, 3 + ! Approximated solution + DO fillMa1 = 1, 3 + DO fillMa2 = fillMa1, 3 IF (fillMa1.EQ.fillMa2) THEN KronDelta = 1.0 ELSE KronDelta = 0.0 END IF - SMat(fillMa1, fillMa2)= KronDelta - (1.-Prandtl)/(2.*Prandtl) & + ! Fill symmetric transformation matrix SMat with anisotopic matrix A = SS + SMat(fillMa1, fillMa2) = KronDelta - (1.-Prandtl)/(2.*Prandtl) & *(3./u2*u0ij(fillMa1, fillMa2)-KronDelta) END DO END DO SMat(2,1)=SMat(1,2) SMat(3,1)=SMat(1,3) SMat(3,2)=SMat(2,3) + ! Generate random normals for the sampling of new velocities of all relaxing particles CALL BGK_BuildTransGaussNums(nRelax, iRanPart) ELSE - !! Exact Solution - DO fillMa1 =1, 3 - DO fillMa2 =fillMa1, 3 + DO fillMa1 = 1, 3 + DO fillMa2 = fillMa1, 3 IF (fillMa1.EQ.fillMa2) THEN KronDelta = 1.0 ELSE KronDelta = 0.0 END IF - A(fillMa1, fillMa2) = KronDelta - (1.-Prandtl)/Prandtl*(3.*u0ij(fillMa1, fillMa2)/u2 - KronDelta) + ! Fill anisotopic matrix A + A(fillMa1, fillMa2) = KronDelta*CellTempRel*BoltzmannConst/MassIC_Mixture - (1.-Prandtl)/Prandtl * & + (CellTempRel/CellTemp*u0ij(fillMa1, fillMa2) - KronDelta*CellTempRel*BoltzmannConst/MassIC_Mixture) END DO END DO IF (ESBGKModel.EQ.2) THEN + ! Exact solution + ! Compute eigenvalues and eigenvectors of matrix A --> output: W is the array that contains the eigenvalues, A then contains + ! the orthonormal eigenvectors of anisotropic matrix A CALL DSYEV('V','U',3,A,3,W,Work,100,INFO) SMat = 0.0 - IF (W(1).LT.0.0) THEN - W(1) = 0.0 - IF (W(2).LT.0) W(2) = 0.0 - END IF - IF (W(3).LT.0) THEN - W(3) = 0.0 - DO fillMa1 =1, 3 - DO fillMa2 =fillMa1, 3 - IF (fillMa1.EQ.fillMa2) THEN - KronDelta = 1.0 - ELSE - KronDelta = 0.0 - END IF - SMat(fillMa1, fillMa2)= KronDelta - (1.-Prandtl)/(2.*Prandtl) & - *(3./u2*u0ij(fillMa1, fillMa2)-KronDelta) - END DO - END DO - SMat(2,1)=SMat(1,2) - SMat(3,1)=SMat(1,3) - SMat(3,2)=SMat(2,3) + IF (W(3).LT.0.0) THEN + ! Due to ascending order of eigenvalues, all three eigenvalues are lower than zero here + ! Fallback to Maxwell BGK + SMat(1,1) = SQRT(BoltzmannConst*CellTempRel/MassIC_Mixture) + SMat(2,2) = SQRT(BoltzmannConst*CellTempRel/MassIC_Mixture) + SMat(3,3) = SQRT(BoltzmannConst*CellTempRel/MassIC_Mixture) ELSE + ! At least W(3) is not negative + ! Set negative eigenvalues to zero to get positive semidefinite matrix + IF (W(1).LT.0.0) THEN + W(1) = 0.0 + IF (W(2).LT.0.0) W(2) = 0.0 + END IF + ! SMat with square roots of the eigenvalues as diagonal elements SMat(1,1) = SQRT(W(1)) SMat(2,2) = SQRT(W(2)) SMat(3,3) = SQRT(W(3)) + ! Diagonalisation of anisotropic matrix, SMat is square root of anisotropic matrix SMat = MATMUL(A, SMat) SMat = MATMUL(SMat, TRANSPOSE(A)) END IF + ! Generate random normals for the sampling of new velocities of all relaxing particles CALL BGK_BuildTransGaussNums(nRelax, iRanPart) ELSE IF (ESBGKModel.EQ.3) THEN + ! Metropolis-Hastings A(2,1)=A(1,2) A(3,1)=A(1,3) A(3,2)=A(2,3) - CALL MetropolisES(nRelax, iRanPart, A) + CALL MetropolisES(nRelax, iRanPart, A*MassIC_Mixture, CellTempRel) END IF END IF - CASE (2) ! Shakov -! CALL MetropolisShakhov(nRelax, iRanPart, u2/3., u2i, Prandtl) + + CASE (2) ! Shakov BGK + ! Acceptance-rejection method CALL ARShakhov(nRelax, iRanPart, u2/3., u2i, Prandtl) + CASE (3) ! Standard BGK (Maxwell target distribution) + ! Generate random normals for the sampling of new velocities of all relaxing particles CALL BGK_BuildTransGaussNums(nRelax, iRanPart) END SELECT + + ! Loop over all particles undergoing a relaxation towards the target distribution function DO iLoop = 1, nRelax iPart = iPartIndx_NodeRelax(iLoop) iSpec = PartSpecies(iPart) - IF ((BGKCollModel.EQ.1).AND.(ESBGKModel.NE.3)) THEN - tempVelo(1:3) = SQRT(BoltzmannConst*CellTemp/Species(iSpec)%MassIC)*iRanPart(1:3,iLoop) + ! Calculation of new velocities of all particles + IF ((BGKCollModel.EQ.1).AND.(ESBGKModel.EQ.1)) THEN + ! Transformation of normalized thermal velocity vector tempVelo (sampled from a Maxwellian distribution) to a thermal velocity + ! vector sampled from the ESBGK target distribution function (anisotropic Gaussian distribution) + tempVelo(1:3) = SQRT(BoltzmannConst*CellTempRel/Species(iSpec)%MassIC)*iRanPart(1:3,iLoop) PartState(4:6,iPart) = vBulkAll(1:3) + MATMUL(SMat,tempVelo) + ELSE IF ((BGKCollModel.EQ.1).AND.(ESBGKModel.EQ.2)) THEN + tempVelo(1:3) = SQRT(MassIC_Mixture/Species(iSpec)%MassIC)*iRanPart(1:3,iLoop) + PartState(4:6,iPart) = vBulkAll(1:3) + MATMUL(SMat,tempVelo) + ELSE IF ((BGKCollModel.EQ.1).AND.(ESBGKModel.EQ.3)) THEN + ! New thermal velocity (in x,y,z) of particle with mass scaling multiplied by normal distributed random vector + PartState(4:6,iPart) = vBulkAll(1:3) + SQRT(1./Species(iSpec)%MassIC)*iRanPart(1:3,iLoop) ELSE - PartState(4:6,iPart) = vBulkAll(1:3) + SQRT(BoltzmannConst*CellTemp/Species(iSpec)%MassIC)*iRanPart(1:3,iLoop) + ! New thermal velocity (in x,y,z) of particle is sqrt(k_B*T/m) multiplied by normal distributed random vector + PartState(4:6,iPart) = vBulkAll(1:3) + SQRT(BoltzmannConst*CellTempRel/Species(iSpec)%MassIC)*iRanPart(1:3,iLoop) END IF partWeight = GetParticleWeight(iPart) + ! Sum up new velocities of relaxing particles for bulk velocity, velocities of non-relaxing particles already calculated in + ! subroutine DetermineRelaxPart vBulk(1:3) = vBulk(1:3) + PartState(4:6,iPart)*Species(iSpec)%MassIC*partWeight END DO END IF ! nRelax.GT.0 @@ -938,7 +1321,8 @@ SUBROUTINE SampleFromTargetDistr(nRelax, iPartIndx_NodeRelax, Prandtl, u2, u0ij, END SUBROUTINE SampleFromTargetDistr -SUBROUTINE EnergyConsVib(nPart, nVibRelax, nVibRelaxSpec, iPartIndx_NodeRelaxVib, NewEnVib, OldEn, nXiVibDOF, Xi_VibSpec, VibEnergyDOF, TEqui) +SUBROUTINE EnergyConsVib(nPart, totalWeight, nVibRelax, VibRelaxWeightSpec, iPartIndx_NodeRelaxVib, NewEnVib, OldEn, nXiVibDOF, & + VibEnergyDOF, Xi_VibSpec, TEqui) !=================================================================================================================================== !> Routine to ensure energy conservation when including vibrational degrees of freedom (continuous and quantized) !=================================================================================================================================== @@ -952,93 +1336,132 @@ SUBROUTINE EnergyConsVib(nPart, nVibRelax, nVibRelaxSpec, iPartIndx_NodeRelaxVib IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES -INTEGER, INTENT(IN) :: nPart,nXiVibDOF -INTEGER, INTENT(IN) :: nVibRelax, iPartIndx_NodeRelaxVib(nPart), nVibRelaxSpec(nSpecies) -REAL, INTENT(IN) :: NewEnVib, VibEnergyDOF(nVibRelax,nXiVibDOF), Xi_VibSpec(nSpecies), TEqui +INTEGER, INTENT(IN) :: nPart, nXiVibDOF +INTEGER, INTENT(IN) :: nVibRelax, iPartIndx_NodeRelaxVib(nPart) +REAL, INTENT(IN) :: VibRelaxWeightSpec(nSpecies), Xi_VibSpec(nSpecies), totalWeight +REAL, INTENT(IN) :: NewEnVib(nSpecies), VibEnergyDOF(nVibRelax,nXiVibDOF), TEqui!, EVibTtransSpecMean(nSpecies) REAL, INTENT(INOUT) :: OldEn !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES INTEGER :: iPart, iLoop, iDOF, iSpec, iQuant, iQuaMax, iPolyatMole -REAL :: alpha, partWeight, betaV, iRan, MaxColQua, Xi_VibTotal +REAL :: Xi_VibTotal, alpha(nSpecies), partWeight, betaV, iRan, MaxColQua !=================================================================================================================================== +! According to F. Hild, M. Pfeiffer, "Multi-species modeling in the particle-based ellipsoidal statistical Bhatnagar-Gross-Krook +! method including internal degrees of freedom", subitted to Phys. Fluids, August 2023 IF(BGKDoVibRelaxation) THEN - IF ((NewEnVib.GT.0.0).AND.(nVibRelax.GT.0)) THEN + ! Vibrational energy is positive for at least one species + there are vibrational relaxations + IF (ANY(NewEnVib.GT.0.0).AND.(nVibRelax.GT.0)) THEN + Xi_VibTotal = 0.0 + DO iSpec = 1, nSpecies + ! Sum of relaxing vibrational degrees of freedom + Xi_VibTotal = Xi_VibTotal + Xi_VibSpec(iSpec)*VibRelaxWeightSpec(iSpec) + END DO + ! Calculate scaling factor alpha per species + ! EVibTtransSpecMean(iSpec)*VibRelaxWeightSpec(iSpec) is energy that should be in vibration + DO iSpec = 1, nSpecies + IF (NewEnVib(iSpec).GT.0.0) THEN + alpha(iSpec) = OldEn/NewEnVib(iSpec)*(Xi_VibSpec(iSpec)*VibRelaxWeightSpec(iSpec)/(3.*(totalWeight-1.)+Xi_VibTotal)) + ELSE + alpha(iSpec) = 0. + END IF + END DO + ! Quantized vibrational energy IF (BGKUseQuantVibEn) THEN - alpha = OldEn/NewEnVib*(Xi_VibSpec(1)*nVibRelax/(3.*(nPart-1.)+Xi_VibSpec(1)*nVibRelax)) DO iLoop = 1, nVibRelax iPart = iPartIndx_NodeRelaxVib(iLoop) partWeight = GetParticleWeight(iPart) iSpec = PartSpecies(iPart) + ! Polyatomic molecules IF(SpecDSMC(iSpec)%PolyatomicMol) THEN PartStateIntEn(1,iPart) = 0.0 iPolyatMole = SpecDSMC(iSpec)%SpecToPolyArray + ! Loop over all vibrational DOF DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF - betaV = alpha*VibEnergyDOF(iLoop,iDOF)/(PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)*BoltzmannConst) + ! Energy per vibrational mode alpha*VibEnergyDOF is reformulated to a quantum number iQuant + betaV = alpha(iSpec)*VibEnergyDOF(iLoop,iDOF)/(PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)*BoltzmannConst) CALL RANDOM_NUMBER(iRan) iQuant = INT(betaV+iRan) + ! Check maximum vibrational quantum number IF(iQuant.GT.PolyatomMolDSMC(iPolyatMole)%MaxVibQuantDOF(iDOF)) iQuant=PolyatomMolDSMC(iPolyatMole)%MaxVibQuantDOF(iDOF) + ! Remaining energy negative, new quantum number needs to be calculated IF ((OldEn - iQuant*PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)*BoltzmannConst*partWeight).LT.0.0) THEN + ! Maximum quantum number MaxColQua = OldEn/(BoltzmannConst*PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)*partWeight) + ! OldEn < k_B*CharaTVibDOF --> iQuant < 1 IF (INT(MaxColQua).EQ.0) THEN iQuant = 0 ELSE CALL RANDOM_NUMBER(iRan) + ! Calculation of new iQuant iQuant = INT(-LOG(iRan)*TEqui/PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)) + ! Determine maximum quantum number iQuaMax = MIN(INT(MaxColQua)+1, PolyatomMolDSMC(iPolyatMole)%MaxVibQuantDOF(iDOF)) + ! Calculation of new iQuant as long as iQuant > maximum quantum number DO WHILE (iQuant.GE.iQuaMax) CALL RANDOM_NUMBER(iRan) iQuant = INT(-LOG(iRan)*TEqui/PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)) END DO END IF END IF + ! Sum up the vibrational energy over all vibrational DOF PartStateIntEn( 1,iPart) = PartStateIntEn( 1,iPart) & + iQuant*PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)*BoltzmannConst VibQuantsPar(iPart)%Quants(iDOF) = iQuant + ! Remaining OldEn for remaining particles OldEn = OldEn - iQuant*PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)*BoltzmannConst*partWeight END DO - PartStateIntEn( 1,iPart) = PartStateIntEn( 1,iPart) & - + SpecDSMC(iSpec)%EZeroPoint + ! Add zero-point energy + PartStateIntEn( 1,iPart) = PartStateIntEn( 1,iPart) + SpecDSMC(iSpec)%EZeroPoint ELSE ! Diatomic molecules - betaV = alpha*PartStateIntEn( 1,iPart)/(SpecDSMC(iSpec)%CharaTVib*BoltzmannConst) + ! Vibrational energy is reformulated to a quantum number iQuant + betaV = alpha(iSpec)*PartStateIntEn( 1,iPart)/(SpecDSMC(iSpec)%CharaTVib*BoltzmannConst) CALL RANDOM_NUMBER(iRan) iQuant = INT(betaV+iRan) + ! Check maximum vibrational quantum number IF (iQuant.GT.SpecDSMC(iSpec)%MaxVibQuant) iQuant = SpecDSMC(iSpec)%MaxVibQuant PartStateIntEn( 1,iPart) = (iQuant + DSMC%GammaQuant)*SpecDSMC(iSpec)%CharaTVib*BoltzmannConst + ! Remaining energy negative, new quantum number needs to be calculated IF ((OldEn - (PartStateIntEn( 1,iPart) - SpecDSMC(iSpec)%EZeroPoint)*partWeight).LT.0.0) THEN + ! Maximum quantum number MaxColQua = OldEn/(BoltzmannConst*SpecDSMC(iSpec)%CharaTVib*partWeight) + ! OldEn < k_B*CharaTVib --> iQuant < 1 IF (INT(MaxColQua).EQ.0) THEN iQuant = 0 ELSE CALL RANDOM_NUMBER(iRan) + ! Calculation of new iQuant iQuant = INT(-LOG(iRan)*TEqui/SpecDSMC(iSpec)%CharaTVib) + ! Determine maximum quantum number iQuaMax = MIN(INT(MaxColQua)+1, SpecDSMC(iSpec)%MaxVibQuant) + ! Calculation of new iQuant as long as iQuant > maximum quantum number DO WHILE (iQuant.GE.iQuaMax) CALL RANDOM_NUMBER(iRan) iQuant = INT(-LOG(iRan)*TEqui/SpecDSMC(iSpec)%CharaTVib) END DO END IF + ! Calculate vibrational energy including zero-point energy PartStateIntEn( 1,iPart) = (iQuant + DSMC%GammaQuant)*SpecDSMC(iSpec)%CharaTVib*BoltzmannConst END IF + ! Remaining OldEn for remaining particles OldEn = OldEn - (PartStateIntEn( 1,iPart) - SpecDSMC(iSpec)%EZeroPoint)*partWeight - END IF ! SpecDSMC(1)%PolyatomicMol + END IF END DO ELSE ! Continuous treatment of vibrational energy - Xi_VibTotal = 0.0 - DO iSpec = 1, nSpecies - Xi_VibTotal = Xi_VibTotal + Xi_VibSpec(iSpec)*nVibRelaxSpec(iSpec) - END DO - alpha = OldEn/NewEnVib*(Xi_VibTotal/(3.*(nPart-1.)+Xi_VibTotal)) DO iLoop = 1, nVibRelax iPart = iPartIndx_NodeRelaxVib(iLoop) iSpec = PartSpecies(iPart) partWeight = GetParticleWeight(iPart) - PartStateIntEn( 1,iPart) = alpha*PartStateIntEn( 1,iPart) + SpecDSMC(iSpec)%EZeroPoint + ! Scaling of vibrational energy with factor alpha + zero-point energy + PartStateIntEn( 1,iPart) = alpha(iSpec)*PartStateIntEn( 1,iPart) + SpecDSMC(iSpec)%EZeroPoint + ! Remaining OldEn for remaining particles OldEn = OldEn - (PartStateIntEn( 1,iPart) - SpecDSMC(iSpec)%EZeroPoint)*partWeight END DO END IF ! BGKUseQuantVibEn - ELSE IF (nVibRelax.GT.0) THEN ! Relaxation towards the vibrational ground-state (new state is simply the zero-point energy) + ! NewEnVib = 0 for all species, relaxation towards the vibrational ground-state (new state is simply the zero-point energy) + ELSE IF (nVibRelax.GT.0) THEN + ! Set zero-point energy as vibrational energy for all particles with vibrational relaxation DO iLoop = 1, nVibRelax iPart = iPartIndx_NodeRelaxVib(iLoop) iSpec = PartSpecies(iPart) @@ -1049,204 +1472,81 @@ SUBROUTINE EnergyConsVib(nPart, nVibRelax, nVibRelaxSpec, iPartIndx_NodeRelaxVib END SUBROUTINE EnergyConsVib -#ifdef WIP -SUBROUTINE ARGrads13(nPart, iRanPart, Vtherm, HeatVec, PressTens) -!=================================================================================================================================== -!> description -!=================================================================================================================================== -! MODULES -USE Ziggurat -! IMPLICIT VARIABLE HANDLING -IMPLICIT NONE -!----------------------------------------------------------------------------------------------------------------------------------- -! INPUT VARIABLES -INTEGER, INTENT(IN) :: nPart -REAL, INTENT(IN) :: HeatVec(3), Vtherm, PressTens(3,3) -!----------------------------------------------------------------------------------------------------------------------------------- -! OUTPUT VARIABLES -REAL, INTENT(OUT) :: iRanPart(:,:) -!----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES -REAL :: Vheat, V2, iRan, OldProb, Envelope, Envelope2, cMat, KronDelta -INTEGER :: iPart, fillMa1, fillMa2 -!=================================================================================================================================== -Envelope = MAX(ABS(HeatVec(1)),ABS(HeatVec(2)),ABS(HeatVec(3)))/Vtherm**(3./2.) -Envelope2 = MAX(ABS(PressTens(1,2)),ABS(PressTens(1,3)),ABS(PressTens(2,3)))/Vtherm -Envelope = 1.+3.*MAX(Envelope, Envelope2) - -DO iPart = 1, nPart - iRanPart(1,iPart) = rnor() - iRanPart(2,iPart) = rnor() - iRanPart(3,iPart) = rnor() - cMat = 0.0 - DO fillMa1 =1, 3 - DO fillMa2 =1, 3 - IF (fillMa1.EQ.fillMa2) THEN - KronDelta = 1.0 - ELSE - KronDelta = 0.0 - END IF - cMat = cMat + iRanPart(fillMa1,iPart)*iRanPart(fillMa2,iPart)*(PressTens(fillMa1,fillMa2)-KronDelta*Vtherm) - END DO - END DO -! cMat=cMat + iRanPart(1,iPart)*iRanPart(2,iPart)*PressTens(1,2) -! cMat=cMat + iRanPart(1,iPart)*iRanPart(3,iPart)*PressTens(1,3) -! cMat=cMat + iRanPart(2,iPart)*iRanPart(3,iPart)*PressTens(2,3) - V2 = iRanPart(1,iPart)*iRanPart(1,iPart) + iRanPart(2,iPart)*iRanPart(2,iPart) + iRanPart(3,iPart)*iRanPart(3,iPart) - Vheat = iRanPart(1,iPart)*HeatVec(1) + iRanPart(2,iPart)*HeatVec(2) + iRanPart(3,iPart)*HeatVec(3) - OldProb = (1. + cMat/(2.*Vtherm) + VHeat/(Vtherm**(3./2.))*(V2/5.-1.)) - CALL RANDOM_NUMBER(iRan) - DO WHILE (Envelope*iRan.GT.OldProb) - iRanPart(1,iPart) = rnor() - iRanPart(2,iPart) = rnor() - iRanPart(3,iPart) = rnor() - cMat = 0.0 - DO fillMa1 =1, 3 - DO fillMa2 =1, 3 - IF (fillMa1.EQ.fillMa2) THEN - KronDelta = 1.0 - ELSE - KronDelta = 0.0 - END IF - cMat = cMat + iRanPart(fillMa1,iPart)*iRanPart(fillMa2,iPart)*(PressTens(fillMa1,fillMa2)-KronDelta*Vtherm) - END DO - END DO -! cMat=cMat + iRanPart(1,iPart)*iRanPart(2,iPart)*PressTens(1,2) -! cMat=cMat + iRanPart(1,iPart)*iRanPart(3,iPart)*PressTens(1,3) -! cMat=cMat + iRanPart(2,iPart)*iRanPart(3,iPart)*PressTens(2,3) - V2 = iRanPart(1,iPart)*iRanPart(1,iPart) + iRanPart(2,iPart)*iRanPart(2,iPart) + iRanPart(3,iPart)*iRanPart(3,iPart) - Vheat = iRanPart(1,iPart)*HeatVec(1) + iRanPart(2,iPart)*HeatVec(2) + iRanPart(3,iPart)*HeatVec(3) - OldProb = (1. + cMat/(2.*Vtherm) + VHeat/(Vtherm**(3./2.))*(V2/5.-1.)) - CALL RANDOM_NUMBER(iRan) - END DO -END DO - -END SUBROUTINE ARGrads13 -SUBROUTINE ARChapEnsk(nPart, iRanPart, Vtherm, HeatVec, PressTens) +SUBROUTINE MetropolisES(nPart, iRanPart, A, CellTempRel) !=================================================================================================================================== -!> description +!> Sampling from ESBGK target distribution function by using a Metropolis-Hastings method !=================================================================================================================================== ! MODULES USE Ziggurat +USE MOD_Basis ,ONLY: INV33 +USE MOD_Globals_Vars ,ONLY: BoltzmannConst ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES INTEGER, INTENT(IN) :: nPart -REAL, INTENT(IN) :: HeatVec(3), Vtherm, PressTens(3,3) -!----------------------------------------------------------------------------------------------------------------------------------- -! OUTPUT VARIABLES -REAL, INTENT(OUT) :: iRanPart(:,:) -!----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES -REAL :: Vheat, V2, iRan, OldProb, Envelope, Envelope2, cMat, cPress -INTEGER :: iPart -!=================================================================================================================================== -Envelope = MAX(ABS(HeatVec(1)),ABS(HeatVec(2)),ABS(HeatVec(3)))/Vtherm**(3./2.) -Envelope2 = MAX(ABS(PressTens(1,2)),ABS(PressTens(1,3)),ABS(PressTens(2,3)))/Vtherm -Envelope = 1.+4.*MAX(Envelope, Envelope2) - -DO iPart = 1, nPart - iRanPart(1,iPart) = rnor() - iRanPart(2,iPart) = rnor() - iRanPart(3,iPart) = rnor() - cMat = 0.0 - cPress = 0.0 - cMat=cMat + iRanPart(1,iPart)*iRanPart(2,iPart)*PressTens(1,2) - cMat=cMat + iRanPart(1,iPart)*iRanPart(3,iPart)*PressTens(1,3) - cMat=cMat + iRanPart(2,iPart)*iRanPart(3,iPart)*PressTens(2,3) - cPress=cPress + (PressTens(1,1)-Vtherm)*(iRanPart(1,iPart)*iRanPart(1,iPart)-iRanPart(3,iPart)*iRanPart(3,iPart)) - cPress=cPress + (PressTens(2,2)-Vtherm)*(iRanPart(2,iPart)*iRanPart(2,iPart)-iRanPart(3,iPart)*iRanPart(3,iPart)) - V2 = iRanPart(1,iPart)*iRanPart(1,iPart) + iRanPart(2,iPart)*iRanPart(2,iPart) + iRanPart(3,iPart)*iRanPart(3,iPart) - Vheat = iRanPart(1,iPart)*HeatVec(1) + iRanPart(2,iPart)*HeatVec(2) + iRanPart(3,iPart)*HeatVec(3) - OldProb = (1. + cMat/Vtherm + cPress/(2.*Vtherm) + VHeat/(2.*Vtherm**(3./2.))*(V2/5.-1.)) - CALL RANDOM_NUMBER(iRan) - DO WHILE (Envelope*iRan.GT.OldProb) - iRanPart(1,iPart) = rnor() - iRanPart(2,iPart) = rnor() - iRanPart(3,iPart) = rnor() - cMat = 0.0 - cPress = 0.0 - cMat=cMat + iRanPart(1,iPart)*iRanPart(2,iPart)*PressTens(1,2) - cMat=cMat + iRanPart(1,iPart)*iRanPart(3,iPart)*PressTens(1,3) - cMat=cMat + iRanPart(2,iPart)*iRanPart(3,iPart)*PressTens(2,3) - cPress=cPress + (PressTens(1,1)-Vtherm)*(iRanPart(1,iPart)*iRanPart(1,iPart)-iRanPart(3,iPart)*iRanPart(3,iPart)) - cPress=cPress + (PressTens(2,2)-Vtherm)*(iRanPart(2,iPart)*iRanPart(2,iPart)-iRanPart(3,iPart)*iRanPart(3,iPart)) - V2 = iRanPart(1,iPart)*iRanPart(1,iPart) + iRanPart(2,iPart)*iRanPart(2,iPart) + iRanPart(3,iPart)*iRanPart(3,iPart) - Vheat = iRanPart(1,iPart)*HeatVec(1) + iRanPart(2,iPart)*HeatVec(2) + iRanPart(3,iPart)*HeatVec(3) - OldProb = (1. + cMat/Vtherm + cPress/(2.*Vtherm) + VHeat/(2.*Vtherm**(3./2.))*(V2/5.-1.)) - CALL RANDOM_NUMBER(iRan) - END DO -END DO - -END SUBROUTINE ARChapEnsk -#endif /*WIP*/ - -SUBROUTINE MetropolisES(nPart, iRanPart, A) -!=================================================================================================================================== -!> description -!=================================================================================================================================== -! MODULES -USE Ziggurat -USE MOD_Basis ,ONLY: INV33 -! IMPLICIT VARIABLE HANDLING -IMPLICIT NONE -!----------------------------------------------------------------------------------------------------------------------------------- -! INPUT VARIABLES -INTEGER, INTENT(IN) :: nPart -REAL, INTENT(IN) :: A(3,3) +REAL, INTENT(IN) :: A(3,3), CellTempRel !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES REAL, INTENT(OUT) :: iRanPart(:,:) !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -REAL :: iRanPartTemp(3), V2, iRan, NewProb, OldProb, NormProb +REAL :: iRanPartTemp(3), V2, iRan, NewProb, OldProb, NormProb, prefacor INTEGER :: iLoop, iPart, iRun LOGICAL :: Changed REAL :: AC(3), AInvers(3,3), detA !=================================================================================================================================== -iRanPart(1,1) = rnor() -iRanPart(2,1) = rnor() -iRanPart(3,1) = rnor() -CALL INV33(A,AInvers, detA) +! Generate normal distributed random vector as start vector for the thermal velocity +prefacor = SQRT(BoltzmannConst*CellTempRel) +iRanPart(1,1) = rnor()*prefacor +iRanPart(2,1) = rnor()*prefacor +iRanPart(3,1) = rnor()*prefacor +! Inverse matrix of A +CALL INV33(A, AInvers, detA) AC(1:3) = MATMUL(AInvers, iRanPart(1:3,1)) V2 = iRanPart(1,1)*AC(1) + iRanPart(2,1)*AC(2) + iRanPart(3,1)*AC(3) OldProb = EXP(-0.5*V2) -!Burn in -DO iLoop = 1, 35 !50 - iRanPartTemp(1) = rnor() - iRanPartTemp(2) = rnor() - iRanPartTemp(3) = rnor() +! Burn-in phase, 35 initial steps +DO iLoop = 1, 35 + ! Generate normal distributed random vector for the thermal velocity + iRanPartTemp(1) = rnor()*prefacor + iRanPartTemp(2) = rnor()*prefacor + iRanPartTemp(3) = rnor()*prefacor AC(1:3) = MATMUL(AInvers, iRanPartTemp(1:3)) V2 = iRanPartTemp(1)*AC(1) + iRanPartTemp(2)*AC(2) + iRanPartTemp(3)*AC(3) NewProb = EXP(-0.5*V2) NormProb = MIN(1.,NewProb/OldProb) CALL RANDOM_NUMBER(iRan) + ! Acceptance of new sample with probability NormProb IF (NormProb.GT.iRan) THEN iRanPart(1:3,1) = iRanPartTemp(1:3) OldProb = NewProb END IF END DO -! All the others +! Main phase, for all following particles DO iPart = 2, nPart + ! Normal distributed random vector from previous particle iRanPart(1,iPart) = iRanPart(1,iPart-1) iRanPart(2,iPart) = iRanPart(2,iPart-1) iRanPart(3,iPart) = iRanPart(3,iPart-1) iRun = 0 Changed = .FALSE. + ! For acception: velocity should be changed at least once and at least ten steps in the Markov chain should be taken DO WHILE ((iRun.LT.10).OR.(.NOT.Changed)) iRun = iRun + 1 - iRanPartTemp(1) = rnor() - iRanPartTemp(2) = rnor() - iRanPartTemp(3) = rnor() + ! Generate normal distributed random vector for the thermal velocity + iRanPartTemp(1) = rnor()*prefacor + iRanPartTemp(2) = rnor()*prefacor + iRanPartTemp(3) = rnor()*prefacor AC(1:3) = MATMUL(AInvers, iRanPartTemp(1:3)) V2 = iRanPartTemp(1)*AC(1) + iRanPartTemp(2)*AC(2) + iRanPartTemp(3)*AC(3) NewProb = EXP(-0.5*V2) NormProb = MIN(1.,NewProb/OldProb) CALL RANDOM_NUMBER(iRan) + ! Acceptance of new sample with probability NormProb, velocity is changed IF (NormProb.GT.iRan) THEN - Changed = .TRUE. + Changed = .TRUE. iRanPart(1:3,iPart) = iRanPartTemp(1:3) OldProb = NewProb END IF @@ -1255,9 +1555,10 @@ SUBROUTINE MetropolisES(nPart, iRanPart, A) END SUBROUTINE MetropolisES + SUBROUTINE ARShakhov(nPart, iRanPart, Vtherm, HeatVec, Prandtl) !=================================================================================================================================== -!> description +!> Acceptance-rejection method for sampling from the Shakhov distribution function !=================================================================================================================================== ! MODULES USE Ziggurat @@ -1275,33 +1576,38 @@ SUBROUTINE ARShakhov(nPart, iRanPart, Vtherm, HeatVec, Prandtl) REAL :: Vheat, V2, iRan, OldProb, Envelope INTEGER :: iPart !=================================================================================================================================== +! Calculate envelope function Envelope = MAX(ABS(HeatVec(1)),ABS(HeatVec(2)),ABS(HeatVec(3)))/Vtherm**(3./2.) Envelope = 1.+4.*Envelope +! Loop over all relaxing particles DO iPart = 1, nPart + ! Generate random normals iRanPart(1,iPart) = rnor() iRanPart(2,iPart) = rnor() iRanPart(3,iPart) = rnor() V2 = iRanPart(1,iPart)*iRanPart(1,iPart) + iRanPart(2,iPart)*iRanPart(2,iPart) + iRanPart(3,iPart)*iRanPart(3,iPart) Vheat = iRanPart(1,iPart)*HeatVec(1) + iRanPart(2,iPart)*HeatVec(2) + iRanPart(3,iPart)*HeatVec(3) - OldProb = (1. + (1.-Prandtl)*VHeat/(5.*Vtherm**(3./2.))*(V2/2.-5./2.)) + OldProb = (1. + (1.-Prandtl)*Vheat/(5.*Vtherm**(3./2.))*(V2/2.-5./2.)) CALL RANDOM_NUMBER(iRan) + ! Acception if Envelope*iRan < OldProb DO WHILE (Envelope*iRan.GT.OldProb) iRanPart(1,iPart) = rnor() iRanPart(2,iPart) = rnor() iRanPart(3,iPart) = rnor() V2 = iRanPart(1,iPart)*iRanPart(1,iPart) + iRanPart(2,iPart)*iRanPart(2,iPart) + iRanPart(3,iPart)*iRanPart(3,iPart) Vheat = iRanPart(1,iPart)*HeatVec(1) + iRanPart(2,iPart)*HeatVec(2) + iRanPart(3,iPart)*HeatVec(3) - OldProb = (1. + (1.-Prandtl)*VHeat/(5.*Vtherm**(3./2.))*(V2/2.-5./2.)) + OldProb = (1. + (1.-Prandtl)*Vheat/(5.*Vtherm**(3./2.))*(V2/2.-5./2.)) CALL RANDOM_NUMBER(iRan) END DO END DO END SUBROUTINE ARShakhov + SUBROUTINE BGK_BuildTransGaussNums(nPart, iRanPart) !=================================================================================================================================== -!> description +!> Generate normal distributed random vector for sampling of new velocities of all relaxing particles relaxing !=================================================================================================================================== ! MODULES USE Ziggurat @@ -1317,6 +1623,7 @@ SUBROUTINE BGK_BuildTransGaussNums(nPart, iRanPart) ! LOCAL VARIABLES INTEGER :: iLoop !=================================================================================================================================== +! Generate three normal distributed random values for all relaxing simulation particles DO iLoop = 1, nPart iRanPart(1,iLoop) = rnor() iRanPart(2,iLoop) = rnor() @@ -1325,401 +1632,20 @@ SUBROUTINE BGK_BuildTransGaussNums(nPart, iRanPart) END SUBROUTINE BGK_BuildTransGaussNums -SUBROUTINE CalcTEqui(nPart, CellTemp, TRot, TVib, Xi_Vib, Xi_Vib_old, RotExp, VibExp, & - TEqui, rotrelaxfreq, vibrelaxfreq, dtCell, DoVibRelaxIn) -!=================================================================================================================================== -! Calculation of the vibrational temperature (zero-point search) for polyatomic molecules -!=================================================================================================================================== -! MODULES -USE MOD_DSMC_Vars, ONLY: SpecDSMC -USE MOD_BGK_Vars, ONLY: BGKDoVibRelaxation -! IMPLICIT VARIABLE HANDLING -IMPLICIT NONE -!----------------------------------------------------------------------------------------------------------------------------------- -! INPUT VARIABLES -REAL, INTENT(IN) :: CellTemp, TRot, TVib, Xi_Vib_old, rotrelaxfreq, vibrelaxfreq, dtCell -INTEGER, INTENT(IN) :: nPart -LOGICAL, OPTIONAL, INTENT(IN) :: DoVibRelaxIn -!----------------------------------------------------------------------------------------------------------------------------------- -! OUTPUT VARIABLES -REAL, INTENT(OUT) :: Xi_vib, TEqui, RotExp, VibExp -!----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES -!----------------------------------------------------------------------------------------------------------------------------------- -REAL :: TEqui_Old, betaR, betaV, RotFrac, VibFrac, TEqui_Old2 -REAL :: eps_prec=1.0E-0 -REAL :: correctFac, correctFacRot, maxexp !, Xi_rel -LOGICAL :: DoVibRelax -!=================================================================================================================================== -IF (PRESENT(DoVibRelaxIn)) THEN - DoVibRelax = DoVibRelaxIn -ELSE - DoVibRelax = BGKDoVibRelaxation -END IF -maxexp = LOG(HUGE(maxexp)) -! Xi_rel = 2.*(2. - CollInf%omega(1,1)) -! correctFac = 1. + (2.*SpecDSMC(1)%CharaTVib / (CellTemp*(EXP(SpecDSMC(1)%CharaTVib / CellTemp)-1.)))**(2.) & -! * EXP(SpecDSMC(1)%CharaTVib /CellTemp) / (2.*Xi_rel) -! correctFacRot = 1. + 2./Xi_rel - -correctFac = 1. -correctFacRot = 1. -RotExp = exp(-rotrelaxfreq*dtCell/correctFacRot) -RotFrac = nPart*(1.-RotExp) -IF(DoVibRelax) THEN - VibExp = exp(-vibrelaxfreq*dtCell/correctFac) - VibFrac = nPart*(1.-VibExp) -ELSE - VibExp = 0.0 - VibFrac = 0.0 - Xi_vib = 0.0 -END IF -TEqui_Old = 0.0 -TEqui = (3.*(nPart-1.)*CellTemp+2.*RotFrac*TRot+Xi_Vib_old*VibFrac*TVib)/(3.*(nPart-1.)+2.*RotFrac+Xi_Vib_old*VibFrac) -DO WHILE ( ABS( TEqui - TEqui_Old ) .GT. eps_prec ) - IF (ABS(TRot-TEqui).LT.1E-3) THEN - RotExp = exp(-rotrelaxfreq*dtCell/correctFacRot) - ELSE - betaR = ((TRot-CellTemp)/(TRot-TEqui))*rotrelaxfreq*dtCell/correctFacRot - IF (-betaR.GT.0.0) THEN - RotExp = 0. - ELSE IF (betaR.GT.maxexp) THEN - RotExp = 0. - ELSE - RotExp = exp(-betaR) - END IF - END IF - RotFrac = nPart*(1.-RotExp) - IF(DoVibRelax) THEN - IF (ABS(TVib-TEqui).LT.1E-3) THEN - VibExp = exp(-vibrelaxfreq*dtCell/correctFac) - ELSE - betaV = ((TVib-CellTemp)/(TVib-TEqui))*vibrelaxfreq*dtCell/correctFac - IF (-betaV.GT.0.0) THEN - VibExp = 0. - ELSE IF (betaV.GT.maxexp) THEN - VibExp = 0. - ELSE - VibExp = exp(-betaV) - END IF - END IF - IF ((SpecDSMC(1)%CharaTVib/TEqui).GT.maxexp) THEN - Xi_Vib = 0.0 - ELSE - Xi_vib = 2.*SpecDSMC(1)%CharaTVib/TEqui/(EXP(SpecDSMC(1)%CharaTVib/TEqui)-1.) - END IF - VibFrac = nPart*(1.-VibExp) - END IF - TEqui_Old = TEqui - TEqui_Old2 = TEqui - TEqui = (3.*(nPart-1.)*CellTemp+2.*RotFrac*TRot+Xi_Vib_old*VibFrac*TVib)/(3.*(nPart-1.)+2.*RotFrac+Xi_Vib*VibFrac) - IF(DoVibRelax) THEN - DO WHILE( ABS( TEqui - TEqui_Old2 ) .GT. eps_prec ) - TEqui =(TEqui + TEqui_Old2)*0.5 - IF ((SpecDSMC(1)%CharaTVib/TEqui).GT.maxexp) THEN - Xi_Vib = 0.0 - ELSE - Xi_vib = 2.*SpecDSMC(1)%CharaTVib/TEqui/(EXP(SpecDSMC(1)%CharaTVib/TEqui)-1.) - END IF - TEqui_Old2 = TEqui - TEqui = (3.*(nPart-1.)*CellTemp+2.*RotFrac*TRot+Xi_Vib_old*VibFrac*TVib) / (3.*(nPart-1.)+2.*RotFrac+Xi_vib*VibFrac) - END DO - END IF -END DO - -END SUBROUTINE CalcTEqui - -SUBROUTINE CalcTEquiMulti(nPart, nSpec, CellTemp, TRotSpec, TVibSpec, Xi_VibSpec, Xi_Vib_oldSpec, RotExpSpec, VibExpSpec, & - TEqui, rotrelaxfreqSpec, vibrelaxfreqSpec, dtCell, DoVibRelaxIn) -!=================================================================================================================================== -! Calculation of the vibrational temperature (zero-point search) for polyatomic molecules -!=================================================================================================================================== -! MODULES -USE MOD_DSMC_Vars, ONLY: SpecDSMC -USE MOD_BGK_Vars, ONLY: BGKDoVibRelaxation -USE MOD_Particle_Vars, ONLY: nSpecies -! IMPLICIT VARIABLE HANDLING -IMPLICIT NONE -!----------------------------------------------------------------------------------------------------------------------------------- -! INPUT VARIABLES -REAL, INTENT(IN) :: CellTemp, TRotSpec(nSpecies), TVibSpec(nSpecies), Xi_Vib_oldSpec(nSpecies) -REAL, INTENT(IN) :: rotrelaxfreqSpec(nSpecies), vibrelaxfreqSpec(nSpecies), dtCell -INTEGER, INTENT(IN) :: nPart, nSpec(nSpecies) -LOGICAL, OPTIONAL, INTENT(IN) :: DoVibRelaxIn -!----------------------------------------------------------------------------------------------------------------------------------- -! OUTPUT VARIABLES -REAL, INTENT(OUT) :: Xi_VibSpec(nSpecies), TEqui, RotExpSpec(nSpecies), VibExpSpec(nSpecies) -!----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES -!----------------------------------------------------------------------------------------------------------------------------------- -REAL :: TEqui_Old, betaR, betaV, RotFracSpec(nSpecies), VibFracSpec(nSpecies), TEqui_Old2 -REAL :: eps_prec=1.0E-0 -REAL :: correctFac, correctFacRot, exparg, TEquiNumDof !, Xi_rel, -LOGICAL :: DoVibRelax -INTEGER :: iSpec -!=================================================================================================================================== -IF (PRESENT(DoVibRelaxIn)) THEN - DoVibRelax = DoVibRelaxIn -ELSE - DoVibRelax = BGKDoVibRelaxation -END IF -! Xi_rel = 2.*(2. - CollInf%omega(1,1)) -! correctFac = 1. + (2.*SpecDSMC(1)%CharaTVib / (CellTemp*(EXP(SpecDSMC(1)%CharaTVib / CellTemp)-1.)))**(2.) & -! * EXP(SpecDSMC(1)%CharaTVib /CellTemp) / (2.*Xi_rel) -! correctFacRot = 1. + 2./Xi_rel - -correctFac = 1. -correctFacRot = 1. -RotFracSpec = 0.0 -VibFracSpec = 0.0 -DO iSpec=1, nSpecies - IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN - RotExpSpec(iSpec) = exp(-rotrelaxfreqSpec(iSpec)*dtCell/correctFacRot) - RotFracSpec(iSpec) = nSpec(iSpec)*(1.-RotExpSpec(iSpec)) - IF(DoVibRelax) THEN - VibExpSpec(iSpec) = exp(-vibrelaxfreqSpec(iSpec)*dtCell/correctFac) - VibFracSpec(iSpec) = nSpec(iSpec)*(1.-VibExpSpec(iSpec)) - ELSE - VibExpSpec(iSpec) = 0.0 - VibFracSpec(iSpec) = 0.0 - Xi_VibSpec(iSpec) = 0.0 - END IF - END IF -END DO -TEqui_Old = 0.0 -TEqui = 3.*(nPart-1.)*CellTemp -TEquiNumDof = 3.*(nPart-1.) -DO iSpec=1, nSpecies - IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN - TEqui = TEqui + 2.*RotFracSpec(iSpec)*TRotSpec(iSpec)+Xi_Vib_oldSpec(iSpec)*VibFracSpec(iSpec)*TVibSpec(iSpec) - TEquiNumDof = TEquiNumDof + 2.*RotFracSpec(iSpec) + Xi_Vib_oldSpec(iSpec)*VibFracSpec(iSpec) - END IF -END DO -TEqui = TEqui / TEquiNumDof -DO WHILE ( ABS( TEqui - TEqui_Old ) .GT. eps_prec ) - DO iSpec = 1, nSpecies - IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN - IF (ABS(TRotSpec(iSpec)-TEqui).LT.1E-3) THEN - RotExpSpec(iSpec) = exp(-rotrelaxfreqSpec(iSpec)*dtCell/correctFacRot) - ELSE - betaR = ((TRotSpec(iSpec)-CellTemp)/(TRotSpec(iSpec)-TEqui))*rotrelaxfreqSpec(iSpec)*dtCell/correctFacRot - IF (-betaR.GT.0.0) THEN - RotExpSpec(iSpec) = 0. - ELSE IF (CHECKEXP(betaR)) THEN - RotExpSpec(iSpec) = exp(-betaR) - ELSE - RotExpSpec(iSpec) = 0. - END IF - END IF - RotFracSpec(iSpec) = nSpec(iSpec)*(1.-RotExpSpec(iSpec)) - IF(DoVibRelax) THEN - IF (ABS(TVibSpec(iSpec)-TEqui).LT.1E-3) THEN - VibExpSpec(iSpec) = exp(-vibrelaxfreqSpec(iSpec)*dtCell/correctFac) - ELSE - betaV = ((TVibSpec(iSpec)-CellTemp)/(TVibSpec(iSpec)-TEqui))*vibrelaxfreqSpec(iSpec)*dtCell/correctFac - IF (-betaV.GT.0.0) THEN - VibExpSpec(iSpec) = 0. - ELSE IF (CHECKEXP(betaV)) THEN - VibExpSpec(iSpec) = exp(-betaV) - ELSE - VibExpSpec(iSpec) = 0. - END IF - END IF - exparg = SpecDSMC(iSpec)%CharaTVib/TEqui - IF(CHECKEXP(exparg))THEN - Xi_VibSpec(iSpec) = 2.*SpecDSMC(iSpec)%CharaTVib/TEqui/(EXP(exparg)-1.) - ELSE - Xi_VibSpec(iSpec) = 0.0 - END IF ! CHECKEXP(exparg) - VibFracSpec(iSpec) = nSpec(iSpec)*(1.-VibExpSpec(iSpec)) - END IF - END IF - END DO - TEqui_Old = TEqui - TEqui_Old2 = TEqui - - TEqui = 3.*(nPart-1.)*CellTemp - TEquiNumDof = 3.*(nPart-1.) - DO iSpec=1, nSpecies - IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN - TEqui = TEqui + 2.*RotFracSpec(iSpec)*TRotSpec(iSpec)+Xi_Vib_oldSpec(iSpec)*VibFracSpec(iSpec)*TVibSpec(iSpec) - TEquiNumDof = TEquiNumDof + 2.*RotFracSpec(iSpec) + Xi_VibSpec(iSpec)*VibFracSpec(iSpec) - END IF - END DO - TEqui = TEqui / TEquiNumDof - IF(DoVibRelax) THEN - DO WHILE( ABS( TEqui - TEqui_Old2 ) .GT. eps_prec ) - TEqui =(TEqui + TEqui_Old2)*0.5 - DO iSpec=1, nSpecies - IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN - exparg = SpecDSMC(iSpec)%CharaTVib/TEqui - IF(CHECKEXP(exparg))THEN - Xi_VibSpec(iSpec) = 2.*SpecDSMC(iSpec)%CharaTVib/TEqui/(EXP(exparg)-1.) - ELSE - Xi_VibSpec(iSpec) = 0.0 - END IF ! CHECKEXP(exparg) - END IF - END DO - TEqui_Old2 = TEqui - TEqui = 3.*(nPart-1.)*CellTemp - TEquiNumDof = 3.*(nPart-1.) - DO iSpec=1, nSpecies - IF((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN - TEqui = TEqui + 2.*RotFracSpec(iSpec)*TRotSpec(iSpec)+Xi_Vib_oldSpec(iSpec)*VibFracSpec(iSpec)*TVibSpec(iSpec) - TEquiNumDof = TEquiNumDof + 2.*RotFracSpec(iSpec) + Xi_VibSpec(iSpec)*VibFracSpec(iSpec) - END IF - END DO - TEqui = TEqui / TEquiNumDof - END DO - END IF -END DO -END SUBROUTINE CalcTEquiMulti - - -SUBROUTINE CalcTEquiPoly(nPart, CellTemp, TRot, TVib, nXiVibDOF, Xi_Vib_DOF, Xi_Vib_old, RotExp, VibExp, TEqui, rotrelaxfreq, vibrelaxfreq, & - dtCell, DoVibRelaxIn) -!=================================================================================================================================== -! Calculation of the vibrational temperature (zero-point search) for polyatomic molecules -!=================================================================================================================================== -! MODULES -USE MOD_DSMC_Vars, ONLY: SpecDSMC, PolyatomMolDSMC -USE MOD_BGK_Vars, ONLY: BGKDoVibRelaxation -! IMPLICIT VARIABLE HANDLING -IMPLICIT NONE -!----------------------------------------------------------------------------------------------------------------------------------- -! INPUT VARIABLES -REAL, INTENT(IN) :: CellTemp, TRot, TVib, Xi_Vib_old, rotrelaxfreq, vibrelaxfreq -INTEGER, INTENT(IN) :: nPart,nXiVibDOF -REAL, INTENT(IN) :: dtCell -LOGICAL, OPTIONAL, INTENT(IN) :: DoVibRelaxIn -!----------------------------------------------------------------------------------------------------------------------------------- -! OUTPUT VARIABLES -REAL, INTENT(OUT) :: Xi_vib_DOF(nXiVibDOF), TEqui, RotExp, VibExp -!----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES -!----------------------------------------------------------------------------------------------------------------------------------- -REAL :: TEqui_Old, betaR, betaV, RotFrac, VibFrac, Xi_Rot, TEqui_Old2, exparg -REAL :: eps_prec=1.0 -REAL :: correctFac, correctFacRot -INTEGER :: iDOF, iPolyatMole -LOGICAL :: DoVibRelax -!=================================================================================================================================== -IF (PRESENT(DoVibRelaxIn)) THEN - DoVibRelax = DoVibRelaxIn -ELSE - DoVibRelax = BGKDoVibRelaxation -END IF - -Xi_Rot = SpecDSMC(1)%Xi_Rot -iPolyatMole = SpecDSMC(1)%SpecToPolyArray -! Xi_rel = 2.*(2. - CollInf%omega(1,1)) -! correctFac = 0.0 -! DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF -! correctFac = correctFac & -! + (2.*PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF) / (CellTemp & -! *(EXP(PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF) / CellTemp)-1.)))**(2.) & -! * EXP(PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF) / CellTemp) / 2. -! END DO -! correctFac = 1. + correctFac/Xi_rel -! correctFacRot = 1. + Xi_Rot/Xi_rel - -correctFac = 1. -correctFacRot = 1. -RotExp = exp(-rotrelaxfreq*dtCell/correctFacRot) -RotFrac = nPart*(1.-RotExp) -IF(DoVibRelax) THEN - VibExp = exp(-vibrelaxfreq*dtCell/correctFac) - VibFrac = nPart*(1.-VibExp) -ELSE - VibExp = 0.0 - VibFrac = 0.0 - Xi_vib_DOF = 0.0 -END IF -TEqui_Old = 0.0 -TEqui = (3.*(nPart-1.)*CellTemp+2.*RotFrac*TRot+Xi_Vib_old*VibFrac*TVib)/(3.*(nPart-1.)+2.*RotFrac+Xi_Vib_old*VibFrac) -DO WHILE ( ABS( TEqui - TEqui_Old ) .GT. eps_prec ) - IF (ABS(TRot-TEqui).LT.1E-3) THEN - RotExp = exp(-rotrelaxfreq*dtCell/correctFacRot) - ELSE - betaR = ((TRot-CellTemp)/(TRot-TEqui))*rotrelaxfreq*dtCell/correctFacRot - IF (-betaR.GT.0.0) THEN - RotExp = 0. - ELSE IF (CHECKEXP(betaR)) THEN - RotExp = exp(-betaR) - ELSE - RotExp = 0. - END IF - END IF - RotFrac = nPart*(1.-RotExp) - IF(DoVibRelax) THEN - IF (ABS(TVib-TEqui).LT.1E-3) THEN - VibExp = exp(-vibrelaxfreq*dtCell/correctFac) - ELSE - betaV = ((TVib-CellTemp)/(TVib-TEqui))*vibrelaxfreq*dtCell/correctFac - IF (-betaV.GT.0.0) THEN - VibExp = 0. - ELSEIF(CHECKEXP(betaV))THEN - VibExp = exp(-betaV) - ELSE - VibExp = 0. - END IF - END IF - DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF - exparg = PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)/TEqui - IF(CHECKEXP(exparg))THEN - IF(exparg.gt.0.)THEN ! positive overflow: exp -> inf - Xi_vib_DOF(iDOF) = 2.*exparg/(EXP(exparg)-1.) - ELSE ! negative overflow: exp -> 0 - Xi_vib_DOF(iDOF) = 2.*exparg/(-1.) - END IF ! exparg.gt.0. - ELSE - Xi_vib_DOF(iDOF) = 0.0 - END IF ! CHECKEXP(exparg) - END DO - VibFrac = nPart*(1.-VibExp) - END IF - TEqui_Old = TEqui - TEqui_Old2 = TEqui - TEqui = (3.*(nPart-1.)*CellTemp+2.*RotFrac*TRot+Xi_Vib_old*VibFrac*TVib) & - / (3.*(nPart-1.)+2.*RotFrac+SUM(Xi_vib_DOF(1:PolyatomMolDSMC(iPolyatMole)%VibDOF))*VibFrac) - IF(DoVibRelax) THEN - DO WHILE( ABS( TEqui - TEqui_Old2 ) .GT. eps_prec ) - TEqui =(TEqui + TEqui_Old2)*0.5 - DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF - exparg = PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)/TEqui - IF(CHECKEXP(exparg))THEN - IF(exparg.gt.0.)THEN ! positive overflow: exp -> inf - Xi_vib_DOF(iDOF) = 2.*exparg/(EXP(exparg)-1.) - ELSE ! negative overflow: exp -> 0 - Xi_vib_DOF(iDOF) = 2.*exparg/(-1.) - END IF ! exparg.gt.0. - ELSE - Xi_vib_DOF(iDOF) = 0.0 - END IF ! CHECKEXP(exparg) - END DO - TEqui_Old2 = TEqui - TEqui = (3.*(nPart-1.)*CellTemp+2.*RotFrac*TRot+Xi_Vib_old*VibFrac*TVib) & - / (3.*(nPart-1.)+2.*RotFrac+SUM(Xi_vib_DOF(1:PolyatomMolDSMC(iPolyatMole)%VibDOF))*VibFrac) - END DO - END IF -END DO - -END SUBROUTINE CalcTEquiPoly - -SUBROUTINE CalcViscosityThermalCondColIntVHS(CellTemp, Xi, dens, Visc, ThermalCond) +SUBROUTINE CalcViscosityThermalCondColIntVHS(CellTemp, Xi, dens, Xi_RotSpec, Xi_VibSpec, Visc, ThermalCond) !=================================================================================================================================== !> Determination of the mixture viscosity and thermal conductivity using collision integrals (derived for the Variable Hard !> Sphere model). Solving an equation system depending on the number of species. !=================================================================================================================================== ! MODULES -USE MOD_DSMC_Vars, ONLY : CollInf +USE MOD_DSMC_Vars, ONLY : CollInf, SpecDSMC USE MOD_Globals_Vars, ONLY : BoltzmannConst USE MOD_Particle_Vars, ONLY : Species, nSpecies IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES -REAL, INTENT(IN) :: CellTemp(nSpecies+1), Xi(nSpecies), dens +REAL, INTENT(IN) :: CellTemp(nSpecies+1), Xi(nSpecies), dens, Xi_RotSpec(nSpecies), Xi_VibSpec(nSpecies) !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES REAL, INTENT(OUT) :: Visc,ThermalCond @@ -1728,29 +1654,55 @@ SUBROUTINE CalcViscosityThermalCondColIntVHS(CellTemp, Xi, dens, Visc, ThermalCo !----------------------------------------------------------------------------------------------------------------------------------- REAL :: Sigma_11, Sigma_22, B_12(nSpecies,nSpecies), A_12(nSpecies,nSpecies), InteractDiam, cv, DiffCoef(nSpecies, nSpecies) REAL :: Mass, ViscSpec(nSpecies), ThermalCondSpec(nSpecies), TVHS, omegaVHS, E_12, CellTemptmp +REAL :: ThermalCondSpec_Vib(nSpecies), ThermalCondSpec_Rot(nSpecies), cv_rot, cv_vib, rhoSpec +REAL :: Xj_Dij(nSpecies,nSpecies), Xi_Dij_tot REAL :: ViscMat(nSpecies, nSpecies), RHSSolve(nSpecies), m0, pressure INTEGER :: iSpec, jSpec, kSpec, IPIV(nSpecies), info_dgesv !=================================================================================================================================== -ViscSpec = 0.; ThermalCondSpec = 0.; DiffCoef =0.; A_12 = 0.; B_12 = 0. +ViscSpec = 0.; ThermalCondSpec = 0.; ThermalCondSpec_Vib = 0.; ThermalCondSpec_Rot = 0.; DiffCoef =0.; A_12 = 0.; B_12 = 0. +Xj_Dij = 0.; cv_rot = 0.; cv_vib = 0.; E_12 = 0. +! Loop over all species combinations DO iSpec = 1, nSpecies + IF (Xi(iSpec).LE.0.0) CYCLE + ! Calculate cv with rotational and vibrational degrees of freedom + IF ((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN + cv_rot = (Xi_RotSpec(iSpec)*BoltzmannConst)/(2.*Species(iSpec)%MassIC) + cv_vib = (Xi_VibSpec(iSpec)*BoltzmannConst)/(2.*Species(iSpec)%MassIC) + END IF DO jSpec = iSpec, nSpecies + IF (Xi(jSpec).LE.0.0) CYCLE + ! Interaction parameters InteractDiam = CollInf%dref(iSpec,jSpec) - Mass = Species(iSpec)%MassIC*Species(jSpec)%MassIC/(Species(iSpec)%MassIC + Species(jSpec)%MassIC) + Mass = Species(iSpec)%MassIC*Species(jSpec)%MassIC/(Species(iSpec)%MassIC + Species(jSpec)%MassIC) ! reduced mass TVHS = CollInf%Tref(iSpec,jSpec) omegaVHS = CollInf%omega(iSpec,jSpec) IF (iSpec.EQ.jSpec) THEN - CellTemptmp = CellTemp(iSpec) + CellTemptmp = CellTemp(iSpec) ! Species temperature or cell temperature for nSpec<2 or u2spec=0 ELSE - CellTemptmp = CellTemp(nSpecies+1) + CellTemptmp = CellTemp(nSpecies+1) ! Cell temperature END IF - Sigma_22 = CalcSigma_22VHS(CellTemptmp,InteractDiam,Mass,TVHS, omegaVHS) + ! Calculation of collision integral Sigma_22 + CALL CalcSigma_22VHS(CellTemptmp, InteractDiam, Mass, TVHS, omegaVHS, Sigma_22) IF (iSpec.EQ.jSpec) THEN - cv= 3./2.*BoltzmannConst/(2.*Mass) + cv= 3./2.*BoltzmannConst/(2.*Mass) ! DOF = 3, translational part + ! Calculation of the viscosity and thermal conductivity + ! S. Chapman and T.G. Cowling, "The mathematical Theory of Non-Uniform Gases", Cambridge University Press, 1970, S. 160 ViscSpec(iSpec) = (5./8.)*(BoltzmannConst*CellTemp(iSpec))/Sigma_22 ThermalCondSpec(iSpec) = (25./16.)*(cv*BoltzmannConst*CellTemp(iSpec))/Sigma_22 - !ThermalCondSpec(iSpec) = (15./4.)*BoltzmannConst/(2.*Mass)*ViscSpec(iSpec) + ! Results in the same as ThermalCondSpec(iSpec) = (15./4.)*BoltzmannConst/(2.*Mass)*ViscSpec(iSpec) + ! Additional calculation of Sigma_11VHS and the diffusion coefficient for molecular species + IF ((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN + CALL CalcSigma_11VHS(CellTemp(nSpecies+1), InteractDiam, Mass, TVHS, omegaVHS, Sigma_11) + E_12 = BoltzmannConst*CellTemp(nSpecies+1)/(8.*Species(iSpec)%MassIC*Species(jSpec)%MassIC & + /(Species(iSpec)%MassIC+Species(jSpec)%MassIC)**2.*Sigma_11) + DiffCoef(iSpec,jSpec) = 3.*E_12/(2.*(Species(iSpec)%MassIC+Species(jSpec)%MassIC)*dens) + END IF ELSE - CALL CalcSigma_11VHS(CellTemp(nSpecies+1),InteractDiam,Mass,TVHS, omegaVHS, Sigma_11) + ! Calculation of collision integral Sigma_11 + CALL CalcSigma_11VHS(CellTemp(nSpecies+1), InteractDiam, Mass, TVHS, omegaVHS, Sigma_11) + ! Parameters for calculation of contribution of species to mixture transport coefficients + ! Pfeiffer et. al., Physics of Fluids 33, 036106 (2021), "Multi-species modeling in the particle-based ellipsoidal + ! statistical Bhatnagar-Gross-Krook method for monatomic gas species" B_12(iSpec,jSpec) = (5.*GAMMA(4.-omegaVHS)-GAMMA(5.-omegaVHS))/(5.*GAMMA(3.-omegaVHS)) B_12(jSpec,iSpec) = B_12(iSpec,jSpec) A_12(iSpec,jSpec) = Sigma_22 / (5.*Sigma_11) @@ -1760,37 +1712,69 @@ SUBROUTINE CalcViscosityThermalCondColIntVHS(CellTemp, Xi, dens, Visc, ThermalCo DiffCoef(iSpec,jSpec) = 3.*E_12/(2.*(Species(iSpec)%MassIC+Species(jSpec)%MassIC)*dens) DiffCoef(jSpec,iSpec) = DiffCoef(iSpec,jSpec) END IF + Xj_Dij(iSpec,jSpec) = Xi(jSpec)/DiffCoef(iSpec,jSpec) + Xj_Dij(jSpec,iSpec) = Xj_Dij(iSpec,jSpec) END DO + IF ((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) THEN + ! Calculation of thermal conductivity of rotation and vibration for each molecular species + ! S. Chapman and T.G. Cowling, "The mathematical Theory of Non-Uniform Gases", Cambridge University Press, 1970, S. 254f + ! F. Hild, M. Pfeiffer, "Multi-species modeling in the particle-based ellipsoidal statistical Bhatnagar-Gross-Krook method + ! including internal degrees of freedom", subitted to Phys. Fluids, August 2023 + Xi_Dij_tot = SUM(Xj_Dij(iSpec,:)) + rhoSpec = dens * Species(iSpec)%MassIC * Xi(iSpec) + ThermalCondSpec_Rot(iSpec) = (rhoSpec*cv_rot/Xi_Dij_tot) + ThermalCondSpec_Vib(iSpec) = (rhoSpec*cv_vib/Xi_Dij_tot) + END IF END DO +! Calculate mixture viscosity by solving a system of linear equations with matrices +! Pfeiffer et. al., Physics of Fluids 33, 036106 (2021), +! "Multi-species modeling in the particle-based ellipsoidal statistical Bhatnagar-Gross-Krook method for monatomic gas species" +! S. Chapman and T.G. Cowling, "The mathematical Theory of Non-Uniform Gases", Cambridge University Press, 1970, S. 352 ViscMat = 0.0 DO iSpec = 1, nSpecies + IF (Xi(iSpec).LE.0.0) THEN + ViscMat(iSpec,iSpec) = 1. ! Ensure invertibility of ViscMat + CYCLE + END IF DO jSpec = 1, nSpecies + IF (Xi(jSpec).LE.0.0) CYCLE IF (iSpec.EQ.jSpec) THEN - ViscMat(iSpec, jSpec) = ViscMat(iSpec, jSpec) + Xi(iSpec)/ViscSpec(iSpec) + ViscMat(iSpec, jSpec) = Xi(iSpec)/ViscSpec(iSpec) DO kSpec = 1, nSpecies - IF(kSpec.EQ.iSpec) CYCLE + IF (Xi(kSpec).LE.0.0) CYCLE + IF (kSpec.EQ.iSpec) CYCLE ViscMat(iSpec, jSpec) = ViscMat(iSpec, jSpec) + 3.*Xi(kSpec) / ((Species(iSpec)%MassIC*dens & - + Species(kSpec)%MassIC*dens)*DiffCoef(iSpec,kSpec))*(2./3.+Species(kSpec)%MassIC/Species(iSpec)%MassIC*A_12(iSpec,kSpec)) + + Species(kSpec)%MassIC*dens)*DiffCoef(iSpec,kSpec))*(2./3.+Species(kSpec)%MassIC/Species(iSpec)%MassIC*A_12(iSpec,kSpec)) END DO ELSE - ViscMat(iSpec, jSpec) = ViscMat(iSpec, jSpec) - Xi(iSpec)*3. / ((Species(iSpec)%MassIC*dens & + ViscMat(iSpec, jSpec) = -Xi(iSpec)*3. / ((Species(iSpec)%MassIC*dens & + Species(jSpec)%MassIC*dens)*DiffCoef(iSpec,jSpec))*(2./3.-A_12(iSpec,jSpec)) END IF END DO - RHSSolve(iSpec) = Xi(iSpec) END DO +RHSSolve(:) = Xi(:) CALL DGESV(nSpecies, 1, ViscMat, nSpecies, IPIV, RHSSolve, nSpecies, info_dgesv) Visc = SUM(RHSSolve) +! Calculate mixture thermal conductivity by solving a system of linear equations with matrices +! Pfeiffer et. al., Physics of Fluids 33, 036106 (2021), +! "Multi-species modeling in the particle-based ellipsoidal statistical Bhatnagar-Gross-Krook method for monatomic gas species" +! S. Chapman and T.G. Cowling, "The mathematical Theory of Non-Uniform Gases", Cambridge University Press, 1970, S. 350f pressure = BoltzmannConst*dens*CellTemp(nSpecies+1) ViscMat = 0.0 DO iSpec = 1, nSpecies + IF (Xi(iSpec).LE.0.0) THEN + ViscMat(iSpec,iSpec) = 1. ! Ensure invertibility of ViscMat + CYCLE + END IF DO jSpec = 1, nSpecies + IF (Xi(jSpec).LE.0.0) CYCLE IF (iSpec.EQ.jSpec) THEN - ViscMat(iSpec, jSpec) = ViscMat(iSpec, jSpec) + Xi(iSpec)/ThermalCondSpec(iSpec) + ViscMat(iSpec, jSpec) = Xi(iSpec)/ThermalCondSpec(iSpec) DO kSpec = 1, nSpecies - IF(kSpec.EQ.iSpec) CYCLE + IF (Xi(kSpec).LE.0.0) CYCLE + IF (kSpec.EQ.iSpec) CYCLE m0 = Species(iSpec)%MassIC+Species(kSpec)%MassIC ViscMat(iSpec, jSpec) = ViscMat(iSpec, jSpec) + CellTemp(nSpecies+1)*Xi(kSpec)/(5.*pressure*DiffCoef(iSpec,kSpec)) & * (6.*Species(iSpec)%MassIC**2./m0**2.+(5.-4.*B_12(iSpec,kSpec))*Species(kSpec)%MassIC**2./m0**2. & @@ -1798,62 +1782,68 @@ SUBROUTINE CalcViscosityThermalCondColIntVHS(CellTemp, Xi, dens, Visc, ThermalCo END DO ELSE m0 = Species(iSpec)%MassIC+Species(jSpec)%MassIC - ViscMat(iSpec, jSpec) = ViscMat(iSpec, jSpec) - Xi(iSpec)*CellTemp(nSpecies+1) & - *(Species(iSpec)%MassIC*Species(jSpec)%MassIC/m0**2.)/(5.*pressure*DiffCoef(iSpec,jSpec)) & - *(11.-4.*B_12(iSpec,jSpec)-8.*A_12(iSpec,jSpec)) + ViscMat(iSpec, jSpec) = -Xi(iSpec)*CellTemp(nSpecies+1) * (Species(iSpec)%MassIC*Species(jSpec)%MassIC/m0**2.) & + /(5.*pressure*DiffCoef(iSpec,jSpec)) *(11.-4.*B_12(iSpec,jSpec)-8.*A_12(iSpec,jSpec)) END IF END DO - RHSSolve(iSpec) = Xi(iSpec) END DO +RHSSolve(:) = Xi(:) CALL DGESV(nSpecies, 1, ViscMat, nSpecies, IPIV, RHSSolve, nSpecies, info_dgesv) -ThermalCond = SUM(RHSSolve) +! Thermal conductivity from translation, rotation and vibration +ThermalCond = SUM(RHSSolve) + SUM(ThermalCondSpec_Rot) + SUM(ThermalCondSpec_Vib) END SUBROUTINE CalcViscosityThermalCondColIntVHS -SUBROUTINE CalcSigma_11VHS(CellTemp,Dref,Mass,Tref, omegaVHS, Sigma_11) +SUBROUTINE CalcSigma_11VHS(CellTemp, Dref, Mass, Tref, omegaVHS, Sigma_11) !=================================================================================================================================== !> !=================================================================================================================================== ! MODULES -USE MOD_Globals_Vars ,ONLY: Pi, BoltzmannConst +USE MOD_Globals_Vars ,ONLY: Pi, BoltzmannConst ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES -REAL, INTENT(IN) :: CellTemp,Dref,Mass,Tref, omegaVHS +REAL, INTENT(IN) :: CellTemp, Dref, Mass, Tref, omegaVHS REAL, INTENT(OUT) :: Sigma_11 !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -REAL :: Prefactor +REAL :: Prefactor !=================================================================================================================================== + ! See Stephani et. al., Physics of Fluids 24, 077101 (2012), + ! “Consistent treatment of transport properties for five-species air direct simulation Monte Carlo/Navier-Stokes applications” Prefactor = Pi/2.*Dref*Dref*SQRT(BoltzmannConst/(2.*Pi*Mass))*Tref**omegaVHS*GAMMA(3.-omegaVHS)/GAMMA(2.-omegaVHS) Sigma_11 = Prefactor*CellTemp**(0.5-omegaVHS) END SUBROUTINE CalcSigma_11VHS -REAL FUNCTION CalcSigma_22VHS(CellTemp,Dref,Mass,Tref, omegaVHS) + +SUBROUTINE CalcSigma_22VHS(CellTemp, Dref, Mass, Tref, omegaVHS, Sigma_22) !=================================================================================================================================== !> !=================================================================================================================================== ! MODULES -USE MOD_Globals_Vars ,ONLY: Pi, BoltzmannConst +USE MOD_Globals_Vars ,ONLY: Pi, BoltzmannConst ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES -REAL, INTENT(IN) :: CellTemp,Dref,Mass,Tref, omegaVHS +REAL, INTENT(IN) :: CellTemp, Dref, Mass, Tref, omegaVHS +REAL, INTENT(OUT) :: Sigma_22 !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -REAL :: Prefactor +REAL :: Prefactor !=================================================================================================================================== + ! See Stephani et. al., Physics of Fluids 24, 077101 (2012), + ! “Consistent treatment of transport properties for five-species air direct simulation Monte Carlo/Navier-Stokes applications” Prefactor = Pi/3.*Dref*Dref*SQRT(BoltzmannConst/(2.*Pi*Mass))*Tref**omegaVHS*GAMMA(4.-omegaVHS)/GAMMA(2.-omegaVHS) - CalcSigma_22VHS = Prefactor*CellTemp**(0.5-omegaVHS) + Sigma_22 = Prefactor*CellTemp**(0.5-omegaVHS) -END FUNCTION CalcSigma_22VHS +END SUBROUTINE CalcSigma_22VHS END MODULE MOD_BGK_CollOperator diff --git a/src/particles/bgk/bgk_init.f90 b/src/particles/bgk/bgk_init.f90 index 67938a2c6..1d8fcae12 100644 --- a/src/particles/bgk/bgk_init.f90 +++ b/src/particles/bgk/bgk_init.f90 @@ -68,7 +68,7 @@ SUBROUTINE DefineParametersBGK() 'cell refinement') CALL prms%CreateLogicalOption('Particles-BGK-MovingAverage', 'Enable a moving average of variables for the calculation '//& 'of the cell temperature for relaxation frequencies','.FALSE.') -CALL prms%CreateRealOption( 'Particles-BGK-MovingAverageFac', 'Use the moving average of moments M with '//& +CALL prms%CreateRealOption( 'Particles-BGK-MovingAverageFac', 'Use the moving average of moments M with '//& 'M^n+1=AverageFac*M+(1-AverageFac)*M^n','0.01') CALL prms%CreateRealOption( 'Particles-BGK-SplittingDens', 'Octree-refinement will only be performed above this number '//& 'density', '0.0') @@ -109,7 +109,6 @@ SUBROUTINE InitBGK() !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES INTEGER :: iSpec, iSpec2 -REAL :: delta_ij LOGICAL :: MoleculePresent !=================================================================================================================================== LBWRITE(UNIT_stdOut,'(A)') ' INIT BGK Solver...' @@ -118,49 +117,46 @@ SUBROUTINE InitBGK() DO iSpec=1, nSpecies IF ((SpecDSMC(iSpec)%InterID.EQ.2).OR.(SpecDSMC(iSpec)%InterID.EQ.20)) MoleculePresent = .TRUE. ALLOCATE(SpecBGK(iSpec)%CollFreqPreFactor(nSpecies)) + ! Calculation of the prefacor of the collision frequency per species + ! S. Chapman and T.G. Cowling, "The mathematical Theory of Non-Uniform Gases", Cambridge University Press, 1970, S. 87f DO iSpec2=1, nSpecies - IF (iSpec.EQ.iSpec2) THEN - delta_ij = 1.0 - ELSE - delta_ij = 0.0 - END IF - SpecBGK(iSpec)%CollFreqPreFactor(iSpec2)= 4.*(2.-delta_ij)*CollInf%dref(iSpec,iSpec2)**2.0 & + SpecBGK(iSpec)%CollFreqPreFactor(iSpec2)= 4.*CollInf%dref(iSpec,iSpec2)**2.0 & * SQRT(Pi*BoltzmannConst*CollInf%Tref(iSpec,iSpec2)*(Species(iSpec)%MassIC + Species(iSpec2)%MassIC) & /(2.*(Species(iSpec)%MassIC * Species(iSpec2)%MassIC)))/CollInf%Tref(iSpec,iSpec2)**(-CollInf%omega(iSpec,iSpec2) +0.5) END DO END DO -IF ((nSpecies.GT.1).AND.(ANY(SpecDSMC(:)%PolyatomicMol))) THEN - CALL abort(__STAMP__,' ERROR Multispec not implemented with polyatomic molecules!') -END IF BGKCollModel = GETINT('Particles-BGK-CollModel') IF ((nSpecies.GT.1).AND.(BGKCollModel.GT.1)) THEN - CALL abort(__STAMP__,' ERROR Multispec only with ESBGK model!') + CALL abort(__STAMP__,'ERROR Multispec only with ESBGK model!') END IF BGKMixtureModel = GETINT('Particles-BGK-MixtureModel') -! ESBGK options -ESBGKModel = GETINT('Particles-ESBGK-Model') ! 1: Approximative, 2: Exact, 3: MetropolisHastings +! ESBGK options for sampling: 1: Approximative, 2: Exact, 3: MetropolisHastings +ESBGKModel = GETINT('Particles-ESBGK-Model') + ! Coupled BGK with DSMC, use a number density as limit above which BGK is used, and below which DSMC is used CoupledBGKDSMC = GETLOGICAL('Particles-CoupledBGKDSMC') IF(CoupledBGKDSMC) THEN IF (DoVirtualCellMerge) THEN - CALL abort(__STAMP__,' Virtual cell merge not implemented for coupled DSMC-BGK simulations!') + CALL abort(__STAMP__,'Virtual cell merge not implemented for coupled DSMC-BGK simulations!') END IF BGKDSMCSwitchDens = GETREAL('Particles-BGK-DSMC-SwitchDens') ELSE IF(RadialWeighting%DoRadialWeighting) RadialWeighting%PerformCloning = .TRUE. END IF + ! Octree-based cell refinement, up to a certain number of particles DoBGKCellAdaptation = GETLOGICAL('Particles-BGK-DoCellAdaptation') IF(DoBGKCellAdaptation) THEN BGKMinPartPerCell = GETINT('Particles-BGK-MinPartsPerCell') IF(.NOT.DSMC%UseOctree) THEN DSMC%UseOctree = .TRUE. - IF(NGeo.GT.PP_N) CALL abort(__STAMP__,' Set PP_N to NGeo, otherwise the volume is not computed correctly.') + IF(NGeo.GT.PP_N) CALL abort(__STAMP__,'Set PP_N to NGeo, otherwise the volume is not computed correctly.') CALL DSMC_init_octree() END IF END IF BGKSplittingDens = GETREAL('Particles-BGK-SplittingDens') + ! Moving Average BGKMovingAverage = GETLOGICAL('Particles-BGK-MovingAverage') IF(BGKMovingAverage) THEN @@ -169,20 +165,16 @@ SUBROUTINE InitBGK() CALL BGK_init_MovingAverage() IF(nSpecies.GT.1) CALL abort(__STAMP__,'nSpecies >1 and molecules not implemented for BGK averaging!') END IF + IF(MoleculePresent) THEN ! Vibrational modelling BGKDoVibRelaxation = GETLOGICAL('Particles-BGK-DoVibRelaxation') BGKUseQuantVibEn = GETLOGICAL('Particles-BGK-UseQuantVibEn') - IF ((nSpecies.GT.1).AND.(BGKUseQuantVibEn)) THEN - CALL abort(& - __STAMP__& - ,' ERROR Multispec not implemented for quantized vibrational energy!') - END IF END IF IF(DSMC%CalcQualityFactors) THEN - ALLOCATE(BGK_QualityFacSamp(1:7,nElems)) - BGK_QualityFacSamp(1:7,1:nElems) = 0.0 + ALLOCATE(BGK_QualityFacSamp(1:9,nElems)) + BGK_QualityFacSamp(1:9,1:nElems) = 0.0 END IF BGKInitDone = .TRUE. @@ -249,6 +241,7 @@ SUBROUTINE FinalizeBGK() END SUBROUTINE FinalizeBGK + SUBROUTINE DeleteElemNodeAverage() !----------------------------------------------------------------------------------------------------------------------------------! ! Delete the pointer tree ElemNodeVol diff --git a/src/particles/bgk/bgk_main.f90 b/src/particles/bgk/bgk_main.f90 index 8d640f43c..d856dd6dc 100644 --- a/src/particles/bgk/bgk_main.f90 +++ b/src/particles/bgk/bgk_main.f90 @@ -48,6 +48,7 @@ SUBROUTINE BGK_DSMC_main(stage_opt) USE MOD_BGK_Vars ,ONLY: BGKMovingAverage,ElemNodeAveraging USE MOD_BGK_Vars ,ONLY: BGK_MeanRelaxFactor,BGK_MeanRelaxFactorCounter,BGK_MaxRelaxFactor,BGK_QualityFacSamp USE MOD_BGK_Vars ,ONLY: BGK_MaxRotRelaxFactor, BGK_PrandtlNumber, BGK_ExpectedPrandtlNumber +USE MOD_BGK_Vars ,ONLY: BGK_Viscosity, BGK_ThermalConductivity USE MOD_BGK_CollOperator ,ONLY: BGK_CollisionOperator USE MOD_DSMC ,ONLY: DSMC_main USE MOD_DSMC_Vars ,ONLY: DSMC, RadialWeighting @@ -128,10 +129,10 @@ SUBROUTINE BGK_DSMC_main(stage_opt) IF(DSMC%CalcQualityFactors) THEN BGK_MeanRelaxFactorCounter = 0; BGK_MeanRelaxFactor = 0.; BGK_MaxRelaxFactor = 0.; BGK_MaxRotRelaxFactor = 0. - BGK_PrandtlNumber=0.; BGK_ExpectedPrandtlNumber=0. + BGK_PrandtlNumber=0.; BGK_ExpectedPrandtlNumber=0.; BGK_Viscosity=0.; BGK_ThermalConductivity=0. END IF IF (BGKMovingAverage) THEN - CALL BGK_CollisionOperator(iPartIndx_Node, nPart, ElemVolume_Shared(CNElemID), ElemNodeAveraging(iElem)%Root%AverageValues(:)) + CALL BGK_CollisionOperator(iPartIndx_Node, nPart, ElemVolume_Shared(CNElemID), ElemNodeAveraging(iElem)%Root%AverageValues(:)) ELSE CALL BGK_CollisionOperator(iPartIndx_Node, nPart, ElemVolume_Shared(CNElemID)) END IF @@ -145,6 +146,8 @@ SUBROUTINE BGK_DSMC_main(stage_opt) BGK_QualityFacSamp(5,iElem) = BGK_QualityFacSamp(5,iElem) + BGK_MaxRotRelaxFactor BGK_QualityFacSamp(6,iElem) = BGK_QualityFacSamp(6,iElem) + BGK_PrandtlNumber BGK_QualityFacSamp(7,iElem) = BGK_QualityFacSamp(7,iElem) + BGK_ExpectedPrandtlNumber + BGK_QualityFacSamp(8,iElem) = BGK_QualityFacSamp(8,iElem) + BGK_Viscosity + BGK_QualityFacSamp(9,iElem) = BGK_QualityFacSamp(9,iElem) + BGK_ThermalConductivity END IF END IF END IF @@ -170,6 +173,7 @@ SUBROUTINE BGK_main(stage_opt) USE MOD_BGK_Vars ,ONLY: DoBGKCellAdaptation, BGKMovingAverage, ElemNodeAveraging USE MOD_BGK_Vars ,ONLY: BGK_MeanRelaxFactor,BGK_MeanRelaxFactorCounter,BGK_MaxRelaxFactor,BGK_QualityFacSamp USE MOD_BGK_Vars ,ONLY: BGK_MaxRotRelaxFactor, BGK_PrandtlNumber, BGK_ExpectedPrandtlNumber +USE MOD_BGK_Vars ,ONLY: BGK_Viscosity, BGK_ThermalConductivity USE MOD_BGK_CollOperator ,ONLY: BGK_CollisionOperator USE MOD_DSMC_Analyze ,ONLY: DSMCMacroSampling USE MOD_Particle_Mesh_Vars ,ONLY: ElemVolume_Shared @@ -266,7 +270,7 @@ SUBROUTINE BGK_main(stage_opt) IF(DSMC%CalcQualityFactors) THEN BGK_MeanRelaxFactorCounter = 0; BGK_MeanRelaxFactor = 0.; BGK_MaxRelaxFactor = 0.; BGK_MaxRotRelaxFactor = 0. - BGK_PrandtlNumber=0.; BGK_ExpectedPrandtlNumber=0. + BGK_PrandtlNumber=0.; BGK_ExpectedPrandtlNumber=0.; BGK_Viscosity=0.; BGK_ThermalConductivity=0. END IF IF (BGKMovingAverage) THEN @@ -284,6 +288,8 @@ SUBROUTINE BGK_main(stage_opt) BGK_QualityFacSamp(5,iElem) = BGK_QualityFacSamp(5,iElem) + BGK_MaxRotRelaxFactor BGK_QualityFacSamp(6,iElem) = BGK_QualityFacSamp(6,iElem) + BGK_PrandtlNumber BGK_QualityFacSamp(7,iElem) = BGK_QualityFacSamp(7,iElem) + BGK_ExpectedPrandtlNumber + BGK_QualityFacSamp(8,iElem) = BGK_QualityFacSamp(8,iElem) + BGK_Viscosity + BGK_QualityFacSamp(9,iElem) = BGK_QualityFacSamp(9,iElem) + BGK_ThermalConductivity END IF END IF END DO diff --git a/src/particles/bgk/bgk_vars.f90 b/src/particles/bgk/bgk_vars.f90 index f2312579e..b1b63a36b 100644 --- a/src/particles/bgk/bgk_vars.f90 +++ b/src/particles/bgk/bgk_vars.f90 @@ -48,6 +48,8 @@ MODULE MOD_BGK_Vars REAL :: BGK_MaxRotRelaxFactor REAL :: BGK_PrandtlNumber REAL :: BGK_ExpectedPrandtlNumber +REAL :: BGK_Viscosity +REAL :: BGK_ThermalConductivity TYPE tElemNodeAveraging TYPE (tNodeAverage), POINTER :: Root => null() diff --git a/src/particles/dsmc/dsmc_analyze.f90 b/src/particles/dsmc/dsmc_analyze.f90 index 9be9fc49e..b8972d88f 100644 --- a/src/particles/dsmc/dsmc_analyze.f90 +++ b/src/particles/dsmc/dsmc_analyze.f90 @@ -722,16 +722,20 @@ SUBROUTINE DSMC_output_calc(nVar,nVar_quality,nVarloc,DSMC_MacroVal) DSMC_MacroVal(nVarCount+2,iElem) = BGK_QualityFacSamp(6,iElem) / BGK_QualityFacSamp(2,iElem) ! Mean expected Prandtl number DSMC_MacroVal(nVarCount+3,iElem) = BGK_QualityFacSamp(7,iElem) / BGK_QualityFacSamp(2,iElem) + ! Mean viscosity + DSMC_MacroVal(nVarCount+4,iElem) = BGK_QualityFacSamp(8,iElem) / BGK_QualityFacSamp(2,iElem) + ! Mean thermal conductivity + DSMC_MacroVal(nVarCount+5,iElem) = BGK_QualityFacSamp(9,iElem) / BGK_QualityFacSamp(2,iElem) END IF IF(BGK_QualityFacSamp(4,iElem).GT.0) THEN ! Max relaxation factor (maximal value of all octree subcells) - DSMC_MacroVal(nVarCount+4,iElem) = BGK_QualityFacSamp(3,iElem) / BGK_QualityFacSamp(4,iElem) + DSMC_MacroVal(nVarCount+6,iElem) = BGK_QualityFacSamp(3,iElem) / BGK_QualityFacSamp(4,iElem) ! Max rotational relaxation factor - DSMC_MacroVal(nVarCount+5,iElem) = BGK_QualityFacSamp(5,iElem) / BGK_QualityFacSamp(4,iElem) + DSMC_MacroVal(nVarCount+7,iElem) = BGK_QualityFacSamp(5,iElem) / BGK_QualityFacSamp(4,iElem) END IF ! Ratio between BGK and DSMC usage per cell - DSMC_MacroVal(nVarCount+6,iElem) = BGK_QualityFacSamp(4,iElem) / iter_loc - nVarCount = nVarCount + 6 + DSMC_MacroVal(nVarCount+8,iElem) = BGK_QualityFacSamp(4,iElem) / iter_loc + nVarCount = nVarCount + 8 END IF ! variable rotation and vibration relaxation IF(Collismode.GT.1) THEN @@ -884,7 +888,7 @@ SUBROUTINE WriteDSMCToHDF5(MeshFileName,OutputTime) IF(UseVarTimeStep) nVar_quality = nVar_quality + 1 IF(DoVirtualCellMerge) nVar_quality = nVar_quality + 1 IF(RadialWeighting%PerformCloning) nVar_quality = nVar_quality + 2 - IF(BGKInitDone) nVar_quality = nVar_quality + 6 + IF(BGKInitDone) nVar_quality = nVar_quality + 8 IF(FPInitDone) nVar_quality = nVar_quality + 5 ELSE nVar_quality=0 @@ -970,10 +974,12 @@ SUBROUTINE WriteDSMCToHDF5(MeshFileName,OutputTime) StrVarNames(nVarCount+1) ='BGK_MeanRelaxationFactor' StrVarNames(nVarCount+2) ='BGK_MeanPrandtlNumber' StrVarNames(nVarCount+3) ='BGK_ExpectedPrandtlNumber' - StrVarNames(nVarCount+4) ='BGK_MaxRelaxationFactor' - StrVarNames(nVarCount+5) ='BGK_MaxRotationRelaxFactor' - StrVarNames(nVarCount+6) ='BGK_DSMC_Ratio' - nVarCount=nVarCount+6 + StrVarNames(nVarCount+4) ='BGK_Viscosity' + StrVarNames(nVarCount+5) ='BGK_ThermalConductivity' + StrVarNames(nVarCount+6) ='BGK_MaxRelaxationFactor' + StrVarNames(nVarCount+7) ='BGK_MaxRotationRelaxFactor' + StrVarNames(nVarCount+8) ='BGK_DSMC_Ratio' + nVarCount=nVarCount+8 END IF IF(FPInitDone) THEN StrVarNames(nVarCount+1) ='FP_MeanRelaxationFactor' diff --git a/src/particles/fp_flow/fpflow_colloperator.f90 b/src/particles/fp_flow/fpflow_colloperator.f90 index 31557eea9..3918580c3 100644 --- a/src/particles/fp_flow/fpflow_colloperator.f90 +++ b/src/particles/fp_flow/fpflow_colloperator.f90 @@ -51,7 +51,6 @@ SUBROUTINE FP_CollisionOperator(iPartIndx_Node, nPart, NodeVolume) USE MOD_DSMC_Vars ,ONLY: CollInf, RadialWeighting USE Ziggurat USE MOD_Particle_Analyze_Tools ,ONLY: CalcTVibPoly -USE MOD_BGK_CollOperator ,ONLY: CalcTEquiPoly, CalcTEqui USE MOD_part_tools ,ONLY: GetParticleWeight ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -186,7 +185,7 @@ SUBROUTINE FP_CollisionOperator(iPartIndx_Node, nPart, NodeVolume) nXiVibDOF = PolyatomMolDSMC(iPolyatMole)%VibDOF ALLOCATE(Xi_vib_DOF(nXiVibDOF)) Xi_vib_DOF(:) = 0. - TVib = CalcTVibPoly(Evib/totalWeight, 1) + TVib = CalcTVibPoly(Evib/totalWeight + SpecDSMC(1)%EZeroPoint, 1) IF (TVib.GT.0.0) THEN DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF Xi_vib = Xi_vib + 2.*PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)/TVib & @@ -703,4 +702,275 @@ SUBROUTINE FP_CollisionOperator(iPartIndx_Node, nPart, NodeVolume) END SUBROUTINE FP_CollisionOperator +SUBROUTINE CalcTEqui(nPart, CellTemp, TRot, TVib, Xi_Vib, Xi_Vib_old, RotExp, VibExp, & + TEqui, rotrelaxfreq, vibrelaxfreq, dtCell, DoVibRelaxIn) +!=================================================================================================================================== +! Calculation of the vibrational temperature (zero-point search) for non-polyatomic molecules for Fokker-Planck +!=================================================================================================================================== +! MODULES +USE MOD_DSMC_Vars, ONLY: SpecDSMC +USE MOD_BGK_Vars, ONLY: BGKDoVibRelaxation +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +REAL, INTENT(IN) :: CellTemp, TRot, TVib, Xi_Vib_old, rotrelaxfreq, vibrelaxfreq, dtCell +INTEGER, INTENT(IN) :: nPart +LOGICAL, OPTIONAL, INTENT(IN) :: DoVibRelaxIn +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +REAL, INTENT(OUT) :: Xi_vib, TEqui, RotExp, VibExp +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +REAL :: TEqui_Old, betaR, betaV, RotFrac, VibFrac, TEqui_Old2 +REAL :: eps_prec=1.0E-0 +REAL :: correctFac, correctFacRot, maxexp !, Xi_rel +LOGICAL :: DoVibRelax +!=================================================================================================================================== +IF (PRESENT(DoVibRelaxIn)) THEN + DoVibRelax = DoVibRelaxIn +ELSE + DoVibRelax = BGKDoVibRelaxation +END IF +maxexp = LOG(HUGE(maxexp)) +! Xi_rel = 2.*(2. - CollInf%omega(1,1)) +! correctFac = 1. + (2.*SpecDSMC(1)%CharaTVib / (CellTemp*(EXP(SpecDSMC(1)%CharaTVib / CellTemp)-1.)))**(2.) & +! * EXP(SpecDSMC(1)%CharaTVib /CellTemp) / (2.*Xi_rel) +! correctFacRot = 1. + 2./Xi_rel + +correctFac = 1. +correctFacRot = 1. +RotExp = exp(-rotrelaxfreq*dtCell/correctFacRot) +RotFrac = nPart*(1.-RotExp) +IF(DoVibRelax) THEN + VibExp = exp(-vibrelaxfreq*dtCell/correctFac) + VibFrac = nPart*(1.-VibExp) +ELSE + VibExp = 0.0 + VibFrac = 0.0 + Xi_vib = 0.0 +END IF +TEqui_Old = 0.0 +TEqui = (3.*(nPart-1.)*CellTemp+2.*RotFrac*TRot+Xi_Vib_old*VibFrac*TVib)/(3.*(nPart-1.)+2.*RotFrac+Xi_Vib_old*VibFrac) +DO WHILE ( ABS( TEqui - TEqui_Old ) .GT. eps_prec ) + IF (ABS(TRot-TEqui).LT.1E-3) THEN + RotExp = exp(-rotrelaxfreq*dtCell/correctFacRot) + ELSE + betaR = ((TRot-CellTemp)/(TRot-TEqui))*rotrelaxfreq*dtCell/correctFacRot + IF (-betaR.GT.0.0) THEN + RotExp = 0. + ELSE IF (betaR.GT.maxexp) THEN + RotExp = 0. + ELSE + RotExp = exp(-betaR) + END IF + END IF + RotFrac = nPart*(1.-RotExp) + IF(DoVibRelax) THEN + IF (ABS(TVib-TEqui).LT.1E-3) THEN + VibExp = exp(-vibrelaxfreq*dtCell/correctFac) + ELSE + betaV = ((TVib-CellTemp)/(TVib-TEqui))*vibrelaxfreq*dtCell/correctFac + IF (-betaV.GT.0.0) THEN + VibExp = 0. + ELSE IF (betaV.GT.maxexp) THEN + VibExp = 0. + ELSE + VibExp = exp(-betaV) + END IF + END IF + IF ((SpecDSMC(1)%CharaTVib/TEqui).GT.maxexp) THEN + Xi_Vib = 0.0 + ELSE + Xi_vib = 2.*SpecDSMC(1)%CharaTVib/TEqui/(EXP(SpecDSMC(1)%CharaTVib/TEqui)-1.) + END IF + VibFrac = nPart*(1.-VibExp) + END IF + TEqui_Old = TEqui + TEqui_Old2 = TEqui + TEqui = (3.*(nPart-1.)*CellTemp+2.*RotFrac*TRot+Xi_Vib_old*VibFrac*TVib)/(3.*(nPart-1.)+2.*RotFrac+Xi_Vib*VibFrac) + IF(DoVibRelax) THEN + DO WHILE( ABS( TEqui - TEqui_Old2 ) .GT. eps_prec ) + TEqui =(TEqui + TEqui_Old2)*0.5 + IF ((SpecDSMC(1)%CharaTVib/TEqui).GT.maxexp) THEN + Xi_Vib = 0.0 + ELSE + Xi_vib = 2.*SpecDSMC(1)%CharaTVib/TEqui/(EXP(SpecDSMC(1)%CharaTVib/TEqui)-1.) + END IF + TEqui_Old2 = TEqui + TEqui = (3.*(nPart-1.)*CellTemp+2.*RotFrac*TRot+Xi_Vib_old*VibFrac*TVib) / (3.*(nPart-1.)+2.*RotFrac+Xi_vib*VibFrac) + END DO + END IF +END DO + +END SUBROUTINE CalcTEqui + + +SUBROUTINE CalcTEquiPoly(nPart, CellTemp, TRot, TVib, nXiVibDOF, Xi_Vib_DOF, Xi_Vib_old, RotExp, VibExp, TEqui, rotrelaxfreq, vibrelaxfreq, & + dtCell, DoVibRelaxIn) +!=================================================================================================================================== +! Calculation of the vibrational temperature (zero-point search) for polyatomic molecules +!=================================================================================================================================== +! MODULES +USE MOD_DSMC_Vars, ONLY: SpecDSMC, PolyatomMolDSMC +USE MOD_BGK_Vars, ONLY: BGKDoVibRelaxation +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +REAL, INTENT(IN) :: CellTemp, TRot, TVib, Xi_Vib_old, rotrelaxfreq, vibrelaxfreq +INTEGER, INTENT(IN) :: nPart,nXiVibDOF +REAL, INTENT(IN) :: dtCell +LOGICAL, OPTIONAL, INTENT(IN) :: DoVibRelaxIn +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +REAL, INTENT(OUT) :: Xi_vib_DOF(nXiVibDOF), TEqui, RotExp, VibExp +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +REAL :: TEqui_Old, betaR, betaV, RotFrac, VibFrac, Xi_Rot, TEqui_Old2, exparg +REAL :: eps_prec=1.0 +REAL :: correctFac, correctFacRot +INTEGER :: iDOF, iPolyatMole +LOGICAL :: DoVibRelax +!=================================================================================================================================== +IF (PRESENT(DoVibRelaxIn)) THEN + DoVibRelax = DoVibRelaxIn +ELSE + DoVibRelax = BGKDoVibRelaxation +END IF + +! rotational degrees of freedom of polyatomic molecule +Xi_Rot = SpecDSMC(1)%Xi_Rot +iPolyatMole = SpecDSMC(1)%SpecToPolyArray + +! Xi_rel = 2.*(2. - CollInf%omega(1,1)) +! correctFac = 0.0 +! DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF +! correctFac = correctFac & +! + (2.*PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF) / (CellTemp & +! *(EXP(PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF) / CellTemp)-1.)))**(2.) & +! * EXP(PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF) / CellTemp) / 2. +! END DO +! correctFac = 1. + correctFac/Xi_rel +! correctFacRot = 1. + Xi_Rot/Xi_rel + +correctFac = 1. +correctFacRot = 1. + +! Calculate number of rotational relaxing molecules with number of molecules * probability of relaxation +! P = 1 - exp(-nu*dt) with relaxation frequency nu and timestep dt +RotExp = exp(-rotrelaxfreq*dtCell/correctFacRot) +RotFrac = nPart*(1.-RotExp) +! Calculate number of vibrational relaxing molecules if enabled with number of molecules * probability of relaxation +! P = 1 - exp(-nu*dt) with relaxation frequency nu and timestep dt +IF(DoVibRelax) THEN + VibExp = exp(-vibrelaxfreq*dtCell/correctFac) + VibFrac = nPart*(1.-VibExp) +ELSE + VibExp = 0.0 + VibFrac = 0.0 + Xi_vib_DOF = 0.0 +END IF +TEqui_Old = 0.0 +! M. Pfeiffer et. al., "Extension of Particle-based BGK Models to Polyatomic Species in Hypersonic Flow around a Flat-faced +! Cylinder", AIP Conference Proceedings 2132, 100001 (2019) +! Solving of equation system for TEqui and betaR and betaV +TEqui = (3.*(nPart-1.)*CellTemp+2.*RotFrac*TRot+Xi_Vib_old*VibFrac*TVib)/(3.*(nPart-1.)+2.*RotFrac+Xi_Vib_old*VibFrac) +! Required condition of Landau-Teller relaxation not fulfilled --> relaxation probabilities of rotation and vibration are +! corrected with a parameter beta for rotation and vibration as suggested by Burt: +! J. Burt and I. Boyd, “Evaluation of a particle method for the ellipsoidal statistical Bhatnagar-Gross-Krook equation”, +! 44th AIAA Aerospace Sciences Meeting and Exhibit (AIAA, 2006), p. 989 +! Solving of equation system until accuracy eps_prec is reached +DO WHILE ( ABS( TEqui - TEqui_Old ) .GT. eps_prec ) + ! if difference too small: beta is not taken into account + IF (ABS(TRot-TEqui).LT.1E-3) THEN + RotExp = exp(-rotrelaxfreq*dtCell/correctFacRot) + ELSE + ! betaR = beta*nu*dt (= correction parameter rotation * relaxation frequency * time step) + betaR = ((TRot-CellTemp)/(TRot-TEqui))*rotrelaxfreq*dtCell/correctFacRot + ! negative betaR would leed to negative relaxation probability! + IF (-betaR.GT.0.0) THEN + RotExp = 0. + ! Check if the exponent is within the range of machine precision + ELSE IF (CHECKEXP(betaR)) THEN + RotExp = exp(-betaR) + ELSE + RotExp = 0. + END IF + END IF + ! new calculation of number of rotational relaxing molecules + RotFrac = nPart*(1.-RotExp) + + IF(DoVibRelax) THEN + ! if difference too small: beta is not taken into account + IF (ABS(TVib-TEqui).LT.1E-3) THEN + VibExp = exp(-vibrelaxfreq*dtCell/correctFac) + ELSE + ! betaV = beta*nu*dt (= correction parameter vibration * relaxation frequency * time step) + betaV = ((TVib-CellTemp)/(TVib-TEqui))*vibrelaxfreq*dtCell/correctFac + ! negative betaV would leed to negative relaxation probability! + IF (-betaV.GT.0.0) THEN + VibExp = 0. + ! Check if the exponent is within the range of machine precision + ELSEIF(CHECKEXP(betaV))THEN + VibExp = exp(-betaV) + ELSE + VibExp = 0. + END IF + END IF + ! new calculation of number of vibrational relaxing molecules + VibFrac = nPart*(1.-VibExp) + + ! Loop over all vibrational degrees of freedom to calculate them using TEqui + DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF + exparg = PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)/TEqui + ! Check if the exponent is within the range of machine precision for calculation of vibrational degrees of freedom + IF(CHECKEXP(exparg))THEN + IF(exparg.gt.0.)THEN ! positive overflow: exp -> inf + Xi_vib_DOF(iDOF) = 2.*exparg/(EXP(exparg)-1.) + ELSE ! negative overflow: exp -> 0 + Xi_vib_DOF(iDOF) = 2.*exparg/(-1.) + END IF ! exparg.gt.0. + ELSE + Xi_vib_DOF(iDOF) = 0.0 + END IF ! CHECKEXP(exparg) + END DO + END IF + TEqui_Old = TEqui + TEqui_Old2 = TEqui + + ! new calculation of equilibrium temperature with new RotFrac, new VibFrac new Xi_vib_DOF(TEqui) in denominator + TEqui = (3.*(nPart-1.)*CellTemp+2.*RotFrac*TRot+Xi_Vib_old*VibFrac*TVib) & + / (3.*(nPart-1.)+2.*RotFrac+SUM(Xi_vib_DOF(1:PolyatomMolDSMC(iPolyatMole)%VibDOF))*VibFrac) + IF(DoVibRelax) THEN + ! accuracy eps_prec not reached yet + DO WHILE( ABS( TEqui - TEqui_Old2 ) .GT. eps_prec ) + ! mean value of old and new equilibrium temperature + TEqui =(TEqui + TEqui_Old2)*0.5 + ! Loop over all vibrational degrees of freedom + DO iDOF = 1, PolyatomMolDSMC(iPolyatMole)%VibDOF + exparg = PolyatomMolDSMC(iPolyatMole)%CharaTVibDOF(iDOF)/TEqui + ! Check if the exponent is within the range of machine precision for calculation of vibrational degrees of freedom + IF(CHECKEXP(exparg))THEN + IF(exparg.gt.0.)THEN ! positive overflow: exp -> inf + Xi_vib_DOF(iDOF) = 2.*exparg/(EXP(exparg)-1.) + ELSE ! negative overflow: exp -> 0 + Xi_vib_DOF(iDOF) = 2.*exparg/(-1.) + END IF ! exparg.gt.0. + ELSE + Xi_vib_DOF(iDOF) = 0.0 + END IF ! CHECKEXP(exparg) + END DO + TEqui_Old2 = TEqui + ! new calculation of equilibrium temperature with corrected vibrational degrees of freedom in denominator + TEqui = (3.*(nPart-1.)*CellTemp+2.*RotFrac*TRot+Xi_Vib_old*VibFrac*TVib) & + / (3.*(nPart-1.)+2.*RotFrac+SUM(Xi_vib_DOF(1:PolyatomMolDSMC(iPolyatMole)%VibDOF))*VibFrac) + END DO + END IF +END DO + +END SUBROUTINE CalcTEquiPoly + END MODULE MOD_FP_CollOperator