diff --git a/.github/workflows/auto-update-gh-pages.yml b/.github/workflows/auto-update-gh-pages.yml index c1a75269fa..333b7983bc 100644 --- a/.github/workflows/auto-update-gh-pages.yml +++ b/.github/workflows/auto-update-gh-pages.yml @@ -39,7 +39,7 @@ jobs: conda deactivate conda activate moose-env git submodule init && git submodule update - make NAVIER_STOKES:='no' PHASE_FIELD:='no' REACTOR:='no' -j 2 + make -j 2 cd doc ./moosedocs.py build --destination html diff --git a/doc/config.yml b/doc/config.yml index 51ff535e5b..1780ef96ab 100644 --- a/doc/config.yml +++ b/doc/config.yml @@ -3,6 +3,11 @@ Content: root_dir: ${ROOT_DIR}/doc/content moose: root_dir: ${MOOSE_DIR}/framework/doc/content + content: + - contrib/** + - css/** + - js/** + - media/** squirrel: root_dir: ${ROOT_DIR}/squirrel/doc/content @@ -20,6 +25,7 @@ Extensions: menu: Getting Started: Install Moltres: getting_started/installation.md + Theory: getting_started/theory.md Tutorials: getting_started/tutorials.md Documentation: Moltres Syntax: syntax/index.md diff --git a/doc/content/getting_started/theory.md b/doc/content/getting_started/theory.md new file mode 100644 index 0000000000..154adff0e4 --- /dev/null +++ b/doc/content/getting_started/theory.md @@ -0,0 +1,174 @@ +# Theory and Methodology + +The *two-step procedure* is a common approach for multiphysics full-core nuclear reactor analysis. +These steps are: + ++Step 1:+ Generate neutron group constant data with a lattice or full-core reactor model on a +high-fidelity, continuous/fine-energy neutronics software at various reactor states. + ++Step 2:+ Use an intermediate-fidelity, computationally cheaper neutronics software with the +neutron group constant data to perform multiphysics reactor analysis. + +Moltres falls under Step 2 of the two-step procedure. Moltres solves the conventional multigroup +neutron diffusion and delayed neutron precursor equations. The precursor equations include +precursor drift terms to account for precursor movement under molten salt flow. Moltres is built on +the MOOSE finite element framework, enabling highly flexible and scalable reactor simulations. + +## Multigroup Neutron Diffusion + +More information to be added in future. + +## Heat Transfer and Fluid Flow + +Moltres compiles with the MOOSE +[Navier-Stokes](https://mooseframework.inl.gov/modules/navier_stokes/index.html) and +[Heat Transfer](https://mooseframework.inl.gov/modules/heat_transfer/index.html) modules for flow +modeling capabilities by default. Moltres couples with these modules natively because they are all +built on the MOOSE framework. + +Past work with Moltres have modeled incompressible salt flow and temperature advection-diffusion +with the +[INS or INSAD implementations](https://mooseframework.inl.gov/modules/navier_stokes/cgfe.html) of +the Navier-Stokes equations [!citep](peterson_overview_2018). Coupling with the compressible or +finite volume implementations within the Navier-Stokes module are likely possible, but have not +been demonstrated yet. + +For turbulence modeling, Moltres has a Spalart-Allmaras (SA) turbulence model +[!citep](spalart_one-equation_1992) implementation with streamline-upwind Petrov-Galerkin (SUPG) +stabilization. On balance the SA model is a complete (does not require prior knowledge of the +actual turbulence behavior) and computationally efficient turbulence model for approximating +wall-bounded turbulent flows. The SA model implementation in Moltres is +compatible with the INSAD implementation only. Additionally, Moltres contains turbulent diffusion +physics kernels for temperature and the delayed neutron precursors. + +The SA model in Moltres follows the (SA-noft2-R) implementation as described on the NASA +[Turbulence Modeling Resource](https://turbmodels.larc.nasa.gov/spalart.html) website. The $f_{t2}$ +term is togglable using the `use_ft2_term` input parameter (false by default). Moltres' +implementation solves for the modified dynamic viscosity $\tilde{\mu}$ which can be easily obtained +using $\nu=mu/\rho$ to convert $\nu$ to $\mu$ as follows: + +!equation id=sa-equation +\rho \frac{\partial\tilde{\mu}}{\partial t} + \rho \mathbf{u}\cdot\nabla\tilde{\mu} = \rho c_{b1} +\left(1-f_{t2}\right)\tilde{S}\tilde{\mu} + \frac{1}{\sigma}\{\nabla\cdot\left[\left(\mu+ +\tilde{\mu}\right)\nabla\tilde{\mu}\right] + c_{b2}|\nabla\tilde{\mu}|^2\} - \left(c_{w1}f_w - +\frac{c_{b1}}{\kappa^2}f_{t2}\right)\left( +\frac{\tilde{\mu}}{d}\right)^2 + +where + +!equation +\mu_t = \tilde{\mu}f_{v1} = \text{turbulent eddy viscosity}, \hspace{4cm} +f_{v1} = \frac{\chi^3}{\chi^3 + c_{v1}^3}, \hspace{4cm} +\chi = \frac{\tilde{\mu}}{\mu}, \hspace{2cm} \\ \hspace{1cm} +\tilde{S} = \Omega + \frac{\tilde{\nu}}{\kappa^2 d^2} f_{v2}, \hspace{6cm} +f_{v2} = 1 - \frac{\chi}{1+\chi f_{v1}}, \hspace{2cm} +\Omega = \sqrt{2W_{ij}W_{ij}} = \text{vorticity magnitude}, \\ +W_{ij} = \frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} - \frac{\partial u_j}{\partial x_i} +\right), \hspace{5cm} +f_w = g\left(\frac{1 + c_{w3}^6}{g^6 + c_{w3}^6}\right)^{1/6}, \hspace{3cm} +g = r + c_{w2}\left(r^6 - r\right), \\ +r = \text{min}\left(\frac{\tilde{\nu}}{\tilde{S}\kappa^2d^2}, 10\right), \hspace{5cm} +f_{t2} = c_{t3} \exp{\left(-c_{t4}\chi^2\right)}, \hspace{6cm} + +and the constants are: + +!equation +\sigma = \frac{2}{3}, \ c_{b1} = 0.1355, \ c_{b2} = 0.622, \ \kappa = 0.41, \ +c_{w1} = \frac{c_{b1}}{\kappa^2} + \frac{1+c_{b2}}{\sigma}, \ c_{w2} = 0.3, \ c_{w3} = 2, \ +c_{v1} = 7.1, \ c_{t3} = 1.2, \ c_{t4} = 0.5 \ \text{.} + +We performed simple verification and validation tests of the SA model in Moltres using reference +problems for channel, pipe, and backward-facing step flow. The Moltres test dataset is available on +[Zenodo](https://doi.org/10.5281/zenodo.10059649). The following subsections show plots comparing +results from Moltres to reference data from literature. + +#### Turbulent Channel Flow Verification Test + +The channel flow test is based on the direct numerical simulation (DNS) of turbulent channel flow +with $Re_\tau\approx395$ +by [!cite](moser_direct_1999). The Moltres input file for this problem may be found in +`moltres/problems/2023-basic-turbulence-cases/channel-flow/`. + +!media channel_vel.png + id=channel-vel + caption=Normalized velocity distribution across the channel. + style=width:33%;padding:10px;float:left; + +!media channel_nondim.png + id=channel-nondim + caption=Dimensionless velocity vs dimensionless wall distance. + style=width:33%;padding:10px;float:left; + +!media channel_stress.png + id=channel-stress + caption=Normalized stress distribution across the channel. + style=width:33%;padding:10px;float:left; + +#### Turbulent Pipe Flow Validation Test + +The pipe flow test is based on the turbulent pipe flow experiment with $Re\approx 40,000$ by +[!cite](laufer_structure_1954). The Moltres input file for this problem may be found in +`moltres/problems/2023-basic-turbulence-cases/pipe-flow/`. + +!media pipe_vel.png + id=pipe-vel + caption=Normalized velocity distribution across the pipe. + style=width:33%;padding:10px;float:left; + +!media pipe_nondim.png + id=pipe-nondim + caption=Dimensionless velocity vs dimensionless wall distance. + style=width:33%;padding:10px;float:left; + +!media pipe_stress.png + id=pipe-stress + caption=Normalized turbulent shear stress distribution across the channel. + style=width:33%;padding:10px;float:left; + +#### Backward-Facing Step Flow Validation Test + +The backward-facing step (BFS) flow test is based on the BFS flow experiment with $Re\approx +36,000$ by +[!cite](driver_features_1985). The Moltres input file for this problem may be found in +`moltres/problems/2023-basic-turbulence-cases/pipe-flow/`. +The SA model implementation in Moltres performs largely similarly to the reference SA model results +provided on the [Turbulence Modeling Resource](https://turbmodels.larc.nasa.gov/spalart.html) +website. + +!table caption=Flow reattachment length estimates (normalized by step height $H$) id=re-table +| Reference experimental data | Reference SA model | Moltres | +| :- | :- | :- | +| $6.26 \pm 0.10$ | $6.1$ | $6.36$ | + +!media bfs_upstream_vel.png + id=bfs-upstream-vel + caption=Normalized velocity distribution at $x/H$ equal -4 (upstream of step). + style=width:33%;padding:10px;float:left; + +!media bfs_downstream_vel.png + id=bfs-downstream-vel + caption=Normalized velocity distribution at various $x/H$ locations (downstream of step). + style=width:33%;padding:10px;float:left; + +!media bfs_cf.png + id=bfs-cf + caption=Skin friction coefficient along the bottom wall. + style=width:33%;padding:10px;float:left; + +!media bfs_cp.png + id=bfs-cp + caption=Skin pressure coefficient along the bottom wall. + style=width:33%;padding:10px;margin-left:17%;float:left; + +!media bfs_stress.png + id=bfs-stress + caption=Normalized turbulent shear stress distribution at various $x/H$ locations + (downstream of step). + style=width:33%;padding:10px;margin-right:17%;float:left; + +!media bfs.png + id=bfs + caption=Close-up of the velocity magnitude and streamlines around the backward-facing step + style=width:80%;float:left;margin-left:10%;margin-right:10%; + +!bibtex bibliography diff --git a/doc/content/getting_started/turbulence.bib b/doc/content/getting_started/turbulence.bib new file mode 100644 index 0000000000..b7b1ec97bc --- /dev/null +++ b/doc/content/getting_started/turbulence.bib @@ -0,0 +1,77 @@ + +@techreport{laufer_structure_1954, + title = {The structure of turbulence in fully developed pipe flow}, + url = {https://ntrs.nasa.gov/citations/19930092199}, + abstract = {Measurements, principally with a hot-wire anemometer, were made in fully developed turbulent flow in a 10-inch pipe at speeds of approximately 10 and 100 feet per second. Emphasis was placed on turbulence and conditions near the wall. The results include relevant mean and statistical quantities, such as Reynolds stresses, triple correlations, turbulent dissipation, and energy spectra. It is shown that rates of turbulent-energy production, dissipation, and diffusion have sharp maximums near the edge of the laminar sublayer and that there exist a strong movement of kinetic energy away from this point and an equally strong movement of pressure energy toward it.}, + urldate = {2023-10-30}, + institution = {National Bureau of Standards}, + author = {Laufer, John}, + month = jan, + year = {1954}, + note = {NTRS Author Affiliations: +NTRS Report/Patent Number: NACA-TR-1174 +NTRS Document ID: 19930092199 +NTRS Research Center: Legacy CDMS (CDMS)}, + file = {19930092199.pdf:C\:\\Users\\Sun Myung\\Zotero\\storage\\L8SHCFP9\\Laufer - 1954 - The structure of turbulence in fully developed pip.pdf:application/pdf;Snapshot:C\:\\Users\\Sun Myung\\Zotero\\storage\\7U6SXNMA\\19930092199.html:text/html}, +} + +@article{moser_direct_1999, + title = {Direct numerical simulation of turbulent channel flow up to {Reτ}=590}, + volume = {11}, + issn = {1070-6631}, + url = {https://doi.org/10.1063/1.869966}, + doi = {10.1063/1.869966}, + abstract = {Numerical simulations of fully developed turbulent channel flow at three Reynolds numbers up to Reτ=590 are reported. It is noted that the higher Reynolds number simulations exhibit fewer low Reynolds number effects than previous simulations at Reτ=180. A comprehensive set of statistics gathered from the simulations is available on the web at http://www.tam.uiuc.edu/Faculty/Moser/channel.}, + number = {4}, + urldate = {2023-10-30}, + journal = {Physics of Fluids}, + author = {Moser, Robert D. and Kim, John and Mansour, Nagi N.}, + month = apr, + year = {1999}, + pages = {943--945}, + file = {Full Text PDF:C\:\\Users\\Sun Myung\\Zotero\\storage\\R6IZTYHZ\\Moser et al. - 1999 - Direct numerical simulation of turbulent channel f.pdf:application/pdf;Snapshot:C\:\\Users\\Sun Myung\\Zotero\\storage\\NN9Z4VLF\\Direct-numerical-simulation-of-turbulent-channel.html:text/html}, +} + +@article{driver_features_1985, + title = {Features of a reattaching turbulent shear layer in divergent channelflow}, + volume = {23}, + issn = {0001-1452, 1533-385X}, + url = {https://arc.aiaa.org/doi/10.2514/3.8890}, + doi = {10.2514/3.8890}, + language = {en}, + number = {2}, + urldate = {2022-01-17}, + journal = {AIAA Journal}, + author = {Driver, David M. and Seegmiller, H. Lee}, + month = feb, + year = {1985}, + pages = {163--171}, + file = {Driver and Seegmiller - 1985 - Features of a reattaching turbulent shear layer in.pdf:C\:\\Users\\Sun Myung\\Zotero\\storage\\PHD7NPMR\\Driver and Seegmiller - 1985 - Features of a reattaching turbulent shear layer in.pdf:application/pdf}, +} + +@inproceedings{spalart_one-equation_1992, + title = {A one-equation turbulence model for aerodynamic flows}, + url = {https://arc.aiaa.org/doi/abs/10.2514/6.1992-439}, + urldate = {2023-11-02}, + booktitle = {30th {Aerospace} {Sciences} {Meeting} and {Exhibit}}, + publisher = {American Institute of Aeronautics and Astronautics}, + author = {Spalart, P. and Allmaras, S.}, + month = jan, + year = {1992}, + doi = {10.2514/6.1992-439}, + note = {\_eprint: https://arc.aiaa.org/doi/pdf/10.2514/6.1992-439}, +} + +@article{peterson_overview_2018, + title = {Overview of the {Incompressible} {Navier}-{Stokes} simulation capabilities in the {MOOSE} {Framework}}, + volume = {119}, + doi = {10.1016/j.advengsoft.2018.02.004}, + abstract = {The Multiphysics Object Oriented Simulation Environment (MOOSE) framework is a high-performance, open source, C++ finite element toolkit developed at Idaho National Laboratory. MOOSE was created with the aim of assisting domain scientists and engineers in creating customizable, high-quality tools for multiphysics simulations. While the core MOOSE framework itself does not contain code for simulating any particular physical application, it is distributed with a number of physics "modules" which are tailored to solving e.g. heat conduction, phase field, and solid/fluid mechanics problems. In this report, we describe the basic equations, finite element formulations, software implementation, and regression/verification tests currently available in MOOSE's navier\_stokes module for solving the Incompressible Navier-Stokes (INS) equations.}, + journal = {Advances in Engineering Software}, + author = {Peterson, John and Lindsay, Alexander and Kong, Fande}, + month = may, + year = {2018}, + keywords = {65N30, Mathematics - Numerical Analysis}, + pages = {68--92}, + file = {arXiv\:1710.08898 PDF:C\:\\Users\\Sun Myung\\Zotero\\storage\\PBHHJKCD\\Peterson et al. - 2017 - Overview of the Incompressible Navier-Stokes simul.pdf:application/pdf;arXiv.org Snapshot:C\:\\Users\\Sun Myung\\Zotero\\storage\\3RRL7ED2\\1710.html:text/html;Full Text PDF:C\:\\Users\\Sun Myung\\Zotero\\storage\\YGQ9TDMU\\Peterson et al. - 2018 - Overview of the Incompressible Navier-Stokes simul.pdf:application/pdf}, +} diff --git a/doc/content/getting_started/tutorials.md b/doc/content/getting_started/tutorials.md index ad5958a547..a912dbb0f6 100644 --- a/doc/content/getting_started/tutorials.md +++ b/doc/content/getting_started/tutorials.md @@ -1,15 +1,6 @@ # Tutorials -The *two-step procedure* is a common approach for multiphysics full-core nuclear reactor analysis. -These steps are: - -+Step 1:+ Generate neutron group constant data with a lattice or full-core reactor model on a -high-fidelity, continuous/fine-energy neutronics software at various reactor states. - -+Step 2:+ Use an intermediate-fidelity, computationally cheaper neutronics software with the -neutron group constant data to perform multiphysics reactor analysis. - -Moltres falls under Step 2 of the two-step procedure. The following sections cover +The following sections cover various tutorials for group constant file generation and running various types of reactor simulations with Moltres. The tutorials assume that readers have a basic understanding of reactor analysis, molten salt reactors, and the [MOOSE framework](https://mooseframework.inl.gov/). @@ -44,7 +35,7 @@ the existing tutorials below. +2b.+ [Eigenvalue Calculation](getting_started/eigenvalue.md) -+2c.+ [Time-dependent Simulation with Thermal-Hydraulic Coupling and Precursor Looping](getting_started/transient.md) ++2c.+ [Time-Dependent Simulation with Thermal-Hydraulic Coupling and Precursor Looping](getting_started/transient.md) +Tip:+ [Recommended Executioner and Preconditioning Settings](getting_started/recommended.md) diff --git a/doc/content/index.md b/doc/content/index.md index b01a6ebf51..b5d649b752 100644 --- a/doc/content/index.md +++ b/doc/content/index.md @@ -12,10 +12,11 @@ Moltres is a finite element-based multiphysics reactor simulation tool for molte (MSRs) and other advanced reactors in various dimensions and geometries. Moltres is built on the MOOSE framework and provides multiphysics capabilities to accurately model delayed neutron precursor drift and the strong temperature -feedback effects inherent in MSRs. Users have the flexibility of setting up fully coupled solves -in which all variables are solved in a single monolithic system or tightly coupled solves by -arbitrarily segregating systems of equations corresponding to separate physics. As a MOOSE-based -application, Moltres can also scale efficiently on high-performance computers. +feedback effects inherent in MSRs. Users have the flexibility of setting up fully coupled +solution schemes in which all variables are solved in a single monolithic system, or tightly +coupled solution schemes by arbitrarily segregating systems of equations corresponding to separate +physics. As a MOOSE-based application, Moltres can also scale efficiently on high-performance +computers. !row! !col! small=12 medium=4 large=6 diff --git a/doc/content/media/bfs.png b/doc/content/media/bfs.png new file mode 100644 index 0000000000..b84bcff364 Binary files /dev/null and b/doc/content/media/bfs.png differ diff --git a/doc/content/media/bfs_cf.png b/doc/content/media/bfs_cf.png new file mode 100644 index 0000000000..7448c37838 Binary files /dev/null and b/doc/content/media/bfs_cf.png differ diff --git a/doc/content/media/bfs_cp.png b/doc/content/media/bfs_cp.png new file mode 100644 index 0000000000..10513ea2b3 Binary files /dev/null and b/doc/content/media/bfs_cp.png differ diff --git a/doc/content/media/bfs_downstream_vel.png b/doc/content/media/bfs_downstream_vel.png new file mode 100644 index 0000000000..5adbfd0f24 Binary files /dev/null and b/doc/content/media/bfs_downstream_vel.png differ diff --git a/doc/content/media/bfs_stress.png b/doc/content/media/bfs_stress.png new file mode 100644 index 0000000000..1e18f4862d Binary files /dev/null and b/doc/content/media/bfs_stress.png differ diff --git a/doc/content/media/bfs_upstream_vel.png b/doc/content/media/bfs_upstream_vel.png new file mode 100644 index 0000000000..d39e2965bf Binary files /dev/null and b/doc/content/media/bfs_upstream_vel.png differ diff --git a/doc/content/media/channel_nondim.png b/doc/content/media/channel_nondim.png new file mode 100644 index 0000000000..7aab83d00e Binary files /dev/null and b/doc/content/media/channel_nondim.png differ diff --git a/doc/content/media/channel_stress.png b/doc/content/media/channel_stress.png new file mode 100644 index 0000000000..48d37c04e1 Binary files /dev/null and b/doc/content/media/channel_stress.png differ diff --git a/doc/content/media/channel_vel.png b/doc/content/media/channel_vel.png new file mode 100644 index 0000000000..00aa43d532 Binary files /dev/null and b/doc/content/media/channel_vel.png differ diff --git a/doc/content/media/pipe_nondim.png b/doc/content/media/pipe_nondim.png new file mode 100644 index 0000000000..177ddab57a Binary files /dev/null and b/doc/content/media/pipe_nondim.png differ diff --git a/doc/content/media/pipe_stress.png b/doc/content/media/pipe_stress.png new file mode 100644 index 0000000000..eb335757cc Binary files /dev/null and b/doc/content/media/pipe_stress.png differ diff --git a/doc/content/media/pipe_vel.png b/doc/content/media/pipe_vel.png new file mode 100644 index 0000000000..6e7529c703 Binary files /dev/null and b/doc/content/media/pipe_vel.png differ diff --git a/doc/content/source/auxkernels/SATurbulentViscosityAux.md b/doc/content/source/auxkernels/SATurbulentViscosityAux.md new file mode 100644 index 0000000000..e114e4fe04 --- /dev/null +++ b/doc/content/source/auxkernels/SATurbulentViscosityAux.md @@ -0,0 +1,18 @@ +# SATurbulentViscosityAux + +!syntax description /AuxKernels/SATurbulentViscosityAux + +## Overview + +This object computes the turbulent viscosity $\mu_t=\tilde{\mu}f_{v1}$ as described +[here](theory.md). + +## Example Input File Syntax + +!! Describe and include an example of how to use the SATurbulentViscosityAux object. + +!syntax parameters /AuxKernels/SATurbulentViscosityAux + +!syntax inputs /AuxKernels/SATurbulentViscosityAux + +!syntax children /AuxKernels/SATurbulentViscosityAux diff --git a/doc/content/source/auxkernels/TurbulentStressAux.md b/doc/content/source/auxkernels/TurbulentStressAux.md new file mode 100644 index 0000000000..f7b7047347 --- /dev/null +++ b/doc/content/source/auxkernels/TurbulentStressAux.md @@ -0,0 +1,18 @@ +# TurbulentStressAux + +!syntax description /AuxKernels/TurbulentStressAux + +## Overview + +This object computes the turbulent stress $\nu_t\sqrt{2S_{ij}S_{ij}}$ under the Spalart-Allmaras +model. + +## Example Input File Syntax + +!! Describe and include an example of how to use the TurbulentStressAux object. + +!syntax parameters /AuxKernels/TurbulentStressAux + +!syntax inputs /AuxKernels/TurbulentStressAux + +!syntax children /AuxKernels/TurbulentStressAux diff --git a/doc/content/source/auxkernels/WallDistanceAux.md b/doc/content/source/auxkernels/WallDistanceAux.md new file mode 100644 index 0000000000..a567f8ff30 --- /dev/null +++ b/doc/content/source/auxkernels/WallDistanceAux.md @@ -0,0 +1,20 @@ +# WallDistanceAux + +!syntax description /AuxKernels/WallDistanceAux + +## Overview + +This object computes the minimum distance of each mesh node to the nearest wall node. The wall +distance calculation may be slightly inaccurate for unstructured or skewed structured meshes as +illustrated on the [Turbulence Modeling Resource](https://turbmodels.larc.nasa.gov/spalart.html) +webpage. + +## Example Input File Syntax + +!! Describe and include an example of how to use the WallDistanceAux object. + +!syntax parameters /AuxKernels/WallDistanceAux + +!syntax inputs /AuxKernels/WallDistanceAux + +!syntax children /AuxKernels/WallDistanceAux diff --git a/doc/content/source/auxkernels/WallShearStressAux.md b/doc/content/source/auxkernels/WallShearStressAux.md new file mode 100644 index 0000000000..a7b9fd42f1 --- /dev/null +++ b/doc/content/source/auxkernels/WallShearStressAux.md @@ -0,0 +1,21 @@ +# WallShearStressAux + +!syntax description /AuxKernels/WallShearStressAux + +## Overview + +This object computes the wall shear stress +$\tau_w=\mu\left(\frac{\partial u_x}{\partial y}\right)_{y=0}$ +where $x$ is the direction of fluid flow and $y$ is the direction perpendicular to the wall. +Moltres automatically detects and sets the appropriate directions. This object is valid for +straight walls and non-separated flows only (i.e. no reversed flow along the wall). + +## Example Input File Syntax + +!! Describe and include an example of how to use the WallShearStressAux object. + +!syntax parameters /AuxKernels/WallShearStressAux + +!syntax inputs /AuxKernels/WallShearStressAux + +!syntax children /AuxKernels/WallShearStressAux diff --git a/doc/content/source/bcs/INSADMassPSPGBC.md b/doc/content/source/bcs/INSADMassPSPGBC.md new file mode 100644 index 0000000000..944d1b29d0 --- /dev/null +++ b/doc/content/source/bcs/INSADMassPSPGBC.md @@ -0,0 +1,18 @@ +# INSADMassPSPGBC + +!syntax description /BCs/INSADMassPSPGBC + +## Overview + +This object removes residual contributions from PSPG stabilization in the INSAD mass equation along +the boundaries. + +## Example Input File Syntax + +!! Describe and include an example of how to use the INSADMassPSPGBC object. + +!syntax parameters /BCs/INSADMassPSPGBC + +!syntax inputs /BCs/INSADMassPSPGBC + +!syntax children /BCs/INSADMassPSPGBC diff --git a/doc/content/source/bcs/INSADMomentumSUPGBC.md b/doc/content/source/bcs/INSADMomentumSUPGBC.md new file mode 100644 index 0000000000..521840c5f2 --- /dev/null +++ b/doc/content/source/bcs/INSADMomentumSUPGBC.md @@ -0,0 +1,18 @@ +# INSADMomentumSUPGBC + +!syntax description /BCs/INSADMomentumSUPGBC + +## Overview + +This object removes residual contributions from SUPG stabilization in the INSAD momentum equation +along the boundaries. + +## Example Input File Syntax + +!! Describe and include an example of how to use the INSADMomentumSUPGBC object. + +!syntax parameters /BCs/INSADMomentumSUPGBC + +!syntax inputs /BCs/INSADMomentumSUPGBC + +!syntax children /BCs/INSADMomentumSUPGBC diff --git a/doc/content/source/bcs/INSADMomentumTurbulentNoBCBC.md b/doc/content/source/bcs/INSADMomentumTurbulentNoBCBC.md new file mode 100644 index 0000000000..895d6c1d17 --- /dev/null +++ b/doc/content/source/bcs/INSADMomentumTurbulentNoBCBC.md @@ -0,0 +1,18 @@ +# INSADMomentumTurbulentNoBCBC + +!syntax description /BCs/INSADMomentumTurbulentNoBCBC + +## Overview + +This object applies the "no BC" boundary condition for the INSAD momentum equations with the +Spalart-Allmaras turbulence modeling. + +## Example Input File Syntax + +!! Describe and include an example of how to use the INSADMomentumTurbulentNoBCBC object. + +!syntax parameters /BCs/INSADMomentumTurbulentNoBCBC + +!syntax inputs /BCs/INSADMomentumTurbulentNoBCBC + +!syntax children /BCs/INSADMomentumTurbulentNoBCBC diff --git a/doc/content/source/bcs/INSOutflowBC.md b/doc/content/source/bcs/INSOutflowBC.md deleted file mode 100644 index 9827613861..0000000000 --- a/doc/content/source/bcs/INSOutflowBC.md +++ /dev/null @@ -1,23 +0,0 @@ -# INSOutflowBC - -!alert construction title=Undocumented Class -The INSOutflowBC has not been documented. The content listed below should be used as a starting point for -documenting the class, which includes the typical automatic documentation associated with a -MooseObject; however, what is contained is ultimately determined by what is necessary to make the -documentation clear for users. - -!syntax description /BCs/INSOutflowBC - -## Overview - -!! Replace these lines with information regarding the INSOutflowBC object. - -## Example Input File Syntax - -!! Describe and include an example of how to use the INSOutflowBC object. - -!syntax parameters /BCs/INSOutflowBC - -!syntax inputs /BCs/INSOutflowBC - -!syntax children /BCs/INSOutflowBC diff --git a/doc/content/source/bcs/INSSymmetryAxisBC.md b/doc/content/source/bcs/INSSymmetryAxisBC.md deleted file mode 100644 index 3810efbfbe..0000000000 --- a/doc/content/source/bcs/INSSymmetryAxisBC.md +++ /dev/null @@ -1,23 +0,0 @@ -# INSSymmetryAxisBC - -!alert construction title=Undocumented Class -The INSSymmetryAxisBC has not been documented. The content listed below should be used as a starting point for -documenting the class, which includes the typical automatic documentation associated with a -MooseObject; however, what is contained is ultimately determined by what is necessary to make the -documentation clear for users. - -!syntax description /BCs/INSSymmetryAxisBC - -## Overview - -!! Replace these lines with information regarding the INSSymmetryAxisBC object. - -## Example Input File Syntax - -!! Describe and include an example of how to use the INSSymmetryAxisBC object. - -!syntax parameters /BCs/INSSymmetryAxisBC - -!syntax inputs /BCs/INSSymmetryAxisBC - -!syntax children /BCs/INSSymmetryAxisBC diff --git a/doc/content/source/bcs/SATurbulentViscosityNoBCBC.md b/doc/content/source/bcs/SATurbulentViscosityNoBCBC.md new file mode 100644 index 0000000000..53c3577178 --- /dev/null +++ b/doc/content/source/bcs/SATurbulentViscosityNoBCBC.md @@ -0,0 +1,18 @@ +# SATurbulentViscosityNoBCBC + +!syntax description /BCs/SATurbulentViscosityNoBCBC + +## Overview + +This object applies the "no BC" boundary condition for the Spalart-Allmaras turbulence model +equation. + +## Example Input File Syntax + +!! Describe and include an example of how to use the SATurbulentViscosityNoBCBC object. + +!syntax parameters /BCs/SATurbulentViscosityNoBCBC + +!syntax inputs /BCs/SATurbulentViscosityNoBCBC + +!syntax children /BCs/SATurbulentViscosityNoBCBC diff --git a/doc/content/source/bcs/SATurbulentViscositySUPGBC.md b/doc/content/source/bcs/SATurbulentViscositySUPGBC.md new file mode 100644 index 0000000000..9b76ed3762 --- /dev/null +++ b/doc/content/source/bcs/SATurbulentViscositySUPGBC.md @@ -0,0 +1,18 @@ +# INSADMomentumSUPGBC + +!syntax description /BCs/INSADMomentumSUPGBC + +## Overview + +This object removes residual contributions from SUPG stabilization in the Spalart-Allmaras +turbulence model equation along the boundaries. + +## Example Input File Syntax + +!! Describe and include an example of how to use the INSADMomentumSUPGBC object. + +!syntax parameters /BCs/INSADMomentumSUPGBC + +!syntax inputs /BCs/INSADMomentumSUPGBC + +!syntax children /BCs/INSADMomentumSUPGBC diff --git a/doc/content/source/dgkernels/DGTurbulentDiffusion.md b/doc/content/source/dgkernels/DGTurbulentDiffusion.md new file mode 100644 index 0000000000..3b3317fb8f --- /dev/null +++ b/doc/content/source/dgkernels/DGTurbulentDiffusion.md @@ -0,0 +1,20 @@ +# DGTurbulentDiffusion + +!syntax description /DGKernels/DGTurbulentDiffusion + +## Overview + +This object adds the $-\nabla D_t \nabla u$ turbulent diffusion term based on turbulence predicted +by the Spalart-Allmaras turbulence model. The $\epsilon$ and $\sigma$ +values may be adjusted to attain the desired level of artificial diffusion for eliminating +undershoots and overshoots in the solution. + +## Example Input File Syntax + +!! Describe and include an example of how to use the DGTurbulentDiffusion object. + +!syntax parameters /DGKernels/DGTurbulentDiffusion + +!syntax inputs /DGKernels/DGTurbulentDiffusion + +!syntax children /DGKernels/DGTurbulentDiffusion diff --git a/doc/content/source/kernels/INSADEnergyTurbulentDiffusion.md b/doc/content/source/kernels/INSADEnergyTurbulentDiffusion.md new file mode 100644 index 0000000000..da803c6200 --- /dev/null +++ b/doc/content/source/kernels/INSADEnergyTurbulentDiffusion.md @@ -0,0 +1,19 @@ +# INSADEnergyTurbulentDiffusion + +!syntax description /Kernels/INSADEnergyTurbulentDiffusion + +## Overview + +This object adds the $-\nabla D_t \nabla T$ turbulent diffusion term, based on turbulence predicted +by the Spalart-Allmaras turbulence model, of +the temperature advection-diffusion equation under the INSAD implementation. + +## Example Input File Syntax + +!! Describe and include an example of how to use the INSADEnergyTurbulentDiffusion object. + +!syntax parameters /Kernels/INSADEnergyTurbulentDiffusion + +!syntax inputs /Kernels/INSADEnergyTurbulentDiffusion + +!syntax children /Kernels/INSADEnergyTurbulentDiffusion diff --git a/doc/content/source/kernels/INSADMomentumCD.md b/doc/content/source/kernels/INSADMomentumCD.md new file mode 100644 index 0000000000..36714dc9bd --- /dev/null +++ b/doc/content/source/kernels/INSADMomentumCD.md @@ -0,0 +1,18 @@ +# INSADMomentumCD + +!syntax description /Kernels/INSADMomentumCD + +## Overview + +This object adds a crosswind diffusion term to the incompressible Navier-Stokes equation under the +INSAD implementation. + +## Example Input File Syntax + +!! Describe and include an example of how to use the INSADMomentumCD object. + +!syntax parameters /Kernels/INSADMomentumCD + +!syntax inputs /Kernels/INSADMomentumCD + +!syntax children /Kernels/INSADMomentumCD diff --git a/doc/content/source/kernels/INSADMomentumTurbulentViscous.md b/doc/content/source/kernels/INSADMomentumTurbulentViscous.md new file mode 100644 index 0000000000..a461489fed --- /dev/null +++ b/doc/content/source/kernels/INSADMomentumTurbulentViscous.md @@ -0,0 +1,20 @@ +# INSADMomentumTurbulentViscous + +!syntax description /Kernels/INSADMomentumTurbulentViscous + +## Overview + +This object adds the $-\mu_t\left(\nabla\vec{u}+\left(\nabla\vec{u}\right)^T\right)$ turbulent +viscous term of the incompressible Navier-Stokes momentum equation under the INSAD implementation. +The term is integrated by parts. Only the traction form is compatible with the Spalart-Allmaras +turbulence model because $\mu_t$ is generally not constant in space. + +## Example Input File Syntax + +!! Describe and include an example of how to use the INSADMomentumTurbulentViscous object. + +!syntax parameters /Kernels/INSADMomentumTurbulentViscous + +!syntax inputs /Kernels/INSADMomentumTurbulentViscous + +!syntax children /Kernels/INSADMomentumTurbulentViscous diff --git a/doc/content/source/kernels/INSMomentumKEpsilon.md b/doc/content/source/kernels/INSMomentumKEpsilon.md deleted file mode 100644 index 6a2acea8d2..0000000000 --- a/doc/content/source/kernels/INSMomentumKEpsilon.md +++ /dev/null @@ -1,23 +0,0 @@ -# INSMomentumKEpsilon - -!alert construction title=Undocumented Class -The INSMomentumKEpsilon has not been documented. The content listed below should be used as a starting point for -documenting the class, which includes the typical automatic documentation associated with a -MooseObject; however, what is contained is ultimately determined by what is necessary to make the -documentation clear for users. - -!syntax description /Kernels/INSMomentumKEpsilon - -## Overview - -!! Replace these lines with information regarding the INSMomentumKEpsilon object. - -## Example Input File Syntax - -!! Describe and include an example of how to use the INSMomentumKEpsilon object. - -!syntax parameters /Kernels/INSMomentumKEpsilon - -!syntax inputs /Kernels/INSMomentumKEpsilon - -!syntax children /Kernels/INSMomentumKEpsilon diff --git a/doc/content/source/kernels/SATimeDerivative.md b/doc/content/source/kernels/SATimeDerivative.md new file mode 100644 index 0000000000..379a01365c --- /dev/null +++ b/doc/content/source/kernels/SATimeDerivative.md @@ -0,0 +1,18 @@ +# SATimeDerivative + +!syntax description /Kernels/SATimeDerivative + +## Overview + +This object adds the $\rho\frac{\partial\tilde{\mu}}{\partial t}$ time derivative term of the +Spalart-Allmaras turbulence model. + +## Example Input File Syntax + +!! Describe and include an example of how to use the SATimeDerivative object. + +!syntax parameters /Kernels/SATimeDerivative + +!syntax inputs /Kernels/SATimeDerivative + +!syntax children /Kernels/SATimeDerivative diff --git a/doc/content/source/kernels/SATurbulentViscosity.md b/doc/content/source/kernels/SATurbulentViscosity.md new file mode 100644 index 0000000000..b096eb2a9f --- /dev/null +++ b/doc/content/source/kernels/SATurbulentViscosity.md @@ -0,0 +1,17 @@ +# SATurbulentViscosity + +!syntax description /Kernels/SATurbulentViscosity + +## Overview + +This object adds the Spalart-Allmaras turbulent viscosity equation as described [here](theory.md). + +## Example Input File Syntax + +!! Describe and include an example of how to use the SATurbulentViscosity object. + +!syntax parameters /Kernels/SATurbulentViscosity + +!syntax inputs /Kernels/SATurbulentViscosity + +!syntax children /Kernels/SATurbulentViscosity diff --git a/doc/content/source/kernels/SATurbulentViscosityCD.md b/doc/content/source/kernels/SATurbulentViscosityCD.md new file mode 100644 index 0000000000..9b6d21386c --- /dev/null +++ b/doc/content/source/kernels/SATurbulentViscosityCD.md @@ -0,0 +1,18 @@ +# SATurbulentViscosityCD + +!syntax description /Kernels/SATurbulentViscosityCD + +## Overview + +This object adds a crosswind diffusion term to the Spalart-Allmaras model turbulent viscosity +equation. + +## Example Input File Syntax + +!! Describe and include an example of how to use the SATurbulentViscosityCD object. + +!syntax parameters /Kernels/SATurbulentViscosityCD + +!syntax inputs /Kernels/SATurbulentViscosityCD + +!syntax children /Kernels/SATurbulentViscosityCD diff --git a/doc/content/source/materials/MoltresJsonMaterial.md b/doc/content/source/materials/MoltresJsonMaterial.md index ee5c353b24..40cef8a840 100644 --- a/doc/content/source/materials/MoltresJsonMaterial.md +++ b/doc/content/source/materials/MoltresJsonMaterial.md @@ -1,16 +1,28 @@ # MoltresJsonMaterial -!alert construction title=Undocumented Class -The MoltresJsonMaterial has not been documented. The content listed below should be used as a starting point for -documenting the class, which includes the typical automatic documentation associated with a -MooseObject; however, what is contained is ultimately determined by what is necessary to make the -documentation clear for users. - !syntax description /Materials/MoltresJsonMaterial ## Overview -!! Replace these lines with information regarding the MoltresJsonMaterial object. +This material class preloads and interpolates JSON-based group constant data provided at various +temperature. `MoltresJsonMaterial` declares the interpolated group constant data as material +properties which may be accessed by the neutron diffusion kernels to model a reactor and simulate +temperature reactivity feedback. Refer to the "Materials Block" subsection of the +[Tutorial 2a](input_syntax.md) for more information. + +!table id=group_const caption=Relevant group constants for neutron energy group $g$ or delayed neutron precursor group $i$ +| Group constant | Formula | Label | +| - | - | - | +| Macroscopic fission cross section | $\Sigma_{f,g}$ | FISSXS | +| Macroscopic removal cross section | $\Sigma_{r,g}$ | REMXS | +| Macroscopic scattering cross section | $\Sigma_s^{g'\rightarrow g}$ | GTRANSFXS | +| Neutron diffusion coefficient | $D_g$ | DIFFCOEF | +| Inverse neutron velocity | $\frac{1}{v_g}$ | RECIPVEL | +| Fission neutron spectrum (total) | $\chi_{t,g}$ | CHI_T | +| Fission neutron spectrum (prompt) | $\chi_{p,g}$ | CHI_P | +| Fission neutron spectrum (delayed) | $\chi_{d,g}$ | CHI_D | +| Effective delayed neutron fraction | $\beta_{eff}$ | BETA_EFF | +| Delayed neutron precursor decay constant | $\lambda_i$ | DECAY_CONSTANT | ## Example Input File Syntax diff --git a/doc/content/source/materials/SATauMaterial.md b/doc/content/source/materials/SATauMaterial.md new file mode 100644 index 0000000000..d07951f036 --- /dev/null +++ b/doc/content/source/materials/SATauMaterial.md @@ -0,0 +1,23 @@ +# SATauMaterial + +!syntax description /Materials/SATauMaterial + +## Overview + +This material class computes strong residuals for the +[Spalart-Allmaras turbulence model equation](theory.md) +and adjusts the $\tau$ stabilization parameter for use in pressure-stabilized (PSPG) and +streamline-upwind (SUPG) kernels. `SATauMaterial` includes all functionality from +[INSADTauMaterial](https://mooseframework.inl.gov/source/materials/INSADTauMaterial.html) for +incompressible Navier-Stokes flow modeling with INSAD kernels from the MOOSE `Navier-Stokes` +physics module. + +## Example Input File Syntax + +!! Describe and include an example of how to use the SATauMaterial object. + +!syntax parameters /Materials/SATauMaterial + +!syntax inputs /Materials/SATauMaterial + +!syntax children /Materials/SATauMaterial diff --git a/doc/content/source/materials/SATauStabilized3Eqn.md b/doc/content/source/materials/SATauStabilized3Eqn.md new file mode 100644 index 0000000000..c7fe5e06b6 --- /dev/null +++ b/doc/content/source/materials/SATauStabilized3Eqn.md @@ -0,0 +1,22 @@ +# SATauStabilized3Eqn + +!syntax description /Materials/SATauStabilized3Eqn + +## Overview + +This material class computes strong residuals for the temperature advection-diffusion equation +and adjusts the $\tau$ stabilization parameter for use in pressure-stabilized (PSPG) and +streamline-upwind (SUPG) kernels. `SATauStabilized3Eqn` includes all functionality from +[INSADStabilized3Eqn](https://mooseframework.inl.gov/source/materials/INSADStabilized3Eqn.html) for +thermal-hydraulics modeling with INSAD kernels from the MOOSE `Navier-Stokes` +physics module and [SATauMaterial](SATauMaterial.md) for Spalart-Allmaras turbulence modeling. + +## Example Input File Syntax + +!! Describe and include an example of how to use the SATauStabilized3Eqn object. + +!syntax parameters /Materials/SATauStabilized3Eqn + +!syntax inputs /Materials/SATauStabilized3Eqn + +!syntax children /Materials/SATauStabilized3Eqn diff --git a/doc/content/syntax/index.md b/doc/content/syntax/index.md index 24a20ff24f..e208766b66 100644 --- a/doc/content/syntax/index.md +++ b/doc/content/syntax/index.md @@ -1,3 +1,3 @@ # Complete Syntax -!syntax complete groups=MoltresApp SquirrelApp +!content location=source diff --git a/include/auxkernels/SATurbulentViscosityAux.h b/include/auxkernels/SATurbulentViscosityAux.h new file mode 100644 index 0000000000..d5553fcc2d --- /dev/null +++ b/include/auxkernels/SATurbulentViscosityAux.h @@ -0,0 +1,20 @@ +#pragma once + +#include "AuxKernel.h" + +/* + * Computes the turbulent dynamic viscosity for the Spalart-Allmaras turbulence model. + */ +class SATurbulentViscosityAux : public AuxKernel +{ +public: + static InputParameters validParams(); + + SATurbulentViscosityAux(const InputParameters & parameters); + +protected: + virtual Real computeValue(); + + const VariableValue & _mu_tilde; + const ADMaterialProperty & _mu; +}; diff --git a/include/auxkernels/TurbulentStressAux.h b/include/auxkernels/TurbulentStressAux.h new file mode 100644 index 0000000000..098559dfad --- /dev/null +++ b/include/auxkernels/TurbulentStressAux.h @@ -0,0 +1,22 @@ +#pragma once + +#include "AuxKernel.h" + +/** + * Auxiliary kernel for computing the norm of turbulent stress + */ +class TurbulentStressAux : public AuxKernel +{ +public: + static InputParameters validParams(); + + TurbulentStressAux(const InputParameters & parameters); + +protected: + virtual Real computeValue(); + + const ADVariableValue & _mu_tilde; + const ADMaterialProperty & _mu; + const ADMaterialProperty & _rho; + const ADMaterialProperty & _strain_mag; +}; diff --git a/include/auxkernels/WallDistanceAux.h b/include/auxkernels/WallDistanceAux.h new file mode 100644 index 0000000000..3c2d0a2f42 --- /dev/null +++ b/include/auxkernels/WallDistanceAux.h @@ -0,0 +1,19 @@ +#pragma once + +#include "AuxKernel.h" + +/* + * Computes the minimum wall distance of the element from the wall boundaries. + */ +class WallDistanceAux : public AuxKernel +{ +public: + static InputParameters validParams(); + + WallDistanceAux(const InputParameters & parameters); + +protected: + virtual Real computeValue() override; + + std::vector _wall_boundary_names; +}; diff --git a/include/auxkernels/WallShearStressAux.h b/include/auxkernels/WallShearStressAux.h new file mode 100644 index 0000000000..f9800ba6ea --- /dev/null +++ b/include/auxkernels/WallShearStressAux.h @@ -0,0 +1,21 @@ +#pragma once + +#include "AuxKernel.h" + +/** + * Auxiliary kernel for computing the wall shear stress along a given boundary. + */ +class WallShearStressAux : public AuxKernel +{ +public: + static InputParameters validParams(); + + WallShearStressAux(const InputParameters & parameters); + +protected: + virtual Real computeValue(); + + const ADVectorVariableGradient & _grad_velocity; + const ADMaterialProperty & _mu; + const MooseArray & _normals; +}; diff --git a/include/bcs/INSADMassPSPGBC.h b/include/bcs/INSADMassPSPGBC.h new file mode 100644 index 0000000000..0e27e104bd --- /dev/null +++ b/include/bcs/INSADMassPSPGBC.h @@ -0,0 +1,21 @@ +#pragma once + +#include "ADIntegratedBC.h" + +/** + * This class removes INS PSPG contributions along the boundary. + */ +class INSADMassPSPGBC : public ADIntegratedBC +{ +public: + static InputParameters validParams(); + + INSADMassPSPGBC(const InputParameters & parameters); + +protected: + ADReal computeQpResidual() override; + + const ADMaterialProperty & _rho; + const ADMaterialProperty & _tau; + const ADMaterialProperty & _momentum_strong_residual; +}; diff --git a/include/bcs/INSADMomentumSUPGBC.h b/include/bcs/INSADMomentumSUPGBC.h new file mode 100644 index 0000000000..895ea8004d --- /dev/null +++ b/include/bcs/INSADMomentumSUPGBC.h @@ -0,0 +1,25 @@ +#pragma once + +#include "ADVectorIntegratedBC.h" + +/** + * This class removes INS SUPG contributions along the boundary. + */ +class INSADMomentumSUPGBC : public ADVectorIntegratedBC +{ +public: + static InputParameters validParams(); + + INSADMomentumSUPGBC(const InputParameters & parameters); + +protected: + ADRealVectorValue computeQpStabilization(); + ADRealVectorValue precomputeQpStrongResidual(); + void computeResidual() override; + void computeResidualsForJacobian() override; + virtual ADReal computeQpResidual() override final; + + const ADMaterialProperty & _rho; + const ADMaterialProperty & _tau; + const ADMaterialProperty & _momentum_strong_residual; +}; diff --git a/include/bcs/INSADMomentumTurbulentNoBCBC.h b/include/bcs/INSADMomentumTurbulentNoBCBC.h new file mode 100644 index 0000000000..1870e8dd84 --- /dev/null +++ b/include/bcs/INSADMomentumTurbulentNoBCBC.h @@ -0,0 +1,20 @@ +#pragma once + +#include "INSADMomentumNoBCBC.h" + +/** + * This class implements the "No BC" boundary condition with additional contributions from + * turbulent viscosity. + */ +class INSADMomentumTurbulentNoBCBC : public INSADMomentumNoBCBC +{ +public: + static InputParameters validParams(); + + INSADMomentumTurbulentNoBCBC(const InputParameters & parameters); + +protected: + ADReal computeQpResidual() override; + + const ADVariableValue & _mu_tilde; +}; diff --git a/include/bcs/INSOutflowBC.h b/include/bcs/INSOutflowBC.h deleted file mode 100644 index 06ae0bc302..0000000000 --- a/include/bcs/INSOutflowBC.h +++ /dev/null @@ -1,50 +0,0 @@ -#pragma once - -#include "IntegratedBC.h" - -/** - * This class computes momentum equation residual and Jacobian - * contributions for the incompressible Navier-Stokes momentum - * equation with a standard k-epsilon turbulence model. - */ -class INSOutflowBC : public IntegratedBC -{ -public: - INSOutflowBC(const InputParameters & parameters); - - static InputParameters validParams(); - - virtual ~INSOutflowBC() {} - -protected: - virtual Real computeQpResidual() override; - virtual Real computeQpJacobian() override; - virtual Real computeQpOffDiagJacobian(unsigned jvar) override; - - // Coupled variables - const VariableValue & _u_vel; - const VariableValue & _v_vel; - const VariableValue & _w_vel; - const VariableValue & _p; - - // Gradients - const VariableGradient & _grad_u_vel; - const VariableGradient & _grad_v_vel; - const VariableGradient & _grad_w_vel; - const VariableGradient & _grad_p; - - // Variable numberings - unsigned int _u_vel_var_number; - unsigned int _v_vel_var_number; - unsigned int _w_vel_var_number; - unsigned int _p_var_number; - - // Material properties - // MaterialProperty & _dynamic_viscosity; - Real _mu; - Real _rho; - - // Parameters - unsigned _component; - bool _integrate_p_by_parts; -}; diff --git a/include/bcs/INSSymmetryAxisBC.h b/include/bcs/INSSymmetryAxisBC.h deleted file mode 100644 index a9adffcd81..0000000000 --- a/include/bcs/INSSymmetryAxisBC.h +++ /dev/null @@ -1,50 +0,0 @@ -#pragma once - -#include "IntegratedBC.h" - -/** - * This class computes momentum equation residual and Jacobian - * contributions for the incompressible Navier-Stokes momentum - * equation with a standard k-epsilon turbulence model. - */ -class INSSymmetryAxisBC : public IntegratedBC -{ -public: - INSSymmetryAxisBC(const InputParameters & parameters); - - static InputParameters validParams(); - - virtual ~INSSymmetryAxisBC() {} - -protected: - virtual Real computeQpResidual() override; - virtual Real computeQpJacobian() override; - virtual Real computeQpOffDiagJacobian(unsigned jvar) override; - - // Coupled variables - const VariableValue & _u_vel; - const VariableValue & _v_vel; - const VariableValue & _w_vel; - const VariableValue & _p; - - // Gradients - const VariableGradient & _grad_u_vel; - const VariableGradient & _grad_v_vel; - const VariableGradient & _grad_w_vel; - const VariableGradient & _grad_p; - - // Variable numberings - unsigned int _u_vel_var_number; - unsigned int _v_vel_var_number; - unsigned int _w_vel_var_number; - unsigned int _p_var_number; - - // Material properties - // MaterialProperty & _dynamic_viscosity; - Real _mu; - Real _rho; - - // Parameters - unsigned _component; - bool _integrate_p_by_parts; -}; diff --git a/include/bcs/SATurbulentViscosityNoBCBC.h b/include/bcs/SATurbulentViscosityNoBCBC.h new file mode 100644 index 0000000000..3d7a609f7e --- /dev/null +++ b/include/bcs/SATurbulentViscosityNoBCBC.h @@ -0,0 +1,19 @@ +#pragma once + +#include "ADVectorIntegratedBC.h" + +/** + * This class implements the "No BC" boundary condition of the Spalart-Allmaras equation + */ +class SATurbulentViscosityNoBCBC : public ADIntegratedBC +{ +public: + static InputParameters validParams(); + + SATurbulentViscosityNoBCBC(const InputParameters & parameters); + +protected: + ADReal computeQpResidual() override; + + const ADMaterialProperty & _mu; +}; diff --git a/include/bcs/SATurbulentViscositySUPGBC.h b/include/bcs/SATurbulentViscositySUPGBC.h new file mode 100644 index 0000000000..5b6f6d9958 --- /dev/null +++ b/include/bcs/SATurbulentViscositySUPGBC.h @@ -0,0 +1,25 @@ +#pragma once + +#include "ADIntegratedBC.h" + +/** + * This class removes SUPG contributions in the Spalart-Allmaras equation along the boundary. + */ +class SATurbulentViscositySUPGBC : public ADIntegratedBC +{ +public: + static InputParameters validParams(); + + SATurbulentViscositySUPGBC(const InputParameters & parameters); + +protected: + ADRealVectorValue computeQpStabilization(); + ADReal precomputeQpStrongResidual(); + void computeResidual() override; + void computeResidualsForJacobian() override; + virtual ADReal computeQpResidual() override final; + + const ADMaterialProperty & _tau_visc; + const ADVectorVariableValue & _velocity; + const ADMaterialProperty & _visc_strong_residual; +}; diff --git a/include/dgkernels/DGTurbulentDiffusion.h b/include/dgkernels/DGTurbulentDiffusion.h new file mode 100644 index 0000000000..91200ad155 --- /dev/null +++ b/include/dgkernels/DGTurbulentDiffusion.h @@ -0,0 +1,31 @@ +#pragma once + +#include "DGKernel.h" + +/** + * Computes residual contributions of the turbulent diffusion term in the delayed neutron + * precursor equation using the discontinuous Galerkin method + */ +class DGTurbulentDiffusion : public DGKernel +{ +public: + static InputParameters validParams(); + + DGTurbulentDiffusion(const InputParameters & parameters); + +protected: + virtual Real computeQpResidual(Moose::DGResidualType type) override; + virtual Real computeQpJacobian(Moose::DGJacobianType type) override; + virtual Real computeQpOffDiagJacobian(Moose::DGJacobianType type, unsigned jvar) override; + + Real _epsilon; + Real _sigma; + const ADMaterialProperty & _mu; + const ADMaterialProperty & _mu_nb; + const ADMaterialProperty & _rho; + const ADMaterialProperty & _rho_nb; + const ADVariableValue & _mu_tilde; + const ADVariableValue & _mu_tilde_nb; + Real _sc; + unsigned int _mu_tilde_var_number; +}; diff --git a/include/kernels/INSADEnergyTurbulentDiffusion.h b/include/kernels/INSADEnergyTurbulentDiffusion.h new file mode 100644 index 0000000000..ab12cfd5d2 --- /dev/null +++ b/include/kernels/INSADEnergyTurbulentDiffusion.h @@ -0,0 +1,22 @@ +#pragma once + +#include "ADDiffusion.h" + +/** + * Computes residual contributions of the turbulent diffusion term in the INSAD energy equation. + */ +class INSADEnergyTurbulentDiffusion : public ADDiffusion +{ +public: + static InputParameters validParams(); + + INSADEnergyTurbulentDiffusion(const InputParameters & parameters); + +protected: + ADRealVectorValue precomputeQpResidual() override; + + const ADMaterialProperty & _cp; + const ADMaterialProperty & _mu; + const ADVariableValue & _mu_tilde; + const Real _pr; +}; diff --git a/include/kernels/INSADMomentumTurbulentViscous.h b/include/kernels/INSADMomentumTurbulentViscous.h new file mode 100644 index 0000000000..1fe2831625 --- /dev/null +++ b/include/kernels/INSADMomentumTurbulentViscous.h @@ -0,0 +1,32 @@ +#pragma once + +#include "INSADMomentumViscous.h" + +/** + * This class computes the momentum equation residual and Jacobian + * contributions for the turbulent viscous term of the incompressible + * Navier-Stokes momentum equation. + */ +class INSADMomentumTurbulentViscous : public INSADMomentumViscous +{ +public: + static InputParameters validParams(); + + INSADMomentumTurbulentViscous(const InputParameters & parameters); + +protected: + ADReal computeQpResidual() override; + + /** + * Computes the cartesian coordinate system viscous term + */ + ADRealTensorValue qpViscousTerm(); + + /** + * Computes an additional contribution to the viscous term present of RZ coordinate systems + * (assumes axisymmetric) + */ + ADRealVectorValue qpAdditionalRZTerm(); + + const ADVariableValue & _mu_tilde; +}; diff --git a/include/kernels/INSMomentumKEpsilon.h b/include/kernels/INSMomentumKEpsilon.h deleted file mode 100644 index 66f22a744e..0000000000 --- a/include/kernels/INSMomentumKEpsilon.h +++ /dev/null @@ -1,54 +0,0 @@ -#pragma once - -#include "Kernel.h" - -/** - * This class computes momentum equation residual and Jacobian - * contributions for the incompressible Navier-Stokes momentum - * equation with a standard k-epsilon turbulence model. - */ -class INSMomentumKEpsilon : public Kernel -{ -public: - INSMomentumKEpsilon(const InputParameters & parameters); - - static InputParameters validParams(); - - virtual ~INSMomentumKEpsilon() {} - -protected: - virtual Real computeQpResidual() override; - virtual Real computeQpJacobian() override; - virtual Real computeQpOffDiagJacobian(unsigned jvar) override; - - // Coupled variables - const VariableValue & _u_vel; - const VariableValue & _v_vel; - const VariableValue & _w_vel; - const VariableValue & _p; - const VariableValue & _k; - - // Gradients - const VariableGradient & _grad_u_vel; - const VariableGradient & _grad_v_vel; - const VariableGradient & _grad_w_vel; - const VariableGradient & _grad_p; - const VariableGradient & _grad_k; - - // Variable numberings - unsigned int _u_vel_var_number; - unsigned int _v_vel_var_number; - unsigned int _w_vel_var_number; - unsigned int _p_var_number; - unsigned int _k_var_number; - - // Material properties - // MaterialProperty & _dynamic_viscosity; - Real _mu; - Real _rho; - RealVectorValue _gravity; - - // Parameters - unsigned _component; - bool _integrate_p_by_parts; -}; diff --git a/include/kernels/SATimeDerivative.h b/include/kernels/SATimeDerivative.h new file mode 100644 index 0000000000..a1310dd9bf --- /dev/null +++ b/include/kernels/SATimeDerivative.h @@ -0,0 +1,20 @@ +#pragma once + +#include "ADTimeDerivative.h" + +/** + * This object computes the residual and Jacobian contributions for the time derivative term + * in the turbulent viscosity equation of the Spalart-Allmaras model. + */ +class SATimeDerivative : public ADTimeDerivative +{ +public: + SATimeDerivative(const InputParameters & parameters); + + static InputParameters validParams(); + +protected: + virtual ADReal precomputeQpResidual() override; + + const ADMaterialProperty & _time_strong_residual; +}; diff --git a/include/kernels/SATurbulentViscosity.h b/include/kernels/SATurbulentViscosity.h new file mode 100644 index 0000000000..05e04fdb09 --- /dev/null +++ b/include/kernels/SATurbulentViscosity.h @@ -0,0 +1,24 @@ +#pragma once + +#include "ADKernel.h" + +/** + * Computes residual contributions of all terms in the turbulent viscosity + * equation of the Spalart-Allmaras model. + */ +class SATurbulentViscosity : public ADKernel +{ +public: + static InputParameters validParams(); + + SATurbulentViscosity(const InputParameters & parameters); + +protected: + virtual ADReal computeQpResidual() override; + + const ADMaterialProperty & _mu; + const ADMaterialProperty & _convection_strong_residual; + const ADMaterialProperty & _destruction_strong_residual; + const ADMaterialProperty & _diffusion_strong_residual; + const ADMaterialProperty & _source_strong_residual; +}; diff --git a/include/kernels/SATurbulentViscositySUPG.h b/include/kernels/SATurbulentViscositySUPG.h new file mode 100644 index 0000000000..ab8c53afd2 --- /dev/null +++ b/include/kernels/SATurbulentViscositySUPG.h @@ -0,0 +1,20 @@ +#pragma once + +#include "ADKernelSUPG.h" + +/** + * This class computes the residual and Jacobian contributions for + * the Spalart-Allmaras equation stabilization + */ +class SATurbulentViscositySUPG : public ADKernelSUPG +{ +public: + static InputParameters validParams(); + + SATurbulentViscositySUPG(const InputParameters & parameters); + +protected: + ADReal precomputeQpStrongResidual() override; + + const ADMaterialProperty & _visc_strong_residual; +}; diff --git a/include/materials/SATauMaterial.h b/include/materials/SATauMaterial.h new file mode 100644 index 0000000000..5d5cc328af --- /dev/null +++ b/include/materials/SATauMaterial.h @@ -0,0 +1,59 @@ +#pragma once + +#include "INSADTauMaterial.h" + +using MetaPhysicL::raw_value; + +template +class SATauMaterialTempl : public T +{ +public: + static InputParameters validParams(); + + SATauMaterialTempl(const InputParameters & parameters); + + void subdomainSetup() override; + +protected: + virtual void computeQpProperties() override; + + const Real _sigma; + const Real _cb1; + const Real _cb2; + const Real _kappa; + const Real _cw1; + const Real _cw2; + const Real _cw3; + const Real _cv1; + const Real _ct1; + const Real _ct2; + const Real _ct3; + const Real _ct4; + const ADVariableValue & _mu_tilde; + const ADVariableGradient & _grad_mu; + const ADVariableValue * _visc_dot; + ADMaterialProperty & _tau_visc; + const VariableValue & _wall_dist; + const bool _use_ft2_term; + + ADMaterialProperty & _strain_mag; + ADMaterialProperty & _convection_strong_residual; + ADMaterialProperty & _destruction_strong_residual; + ADMaterialProperty & _diffusion_strong_residual; + ADMaterialProperty & _source_strong_residual; + ADMaterialProperty & _time_strong_residual; + ADMaterialProperty & _visc_strong_residual; + + using T::_alpha; + using T::_dt; + using T::_has_transient; + using T::_hmax; + using T::_mu; + using T::_qp; + using T::_rho; + using T::_velocity; + using T::_grad_velocity; + using T::_tau; +}; + +typedef SATauMaterialTempl SATauMaterial; diff --git a/include/materials/SATauStabilized3Eqn.h b/include/materials/SATauStabilized3Eqn.h new file mode 100644 index 0000000000..41a02d5172 --- /dev/null +++ b/include/materials/SATauStabilized3Eqn.h @@ -0,0 +1,34 @@ +#pragma once + +#include "SATauMaterial.h" + +class INSADMaterial; +class INSADStabilized3Eqn; + +#include "SATauMaterial.h" +#include "INSADStabilized3Eqn.h" + +class SATauStabilized3Eqn : public SATauMaterialTempl +{ +public: + static InputParameters validParams(); + + SATauStabilized3Eqn(const InputParameters & parameters); + +protected: + virtual void computeQpProperties() override; + + const Real _pr; + + using SATauMaterialTempl::_mu_tilde; + using SATauMaterialTempl::_grad_mu; + using SATauMaterialTempl::_mu; + using SATauMaterialTempl::_k; + using SATauMaterialTempl::_rho; + using SATauMaterialTempl::_cp; + using SATauMaterialTempl::_alpha; + using SATauMaterialTempl::_temperature_strong_residual; + using SATauMaterialTempl::_has_energy_transient; + using SATauMaterialTempl::_grad_temperature; + using SATauMaterialTempl::_second_temperature; +}; diff --git a/problems/2023-basic-turbulence-cases/backward-facing-step/bfs_flow_step.i b/problems/2023-basic-turbulence-cases/backward-facing-step/bfs_flow_step.i new file mode 100644 index 0000000000..db4306b3cf --- /dev/null +++ b/problems/2023-basic-turbulence-cases/backward-facing-step/bfs_flow_step.i @@ -0,0 +1,415 @@ +Re = 3.6e4 +density = 1 # kg cm-3 +D = 1 +inlet = 1 +viscosity = '${fparse density * inlet * D / Re}' # dynamic viscosity, mu = nu * rho, kg cm-1 s-1 +alpha = 0.33333 # INS SUPG and PSPG stabilization parameter + +[GlobalParams] + integrate_p_by_parts = false + viscous_form = 'traction' + pressure = p + velocity = vel +[] + +[Mesh] + [block] + type = CartesianMeshGenerator + dim = 2 + dx = '4 0.8 0.4 1.6 0.4 6.8 4 2 2 2.5 3 3.5 25' + ix = '64 16 10 50 10 136 64 25 20 20 20 20 125' + dy = '0.025 0.025 0.9 0.025 0.025 0.025 0.025 0.2 1.25 1 2.5 1.5 1.25 0.2 0.025 0.025' + iy = '16 8 36 5 10 20 5 10 25 20 5 30 25 10 5 20' + [] + [rename] + type = SubdomainBoundingBoxGenerator + input = block + bottom_left = '0 0 0' + top_right = '6 1 0' + block_id = 1 + [] + [delete] + type = BlockDeletionGenerator + input = rename + block = 1 + new_boundary = 'step' + [] + [bottom] + type = SideSetsFromNormalsGenerator + input = delete + normals = '0 -1 0' + new_boundary = 'bottom' + fixed_normal = true + replace = true + [] + [transform] + type = TransformGenerator + input = bottom + transform = TRANSLATE + vector_value = '104 0 0' + [] + [corner_node] + type = ExtraNodesetGenerator + input = transform + new_boundary = 'pinned_node' + coord = '160 0 0' + [] +[] + +[Problem] + type = FEProblem +[] + +[Variables] + [vel] + family = LAGRANGE_VEC + order = FIRST + [] + [p] + [] + [mu_tilde] + initial_condition = '${fparse viscosity * 3}' + [] +[] + +[AuxVariables] + [mu_turbulence] + [] + [wall_dist] + [] + [velx] + [] + [vely] + [] + [wall_shear_stress] + family = MONOMIAL + order = CONSTANT + [] + [turbulent_stress] + [] +[] + +[Kernels] + [mass] + type = INSADMass + variable = p + [] + [mass_pspg] + type = INSADMassPSPG + variable = p + [] + [momentum_advection] + type = INSADMomentumAdvection + variable = vel + [] + [momentum_viscous] + type = INSADMomentumViscous + variable = vel + [] + [momentum_turbulent_viscous] + type = INSADMomentumTurbulentViscous + variable = vel + mu_tilde = mu_tilde + [] + [momentum_pressure] + type = INSADMomentumPressure + variable = vel + [] + [momentum_supg] + type = INSADMomentumSUPG + variable = vel + [] + [momentum_time] + type = INSADMomentumTimeDerivative + variable = vel + [] + + [mu_tilde_time] + type = SATimeDerivative + variable = mu_tilde + [] + [mu_tilde_space] + type = SATurbulentViscosity + variable = mu_tilde + [] + [mu_tilde_supg] + type = SATurbulentViscositySUPG + variable = mu_tilde + [] +[] + +[AuxKernels] + [mu_space] + type = SATurbulentViscosityAux + variable = mu_turbulence + mu_tilde = mu_tilde + [] + [wall_distance] + type = WallDistanceAux + variable = wall_dist + walls = 'top bottom step' + execute_on = initial + [] + [velx] + type = VectorVariableComponentAux + variable = velx + vector_variable = vel + component = x + [] + [vely] + type = VectorVariableComponentAux + variable = vely + vector_variable = vel + component = y + [] + [wall_shear_stress] + type = WallShearStressAux + variable = wall_shear_stress + boundary = 'top bottom' + [] + [turbulent_stress] + type = TurbulentStressAux + variable = turbulent_stress + mu_tilde = mu_tilde + [] +[] + +[UserObjects] + [inlet] + type = SolutionUserObject + mesh = bfs_flow_upstream_exodus.e + system_variables = 'velx vely mu_tilde' + timestep = LATEST + execute_on = initial + [] +[] + +[Functions] + [xfunc] + type = ParsedFunction + expression = '${fparse inlet}' + [] + [velx_bc] + type = SolutionFunction + solution = inlet + from_variable = velx + [] + [vely_bc] + type = SolutionFunction + solution = inlet + from_variable = vely + [] + [mu_tilde_bc] + type = SolutionFunction + solution = inlet + from_variable = mu_tilde + [] +[] + +[ICs] + [velocity] + type = VectorFunctionIC + function_x = xfunc + function_y = 0 + variable = vel + [] +[] + +[BCs] + [no_slip] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'top bottom step' + function_x = 0 + function_y = 0 + [] + [inlet] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'left' + function_x = velx_bc + function_y = vely_bc + [] + [outlet] + type = INSADMomentumTurbulentNoBCBC + variable = vel + boundary = 'right' + mu_tilde = mu_tilde + [] + [outlet_supg] + type = INSADMomentumSUPGBC + variable = vel + boundary = 'right' + [] + [outlet_pspg] + type = INSADMassPSPGBC + variable = p + boundary = 'right' + [] + + [pressure_pin] + type = DirichletBC + variable = p + boundary = 'pinned_node' + value = 0 + [] + + [mu_inlet] + type = FunctionDirichletBC + variable = mu_tilde + boundary = 'left' + function = mu_tilde_bc + [] + [mu_wall] + type = DirichletBC + variable = mu_tilde + boundary = 'top bottom step' + value = 0 + [] + [mu_outlet] + type = SATurbulentViscosityNoBCBC + variable = mu_tilde + boundary = 'right' + [] + [mu_outlet_supg] + type = SATurbulentViscositySUPGBC + variable = mu_tilde + boundary = 'right' + [] +[] + +[Materials] + [ad_mat] + type = ADGenericConstantMaterial + prop_names = 'rho mu' + prop_values = '${density} ${viscosity}' + [] + [sa_mat] + type = SATauMaterial + alpha = ${alpha} + mu_tilde = mu_tilde + wall_distance_var = wall_dist + [] +[] + +[VectorPostprocessors] + [vel1] + type = LineValueSampler + variable = 'velx mu_turbulence' + start_point = '106 1 0' + end_point = '106 9 0' + num_points = 25601 + sort_by = y + execute_on = final + [] + [vel2] + type = LineValueSampler + variable = 'velx mu_turbulence' + start_point = '111 0 0' + end_point = '111 9 0' + num_points = 28801 + sort_by = y + execute_on = final + [] + [vel3] + type = LineValueSampler + variable = 'velx mu_turbulence' + start_point = '114 0 0' + end_point = '114 9 0' + num_points = 28801 + sort_by = y + execute_on = final + [] + [vel4] + type = LineValueSampler + variable = 'velx mu_turbulence' + start_point = '116 0 0' + end_point = '116 9 0' + num_points = 28801 + sort_by = y + execute_on = final + [] + [vel5] + type = LineValueSampler + variable = 'velx mu_turbulence' + start_point = '120 0 0' + end_point = '120 9 0' + num_points = 28801 + sort_by = y + execute_on = final + [] + [upstream] + type = LineValueSampler + variable = 'velx p' + start_point = '104 1.00125 0' + end_point = '110 1.00125 0' + num_points = 7201 + sort_by = x + execute_on = final + [] + [downstream] + type = LineValueSampler + variable = 'velx p' + start_point = '110 0.0015625 0' + end_point = '160 0.0015625 0' + num_points = 60001 + sort_by = x + execute_on = final + [] +[] + +[Executioner] + type = Transient + + solve_type = 'NEWTON' + petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor' + petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_package' + petsc_options_value = 'lu NONZERO superlu_dist' + line_search = none + + automatic_scaling = true + compute_scaling_once = false + resid_vs_jac_scaling_param = .1 + scaling_group_variables = 'vel; p; mu_tilde' + off_diagonals_in_auto_scaling = true + + nl_rel_tol = 1e-8 + nl_abs_tol = 1e-12 + nl_max_its = 10 + l_max_its = 200 + + steady_state_detection = true + steady_state_tolerance = 1e-6 + dtmin = 1e-2 + dtmax = 10 + [TimeStepper] + type = IterationAdaptiveDT + dt = 1e-1 + cutback_factor = 0.9 + growth_factor = 1.1 + optimal_iterations = 6 + iteration_window = 1 + linear_iteration_ratio = 1000 + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Outputs] + [exodus] + type = Exodus + [] + [csv] + type = CSV + execute_on = final + [] +[] + +[Debug] + show_var_residual_norms = true +[] diff --git a/problems/2023-basic-turbulence-cases/backward-facing-step/bfs_flow_upstream.i b/problems/2023-basic-turbulence-cases/backward-facing-step/bfs_flow_upstream.i new file mode 100644 index 0000000000..f01a034f91 --- /dev/null +++ b/problems/2023-basic-turbulence-cases/backward-facing-step/bfs_flow_upstream.i @@ -0,0 +1,314 @@ +Re = 3.6e4 +density = 1 # kg cm-3 +D = 1 +inlet = 1 +viscosity = '${fparse density * inlet * D / Re}' # dynamic viscosity, mu = nu * rho, kg cm-1 s-1 +alpha = 0.33333 # INS SUPG and PSPG stabilization parameter + +[GlobalParams] + integrate_p_by_parts = false + viscous_form = 'traction' + pressure = p + velocity = vel +[] + +[Mesh] + [block] + type = CartesianMeshGenerator + dim = 2 + dx = '0.25 0.75 103' + ix = '40 1 103' + dy = '0.025 0.025 0.9 0.025 0.025 0.025 0.025 0.2 1.25 1 2.5 1.5 1.25 0.2 0.025 0.025' + iy = '16 8 36 5 10 20 5 10 25 20 5 30 25 10 5 20' + [] + [rename] + type = SubdomainBoundingBoxGenerator + input = block + bottom_left = '0 0 0' + top_right = '104 1 0' + block_id = 1 + [] + [delete] + type = BlockDeletionGenerator + input = rename + block = 1 + new_boundary = 'bottom' + [] + [corner_node] + type = ExtraNodesetGenerator + input = delete + new_boundary = 'pinned_node' + coord = '104 9' + [] +[] + +[Problem] + type = FEProblem +[] + +[Variables] + [vel] + family = LAGRANGE_VEC + order = FIRST + [] + [p] + [] + [mu_tilde] + initial_condition = '${fparse viscosity * 3}' + [] +[] + +[AuxVariables] + [mu_turbulence] + [] + [wall_dist] + [] + [velx] + [] + [vely] + [] + [wall_shear_stress] + family = MONOMIAL + order = CONSTANT + [] + [turbulent_stress] + [] +[] + +[Kernels] + [mass] + type = INSADMass + variable = p + [] + [mass_pspg] + type = INSADMassPSPG + variable = p + [] + [momentum_advection] + type = INSADMomentumAdvection + variable = vel + [] + [momentum_viscous] + type = INSADMomentumViscous + variable = vel + [] + [momentum_turbulent_viscous] + type = INSADMomentumTurbulentViscous + variable = vel + mu_tilde = mu_tilde + [] + [momentum_pressure] + type = INSADMomentumPressure + variable = vel + [] + [momentum_supg] + type = INSADMomentumSUPG + variable = vel + [] + [momentum_time] + type = INSADMomentumTimeDerivative + variable = vel + [] + + [mu_tilde_time] + type = SATimeDerivative + variable = mu_tilde + [] + [mu_tilde_space] + type = SATurbulentViscosity + variable = mu_tilde + [] + [mu_tilde_supg] + type = SATurbulentViscositySUPG + variable = mu_tilde + [] +[] + +[AuxKernels] + [mu_space] + type = SATurbulentViscosityAux + variable = mu_turbulence + mu_tilde = mu_tilde + [] + [wall_distance] + type = WallDistanceAux + variable = wall_dist + walls = 'top bottom' + execute_on = initial + [] + [velx] + type = VectorVariableComponentAux + variable = velx + vector_variable = vel + component = x + [] + [vely] + type = VectorVariableComponentAux + variable = vely + vector_variable = vel + component = y + [] + [wall_shear_stress] + type = WallShearStressAux + variable = wall_shear_stress + boundary = 'top bottom' + [] + [turbulent_stress] + type = TurbulentStressAux + variable = turbulent_stress + mu_tilde = mu_tilde + [] +[] + +[Functions] + [xfunc] + type = ParsedFunction + expression = '${fparse inlet}' + [] + [mu_func] + type = ParsedFunction + expression = 'if(y=1, 0, if(y=9, 0, ${fparse viscosity * 3}))' + [] +[] + +[ICs] + [velocity] + type = VectorFunctionIC + function_x = xfunc + function_y = 0 + variable = vel + [] +[] + +[BCs] + [no_slip] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'top bottom' + function_x = 0 + function_y = 0 + [] + [inlet] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'left' + function_x = xfunc + function_y = 0 + [] + [outlet] + type = INSADMomentumTurbulentNoBCBC + variable = vel + boundary = 'right' + mu_tilde = mu_tilde + [] + [outlet_supg] + type = INSADMomentumSUPGBC + variable = vel + boundary = 'right' + [] + [outlet_pspg] + type = INSADMassPSPGBC + variable = p + boundary = 'right' + [] + + [pressure_pin] + type = DirichletBC + variable = p + boundary = 'pinned_node' + value = 0 + [] + + [mu_inlet] + type = FunctionDirichletBC + variable = mu_tilde + boundary = 'left' + function = mu_func + [] + [mu_wall] + type = DirichletBC + variable = mu_tilde + boundary = 'top bottom' + value = 0 + [] + [mu_outlet] + type = SATurbulentViscosityNoBCBC + variable = mu_tilde + boundary = 'right' + [] + [mu_outlet_supg] + type = SATurbulentViscositySUPGBC + variable = mu_tilde + boundary = 'right' + [] +[] + +[Materials] + [ad_mat] + type = ADGenericConstantMaterial + prop_names = 'rho mu' + prop_values = '${density} ${viscosity}' + [] + [sa_mat] + type = SATauMaterial + alpha = ${alpha} + mu_tilde = mu_tilde + wall_distance_var = wall_dist + [] +[] + +[Executioner] + type = Transient + + solve_type = 'NEWTON' + petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor' + petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_package' + petsc_options_value = 'lu NONZERO superlu_dist' + line_search = none + + automatic_scaling = true + compute_scaling_once = false + resid_vs_jac_scaling_param = 0.1 + scaling_group_variables = 'vel; p; mu_tilde' + off_diagonals_in_auto_scaling = true + + nl_rel_tol = 1e-8 + nl_abs_tol = 1e-14 + nl_max_its = 10 + l_max_its = 200 + + steady_state_detection = true + steady_state_tolerance = 1e-6 + dtmin = 1e-2 + dtmax = 10 + [TimeStepper] + type = IterationAdaptiveDT + dt = 1e-1 + cutback_factor = 0.9 + growth_factor = 1.1 + optimal_iterations = 6 + iteration_window = 1 + linear_iteration_ratio = 1000 + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Outputs] + [exodus] + type = Exodus + [] + [csv] + type = CSV + execute_on = final + [] +[] + +[Debug] + show_var_residual_norms = true +[] diff --git a/problems/2023-basic-turbulence-cases/channel/channel_flow.i b/problems/2023-basic-turbulence-cases/channel/channel_flow.i new file mode 100644 index 0000000000..98876783cb --- /dev/null +++ b/problems/2023-basic-turbulence-cases/channel/channel_flow.i @@ -0,0 +1,333 @@ +Re = 1.375e4 +density = 1 # kg cm-3 +D = 1 +inlet = 1 +viscosity = '${fparse density * inlet * D / Re}' # dynamic viscosity, mu = nu * rho, kg cm-1 s-1 +alpha = 0.33333 # INS SUPG and PSPG stabilization parameter + +[GlobalParams] + integrate_p_by_parts = false + viscous_form = 'traction' + pressure = p + velocity = vel +[] + +[Mesh] + [channel] + type = CartesianMeshGenerator + dim = 2 + dx = '0.05 0.2 139.75' + dy = '0.3 0.1 0.05 0.05' + ix = '20 1 559' + iy = '20 10 10 40' + [] + [corner_node] + type = ExtraNodesetGenerator + input = channel + new_boundary = 'pinned_node' + coord = '140 0.5' + [] +[] + +[Problem] + type = FEProblem +[] + +[Variables] + [vel] + family = LAGRANGE_VEC + order = FIRST + [] + [p] + [] + [mu_tilde] + initial_condition = '${fparse viscosity * 5}' + [] +[] + +[AuxVariables] + [mu_turbulence] + [] + [wall_dist] + [] + [velx] + [] + [vely] + [] + [wall_shear_stress] + family = MONOMIAL + order = CONSTANT + [] + [turbulent_stress] + family = MONOMIAL + order = CONSTANT + [] +[] + +[Kernels] + [mass] + type = INSADMass + variable = p + [] + [mass_pspg] + type = INSADMassPSPG + variable = p + [] + [momentum_advection] + type = INSADMomentumAdvection + variable = vel + [] + [momentum_viscous] + type = INSADMomentumViscous + variable = vel + [] + [momentum_turbulent_viscous] + type = INSADMomentumTurbulentViscous + variable = vel + mu_tilde = mu_tilde + [] + [momentum_pressure] + type = INSADMomentumPressure + variable = vel + [] + [momentum_supg] + type = INSADMomentumSUPG + variable = vel + [] + [momentum_time] + type = INSADMomentumTimeDerivative + variable = vel + [] + + [mu_tilde_time] + type = SATimeDerivative + variable = mu_tilde + [] + [mu_tilde_space] + type = SATurbulentViscosity + variable = mu_tilde + [] + [mu_tilde_supg] + type = SATurbulentViscositySUPG + variable = mu_tilde + [] +[] + +[AuxKernels] + [mu_space] + type = SATurbulentViscosityAux + variable = mu_turbulence + mu_tilde = mu_tilde + [] + [wall_distance] + type = WallDistanceAux + variable = wall_dist + walls = 'top' + execute_on = initial + [] + [velx] + type = VectorVariableComponentAux + variable = velx + vector_variable = vel + component = x + [] + [vely] + type = VectorVariableComponentAux + variable = vely + vector_variable = vel + component = y + [] + [wall_shear_stress] + type = WallShearStressAux + variable = wall_shear_stress + boundary = top + [] + [turbulent_stress] + type = TurbulentStressAux + variable = turbulent_stress + mu_tilde = mu_tilde + [] +[] + +[Functions] + [xfunc] + type = ParsedFunction + expression = '${fparse inlet}' + [] + [mu_func] + type = ParsedFunction + expression = 'if(y=0.5, 0, ${fparse viscosity * 5})' + [] +[] + +[ICs] + [velocity] + type = VectorFunctionIC + function_x = xfunc + function_y = 0 + variable = vel + [] +[] + +[BCs] + [no_slip] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'top' + function_x = 0 + function_y = 0 + [] + [inlet] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'left' + function_x = xfunc + function_y = 0 + [] + [outlet] + type = INSADMomentumTurbulentNoBCBC + variable = vel + boundary = 'right' + mu_tilde = mu_tilde + [] + [outlet_supg] + type = INSADMomentumSUPGBC + variable = vel + boundary = 'right' + [] + [outlet_pspg] + type = INSADMassPSPGBC + variable = p + boundary = 'right' + [] + [symmetry] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'bottom' + set_x_comp = false + function_y = 0 + [] + + [pressure_pin] + type = DirichletBC + variable = p + boundary = 'pinned_node' + value = 0 + [] + + [mu_inlet] + type = FunctionDirichletBC + variable = mu_tilde + boundary = 'left' + function = mu_func + [] + [mu_wall] + type = DirichletBC + variable = mu_tilde + boundary = 'top' + value = 0 + [] + [mu_outlet] + type = SATurbulentViscosityNoBCBC + variable = mu_tilde + boundary = 'right' + [] + [mu_outlet_supg] + type = SATurbulentViscositySUPGBC + variable = mu_tilde + boundary = 'right' + [] +[] + +[Materials] + [ad_mat] + type = ADGenericConstantMaterial + prop_names = 'rho mu' + prop_values = '${density} ${viscosity}' + [] + [sa_mat] + type = SATauMaterial + alpha = ${alpha} + mu_tilde = mu_tilde + wall_distance_var = wall_dist + [] +[] + +[VectorPostprocessors] + [outlet] + type = NodalValueSampler + variable = velx + boundary = 'right' + sort_by = y + execute_on = final + [] + [prior] + type = NodalValueSampler + variable = velx + boundary = 'right' + sort_by = y + execute_on = final + [] + [stress] + type = SideValueSampler + variable = turbulent_stress + boundary = 'right' + sort_by = y + execute_on = final + [] +[] + +[Executioner] + type = Transient + + solve_type = 'NEWTON' + petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor' + petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_package' + petsc_options_value = 'lu NONZERO superlu_dist' + line_search = none + + automatic_scaling = true + compute_scaling_once = false + resid_vs_jac_scaling_param = 0.1 + scaling_group_variables = 'vel; p; mu_tilde' + off_diagonals_in_auto_scaling = true + + nl_rel_tol = 1e-8 + nl_abs_tol = 1e-12 + nl_max_its = 10 + + steady_state_detection = true + steady_state_tolerance = 1e-8 + dtmin = 1e-2 + dtmax = 10 + [TimeStepper] + type = IterationAdaptiveDT + dt = 1e-1 + cutback_factor = 0.9 + growth_factor = 1.1 + optimal_iterations = 6 + iteration_window = 1 + linear_iteration_ratio = 1000 + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Outputs] + [exodus] + type = Exodus + [] + [csv] + type = CSV + execute_on = final + [] +[] + +[Debug] + show_var_residual_norms = true +[] diff --git a/problems/2023-basic-turbulence-cases/pipe/pipe_flow.i b/problems/2023-basic-turbulence-cases/pipe/pipe_flow.i new file mode 100644 index 0000000000..f402445b1d --- /dev/null +++ b/problems/2023-basic-turbulence-cases/pipe/pipe_flow.i @@ -0,0 +1,320 @@ +Re = 4e4 +density = 1 # kg cm-3 +D = 1 +inlet = 1 +viscosity = '${fparse density * inlet * D / Re}' # dynamic viscosity, mu = nu * rho, kg cm-1 s-1 +alpha = 0.33333 # INS SUPG and PSPG stabilization parameter + +[GlobalParams] + integrate_p_by_parts = false + viscous_form = 'traction' + pressure = p + velocity = vel +[] + +[Mesh] + coord_type = 'RZ' + [pipe] + type = CartesianMeshGenerator + dim = 2 + dx = '0.3 0.1 0.05 0.034 0.016' + dy = '0.025 0.475 149.5' + ix = '20 10 10 17 40' + iy = '20 1 299' + [] + [corner_node] + type = ExtraNodesetGenerator + input = pipe + new_boundary = 'pinned_node' + coord = '0.5 150' + [] +[] + +[Problem] + type = FEProblem +[] + +[Variables] + [vel] + family = LAGRANGE_VEC + order = FIRST + [] + [p] + [] + [mu_tilde] + initial_condition = '${fparse viscosity * 5}' + [] +[] + +[AuxVariables] + [mu_turbulence] + [] + [wall_dist] + [] + [velx] + [] + [vely] + [] + [wall_shear_stress] + family = MONOMIAL + order = CONSTANT + [] + [turbulent_stress] + family = MONOMIAL + order = CONSTANT + [] +[] + +[Kernels] + [mass] + type = INSADMass + variable = p + [] + [mass_pspg] + type = INSADMassPSPG + variable = p + [] + [momentum_advection] + type = INSADMomentumAdvection + variable = vel + [] + [momentum_viscous] + type = INSADMomentumViscous + variable = vel + [] + [momentum_turbulent_viscous] + type = INSADMomentumTurbulentViscous + variable = vel + mu_tilde = mu_tilde + [] + [momentum_pressure] + type = INSADMomentumPressure + variable = vel + [] + [momentum_supg] + type = INSADMomentumSUPG + variable = vel + [] + [momentum_time] + type = INSADMomentumTimeDerivative + variable = vel + [] + + [mu_tilde_time] + type = SATimeDerivative + variable = mu_tilde + [] + [mu_tilde_space] + type = SATurbulentViscosity + variable = mu_tilde + [] + [mu_tilde_supg] + type = SATurbulentViscositySUPG + variable = mu_tilde + [] +[] + +[AuxKernels] + [mu_space] + type = SATurbulentViscosityAux + variable = mu_turbulence + mu_tilde = mu_tilde + [] + [wall_distance] + type = WallDistanceAux + variable = wall_dist + walls = 'right' + execute_on = initial + [] + [velx] + type = VectorVariableComponentAux + variable = velx + vector_variable = vel + component = x + [] + [vely] + type = VectorVariableComponentAux + variable = vely + vector_variable = vel + component = y + [] + [wall_shear_stress] + type = WallShearStressAux + variable = wall_shear_stress + boundary = right + [] + [turbulent_stress] + type = TurbulentStressAux + variable = turbulent_stress + mu_tilde = mu_tilde + [] +[] + +[Functions] + [yfunc] + type = ParsedFunction + expression = '${fparse inlet}' + [] + [mu_func] + type = ParsedFunction + expression = 'if(x=0.5, 0, ${fparse viscosity * 5})' + [] +[] + +[ICs] + [velocity] + type = VectorFunctionIC + function_x = 0 + function_y = yfunc + variable = vel + [] +[] + +[BCs] + [no_slip] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'right' + function_x = 0 + function_y = 0 + [] + [inlet] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'bottom' + function_x = 0 + function_y = yfunc + [] + [outlet] + type = INSADMomentumTurbulentNoBCBC + variable = vel + boundary = 'top' + mu_tilde = mu_tilde + [] + [outlet_supg] + type = INSADMomentumSUPGBC + variable = vel + boundary = 'top' + [] + [outlet_pspg] + type = INSADMassPSPGBC + variable = p + boundary = 'top' + [] + + [pressure_pin] + type = DirichletBC + variable = p + boundary = 'pinned_node' + value = 0 + [] + + [mu_inlet] + type = FunctionDirichletBC + variable = mu_tilde + boundary = 'bottom' + function = mu_func + [] + [mu_wall] + type = DirichletBC + variable = mu_tilde + boundary = 'right' + value = 0 + [] + [mu_outlet] + type = SATurbulentViscosityNoBCBC + variable = mu_tilde + boundary = 'top' + [] + [mu_outlet_supg] + type = SATurbulentViscositySUPGBC + variable = mu_tilde + boundary = 'top' + [] +[] + +[Materials] + [ad_mat] + type = ADGenericConstantMaterial + prop_names = 'rho mu' + prop_values = '${density} ${viscosity}' + [] + [sa_mat] + type = SATauMaterial + alpha = ${alpha} + mu_tilde = mu_tilde + wall_distance_var = wall_dist + [] +[] + +[VectorPostprocessors] + [outlet] + type = NodalValueSampler + variable = vely + boundary = 'top' + sort_by = x + execute_on = final + [] + [stress] + type = SideValueSampler + variable = turbulent_stress + boundary = 'top' + sort_by = x + execute_on = final + [] +[] + +[Executioner] + type = Transient + + solve_type = 'NEWTON' + petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor' + petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_package' + petsc_options_value = 'lu NONZERO superlu_dist' + line_search = none + + automatic_scaling = true + compute_scaling_once = false + resid_vs_jac_scaling_param = 0.1 + scaling_group_variables = 'vel; p; mu_tilde' + off_diagonals_in_auto_scaling = true + + nl_rel_tol = 1e-8 + nl_abs_tol = 1e-12 + nl_max_its = 10 + + steady_state_detection = true + steady_state_tolerance = 1e-8 + dtmin = 1e-2 + dtmax = 10 + [TimeStepper] + type = IterationAdaptiveDT + dt = 1e-2 + cutback_factor = 0.9 + growth_factor = 1.1 + optimal_iterations = 6 + iteration_window = 1 + linear_iteration_ratio = 1000 + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Outputs] + [exodus] + type = Exodus + [] + [csv] + type = CSV + execute_on = final + [] +[] + +[Debug] + show_var_residual_norms = true +[] diff --git a/src/auxkernels/SATurbulentViscosityAux.C b/src/auxkernels/SATurbulentViscosityAux.C new file mode 100644 index 0000000000..b31830fb7e --- /dev/null +++ b/src/auxkernels/SATurbulentViscosityAux.C @@ -0,0 +1,34 @@ +#include "SATurbulentViscosityAux.h" +#include "metaphysicl/raw_type.h" + +registerMooseObject("MoltresApp", SATurbulentViscosityAux); + +InputParameters +SATurbulentViscosityAux::validParams() +{ + InputParameters params = AuxKernel::validParams(); + params.addClassDescription("Computes the turbulent viscosity for the Spalart-Allmaras " + "turbulence model."); + params.addRequiredCoupledVar("mu_tilde", + "Scaled turbulent viscosity variable from the Spalart-Allmaras " + "transport equation."); + params.addParam( + "mu_name", "mu", "The name of the viscosity material property"); + return params; +} + +SATurbulentViscosityAux::SATurbulentViscosityAux(const InputParameters & parameters) + : AuxKernel(parameters), + _mu_tilde(coupledValue("mu_tilde")), + _mu(getADMaterialProperty("mu_name")) +{ +} + +Real +SATurbulentViscosityAux::computeValue() +{ + Real chi = _mu_tilde[_qp] / MetaPhysicL::raw_value(_mu[_qp]); + Real cv1 = 7.1; + Real fv1 = std::pow(chi, 3) / (std::pow(chi, 3) + std::pow(cv1, 3)); + return _mu_tilde[_qp] * fv1; +} diff --git a/src/auxkernels/TurbulentStressAux.C b/src/auxkernels/TurbulentStressAux.C new file mode 100644 index 0000000000..482edbbf46 --- /dev/null +++ b/src/auxkernels/TurbulentStressAux.C @@ -0,0 +1,37 @@ +#include "TurbulentStressAux.h" + +using MetaPhysicL::raw_value; + +registerMooseObject("MoltresApp", TurbulentStressAux); + +InputParameters +TurbulentStressAux::validParams() +{ + InputParameters params = AuxKernel::validParams(); + params.addClassDescription("Computes the turbulent stress introduced by the Spalart-Allmaras " + "model."); + params.addRequiredCoupledVar("mu_tilde", "Spalart-Allmaras turbulence viscosity variable"); + params.addParam( + "mu_name", "mu", "The name of the viscosity material property"); + params.addParam( + "rho_name", "rho", "The name of the density material property"); + return params; +} + +TurbulentStressAux::TurbulentStressAux(const InputParameters & parameters) + : AuxKernel(parameters), + _mu_tilde(adCoupledValue("mu_tilde")), + _mu(getADMaterialProperty("mu_name")), + _rho(getADMaterialProperty("rho_name")), + _strain_mag(getADMaterialProperty("strain_mag")) +{ +} + +Real +TurbulentStressAux::computeValue() +{ + ADReal chi = _mu_tilde[_qp] / _mu[_qp]; + ADReal fv1 = std::pow(chi, 3) / (std::pow(chi, 3) + std::pow(7.1, 3)); + ADReal nu_T = _mu_tilde[_qp] * fv1 / _rho[_qp]; + return raw_value(nu_T * _strain_mag[_qp]); +} diff --git a/src/auxkernels/WallDistanceAux.C b/src/auxkernels/WallDistanceAux.C new file mode 100644 index 0000000000..c08212ebdc --- /dev/null +++ b/src/auxkernels/WallDistanceAux.C @@ -0,0 +1,54 @@ +#include "WallDistanceAux.h" + +registerMooseObject("MoltresApp", WallDistanceAux); + +InputParameters +WallDistanceAux::validParams() +{ + InputParameters params = AuxKernel::validParams(); + params.addClassDescription("Computes the minimum wall distance of each node from wall " + "boundaries"); + params.addRequiredParam>( + "walls", + "Boundaries that correspond to solid walls"); + return params; +} + +WallDistanceAux::WallDistanceAux(const InputParameters & parameters) + : AuxKernel(parameters), + _wall_boundary_names(getParam>("walls")) +{ + const MeshBase & mesh = _subproblem.mesh().getMesh(); + if (!mesh.is_replicated()) + mooseError("WallDistanceAux only supports replicated meshes"); + if (!isNodal()) + mooseError("WallDistanceAux only works on nodal wall distance variable fields " + "(e.g. LAGRANGE)."); +} + +Real +WallDistanceAux::computeValue() +{ + // Get the ids of the wall boundaries + std::vector vec_ids = _mesh.getBoundaryIDs(_wall_boundary_names, true); + + // Loop over all boundary nodes + Real min_dist = 1e9; + for (const auto & bnode : *_mesh.getBoundaryNodeRange()) + { + BoundaryID bnd_id = bnode->_bnd_id; + // Loop over all boundary ids + for (BoundaryID bid : vec_ids) + { + // Check if boundary node falls on wall boundary + if (bid == bnd_id) + { + // Find distance to closest wall boundary node + Node * node = bnode->_node; + const auto dist = (*node - *_current_node).norm(); + min_dist = std::min(min_dist, dist); + } + } + } + return min_dist; +} diff --git a/src/auxkernels/WallShearStressAux.C b/src/auxkernels/WallShearStressAux.C new file mode 100644 index 0000000000..046948785b --- /dev/null +++ b/src/auxkernels/WallShearStressAux.C @@ -0,0 +1,44 @@ +#include "WallShearStressAux.h" +#include "Assembly.h" + +using MetaPhysicL::raw_value; + +registerMooseObject("MoltresApp", WallShearStressAux); + +InputParameters +WallShearStressAux::validParams() +{ + InputParameters params = AuxKernel::validParams(); + params.addClassDescription("Computes the wall shear stress along the given boundary."); + params.addRequiredCoupledVar("velocity", "The velocity"); + params.addParam("mu_name", "mu", "The name of the dynamic viscosity"); + return params; +} + +WallShearStressAux::WallShearStressAux(const InputParameters & parameters) + : AuxKernel(parameters), + _grad_velocity(adCoupledVectorGradient("velocity")), + _mu(getADMaterialProperty("mu_name")), + _normals(_assembly.normals()) +{ + if (!isParamValid("boundary")) + paramError("boundary", + "The wall boundary must be provided to calculate the wall shear stress."); + if (_var.feType().family != MONOMIAL || _var.feType().order != CONSTANT) + mooseError("WallShearStressAux is compatible with CONSTANT MONOMIAL " + "wall shear stress variable only."); +} + +Real +WallShearStressAux::computeValue() +{ + RealVectorValue grad_vel_wall(0., 0., 0.); + // Velocity gradient along normal to wall (pointing in) + for (unsigned int i = 0; i < 3; i++) + for (unsigned int j = 0; j < 3; j++) + grad_vel_wall(i) -= raw_value(_grad_velocity[_qp](i, j)) * _normals[_qp](j); + // Parallel component of prior velocity gradient + grad_vel_wall -= (_normals[_qp] * grad_vel_wall) * _normals[_qp]; + + return raw_value(_mu[_qp]) * std::norm(grad_vel_wall); +} diff --git a/src/bcs/INSADMassPSPGBC.C b/src/bcs/INSADMassPSPGBC.C new file mode 100644 index 0000000000..8ef0b13353 --- /dev/null +++ b/src/bcs/INSADMassPSPGBC.C @@ -0,0 +1,26 @@ +#include "INSADMassPSPGBC.h" + +registerMooseObject("MoltresApp", INSADMassPSPGBC); + +InputParameters +INSADMassPSPGBC::validParams() +{ + InputParameters params = ADIntegratedBC::validParams(); + params.addClassDescription("This class removes INS PSPG contributions along the boundary."); + params.addParam("rho_name", "rho", "The name of the density"); + return params; +} + +INSADMassPSPGBC::INSADMassPSPGBC(const InputParameters & parameters) + : ADIntegratedBC(parameters), + _rho(getADMaterialProperty("rho_name")), + _tau(getADMaterialProperty("tau")), + _momentum_strong_residual(getADMaterialProperty("momentum_strong_residual")) +{ +} + +ADReal +INSADMassPSPGBC::computeQpResidual() +{ + return _tau[_qp] / _rho[_qp] * _momentum_strong_residual[_qp] * _test[_i][_qp] * _normals[_qp]; +} diff --git a/src/bcs/INSADMomentumSUPGBC.C b/src/bcs/INSADMomentumSUPGBC.C new file mode 100644 index 0000000000..ac27c6aad0 --- /dev/null +++ b/src/bcs/INSADMomentumSUPGBC.C @@ -0,0 +1,101 @@ +#include "INSADMomentumSUPGBC.h" +#include "Assembly.h" + +using MetaPhysicL::raw_value; + +registerMooseObject("MoltresApp", INSADMomentumSUPGBC); + +InputParameters +INSADMomentumSUPGBC::validParams() +{ + InputParameters params = ADVectorIntegratedBC::validParams(); + params.addClassDescription("This class removes INS SUPG contributions along the boundary."); + params.addParam( + "rho_name", "rho", "Density"); + params.addParam( + "tau_name", "tau", "The name of the stabilization parameter tau."); + return params; +} + +INSADMomentumSUPGBC::INSADMomentumSUPGBC(const InputParameters & parameters) + : ADVectorIntegratedBC(parameters), + _rho(getADMaterialProperty("rho_name")), + _tau(getADMaterialProperty("tau_name")), + _momentum_strong_residual(getADMaterialProperty("momentum_strong_residual")) +{ +} + +ADRealVectorValue +INSADMomentumSUPGBC::computeQpStabilization() +{ + return _u[_qp] * _tau[_qp]; +} + +ADRealVectorValue +INSADMomentumSUPGBC::precomputeQpStrongResidual() +{ + return -_momentum_strong_residual[_qp]; +} + +void +INSADMomentumSUPGBC::computeResidual() +{ + _residuals.resize(_test.size(), 0); + for (auto & r : _residuals) + r = 0; + + if (_use_displaced_mesh) + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + const auto value = precomputeQpStrongResidual() * _ad_JxW[_qp] * _ad_coord[_qp]; + for (_i = 0; _i < _test.size(); _i++) + _residuals[_i] += raw_value((_test[_i][_qp] * _normals[_qp]) * computeQpStabilization() * + value); + } + else + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + const auto value = precomputeQpStrongResidual() * _JxW[_qp] * _coord[_qp]; + for (_i = 0; _i < _test.size(); _i ++) + _residuals[_i] += + raw_value((_test[_i][_qp] * _normals[_qp]) * computeQpStabilization() * value); + } + + addResiduals(_assembly, _residuals, _var.dofIndices(), _var.scalingFactor()); + + if (_has_save_in) + for (unsigned int i = 0; i < _save_in.size(); i++) + _save_in[i]->sys().solution().add_vector(_residuals.data(), _save_in[i]->dofIndices()); +} + +void +INSADMomentumSUPGBC::computeResidualsForJacobian() +{ + if (_residuals_and_jacobians.size() != _test.size()) + _residuals_and_jacobians.resize(_test.size(), 0); + for (auto & r : _residuals_and_jacobians) + r = 0; + + if (_use_displaced_mesh) + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + const auto value = precomputeQpStrongResidual() * _ad_JxW[_qp] * _ad_coord[_qp]; + for (_i = 0; _i < _test.size(); _i++) + _residuals_and_jacobians[_i] += (_test[_i][_qp] * _normals[_qp]) * + computeQpStabilization() * value; + } + else + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + const auto value = precomputeQpStrongResidual() * _JxW[_qp] * _coord[_qp]; + for (_i = 0; _i < _test.size(); _i ++) + _residuals_and_jacobians[_i] += + (_test[_i][_qp] * _normals[_qp]) * computeQpStabilization() * value; + } +} + +ADReal +INSADMomentumSUPGBC::computeQpResidual() +{ + mooseError("INSADMomentumSUPGBC should not use computeQpResidual."); +} diff --git a/src/bcs/INSADMomentumTurbulentNoBCBC.C b/src/bcs/INSADMomentumTurbulentNoBCBC.C new file mode 100644 index 0000000000..92ca939375 --- /dev/null +++ b/src/bcs/INSADMomentumTurbulentNoBCBC.C @@ -0,0 +1,35 @@ +#include "INSADMomentumTurbulentNoBCBC.h" + +registerMooseObject("MoltresApp", INSADMomentumTurbulentNoBCBC); + +InputParameters +INSADMomentumTurbulentNoBCBC::validParams() +{ + InputParameters params = INSADMomentumNoBCBC::validParams(); + params.addClassDescription("This class implements the 'No BC' boundary condition with " + "additional contributions from turbulent viscosity."); + params.addRequiredCoupledVar("mu_tilde", "Spalart-Allmaras turbulence viscosity variable"); + return params; +} + +INSADMomentumTurbulentNoBCBC::INSADMomentumTurbulentNoBCBC(const InputParameters & parameters) + : INSADMomentumNoBCBC(parameters), + _mu_tilde(adCoupledValue("mu_tilde")) +{ + if (_form == "laplace") + mooseError("The Laplace form of the INS equations is incompatible with the Spalart-Allmaras " + "turbulent viscosity variable."); +} + +ADReal +INSADMomentumTurbulentNoBCBC::computeQpResidual() +{ + ADReal chi = _mu_tilde[_qp] / _mu[_qp]; + ADReal fv1 = std::pow(chi, 3) / (std::pow(chi, 3) + std::pow(7.1, 3)); + ADReal residual; + residual = -(_mu[_qp] + _mu_tilde[_qp] * fv1) * ((_grad_u[_qp] + _grad_u[_qp].transpose()) * \ + _normals[_qp]) * _test[_i][_qp]; + if (_integrate_p_by_parts) + residual += _p[_qp] * _normals[_qp] * _test[_i][_qp]; + return residual; +} diff --git a/src/bcs/INSOutflowBC.C b/src/bcs/INSOutflowBC.C deleted file mode 100644 index 6667319030..0000000000 --- a/src/bcs/INSOutflowBC.C +++ /dev/null @@ -1,115 +0,0 @@ -#include "INSOutflowBC.h" - -registerMooseObject("MoltresApp", INSOutflowBC); - -InputParameters -INSOutflowBC::validParams() -{ - InputParameters params = IntegratedBC::validParams(); - - // Coupled variables - params.addRequiredCoupledVar("u", "x-velocity"); - params.addCoupledVar("v", 0, "y-velocity"); // only required in 2D and 3D - params.addCoupledVar("w", 0, "z-velocity"); // only required in 3D - params.addRequiredCoupledVar("p", "pressure"); - - // Required parameters - params.addRequiredParam("mu", "dynamic viscosity"); - params.addRequiredParam("rho", "density"); - params.addRequiredParam( - "component", - "0,1,2 depending on if we are solving the x,y,z component of the momentum equation"); - params.addParam("integrate_p_by_parts", - true, - "Allows simulations to be run with pressure BC if set to false"); - - return params; -} - -INSOutflowBC::INSOutflowBC(const InputParameters & parameters) - : IntegratedBC(parameters), - - // Coupled variables - _u_vel(coupledValue("u")), - _v_vel(coupledValue("v")), - _w_vel(coupledValue("w")), - _p(coupledValue("p")), - - // Gradients - _grad_u_vel(coupledGradient("u")), - _grad_v_vel(coupledGradient("v")), - _grad_w_vel(coupledGradient("w")), - _grad_p(coupledGradient("p")), - - // Variable numberings - _u_vel_var_number(coupled("u")), - _v_vel_var_number(coupled("v")), - _w_vel_var_number(coupled("w")), - _p_var_number(coupled("p")), - - // Required parameters - _mu(getParam("mu")), - _rho(getParam("rho")), - _component(getParam("component")), - _integrate_p_by_parts(getParam("integrate_p_by_parts")) - -// Material properties -// _dynamic_viscosity(getMaterialProperty("dynamic_viscosity")) -{ -} - -Real -INSOutflowBC::computeQpResidual() -{ - // The pressure part, -p (div v) or (dp/dx_{component}) * test if not integrated by parts. - Real pressure_part = 0.; - if (_integrate_p_by_parts) - pressure_part = _p[_qp] * _normals[_qp](_component) * _test[_i][_qp]; - - // The component'th row (or col, it's symmetric) of the viscous stress tensor - RealVectorValue tau_row; - - // Apply BC that the gradient of each velocity component in the normal direction is zero - switch (_component) - { - case 0: - tau_row(0) = _grad_u_vel[_qp](0); - tau_row(1) = _grad_v_vel[_qp](0); - tau_row(2) = _grad_w_vel[_qp](0); - break; - - case 1: - tau_row(0) = _grad_u_vel[_qp](1); - tau_row(1) = _grad_v_vel[_qp](1); - tau_row(2) = _grad_w_vel[_qp](1); - break; - - case 2: - tau_row(0) = _grad_u_vel[_qp](2); - tau_row(1) = _grad_v_vel[_qp](2); - tau_row(2) = _grad_w_vel[_qp](2); - break; - - default: - mooseError("Unrecognized _component requested."); - } - - // The viscous part, tau : grad(v) - Real viscous_part = -_mu * (tau_row * _normals[_qp] * _test[_i][_qp]); - - return pressure_part + viscous_part; -} - -// Implement Jacobian - -Real -INSOutflowBC::computeQpJacobian() -{ - return 0.; -} - -Real -INSOutflowBC::computeQpOffDiagJacobian(unsigned /*jvar*/) -{ - return 0.; -} diff --git a/src/bcs/INSSymmetryAxisBC.C b/src/bcs/INSSymmetryAxisBC.C deleted file mode 100644 index a5ef52cb4f..0000000000 --- a/src/bcs/INSSymmetryAxisBC.C +++ /dev/null @@ -1,110 +0,0 @@ -#include "INSSymmetryAxisBC.h" - -registerMooseObject("MoltresApp", INSSymmetryAxisBC); - -InputParameters -INSSymmetryAxisBC::validParams() -{ - InputParameters params = IntegratedBC::validParams(); - - // Coupled variables - params.addRequiredCoupledVar("u", "x-velocity"); - params.addCoupledVar("v", 0, "y-velocity"); // only required in 2D and 3D - params.addCoupledVar("w", 0, "z-velocity"); // only required in 3D - params.addRequiredCoupledVar("p", "pressure"); - - // Required parameters - params.addRequiredParam("mu", "dynamic viscosity"); - params.addRequiredParam("rho", "density"); - params.addRequiredParam( - "component", - "0,1,2 depending on if we are solving the x,y,z component of the momentum equation"); - params.addParam("integrate_p_by_parts", - true, - "Allows simulations to be run with pressure BC if set to false"); - - return params; -} - -INSSymmetryAxisBC::INSSymmetryAxisBC(const InputParameters & parameters) - : IntegratedBC(parameters), - - // Coupled variables - _u_vel(coupledValue("u")), - _v_vel(coupledValue("v")), - _w_vel(coupledValue("w")), - _p(coupledValue("p")), - - // Gradients - _grad_u_vel(coupledGradient("u")), - _grad_v_vel(coupledGradient("v")), - _grad_w_vel(coupledGradient("w")), - _grad_p(coupledGradient("p")), - - // Variable numberings - _u_vel_var_number(coupled("u")), - _v_vel_var_number(coupled("v")), - _w_vel_var_number(coupled("w")), - _p_var_number(coupled("p")), - - // Required parameters - _mu(getParam("mu")), - _rho(getParam("rho")), - _component(getParam("component")), - _integrate_p_by_parts(getParam("integrate_p_by_parts")) - -// Material properties -// _dynamic_viscosity(getMaterialProperty("dynamic_viscosity")) -{ -} - -Real -INSSymmetryAxisBC::computeQpResidual() -{ - // The pressure part, -p (div v) or (dp/dx_{component}) * test if not integrated by parts. - Real pressure_part = 0.; - if (_integrate_p_by_parts) - pressure_part = _p[_qp] * _normals[_qp](_component) * _test[_i][_qp]; - - // The component'th row (or col, it's symmetric) of the viscous stress tensor - RealVectorValue tau_row; - - // Apply BC that the derivative of the z-component of the velocity in the radial direction is - // zero. - switch (_component) - { - case 0: - mooseError("This BC is only applicable to the z-component of the velocity"); - - case 1: - tau_row(0) = _grad_u_vel[_qp](1); // dv/dx1 + du/dx2 = 0 + du/dx2 - tau_row(1) = 2. * _grad_v_vel[_qp](1); // 2*dv/dx2 - tau_row(2) = _grad_v_vel[_qp](2) + _grad_w_vel[_qp](1); // dv/dx3 + dw/dx2 - break; - - case 2: - mooseError("This BC is only applicable to the z-component of the velocity"); - - default: - mooseError("Unrecognized _component requested."); - } - - // The viscous part, tau : grad(v) - Real viscous_part = -_mu * (tau_row * _normals[_qp] * _test[_i][_qp]); - - return pressure_part + viscous_part; -} - -// Implement Jacobian - -Real -INSSymmetryAxisBC::computeQpJacobian() -{ - return 0.; -} - -Real -INSSymmetryAxisBC::computeQpOffDiagJacobian(unsigned /*jvar*/) -{ - return 0.; -} diff --git a/src/bcs/SATurbulentViscosityNoBCBC.C b/src/bcs/SATurbulentViscosityNoBCBC.C new file mode 100644 index 0000000000..9c8b35fa10 --- /dev/null +++ b/src/bcs/SATurbulentViscosityNoBCBC.C @@ -0,0 +1,26 @@ +#include "SATurbulentViscosityNoBCBC.h" + +registerMooseObject("MoltresApp", SATurbulentViscosityNoBCBC); + +InputParameters +SATurbulentViscosityNoBCBC::validParams() +{ + InputParameters params = ADIntegratedBC::validParams(); + params.addClassDescription("This class implements the 'No BC' boundary condition for the " + "Spalart-Allmaras equation"); + params.addParam("mu_name", "mu", "The name of the dynamic viscosity"); + return params; +} + +SATurbulentViscosityNoBCBC::SATurbulentViscosityNoBCBC(const InputParameters & parameters) + : ADIntegratedBC(parameters), + _mu(getADMaterialProperty("mu_name")) +{ +} + +ADReal +SATurbulentViscosityNoBCBC::computeQpResidual() +{ + Real sigma = 2.0 / 3.0; + return (-1.0 / sigma) * (_mu[_qp] + _u[_qp]) * (_grad_u[_qp] * _normals[_qp]) * _test[_i][_qp]; +} diff --git a/src/bcs/SATurbulentViscositySUPGBC.C b/src/bcs/SATurbulentViscositySUPGBC.C new file mode 100644 index 0000000000..d451d98447 --- /dev/null +++ b/src/bcs/SATurbulentViscositySUPGBC.C @@ -0,0 +1,98 @@ +#include "SATurbulentViscositySUPGBC.h" + +registerMooseObject("MoltresApp", SATurbulentViscositySUPGBC); + +InputParameters +SATurbulentViscositySUPGBC::validParams() +{ + InputParameters params = ADIntegratedBC::validParams(); + params.addClassDescription("This class removes SUPG contributions in the Spalart-Allmaras" + " equation along the boundary."); + params.addParam("tau_name", "tau_viscosity", + "The name of the dynamic viscosity material property"); + params.addRequiredCoupledVar("velocity", "The velocity variable."); + return params; +} + +SATurbulentViscositySUPGBC::SATurbulentViscositySUPGBC(const InputParameters & parameters) + : ADIntegratedBC(parameters), + _tau_visc(getADMaterialProperty("tau_name")), + _velocity(adCoupledVectorValue("velocity")), + _visc_strong_residual(getADMaterialProperty("visc_strong_residual")) +{ +} + +ADRealVectorValue +SATurbulentViscositySUPGBC::computeQpStabilization() +{ + return _velocity[_qp] * _tau_visc[_qp]; +} + +ADReal +SATurbulentViscositySUPGBC::precomputeQpStrongResidual() +{ + return -_visc_strong_residual[_qp]; +} + +void +SATurbulentViscositySUPGBC::computeResidual() +{ + _residuals.resize(_test.size(), 0); + for (auto & r : _residuals) + r = 0; + + if (_use_displaced_mesh) + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + const auto value = precomputeQpStrongResidual() * _ad_JxW[_qp] * _ad_coord[_qp]; + for (_i = 0; _i < _test.size(); _i++) + _residuals[_i] += raw_value(_test[_i][_qp] * _normals[_qp] * computeQpStabilization() * + value); + } + else + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + const auto value = precomputeQpStrongResidual() * _JxW[_qp] * _coord[_qp]; + for (_i = 0; _i < _test.size(); _i ++) + _residuals[_i] += + raw_value((_test[_i][_qp] * _normals[_qp]) * computeQpStabilization() * value); + } + + addResiduals(_assembly, _residuals, _var.dofIndices(), _var.scalingFactor()); + + if (_has_save_in) + for (unsigned int i = 0; i < _save_in.size(); i++) + _save_in[i]->sys().solution().add_vector(_residuals.data(), _save_in[i]->dofIndices()); +} + +void +SATurbulentViscositySUPGBC::computeResidualsForJacobian() +{ + if (_residuals_and_jacobians.size() != _test.size()) + _residuals_and_jacobians.resize(_test.size(), 0); + for (auto & r : _residuals_and_jacobians) + r = 0; + + if (_use_displaced_mesh) + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + const auto value = precomputeQpStrongResidual() * _ad_JxW[_qp] * _ad_coord[_qp]; + for (_i = 0; _i < _test.size(); _i++) + _residuals_and_jacobians[_i] += _test[_i][_qp] * _normals[_qp] * computeQpStabilization() * + value; + } + else + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + const auto value = precomputeQpStrongResidual() * _JxW[_qp] * _coord[_qp]; + for (_i = 0; _i < _test.size(); _i ++) + _residuals_and_jacobians[_i] += + (_test[_i][_qp] * _normals[_qp]) * computeQpStabilization() * value; + } +} + +ADReal +SATurbulentViscositySUPGBC::computeQpResidual() +{ + mooseError("SATurbulentViscositySUPGBC should not use computeQpResidual."); +} diff --git a/src/dgkernels/DGTurbulentDiffusion.C b/src/dgkernels/DGTurbulentDiffusion.C new file mode 100644 index 0000000000..f1da585277 --- /dev/null +++ b/src/dgkernels/DGTurbulentDiffusion.C @@ -0,0 +1,192 @@ +#include "DGTurbulentDiffusion.h" +#include "metaphysicl/raw_type.h" + +using MetaPhysicL::raw_value; + +registerMooseObject("MoltresApp", DGTurbulentDiffusion); + +InputParameters +DGTurbulentDiffusion::validParams() +{ + InputParameters params = DGKernel::validParams(); + params.addClassDescription( + "Adds the turbulent diffusion term for delayed " + "neutron precursors using the discontinuous Galerkin method."); + params.addParam("sigma", 6, "sigma"); + params.addParam("epsilon", -1, "epsilon"); + params.addParam( + "mu_name", "mu", "The name of the dynamic viscosity material proerty"); + params.addParam( + "rho_name", "rho", "The name of the density material property"); + params.addRequiredCoupledVar("mu_tilde", "Spalart-Allmaras turbulence viscosity variable"); + params.addParam("schmidt_number", 0.85, "Turbulent Schmidt number"); + return params; +} + +DGTurbulentDiffusion::DGTurbulentDiffusion(const InputParameters & parameters) + : DGKernel(parameters), + _epsilon(getParam("epsilon")), + _sigma(getParam("sigma")), + _mu(getADMaterialProperty("mu_name")), + _mu_nb(getNeighborADMaterialProperty("mu_name")), + _rho(getADMaterialProperty("rho_name")), + _rho_nb(getNeighborADMaterialProperty("rho_name")), + _mu_tilde(adCoupledValue("mu_tilde")), + _mu_tilde_nb(adCoupledNeighborValue("mu_tilde")), + _sc(getParam("schmidt_number")), + _mu_tilde_var_number(coupled("mu_tilde")) +{ +} + +Real +DGTurbulentDiffusion::computeQpResidual(Moose::DGResidualType type) +{ + Real chi = raw_value(_mu_tilde[_qp] / _mu[_qp]); + Real cv1 = 7.1; + Real fv1 = Utility::pow<3>(chi) / (Utility::pow<3>(chi) + Utility::pow<3>(cv1)); + // Turbulent diffusion coefficient + Real D = raw_value(_mu_tilde[_qp] / _rho[_qp]) * fv1 / _sc; + + Real chi_nb = raw_value(_mu_tilde_nb[_qp]) / raw_value(_mu_nb[_qp]); + Real cv1_nb = 7.1; + Real fv1_nb = Utility::pow<3>(chi_nb) / (Utility::pow<3>(chi_nb) + Utility::pow<3>(cv1_nb)); + // Neighbor turbulent diffusion coefficient + Real D_nb = raw_value(_mu_tilde_nb[_qp] / _rho_nb[_qp]) * fv1_nb / _sc; + + Real r = 0.0; + + const int elem_b_order = std::max(libMesh::Order(1), _var.order()); + const Real h_elem = + _current_elem_volume / _current_side_volume * 1.0 / Utility::pow<2>(elem_b_order); + + switch (type) + { + case Moose::Element: + r -= 0.5 * + (D * _grad_u[_qp] * _normals[_qp] + + D_nb * _grad_u_neighbor[_qp] * _normals[_qp]) * + _test[_i][_qp]; + r += _epsilon * 0.5 * (_u[_qp] - _u_neighbor[_qp]) * D * _grad_test[_i][_qp] * + _normals[_qp]; + r += _sigma / h_elem * (_u[_qp] - _u_neighbor[_qp]) * _test[_i][_qp]; + break; + + case Moose::Neighbor: + r += 0.5 * + (D * _grad_u[_qp] * _normals[_qp] + + D_nb * _grad_u_neighbor[_qp] * _normals[_qp]) * + _test_neighbor[_i][_qp]; + r += _epsilon * 0.5 * (_u[_qp] - _u_neighbor[_qp]) * D_nb * + _grad_test_neighbor[_i][_qp] * _normals[_qp]; + r -= _sigma / h_elem * (_u[_qp] - _u_neighbor[_qp]) * _test_neighbor[_i][_qp]; + break; + } + + return r; +} + +Real +DGTurbulentDiffusion::computeQpJacobian(Moose::DGJacobianType type) +{ + Real chi = raw_value(_mu_tilde[_qp] / _mu[_qp]); + Real cv1 = 7.1; + Real fv1 = Utility::pow<3>(chi) / (Utility::pow<3>(chi) + Utility::pow<3>(cv1)); + // Turbulent diffusion coefficient + Real D = raw_value(_mu_tilde[_qp] / _rho[_qp]) * fv1 / _sc; + + Real chi_nb = raw_value(_mu_tilde_nb[_qp]) / raw_value(_mu_nb[_qp]); + Real cv1_nb = 7.1; + Real fv1_nb = Utility::pow<3>(chi_nb) / (Utility::pow<3>(chi_nb) + Utility::pow<3>(cv1_nb)); + // Neighbor turbulent diffusion coefficient + Real D_nb = raw_value(_mu_tilde_nb[_qp] / _rho_nb[_qp]) * fv1_nb / _sc; + + Real r = 0.0; + + const int elem_b_order = std::max(libMesh::Order(1), _var.order()); + const Real h_elem = + _current_elem_volume / _current_side_volume * 1.0 / Utility::pow<2>(elem_b_order); + + switch (type) + { + case Moose::ElementElement: + r -= 0.5 * D * _grad_phi[_j][_qp] * _normals[_qp] * _test[_i][_qp]; + r += _epsilon * 0.5 * _phi[_j][_qp] * D * _grad_test[_i][_qp] * _normals[_qp]; + r += _sigma / h_elem * _phi[_j][_qp] * _test[_i][_qp]; + break; + + case Moose::ElementNeighbor: + r -= 0.5 * D_nb * _grad_phi_neighbor[_j][_qp] * _normals[_qp] * _test[_i][_qp]; + r += _epsilon * 0.5 * -_phi_neighbor[_j][_qp] * D * _grad_test[_i][_qp] * + _normals[_qp]; + r += _sigma / h_elem * -_phi_neighbor[_j][_qp] * _test[_i][_qp]; + break; + + case Moose::NeighborElement: + r += 0.5 * D * _grad_phi[_j][_qp] * _normals[_qp] * _test_neighbor[_i][_qp]; + r += _epsilon * 0.5 * _phi[_j][_qp] * D_nb * _grad_test_neighbor[_i][_qp] * + _normals[_qp]; + r -= _sigma / h_elem * _phi[_j][_qp] * _test_neighbor[_i][_qp]; + break; + + case Moose::NeighborNeighbor: + r += 0.5 * D_nb * _grad_phi_neighbor[_j][_qp] * _normals[_qp] * + _test_neighbor[_i][_qp]; + r += _epsilon * 0.5 * -_phi_neighbor[_j][_qp] * D_nb * + _grad_test_neighbor[_i][_qp] * _normals[_qp]; + r -= _sigma / h_elem * -_phi_neighbor[_j][_qp] * _test_neighbor[_i][_qp]; + break; + } + + return r; +} + +Real +DGTurbulentDiffusion::computeQpOffDiagJacobian(Moose::DGJacobianType type, unsigned jvar) +{ + Real r = 0.0; + + if (jvar == _mu_tilde_var_number) + { + Real chi = raw_value(_mu_tilde[_qp] / _mu[_qp]); + Real cv1 = 7.1; + Real fv1 = Utility::pow<3>(chi) / (Utility::pow<3>(chi) + Utility::pow<3>(cv1)); + // Turbulent diffusion coefficient + Real D = _phi[_j][_qp] * fv1 / raw_value(_rho[_qp]) / _sc * + (1 + 3 * Utility::pow<3>(cv1) / (Utility::pow<3>(chi) + Utility::pow<3>(cv1))); + + Real chi_nb = raw_value(_mu_tilde_nb[_qp] / _mu_nb[_qp]); + Real fv1_nb = Utility::pow<3>(chi_nb) / (Utility::pow<3>(chi_nb) + Utility::pow<3>(cv1)); + // Neighbor turbulent diffusion coefficient + Real D_nb = _phi_neighbor[_j][_qp] * fv1_nb / raw_value(_rho_nb[_qp]) / _sc * + (1 + 3 * Utility::pow<3>(cv1) / (Utility::pow<3>(chi_nb) + Utility::pow<3>(cv1))); + + switch (type) + { + case Moose::ElementElement: + r -= 0.5 * D * _grad_u[_qp] * _normals[_qp] * _test[_i][_qp]; + r += _epsilon * 0.5 * _u[_qp] * D * _grad_test[_i][_qp] * _normals[_qp]; + break; + + case Moose::ElementNeighbor: + r -= 0.5 * D_nb * _grad_u_neighbor[_qp] * _normals[_qp] * _test[_i][_qp]; + r += _epsilon * 0.5 * -_u_neighbor[_qp] * D * _grad_test[_i][_qp] * + _normals[_qp]; + break; + + case Moose::NeighborElement: + r += 0.5 * D * _grad_u[_qp] * _normals[_qp] * _test_neighbor[_i][_qp]; + r += _epsilon * 0.5 * _u[_qp] * D_nb * _grad_test_neighbor[_i][_qp] * + _normals[_qp]; + break; + + case Moose::NeighborNeighbor: + r += 0.5 * D_nb * _grad_u_neighbor[_qp] * _normals[_qp] * + _test_neighbor[_i][_qp]; + r += _epsilon * 0.5 * -_u_neighbor[_qp] * D_nb * + _grad_test_neighbor[_i][_qp] * _normals[_qp]; + break; + } + } + + return r; +} diff --git a/src/kernels/INSADEnergyTurbulentDiffusion.C b/src/kernels/INSADEnergyTurbulentDiffusion.C new file mode 100644 index 0000000000..19fddbbdee --- /dev/null +++ b/src/kernels/INSADEnergyTurbulentDiffusion.C @@ -0,0 +1,39 @@ +#include "INSADEnergyTurbulentDiffusion.h" + +#include "metaphysicl/raw_type.h" + +using MetaPhysicL::raw_value; + +registerMooseObject("MoltresApp", INSADEnergyTurbulentDiffusion); + +InputParameters +INSADEnergyTurbulentDiffusion::validParams() +{ + InputParameters params = ADDiffusion::validParams(); + params.addClassDescription("Adds the turbulent diffusion term in the INSAD heat equation"); + params.addParam("cp_name", "cp", "The name of the specific heat capacity"); + params.addParam("mu_name", "mu", + "The name of the dynamic viscosity material property"); + params.addRequiredCoupledVar("mu_tilde", "Spalart-Allmaras turbulence viscosity variable"); + params.addParam("prandtl_number", 0.85, "Prandtl number"); + return params; +} + +INSADEnergyTurbulentDiffusion::INSADEnergyTurbulentDiffusion(const InputParameters & parameters) + : ADDiffusion(parameters), + _cp(getADMaterialProperty("cp_name")), + _mu(getADMaterialProperty("mu_name")), + _mu_tilde(adCoupledValue("mu_tilde")), + _pr(getParam("prandtl_number")) +{ +} + +ADRealVectorValue +INSADEnergyTurbulentDiffusion::precomputeQpResidual() +{ + ADReal chi = _mu_tilde[_qp] / _mu[_qp]; + Real cv1 = 7.1; + ADReal fv1 = Utility::pow<3>(chi) / (Utility::pow<3>(chi) + Utility::pow<3>(cv1)); + + return _cp[_qp] * _mu_tilde[_qp] * fv1 / _pr * ADDiffusion::precomputeQpResidual(); +} diff --git a/src/kernels/INSADMomentumTurbulentViscous.C b/src/kernels/INSADMomentumTurbulentViscous.C new file mode 100644 index 0000000000..6626b17103 --- /dev/null +++ b/src/kernels/INSADMomentumTurbulentViscous.C @@ -0,0 +1,57 @@ +#include "INSADMomentumTurbulentViscous.h" + +#include "metaphysicl/raw_type.h" + +using MetaPhysicL::raw_value; + +registerMooseObject("MoltresApp", INSADMomentumTurbulentViscous); + +InputParameters +INSADMomentumTurbulentViscous::validParams() +{ + InputParameters params = INSADMomentumViscous::validParams(); + params.addClassDescription("Adds the turbulent viscous term to the INS momentum equation"); + params.addRequiredCoupledVar("mu_tilde", "Spalart-Allmaras turbulence viscosity variable"); + return params; +} + +INSADMomentumTurbulentViscous::INSADMomentumTurbulentViscous(const InputParameters & parameters) + : INSADMomentumViscous(parameters), + _mu_tilde(adCoupledValue("mu_tilde")) +{ + if (_form == "laplace") + mooseError("The Laplace form of the INS equations is incompatible with the Spalart-Allmaras " + "turbulent viscosity variable."); +} + +ADRealTensorValue +INSADMomentumTurbulentViscous::qpViscousTerm() +{ + ADReal chi = _mu_tilde[_qp] / _mu[_qp]; + Real cv1 = 7.1; + ADReal fv1 = Utility::pow<3>(chi) / (Utility::pow<3>(chi) + Utility::pow<3>(cv1)); + + return _mu_tilde[_qp] * fv1 * (_grad_u[_qp] + _grad_u[_qp].transpose()); +} + +ADRealVectorValue +INSADMomentumTurbulentViscous::qpAdditionalRZTerm() +{ + ADReal chi = _mu_tilde[_qp] / _mu[_qp]; + Real cv1 = 7.1; + ADReal fv1 = Utility::pow<3>(chi) / (Utility::pow<3>(chi) + Utility::pow<3>(cv1)); + + // Add the u_r / r^2 term. There is an extra factor of 2 for the traction form + ADReal resid = _mu_tilde[_qp] * fv1 * _u[_qp](_rz_radial_coord) * 2.0; + + if (_use_displaced_mesh) + return resid / (_ad_q_point[_qp](_rz_radial_coord) * _ad_q_point[_qp](_rz_radial_coord)); + else + return resid / (_q_point[_qp](_rz_radial_coord) * _q_point[_qp](_rz_radial_coord)); +} + +ADReal +INSADMomentumTurbulentViscous::computeQpResidual() +{ + mooseError("computeQpResidual is not used in the INSADMomentumTurbulentViscous class"); +} diff --git a/src/kernels/INSMomentumKEpsilon.C b/src/kernels/INSMomentumKEpsilon.C deleted file mode 100644 index f229a90d08..0000000000 --- a/src/kernels/INSMomentumKEpsilon.C +++ /dev/null @@ -1,188 +0,0 @@ -#include "INSMomentumKEpsilon.h" - -registerMooseObject("MoltresApp", INSMomentumKEpsilon); - -InputParameters -INSMomentumKEpsilon::validParams() -{ - InputParameters params = Kernel::validParams(); - - // Coupled variables - params.addRequiredCoupledVar("u", "x-velocity"); - params.addCoupledVar("v", 0, "y-velocity"); // only required in 2D and 3D - params.addCoupledVar("w", 0, "z-velocity"); // only required in 3D - params.addRequiredCoupledVar("p", "pressure"); - params.addRequiredCoupledVar("k", "turbulent kinetic energy"); - - // Required parameters - params.addRequiredParam("mu", "dynamic viscosity"); - params.addRequiredParam("rho", "density"); - params.addRequiredParam("gravity", "Direction of the gravity vector"); - params.addRequiredParam( - "component", - "0,1,2 depending on if we are solving the x,y,z component of the momentum equation"); - params.addParam("integrate_p_by_parts", - true, - "Allows simulations to be run with pressure BC if set to false"); - - return params; -} - -INSMomentumKEpsilon::INSMomentumKEpsilon(const InputParameters & parameters) - : Kernel(parameters), - - // Coupled variables - _u_vel(coupledValue("u")), - _v_vel(coupledValue("v")), - _w_vel(coupledValue("w")), - _p(coupledValue("p")), - _k(coupledValue("k")), - - // Gradients - _grad_u_vel(coupledGradient("u")), - _grad_v_vel(coupledGradient("v")), - _grad_w_vel(coupledGradient("w")), - _grad_p(coupledGradient("p")), - _grad_k(coupledGradient("k")), - - // Variable numberings - _u_vel_var_number(coupled("u")), - _v_vel_var_number(coupled("v")), - _w_vel_var_number(coupled("w")), - _p_var_number(coupled("p")), - _k_var_number(coupled("k")), - - // Required parameters - _mu(getParam("mu")), - _rho(getParam("rho")), - _gravity(getParam("gravity")), - _component(getParam("component")), - _integrate_p_by_parts(getParam("integrate_p_by_parts")) - -// Material properties -// _dynamic_viscosity(getMaterialProperty("dynamic_viscosity")) -{ -} - -Real -INSMomentumKEpsilon::computeQpResidual() -{ - // The convection part, rho * (u.grad) * u_component * v. - // Note: _grad_u is the gradient of the _component entry of the velocity vector. - Real convective_part = _rho * (_u_vel[_qp] * _grad_u[_qp](0) + _v_vel[_qp] * _grad_u[_qp](1) + - _w_vel[_qp] * _grad_u[_qp](2)) * - _test[_i][_qp]; - - // The pressure part, -p (div v) or (dp/dx_{component}) * test if not integrated by parts. - Real pressure_part = 0.; - if (_integrate_p_by_parts) - pressure_part = -_p[_qp] * _grad_test[_i][_qp](_component); - else - pressure_part = _grad_p[_qp](_component) * _test[_i][_qp]; - - // The turbulent energy part - Real turbulent_part = _test[_i][_qp] * 2. / 3. * _rho * _grad_k[_qp](_component); - - // The component'th row (or col, it's symmetric) of the viscous stress tensor - RealVectorValue tau_row; - - switch (_component) - { - case 0: - tau_row(0) = 2. * _grad_u_vel[_qp](0); // 2*du/dx1 - tau_row(1) = _grad_u_vel[_qp](1) + _grad_v_vel[_qp](0); // du/dx2 + dv/dx1 - tau_row(2) = _grad_u_vel[_qp](2) + _grad_w_vel[_qp](0); // du/dx3 + dw/dx1 - break; - - case 1: - tau_row(0) = _grad_v_vel[_qp](0) + _grad_u_vel[_qp](1); // dv/dx1 + du/dx2 - tau_row(1) = 2. * _grad_v_vel[_qp](1); // 2*dv/dx2 - tau_row(2) = _grad_v_vel[_qp](2) + _grad_w_vel[_qp](1); // dv/dx3 + dw/dx2 - break; - - case 2: - tau_row(0) = _grad_w_vel[_qp](0) + _grad_u_vel[_qp](2); // dw/dx1 + du/dx3 - tau_row(1) = _grad_w_vel[_qp](1) + _grad_v_vel[_qp](2); // dw/dx2 + dv/dx3 - tau_row(2) = 2. * _grad_w_vel[_qp](2); // 2*dw/dx3 - break; - - default: - mooseError("Unrecognized _component requested."); - } - - // The viscous part, tau : grad(v) - Real viscous_part = _mu * (tau_row * _grad_test[_i][_qp]); - - // Simplified version: mu * Laplacian(u_component) - // Real viscous_part = _mu * (_grad_u[_qp] * _grad_test[_i][_qp]); - - // Body force term. For truly incompressible flow, this term is constant, and - // since it is proportional to g, can be written as the gradient of some scalar - // and absorbed into the pressure definition. - // Real body_force_part = - _rho * _gravity(_component); - - return convective_part + pressure_part + turbulent_part + viscous_part /*+ body_force_part*/; -} - -Real -INSMomentumKEpsilon::computeQpJacobian() -{ - RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]); - - // Convective part - Real convective_part = - _rho * ((U * _grad_phi[_j][_qp]) + _phi[_j][_qp] * _grad_u[_qp](_component)) * _test[_i][_qp]; - - // Viscous part, Stokes/Laplacian version - // Real viscous_part = _mu * (_grad_phi[_j][_qp] * _grad_test[_i][_qp]); - - // Viscous part, full stress tensor. The extra contribution comes from the "2" - // on the diagonal of the viscous stress tensor. - Real viscous_part = _mu * (_grad_phi[_j][_qp] * _grad_test[_i][_qp] + - _grad_phi[_j][_qp](_component) * _grad_test[_i][_qp](_component)); - - return convective_part + viscous_part; -} - -Real -INSMomentumKEpsilon::computeQpOffDiagJacobian(unsigned jvar) -{ - // In Stokes/Laplacian version, off-diag Jacobian entries wrt u,v,w are zero - if (jvar == _u_vel_var_number) - { - Real convective_part = _phi[_j][_qp] * _grad_u[_qp](0) * _test[_i][_qp]; - Real viscous_part = _mu * _grad_phi[_j][_qp](_component) * _grad_test[_i][_qp](0); - - return convective_part + viscous_part; - } - - else if (jvar == _v_vel_var_number) - { - Real convective_part = _phi[_j][_qp] * _grad_u[_qp](1) * _test[_i][_qp]; - Real viscous_part = _mu * _grad_phi[_j][_qp](_component) * _grad_test[_i][_qp](1); - - return convective_part + viscous_part; - } - - else if (jvar == _w_vel_var_number) - { - Real convective_part = _phi[_j][_qp] * _grad_u[_qp](2) * _test[_i][_qp]; - Real viscous_part = _mu * _grad_phi[_j][_qp](_component) * _grad_test[_i][_qp](2); - - return convective_part + viscous_part; - } - - else if (jvar == _p_var_number) - { - if (_integrate_p_by_parts) - return -_phi[_j][_qp] * _grad_test[_i][_qp](_component); - else - return _grad_phi[_j][_qp](_component) * _test[_i][_qp]; - } - - else if (jvar == _k_var_number) - return _test[_i][_qp] * 2. / 3. * _rho * _grad_phi[_j][_qp](_component); - - else - return 0; -} diff --git a/src/kernels/SATimeDerivative.C b/src/kernels/SATimeDerivative.C new file mode 100644 index 0000000000..72b881bd6f --- /dev/null +++ b/src/kernels/SATimeDerivative.C @@ -0,0 +1,24 @@ +#include "SATimeDerivative.h" + +registerMooseObject("MoltresApp", SATimeDerivative); + +InputParameters +SATimeDerivative::validParams() +{ + InputParameters params = ADTimeDerivative::validParams(); + params.addClassDescription("Adds the time derivative term of the Spalart-Allmaras " + "turbulence model."); + return params; +} + +SATimeDerivative::SATimeDerivative(const InputParameters & parameters) + : ADTimeDerivative(parameters), + _time_strong_residual(getADMaterialProperty("time_strong_residual")) +{ +} + +ADReal +SATimeDerivative::precomputeQpResidual() +{ + return _time_strong_residual[_qp]; +} diff --git a/src/kernels/SATurbulentViscosity.C b/src/kernels/SATurbulentViscosity.C new file mode 100644 index 0000000000..41d3cddc01 --- /dev/null +++ b/src/kernels/SATurbulentViscosity.C @@ -0,0 +1,36 @@ +#include "SATurbulentViscosity.h" + +registerMooseObject("MoltresApp", SATurbulentViscosity); + +InputParameters +SATurbulentViscosity::validParams() +{ + InputParameters params = ADKernel::validParams(); + params.addClassDescription("Adds all terms in the turbulent viscosity equation from the " + "Spalart-Allmaras model"); + params.addParam("mu_name", "mu", + "The name of the dynamic viscosity material property"); + return params; +} + +SATurbulentViscosity::SATurbulentViscosity(const InputParameters & parameters) + : ADKernel(parameters), + _mu(getADMaterialProperty("mu_name")), + _convection_strong_residual(getADMaterialProperty("convection_strong_residual")), + _destruction_strong_residual(getADMaterialProperty("destruction_strong_residual")), + _diffusion_strong_residual(getADMaterialProperty("diffusion_strong_residual")), + _source_strong_residual(getADMaterialProperty("source_strong_residual")) +{ +} + +ADReal +SATurbulentViscosity::computeQpResidual() +{ + Real sigma = 2.0 / 3.0; + Real cb2 = 0.622; + ADReal diffusion_part = 1.0 / sigma * ((_mu[_qp] + _u[_qp]) * \ + _grad_u[_qp] * _grad_test[_i][_qp] - cb2 * _grad_u[_qp] * _grad_u[_qp] * \ + _test[_i][_qp]); + return (_convection_strong_residual[_qp] + _destruction_strong_residual[_qp] + \ + _source_strong_residual[_qp]) * _test[_i][_qp] + diffusion_part; +} diff --git a/src/kernels/SATurbulentViscositySUPG.C b/src/kernels/SATurbulentViscositySUPG.C new file mode 100644 index 0000000000..6eaedbdd46 --- /dev/null +++ b/src/kernels/SATurbulentViscositySUPG.C @@ -0,0 +1,25 @@ +#include "SATurbulentViscositySUPG.h" + +registerMooseObject("MoltresApp", SATurbulentViscositySUPG); + +InputParameters +SATurbulentViscositySUPG::validParams() +{ + InputParameters params = ADKernelSUPG::validParams(); + params.addClassDescription("Adds the SUPG stabilization to the Spalart-Allmaras turbulent " + "viscosity equation"); + params.set("tau_name") = "tau_viscosity"; + return params; +} + +SATurbulentViscositySUPG::SATurbulentViscositySUPG(const InputParameters & parameters) + : ADKernelSUPG(parameters), + _visc_strong_residual(getADMaterialProperty("visc_strong_residual")) +{ +} + +ADReal +SATurbulentViscositySUPG::precomputeQpStrongResidual() +{ + return _visc_strong_residual[_qp]; +} diff --git a/src/materials/SATauMaterial.C b/src/materials/SATauMaterial.C new file mode 100644 index 0000000000..54eebf5972 --- /dev/null +++ b/src/materials/SATauMaterial.C @@ -0,0 +1,146 @@ +#include "SATauMaterial.h" +#include "INSADTauMaterial.h" + +registerMooseObject("MoltresApp", SATauMaterial); + +template +InputParameters +SATauMaterialTempl::validParams() +{ + InputParameters params = T::validParams(); + params.addClassDescription( + "This is the material class used to compute the stabilization parameter tau_viscosity " + "for the Spalart-Allmaras turbulent viscosity equation."); + params.addRequiredCoupledVar("mu_tilde", "Spalart-Allmaras turbulence viscosity variable"); + params.addRequiredCoupledVar("wall_distance_var", "Wall distance aux variable name"); + params.addParam("use_ft2_term", false, "Whether to apply the f_t2 term in the equation"); + return params; +} + +template +SATauMaterialTempl::SATauMaterialTempl(const InputParameters & parameters) + : T(parameters), + _sigma(2.0/3.0), + _cb1(0.1355), + _cb2(0.622), + _kappa(0.41), + _cw1(_cb1 / _kappa / _kappa + (1 + _cb2) / _sigma), + _cw2(0.3), + _cw3(2.0), + _cv1(7.1), + _ct1(1.0), + _ct2(2.0), + _ct3(1.2), + _ct4(0.5), + _mu_tilde(this->template adCoupledValue("mu_tilde")), + _grad_mu(this->template adCoupledGradient("mu_tilde")), + _tau_visc(this->template declareADProperty("tau_viscosity")), + _wall_dist(this->template coupledValue("wall_distance_var")), + _use_ft2_term(this->template getParam("use_ft2_term")), + _strain_mag(this->template declareADProperty("strain_mag")), + _convection_strong_residual(this->template + declareADProperty("convection_strong_residual")), + _destruction_strong_residual(this->template + declareADProperty("destruction_strong_residual")), + _diffusion_strong_residual(this->template + declareADProperty("diffusion_strong_residual")), + _source_strong_residual(this->template + declareADProperty("source_strong_residual")), + _time_strong_residual(this->template + declareADProperty("time_strong_residual")), + _visc_strong_residual(this->template + declareADProperty("visc_strong_residual")) +{ +} + +template +void +SATauMaterialTempl::subdomainSetup() +{ + T::subdomainSetup(); + + if (_has_transient) + _visc_dot = & this->template adCoupledDot("mu_tilde"); + else + _visc_dot = nullptr; +} + +template +void +SATauMaterialTempl::computeQpProperties() +{ + T::computeQpProperties(); + + // Compute strain rate and vorticity magnitudes + _strain_mag[_qp] = 2.0 * Utility::pow<2>(_grad_velocity[_qp](0, 0)) + + 2.0 * Utility::pow<2>(_grad_velocity[_qp](1, 1)) + + 2.0 * Utility::pow<2>(_grad_velocity[_qp](2, 2)) + + Utility::pow<2>(_grad_velocity[_qp](0, 2) + _grad_velocity[_qp](2, 0)) + + Utility::pow<2>(_grad_velocity[_qp](0, 1) + _grad_velocity[_qp](1, 0)) + + Utility::pow<2>(_grad_velocity[_qp](1, 2) + _grad_velocity[_qp](2, 1)); + _strain_mag[_qp] = std::sqrt(_strain_mag[_qp] + 1e-16); + ADReal vorticity_mag = + Utility::pow<2>(_grad_velocity[_qp](0, 2) - _grad_velocity[_qp](2, 0)) + + Utility::pow<2>(_grad_velocity[_qp](0, 1) - _grad_velocity[_qp](1, 0)) + + Utility::pow<2>(_grad_velocity[_qp](1, 2) - _grad_velocity[_qp](2, 1)); + vorticity_mag = std::sqrt(vorticity_mag + 1e-16); + + // Compute relevant parameters for the SA equation + const Real d = std::max(_wall_dist[_qp], 1e-16); // Avoid potential division by zero + const ADReal chi = _mu_tilde[_qp] / _mu[_qp]; + const ADReal fv1 = Utility::pow<3>(chi) / (Utility::pow<3>(chi) + Utility::pow<3>(_cv1)); + const ADReal fv2 = 1.0 - chi / (1. + chi * fv1); + const ADReal S_tilde = + vorticity_mag + _mu_tilde[_qp] * fv2 / (_kappa * _kappa * d * d * _rho[_qp]); + const ADReal S = S_tilde + 2 * std::min(0.0, _strain_mag[_qp] - vorticity_mag); + ADReal r; + if (S_tilde <= 0.0) // Avoid potential division by zero + r = 10.; + else + r = std::min(_mu_tilde[_qp] / (S_tilde * _kappa * _kappa * d * d * _rho[_qp]), 10.0); + const ADReal g = r + _cw2 * (Utility::pow<6>(r) - r); + const ADReal fw = g * std::pow((1. + Utility::pow<6>(_cw3)) / + (Utility::pow<6>(g) + Utility::pow<6>(_cw3)), + 1.0 / 6.0); + + // Compute strong forms of the SA equation + if (_use_ft2_term) // Whether to apply the f_t2 term in the SA equation + { + const ADReal ft2 = _ct3 * std::exp(-_ct4 * chi * chi); + _destruction_strong_residual[_qp] = + (_cw1 * fw - _cb1 * ft2 / _kappa / _kappa) * Utility::pow<2>(_mu_tilde[_qp] / d); + _source_strong_residual[_qp] = -(1 - ft2) * _rho[_qp] * _cb1 * S * _mu_tilde[_qp]; + } + else + { + _destruction_strong_residual[_qp] = _cw1 * fw * Utility::pow<2>(_mu_tilde[_qp] / d); + _source_strong_residual[_qp] = -_rho[_qp] * _cb1 * S * _mu_tilde[_qp]; + } + _convection_strong_residual[_qp] = _rho[_qp] * _velocity[_qp] * _grad_mu[_qp]; + _diffusion_strong_residual[_qp] = -1.0 / _sigma * _cb2 * (_grad_mu[_qp] * _grad_mu[_qp]); + if (_has_transient) + _time_strong_residual[_qp] = (*_visc_dot)[_qp] * _rho[_qp]; + _visc_strong_residual[_qp] = _has_transient ? _time_strong_residual[_qp] : 0.0; + _visc_strong_residual[_qp] += (_convection_strong_residual[_qp] + + _destruction_strong_residual[_qp] + + _diffusion_strong_residual[_qp] + + _source_strong_residual[_qp]); + + // Compute the tau stabilization parameter for mu_tilde SUPG stabilization + const auto nu_visc = (_mu[_qp] + _mu_tilde[_qp]) / _rho[_qp] / _sigma; + const auto transient_part = _has_transient ? 4.0 / (_dt * _dt) : 0.0; + const auto speed = NS::computeSpeed(_velocity[_qp]); + _tau_visc[_qp] = _alpha / std::sqrt(transient_part + + (2.0 * speed / _hmax) * (2.0 * speed / _hmax) + + 9.0 * (4.0 * nu_visc / (_hmax * _hmax)) * + 4.0 * (nu_visc / (_hmax * _hmax))); + + // Replace the nu value in the tau stabilization parameter for INS SUPG stabilization + const auto nu = (_mu[_qp] + _mu_tilde[_qp] * fv1) / _rho[_qp]; + _tau[_qp] = _alpha / std::sqrt(transient_part + + (2.0 * speed / _hmax) * (2.0 * speed / _hmax) + + 9.0 * (4.0 * nu / (_hmax * _hmax)) * + 4.0 * (nu / (_hmax * _hmax))); +} + +template class SATauMaterialTempl; diff --git a/src/materials/SATauStabilized3Eqn.C b/src/materials/SATauStabilized3Eqn.C new file mode 100644 index 0000000000..17a4b6d6e4 --- /dev/null +++ b/src/materials/SATauStabilized3Eqn.C @@ -0,0 +1,44 @@ +#include "SATauStabilized3Eqn.h" + +registerMooseObject("MoltresApp", SATauStabilized3Eqn); + +InputParameters +SATauStabilized3Eqn::validParams() +{ + InputParameters params = SATauMaterialTempl::validParams(); + params.addClassDescription("This is the material class used to add the turbulent diffusion " + "strong residual contribution for the stabilization of the " + "energy equation."); + params.addParam("prandtl_number", 0.85, "Turbulent Prandtl number"); + return params; +} + +SATauStabilized3Eqn::SATauStabilized3Eqn(const InputParameters & parameters) + : SATauMaterialTempl(parameters), + _pr(getParam("prandtl_number")) +{ +} + +void +SATauStabilized3Eqn::computeQpProperties() +{ + SATauMaterialTempl::computeQpProperties(); + + const ADReal chi = _mu_tilde[_qp] / _mu[_qp]; + const Real cv1 = 7.1; + const ADReal fv1 = Utility::pow<3>(chi) / (Utility::pow<3>(chi) + Utility::pow<3>(cv1)); + + const auto dissipation_coefficient = + _k[_qp] / (_rho[_qp] * _cp[_qp]) + _mu_tilde[_qp] * fv1 / (_pr * _rho[_qp]); + const auto transient_part = _has_energy_transient ? 4.0 / (_dt * _dt) : 0.0; + const auto speed = NS::computeSpeed(_velocity[_qp]); + _tau_energy[_qp] = + _alpha / std::sqrt(transient_part + (2.0 * speed / _hmax) * (2.0 * speed / _hmax) + + 9.0 * (4.0 * dissipation_coefficient / (_hmax * _hmax)) * + (4.0 * dissipation_coefficient / (_hmax * _hmax))); + + // Turbulent diffusion + _temperature_strong_residual[_qp] -= + _cp[_qp] * fv1 / _pr * (_mu_tilde[_qp] * _second_temperature[_qp].tr() + + _grad_mu[_qp] * _grad_temperature[_qp]); +} diff --git a/tests/sa-model/channel_flow.i b/tests/sa-model/channel_flow.i new file mode 100644 index 0000000000..a71aaade11 --- /dev/null +++ b/tests/sa-model/channel_flow.i @@ -0,0 +1,307 @@ +Re = 1.375e4 +density = 1 # kg cm-3 +D = 1 +inlet = 1 +viscosity = '${fparse density * inlet * D / Re}' # dynamic viscosity, mu = nu * rho, kg cm-1 s-1 +alpha = 0.33333 # INS SUPG and PSPG stabilization parameter + +[GlobalParams] + integrate_p_by_parts = false + viscous_form = 'traction' + pressure = p + velocity = vel +[] + +[Mesh] + [channel] + type = CartesianMeshGenerator + dim = 2 + dx = '0.05 0.2 19.75' + dy = '0.3 0.1 0.05 0.05' + ix = '20 1 79' + iy = '10 5 5 10' + [] + [corner_node] + type = ExtraNodesetGenerator + input = channel + new_boundary = 'pinned_node' + coord = '20 0.5' + [] +[] + +[Problem] + type = FEProblem +[] + +[Variables] + [vel] + family = LAGRANGE_VEC + order = FIRST + [] + [p] + [] + [mu_tilde] + initial_condition = '${fparse viscosity * 5}' + [] +[] + +[AuxVariables] + [mu_turbulence] + [] + [wall_dist] + [] + [velx] + [] + [vely] + [] + [wall_shear_stress] + family = MONOMIAL + order = CONSTANT + [] + [turbulent_stress] + family = MONOMIAL + order = CONSTANT + [] +[] + +[Kernels] + [mass] + type = INSADMass + variable = p + [] + [mass_pspg] + type = INSADMassPSPG + variable = p + [] + [momentum_advection] + type = INSADMomentumAdvection + variable = vel + [] + [momentum_viscous] + type = INSADMomentumViscous + variable = vel + [] + [momentum_turbulent_viscous] + type = INSADMomentumTurbulentViscous + variable = vel + mu_tilde = mu_tilde + [] + [momentum_pressure] + type = INSADMomentumPressure + variable = vel + [] + [momentum_supg] + type = INSADMomentumSUPG + variable = vel + [] + [momentum_time] + type = INSADMomentumTimeDerivative + variable = vel + [] + + [mu_tilde_time] + type = SATimeDerivative + variable = mu_tilde + [] + [mu_tilde_space] + type = SATurbulentViscosity + variable = mu_tilde + [] + [mu_tilde_supg] + type = SATurbulentViscositySUPG + variable = mu_tilde + [] +[] + +[AuxKernels] + [mu_space] + type = SATurbulentViscosityAux + variable = mu_turbulence + mu_tilde = mu_tilde + [] + [wall_distance] + type = WallDistanceAux + variable = wall_dist + walls = 'top' + execute_on = initial + [] + [velx] + type = VectorVariableComponentAux + variable = velx + vector_variable = vel + component = x + [] + [vely] + type = VectorVariableComponentAux + variable = vely + vector_variable = vel + component = y + [] + [wall_shear_stress] + type = WallShearStressAux + variable = wall_shear_stress + boundary = top + [] + [turbulent_stress] + type = TurbulentStressAux + variable = turbulent_stress + mu_tilde = mu_tilde + [] +[] + +[Functions] + [xfunc] + type = ParsedFunction + expression = '${fparse inlet}' + [] + [mu_func] + type = ParsedFunction + expression = 'if(y=0.5, 0, ${fparse viscosity * 5})' + [] +[] + +[ICs] + [velocity] + type = VectorFunctionIC + function_x = xfunc + function_y = 0 + variable = vel + [] +[] + +[BCs] + [no_slip] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'top' + function_x = 0 + function_y = 0 + [] + [inlet] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'left' + function_x = xfunc + function_y = 0 + [] + [outlet] + type = INSADMomentumTurbulentNoBCBC + variable = vel + boundary = 'right' + mu_tilde = mu_tilde + [] + [outlet_supg] + type = INSADMomentumSUPGBC + variable = vel + boundary = 'right' + [] + [outlet_pspg] + type = INSADMassPSPGBC + variable = p + boundary = 'right' + [] + [symmetry] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'bottom' + set_x_comp = false + function_y = 0 + [] + + [pressure_pin] + type = DirichletBC + variable = p + boundary = 'pinned_node' + value = 0 + [] + + [mu_inlet] + type = FunctionDirichletBC + variable = mu_tilde + boundary = 'left' + function = mu_func + [] + [mu_wall] + type = DirichletBC + variable = mu_tilde + boundary = 'top' + value = 0 + [] + [mu_outlet] + type = SATurbulentViscosityNoBCBC + variable = mu_tilde + boundary = 'right' + [] + [mu_outlet_supg] + type = SATurbulentViscositySUPGBC + variable = mu_tilde + boundary = 'right' + [] +[] + +[Materials] + [ad_mat] + type = ADGenericConstantMaterial + prop_names = 'rho mu' + prop_values = '${density} ${viscosity}' + [] + [sa_mat] + type = SATauMaterial + alpha = ${alpha} + mu_tilde = mu_tilde + wall_distance_var = wall_dist + use_ft2_term = true + [] +[] + +[Executioner] + type = Transient + + solve_type = 'NEWTON' + petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor' + petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_package' + petsc_options_value = 'lu NONZERO superlu_dist' + line_search = none + + automatic_scaling = true + compute_scaling_once = false + resid_vs_jac_scaling_param = 0.1 + scaling_group_variables = 'vel; p; mu_tilde' + off_diagonals_in_auto_scaling = true + + nl_rel_tol = 1e-8 + nl_abs_tol = 1e-14 + nl_max_its = 10 + + steady_state_detection = true + steady_state_tolerance = 1e-10 + dtmin = 1e-0 + dtmax = 10 + [TimeStepper] + type = IterationAdaptiveDT + dt = 1e-0 + cutback_factor = 0.9 + growth_factor = 1.1 + optimal_iterations = 6 + iteration_window = 1 + linear_iteration_ratio = 1000 + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Outputs] + [exodus] + type = Exodus + execute_on = final + [] +[] + +[Debug] + show_var_residual_norms = true +[] diff --git a/tests/sa-model/channel_flow_with_heat.i b/tests/sa-model/channel_flow_with_heat.i new file mode 100644 index 0000000000..f98a8bd688 --- /dev/null +++ b/tests/sa-model/channel_flow_with_heat.i @@ -0,0 +1,347 @@ +Re = 1.375e4 +density = 1 # kg cm-3 +D = 1 +inlet = 1 +viscosity = '${fparse density * inlet * D / Re}' # dynamic viscosity, mu = nu * rho, kg cm-1 s-1 +alpha = 0.33333 # INS SUPG and PSPG stabilization parameter +cp = 1 +k = '${fparse viscosity * cp / 0.85}' + +[GlobalParams] + integrate_p_by_parts = false + viscous_form = 'traction' + pressure = p + velocity = vel + temperature = temp +[] + +[Mesh] + [channel] + type = CartesianMeshGenerator + dim = 2 + dx = '0.05 0.2 19.75' + dy = '0.3 0.1 0.05 0.05' + ix = '20 1 79' + iy = '10 5 5 10' + [] + [corner_node] + type = ExtraNodesetGenerator + input = channel + new_boundary = 'pinned_node' + coord = '20 0.5' + [] +[] + +[Problem] + type = FEProblem +[] + +[Variables] + [vel] + family = LAGRANGE_VEC + order = FIRST + [] + [p] + [] + [mu_tilde] + initial_condition = '${fparse viscosity * 5}' + [] + [temp] + initial_condition = 300 + [] +[] + +[AuxVariables] + [mu_turbulence] + [] + [wall_dist] + [] + [velx] + [] + [vely] + [] + [wall_shear_stress] + family = MONOMIAL + order = CONSTANT + [] + [turbulent_stress] + family = MONOMIAL + order = CONSTANT + [] +[] + +[Kernels] + [mass] + type = INSADMass + variable = p + [] + [mass_pspg] + type = INSADMassPSPG + variable = p + [] + + [momentum_advection] + type = INSADMomentumAdvection + variable = vel + [] + [momentum_viscous] + type = INSADMomentumViscous + variable = vel + [] + [momentum_turbulent_viscous] + type = INSADMomentumTurbulentViscous + variable = vel + mu_tilde = mu_tilde + [] + [momentum_pressure] + type = INSADMomentumPressure + variable = vel + [] + [momentum_supg] + type = INSADMomentumSUPG + variable = vel + [] + [momentum_time] + type = INSADMomentumTimeDerivative + variable = vel + [] + + [mu_tilde_time] + type = SATimeDerivative + variable = mu_tilde + [] + [mu_tilde_space] + type = SATurbulentViscosity + variable = mu_tilde + [] + [mu_tilde_supg] + type = SATurbulentViscositySUPG + variable = mu_tilde + [] + + [temperature_advection] + type = INSADEnergyAdvection + variable = temp + [] + [temp_conduction] + type = ADHeatConduction + variable = temp + thermal_conductivity = 'k' + [] + [temp_turbulent_diffusion] + type = INSADEnergyTurbulentDiffusion + variable = temp + mu_tilde = mu_tilde + [] + [temp_source] + type = INSADEnergySource + variable = temp + source_function = '1' + [] + [temp_supg] + type = INSADEnergySUPG + variable = temp + [] + [temp_time] + type = INSADHeatConductionTimeDerivative + variable = temp + [] +[] + +[AuxKernels] + [mu_space] + type = SATurbulentViscosityAux + variable = mu_turbulence + mu_tilde = mu_tilde + [] + [wall_distance] + type = WallDistanceAux + variable = wall_dist + walls = 'top' + execute_on = initial + [] + [velx] + type = VectorVariableComponentAux + variable = velx + vector_variable = vel + component = x + [] + [vely] + type = VectorVariableComponentAux + variable = vely + vector_variable = vel + component = y + [] + [wall_shear_stress] + type = WallShearStressAux + variable = wall_shear_stress + boundary = top + [] + [turbulent_stress] + type = TurbulentStressAux + variable = turbulent_stress + mu_tilde = mu_tilde + [] +[] + +[Functions] + [xfunc] + type = ParsedFunction + expression = '${fparse inlet}' + [] + [mu_func] + type = ParsedFunction + expression = 'if(y=0.5, 0, ${fparse viscosity * 5})' + [] +[] + +[ICs] + [velocity] + type = VectorFunctionIC + function_x = xfunc + function_y = 0 + variable = vel + [] +[] + +[BCs] + [no_slip] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'top' + function_x = 0 + function_y = 0 + [] + [inlet] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'left' + function_x = xfunc + function_y = 0 + [] + [outlet] + type = INSADMomentumTurbulentNoBCBC + variable = vel + boundary = 'right' + mu_tilde = mu_tilde + [] + [outlet_supg] + type = INSADMomentumSUPGBC + variable = vel + boundary = 'right' + [] + [outlet_pspg] + type = INSADMassPSPGBC + variable = p + boundary = 'right' + [] + [symmetry] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'bottom' + set_x_comp = false + function_y = 0 + [] + + [pressure_pin] + type = DirichletBC + variable = p + boundary = 'pinned_node' + value = 0 + [] + + [mu_inlet] + type = FunctionDirichletBC + variable = mu_tilde + boundary = 'left' + function = mu_func + [] + [mu_wall] + type = DirichletBC + variable = mu_tilde + boundary = 'top' + value = 0 + [] + [mu_outlet] + type = SATurbulentViscosityNoBCBC + variable = mu_tilde + boundary = 'right' + [] + [mu_outlet_supg] + type = SATurbulentViscositySUPGBC + variable = mu_tilde + boundary = 'right' + [] + + [temperature_inlet] + type = DirichletBC + variable = temp + boundary = 'left' + value = 300 + [] +[] + +[Materials] + [ad_mat] + type = ADGenericConstantMaterial + prop_names = 'rho mu cp k grad_k' + prop_values = '${density} ${viscosity} ${cp} ${k} 0' + [] + [sa_mat] + type = SATauStabilized3Eqn + alpha = ${alpha} + mu_tilde = mu_tilde + wall_distance_var = wall_dist + use_ft2_term = true + [] +[] + +[Executioner] + type = Transient + end_time = 30 + + solve_type = 'NEWTON' + petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor' + petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_package' + petsc_options_value = 'lu NONZERO superlu_dist' + line_search = none + + automatic_scaling = true + compute_scaling_once = false + resid_vs_jac_scaling_param = 0.1 + off_diagonals_in_auto_scaling = true + + nl_rel_tol = 1e-8 + nl_abs_tol = 1e-13 + nl_max_its = 10 + + dtmin = 1 + dtmax = 10 + [TimeStepper] + type = IterationAdaptiveDT + dt = 1 + cutback_factor = 0.9 + growth_factor = 1.2 + optimal_iterations = 6 + iteration_window = 1 + linear_iteration_ratio = 1000 + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Outputs] + [exodus] + type = Exodus + execute_on = final + [] +[] + +[Debug] + show_var_residual_norms = true +[] diff --git a/tests/sa-model/channel_flow_with_precursors.i b/tests/sa-model/channel_flow_with_precursors.i new file mode 100644 index 0000000000..d55987ab26 --- /dev/null +++ b/tests/sa-model/channel_flow_with_precursors.i @@ -0,0 +1,360 @@ +Re = 1.375e4 +density = 1 # kg cm-3 +D = 1 +inlet = 1 +viscosity = '${fparse density * inlet * D / Re}' # dynamic viscosity, mu = nu * rho, kg cm-1 s-1 +alpha = .33333 # INS SUPG and PSPG stabilization parameter +diff = '${fparse viscosity / density / 0.85}' + +[GlobalParams] + integrate_p_by_parts = false + viscous_form = 'traction' + pressure = p + velocity = vel + use_exp_form = false +[] + +[Mesh] + [channel] + type = CartesianMeshGenerator + dim = 2 + dx = '0.05 0.2 19.75' + dy = '0.3 0.1 0.05 0.05' + ix = '20 1 79' + iy = '10 5 5 10' + [] + [corner_node] + type = ExtraNodesetGenerator + input = channel + new_boundary = 'pinned_node' + coord = '20 0.5' + [] +[] + +[Problem] + type = FEProblem +[] + +[Variables] + [vel] + family = LAGRANGE_VEC + order = FIRST + [] + [p] + [] + [mu_tilde] + initial_condition = '${fparse viscosity * 5}' + [] + [prec] + family = MONOMIAL + order = CONSTANT + initial_condition = 300 + scaling = 1e-3 + [] +[] + +[AuxVariables] + [mu_turbulence] + [] + [wall_dist] + [] + [velx] + [] + [vely] + [] + [wall_shear_stress] + family = MONOMIAL + order = CONSTANT + [] + [turbulent_stress] + family = MONOMIAL + order = CONSTANT + [] +[] + +[Kernels] + [mass] + type = INSADMass + variable = p + [] + [mass_pspg] + type = INSADMassPSPG + variable = p + [] + [momentum_advection] + type = INSADMomentumAdvection + variable = vel + [] + [momentum_viscous] + type = INSADMomentumViscous + variable = vel + [] + [momentum_turbulent_viscous] + type = INSADMomentumTurbulentViscous + variable = vel + mu_tilde = mu_tilde + [] + [momentum_pressure] + type = INSADMomentumPressure + variable = vel + [] + [momentum_supg] + type = INSADMomentumSUPG + variable = vel + [] + [momentum_time] + type = INSADMomentumTimeDerivative + variable = vel + [] + + [mu_tilde_time] + type = SATimeDerivative + variable = mu_tilde + [] + [mu_tilde_space] + type = SATurbulentViscosity + variable = mu_tilde + [] + [mu_tilde_supg] + type = SATurbulentViscositySUPG + variable = mu_tilde + [] + + [prec_source] + type = BodyForce + variable = prec + [] + [prec_time] + type = ScalarTransportTimeDerivative + variable = prec + [] +[] + +[DGKernels] + [diffusion] + type = DGDiffusion + variable = prec + epsilon = -1 + sigma = 1e-3 + diff = '${diff}' + [] + [turbulent_diffusion] + type = DGTurbulentDiffusion + variable = prec + mu_tilde = mu_tilde + epsilon = -1 + sigma = 1e-3 + [] + [advection] + type = DGCoupledAdvection + variable = prec + uvel = velx + vvel = vely + [] +[] + +[AuxKernels] + [mu_space] + type = SATurbulentViscosityAux + variable = mu_turbulence + mu_tilde = mu_tilde + [] + [wall_distance] + type = WallDistanceAux + variable = wall_dist + walls = 'top' + execute_on = initial + [] + [velx] + type = VectorVariableComponentAux + variable = velx + vector_variable = vel + component = x + [] + [vely] + type = VectorVariableComponentAux + variable = vely + vector_variable = vel + component = y + [] + [wall_shear_stress] + type = WallShearStressAux + variable = wall_shear_stress + boundary = top + [] + [turbulent_stress] + type = TurbulentStressAux + variable = turbulent_stress + mu_tilde = mu_tilde + [] +[] + +[Functions] + [xfunc] + type = ParsedFunction + expression = '${fparse inlet}' + [] + [mu_func] + type = ParsedFunction + expression = 'if(y=0.5, 0, ${fparse viscosity * 5})' + [] +[] + +[ICs] + [velocity] + type = VectorFunctionIC + function_x = xfunc + function_y = 0 + variable = vel + [] +[] + +[BCs] + [no_slip] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'top' + function_x = 0 + function_y = 0 + [] + [inlet] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'left' + function_x = xfunc + function_y = 0 + [] + [outlet] + type = INSADMomentumTurbulentNoBCBC + variable = vel + boundary = 'right' + mu_tilde = mu_tilde + [] + [outlet_supg] + type = INSADMomentumSUPGBC + variable = vel + boundary = 'right' + [] + [outlet_pspg] + type = INSADMassPSPGBC + variable = p + boundary = 'right' + [] + [symmetry] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'bottom' + set_x_comp = false + function_y = 0 + [] + + [pressure_pin] + type = DirichletBC + variable = p + boundary = 'pinned_node' + value = 0 + [] + + [mu_inlet] + type = FunctionDirichletBC + variable = mu_tilde + boundary = 'left' + function = mu_func + [] + [mu_wall] + type = DirichletBC + variable = mu_tilde + boundary = 'top' + value = 0 + [] + [mu_outlet] + type = SATurbulentViscosityNoBCBC + variable = mu_tilde + boundary = 'right' + [] + [mu_outlet_supg] + type = SATurbulentViscositySUPGBC + variable = mu_tilde + boundary = 'right' + [] + + [prec_inlet] + type = InflowBC + variable = prec + boundary = 'left' + uu = '${inlet}' + inlet_conc = 300 + [] + [prec_outlet] + type = CoupledScalarAdvectionNoBCBC + variable = prec + boundary = 'right' + u = velx + [] +[] + +[Materials] + [ad_mat] + type = ADGenericConstantMaterial + prop_names = 'rho mu' + prop_values = '${density} ${viscosity}' + [] + [sa_mat] + type = SATauMaterial + alpha = ${alpha} + mu_tilde = mu_tilde + wall_distance_var = wall_dist + use_ft2_term = true + [] +[] + +[Executioner] + type = Transient + end_time = 30 + + solve_type = 'NEWTON' + petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor' + petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_package' + petsc_options_value = 'lu NONZERO superlu_dist' + line_search = none + + automatic_scaling = true + compute_scaling_once = false + resid_vs_jac_scaling_param = .1 + ignore_variables_for_autoscaling = 'prec' + off_diagonals_in_auto_scaling = true + + nl_rel_tol = 1e-8 + nl_abs_tol = 1e-13 + nl_max_its = 10 + + dtmin = 1 + dtmax = 10 + [TimeStepper] + type = IterationAdaptiveDT + dt = 1 + cutback_factor = 0.9 + growth_factor = 1.2 + optimal_iterations = 6 + iteration_window = 1 + linear_iteration_ratio = 1000 + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Outputs] + [exodus] + type = Exodus + execute_on = final + [] +[] + +[Debug] + show_var_residual_norms = true +[] diff --git a/tests/sa-model/gold/channel_flow_exodus.e b/tests/sa-model/gold/channel_flow_exodus.e new file mode 100644 index 0000000000..b8b8f7530b Binary files /dev/null and b/tests/sa-model/gold/channel_flow_exodus.e differ diff --git a/tests/sa-model/gold/channel_flow_with_heat_exodus.e b/tests/sa-model/gold/channel_flow_with_heat_exodus.e new file mode 100644 index 0000000000..25d05c63df Binary files /dev/null and b/tests/sa-model/gold/channel_flow_with_heat_exodus.e differ diff --git a/tests/sa-model/gold/channel_flow_with_precursors_exodus.e b/tests/sa-model/gold/channel_flow_with_precursors_exodus.e new file mode 100644 index 0000000000..810a2828e3 Binary files /dev/null and b/tests/sa-model/gold/channel_flow_with_precursors_exodus.e differ diff --git a/tests/sa-model/gold/heat_turbulent_diffusion_exodus.e b/tests/sa-model/gold/heat_turbulent_diffusion_exodus.e new file mode 100644 index 0000000000..83c56bd5c0 Binary files /dev/null and b/tests/sa-model/gold/heat_turbulent_diffusion_exodus.e differ diff --git a/tests/sa-model/gold/pipe_flow_exodus.e b/tests/sa-model/gold/pipe_flow_exodus.e new file mode 100644 index 0000000000..3b39a05508 Binary files /dev/null and b/tests/sa-model/gold/pipe_flow_exodus.e differ diff --git a/tests/sa-model/gold/precursor_turbulent_diffusion_exodus.e b/tests/sa-model/gold/precursor_turbulent_diffusion_exodus.e new file mode 100644 index 0000000000..a0ec78d138 Binary files /dev/null and b/tests/sa-model/gold/precursor_turbulent_diffusion_exodus.e differ diff --git a/tests/sa-model/heat_turbulent_diffusion.i b/tests/sa-model/heat_turbulent_diffusion.i new file mode 100644 index 0000000000..bcf039d458 --- /dev/null +++ b/tests/sa-model/heat_turbulent_diffusion.i @@ -0,0 +1,87 @@ +[GlobalParams] + pressure = p + velocity = vel +[] + +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 2 + nx = 10 + ny = 10 + [] +[] + +[Variables] + [temp] + [] +[] + +[AuxVariables] + [mu_tilde] + initial_condition = 1 + [] + [vel] + family = LAGRANGE_VEC + order = FIRST + [] + [p] + [] + [wall_dist] + [] +[] + +[Kernels] + [turbulent_diffusion] + type = INSADEnergyTurbulentDiffusion + variable = temp + mu_tilde = mu_tilde + [] +[] + +[BCs] + [left] + type = ADDirichletBC + variable = temp + boundary = 'left' + value = 1 + [] + [right] + type = ADDirichletBC + variable = temp + boundary = 'right' + value = 0 + [] +[] + +[Materials] + [const_mat] + type = ADGenericConstantMaterial + prop_names = 'cp rho mu k grad_k' + prop_values = '1 1 1 1 0' + [] + [sa_mat] + type = SATauStabilized3Eqn + alpha = 0.33333 + mu_tilde = mu_tilde + wall_distance_var = wall_dist + use_ft2_term = true + temperature = temp + [] +[] + +[Executioner] + type = Steady + solve_type = 'NEWTON' +[] + +[Outputs] + [exodus] + type = Exodus + execute_on = final + [] +[] + +[Debug] + show_var_residual_norms = true +[] diff --git a/tests/sa-model/pipe_flow.i b/tests/sa-model/pipe_flow.i new file mode 100644 index 0000000000..77f6df6a5d --- /dev/null +++ b/tests/sa-model/pipe_flow.i @@ -0,0 +1,300 @@ +Re = 4e4 +density = 1 # kg cm-3 +D = 1 +inlet = 1 +viscosity = '${fparse density * inlet * D / Re}' # dynamic viscosity, mu = nu * rho, kg cm-1 s-1 +alpha = 0.33333 # INS SUPG and PSPG stabilization parameter + +[GlobalParams] + integrate_p_by_parts = false + viscous_form = 'traction' + pressure = p + velocity = vel +[] + +[Mesh] + coord_type = 'RZ' + [pipe] + type = CartesianMeshGenerator + dim = 2 + dx = '0.3 0.1 0.05 0.034 0.016' + dy = '0.025 0.475 19.5' + ix = '10 5 5 4 10' + iy = '20 1 39' + [] + [corner_node] + type = ExtraNodesetGenerator + input = pipe + new_boundary = 'pinned_node' + coord = '0.5 20' + [] +[] + +[Problem] + type = FEProblem +[] + +[Variables] + [vel] + family = LAGRANGE_VEC + order = FIRST + [] + [p] + [] + [mu_tilde] + initial_condition = '${fparse viscosity * 5}' + [] +[] + +[AuxVariables] + [mu_turbulence] + [] + [wall_dist] + [] + [velx] + [] + [vely] + [] + [wall_shear_stress] + family = MONOMIAL + order = CONSTANT + [] + [turbulent_stress] + family = MONOMIAL + order = CONSTANT + [] +[] + +[Kernels] + [mass] + type = INSADMass + variable = p + [] + [mass_pspg] + type = INSADMassPSPG + variable = p + [] + [momentum_advection] + type = INSADMomentumAdvection + variable = vel + [] + [momentum_viscous] + type = INSADMomentumViscous + variable = vel + [] + [momentum_turbulent_viscous] + type = INSADMomentumTurbulentViscous + variable = vel + mu_tilde = mu_tilde + [] + [momentum_pressure] + type = INSADMomentumPressure + variable = vel + [] + [momentum_supg] + type = INSADMomentumSUPG + variable = vel + [] + [momentum_time] + type = INSADMomentumTimeDerivative + variable = vel + [] + + [mu_tilde_time] + type = SATimeDerivative + variable = mu_tilde + [] + [mu_tilde_space] + type = SATurbulentViscosity + variable = mu_tilde + [] + [mu_tilde_supg] + type = SATurbulentViscositySUPG + variable = mu_tilde + [] +[] + +[AuxKernels] + [mu_space] + type = SATurbulentViscosityAux + variable = mu_turbulence + mu_tilde = mu_tilde + [] + [wall_distance] + type = WallDistanceAux + variable = wall_dist + walls = 'right' + execute_on = initial + [] + [velx] + type = VectorVariableComponentAux + variable = velx + vector_variable = vel + component = x + [] + [vely] + type = VectorVariableComponentAux + variable = vely + vector_variable = vel + component = y + [] + [wall_shear_stress] + type = WallShearStressAux + variable = wall_shear_stress + boundary = right + [] + [turbulent_stress] + type = TurbulentStressAux + variable = turbulent_stress + mu_tilde = mu_tilde + [] +[] + +[Functions] + [yfunc] + type = ParsedFunction + expression = '${fparse inlet}' + [] + [mu_func] + type = ParsedFunction + expression = 'if(x=0.5, 0, ${fparse viscosity * 5})' + [] +[] + +[ICs] + [velocity] + type = VectorFunctionIC + function_x = 0 + function_y = yfunc + variable = vel + [] +[] + +[BCs] + [no_slip] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'right' + function_x = 0 + function_y = 0 + [] + [inlet] + type = ADVectorFunctionDirichletBC + variable = vel + boundary = 'bottom' + function_x = 0 + function_y = yfunc + [] + [outlet] + type = INSADMomentumTurbulentNoBCBC + variable = vel + boundary = 'top' + mu_tilde = mu_tilde + [] + [outlet_supg] + type = INSADMomentumSUPGBC + variable = vel + boundary = 'top' + [] + [outlet_pspg] + type = INSADMassPSPGBC + variable = p + boundary = 'top' + [] + + [pressure_pin] + type = DirichletBC + variable = p + boundary = 'pinned_node' + value = 0 + [] + + [mu_inlet] + type = FunctionDirichletBC + variable = mu_tilde + boundary = 'bottom' + function = mu_func + [] + [mu_wall] + type = DirichletBC + variable = mu_tilde + boundary = 'right' + value = 0 + [] + [mu_outlet] + type = SATurbulentViscosityNoBCBC + variable = mu_tilde + boundary = 'top' + [] + [mu_outlet_supg] + type = SATurbulentViscositySUPGBC + variable = mu_tilde + boundary = 'top' + [] +[] + +[Materials] + [ad_mat] + type = ADGenericConstantMaterial + prop_names = 'rho mu' + prop_values = '${density} ${viscosity}' + [] + [sa_mat] + type = SATauMaterial + alpha = ${alpha} + mu_tilde = mu_tilde + wall_distance_var = wall_dist + [] +[] + +[Executioner] + type = Transient + + solve_type = 'NEWTON' + petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor' + petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_package' + petsc_options_value = 'lu NONZERO superlu_dist' + line_search = none + + automatic_scaling = true + compute_scaling_once = false + resid_vs_jac_scaling_param = 0.1 + scaling_group_variables = 'vel; p; mu_tilde' + off_diagonals_in_auto_scaling = true + + nl_rel_tol = 1e-8 + nl_abs_tol = 1e-14 + nl_max_its = 10 + + steady_state_detection = true + steady_state_tolerance = 1e-10 + dtmin = 1e-0 + dtmax = 10 + [TimeStepper] + type = IterationAdaptiveDT + dt = 1e-0 + cutback_factor = 0.9 + growth_factor = 1.1 + optimal_iterations = 6 + iteration_window = 1 + linear_iteration_ratio = 1000 + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Outputs] + [exodus] + type = Exodus + execute_on = final + [] +[] + +[Debug] + show_var_residual_norms = true +[] diff --git a/tests/sa-model/precursor_turbulent_diffusion.i b/tests/sa-model/precursor_turbulent_diffusion.i new file mode 100644 index 0000000000..dd5bf52925 --- /dev/null +++ b/tests/sa-model/precursor_turbulent_diffusion.i @@ -0,0 +1,90 @@ +[GlobalParams] + pressure = p + velocity = vel +[] + +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 2 + nx = 10 + ny = 10 + [] +[] + +[Variables] + [prec] + family = MONOMIAL + order = CONSTANT + [] +[] + +[AuxVariables] + [mu_tilde] + initial_condition = 1 + [] + [vel] + family = LAGRANGE_VEC + order = FIRST + [] + [p] + [] + [wall_dist] + [] +[] + +[DGKernels] + [turbulent_diffusion] + type = DGTurbulentDiffusion + variable = prec + mu_tilde = mu_tilde + [] +[] + +[BCs] + [left] + type = PenaltyDirichletBC + variable = prec + boundary = 'left' + value = 1 + penalty = 1e5 + [] + [right] + type = PenaltyDirichletBC + variable = prec + boundary = 'right' + value = 0 + penalty = 1e6 + [] +[] + +[Materials] + [const_mat] + type = ADGenericConstantMaterial + prop_names = 'rho mu' + prop_values = '1 1' + [] + [sa_mat] + type = SATauMaterial + alpha = 0.33333 + mu_tilde = mu_tilde + wall_distance_var = wall_dist + use_ft2_term = true + [] +[] + +[Executioner] + type = Steady + solve_type = 'NEWTON' +[] + +[Outputs] + [exodus] + type = Exodus + execute_on = final + [] +[] + +[Debug] + show_var_residual_norms = true +[] diff --git a/tests/sa-model/tests b/tests/sa-model/tests new file mode 100644 index 0000000000..c6bea0caac --- /dev/null +++ b/tests/sa-model/tests @@ -0,0 +1,40 @@ +[Tests] + [channel_flow] + type = 'Exodiff' + input = 'channel_flow.i' + exodiff = 'channel_flow_exodus.e' + heavy = true + max_time = 300 + [] + [pipe_flow] + type = 'Exodiff' + input = 'pipe_flow.i' + exodiff = 'pipe_flow_exodus.e' + heavy = true + max_time = 300 + [] + [heat_turbulent_diffusion] + type = 'Exodiff' + input = 'heat_turbulent_diffusion.i' + exodiff = 'heat_turbulent_diffusion_exodus.e' + [] + [precursor_turbulent_diffusion] + type = 'Exodiff' + input = 'precursor_turbulent_diffusion.i' + exodiff = 'precursor_turbulent_diffusion_exodus.e' + [] + [channel_flow_heat] + type = 'Exodiff' + input = 'channel_flow_with_heat.i' + exodiff = 'channel_flow_with_heat_exodus.e' + heavy = true + max_time = 300 + [] + [channel_flow_precursors] + type = 'Exodiff' + input = 'channel_flow_with_precursors.i' + exodiff = 'channel_flow_with_precursors_exodus.e' + heavy = true + max_time = 300 + [] +[]