From 965715029561b68b314767d13dad3156f6d14b5b Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Apr 2025 12:10:31 +0200 Subject: [PATCH 01/26] add sgs --- src/schemes/fluid/viscosity.jl | 74 ++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 1a0665a715..acd5618268 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -208,6 +208,25 @@ struct ViscosityAdami{ELTYPE} end end +function adami_viscosity_force(h, pos_diff, distance, grad_kernel, m_a, + m_b, rho_a, rho_b, v_diff, nu_eff; ε=0.01) + # Compute dynamic viscosities: η = ρ * ν_eff + eta_a = rho_a * nu_eff + eta_b = rho_b * nu_eff + # Symmetric averaging (as used in Adami's formulation) + eta_tilde = 2 * (eta_a * eta_b) / (eta_a + eta_b) + + tmp = eta_tilde / (distance^2 + ε * h^2) + + # Compute particle volumes. + volume_a = m_a / rho_a + volume_b = m_b / rho_b + + # The viscous force contribution (per unit mass) is then + visc = (volume_a^2 + volume_b^2) * dot(grad_kernel, pos_diff) * tmp / m_a + return m_b * (visc .* v_diff) +end + @inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, @@ -257,3 +276,58 @@ end @propagate_inbounds function viscous_velocity(v, system, particle) return current_velocity(v, system, particle) end + +struct ViscosityAdamiSGS{ELTYPE} + nu::ELTYPE # Standard (molecular) kinematic viscosity [e.g., 1e-6 m²/s] + C_S::ELTYPE # Smagorinsky constant [e.g., 0.1-0.2] + epsilon::ELTYPE # Epsilon for singularity prevention [e.g., 0.01] +end + +# Convenient constructor with default values for C_S and epsilon. +ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.01) = ViscosityAdamiSGS(nu, C_S, epsilon) + +function compute_SGS_nu(h, v_diff, distance; C_S=0.1, δ=1e-6) + # Estimate the strain rate magnitude |S| + S_mag = norm(v_diff) / (distance + δ) + return (C_S * h)^2 * S_mag +end + +@propagate_inbounds function (viscosity::ViscosityAdamiSGS)(particle_system, + neighbor_system, + v_particle_system, + v_neighbor_system, + particle, neighbor, pos_diff, + distance, sound_speed, m_a, m_b, + rho_a, rho_b, grad_kernel) + # Retrieve the local smoothing length for the reference particle. + h = smoothing_length(particle_system, particle) + + # Retrieve the velocities + v_a = viscous_velocity(v_particle_system, particle_system, particle) + v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) + v_diff = v_a - v_b + + # ------------------------------------------------------------------------------ + # SGS part: Compute the subgrid-scale eddy viscosity. + # ------------------------------------------------------------------------------ + nu_SGS = compute_SGS_nu(h, v_diff, distance; C_S=viscosity.C_S, δ=1e-6) + + # Effective kinematic viscosity is the sum of the standard and SGS parts. + nu_eff = viscosity.nu + nu_SGS + + # ------------------------------------------------------------------------------ + # Common Adami formulation: Compute the viscous acceleration. + # ------------------------------------------------------------------------------ + return adami_viscosity_force(h, pos_diff, distance, grad_kernel, m_a, m_b, + rho_a, rho_b, v_diff, nu_eff; ε=viscosity.epsilon) +end + +# ------------------------------------------------------------------------------ +# A helper to compute the effective kinematic viscosity for diagnostics, if needed. +# ------------------------------------------------------------------------------ +function kinematic_viscosity(particle_system, viscosity::ViscosityAdamiSGS) + # In this diagnostic, we return the baseline (molecular) viscosity. + # The effective viscosity may vary locally. + (; smoothing_length) = particle_system + return viscosity.nu +end From 639be871c9c0e57dcfdb75c7dea120a1f2977763 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Apr 2025 13:29:24 +0200 Subject: [PATCH 02/26] update --- src/TrixiParticles.jl | 3 +- src/schemes/fluid/viscosity.jl | 222 ++++++++++++++++++++++++--------- 2 files changed, 168 insertions(+), 57 deletions(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index e8acf5d54d..8ec0d19929 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -67,7 +67,8 @@ export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel, WendlandC6Kernel, SpikyKernel, Poly6Kernel export StateEquationCole, StateEquationIdealGas -export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris +export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris, ViscosityAdamiSGS, + ViscosityMorrisSGS export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari, DensityDiffusionAntuono export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 03e02ba468..b385d5cfd1 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -212,23 +212,31 @@ struct ViscosityAdami{ELTYPE} end end -function adami_viscosity_force(h, pos_diff, distance, grad_kernel, m_a, - m_b, rho_a, rho_b, v_diff, nu_eff; ε=0.01) - # Compute dynamic viscosities: η = ρ * ν_eff - eta_a = rho_a * nu_eff - eta_b = rho_b * nu_eff - # Symmetric averaging (as used in Adami's formulation) +function adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel, + m_a, + m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon) + eta_a = nu_a * rho_a + eta_b = nu_b * rho_b + eta_tilde = 2 * (eta_a * eta_b) / (eta_a + eta_b) - tmp = eta_tilde / (distance^2 + ε * h^2) + tmp = eta_tilde / (distance^2 + epsilon * smoothing_length_average^2) - # Compute particle volumes. volume_a = m_a / rho_a volume_b = m_b / rho_b - # The viscous force contribution (per unit mass) is then + # This formulation was introduced by Hu and Adams (2006). https://doi.org/10.1016/j.jcp.2005.09.001 + # They argued that the formulation is more flexible because of the possibility to formulate + # different inter-particle averages or to assume different inter-particle distributions. + # Ramachandran (2019) and Adami (2012) use this formulation also for the pressure acceleration. + # + # TODO: Is there a better formulation to discretize the Laplace operator? + # Because when using this formulation for the pressure acceleration, it is not + # energy conserving. + # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394 visc = (volume_a^2 + volume_b^2) * dot(grad_kernel, pos_diff) * tmp / m_a - return m_b * (visc .* v_diff) + + return visc .* v_diff end @inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system, @@ -240,6 +248,7 @@ end smoothing_length_particle = smoothing_length(particle_system, particle) smoothing_length_neighbor = smoothing_length(particle_system, neighbor) + smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2 nu_a = kinematic_viscosity(particle_system, viscosity_model(neighbor_system, particle_system), @@ -252,29 +261,8 @@ end v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) v_diff = v_a - v_b - eta_a = nu_a * rho_a - eta_b = nu_b * rho_b - - eta_tilde = 2 * (eta_a * eta_b) / (eta_a + eta_b) - - smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2 - tmp = eta_tilde / (distance^2 + epsilon * smoothing_length_average^2) - - volume_a = m_a / rho_a - volume_b = m_b / rho_b - - # This formulation was introduced by Hu and Adams (2006). https://doi.org/10.1016/j.jcp.2005.09.001 - # They argued that the formulation is more flexible because of the possibility to formulate - # different inter-particle averages or to assume different inter-particle distributions. - # Ramachandran (2019) and Adami (2012) use this formulation also for the pressure acceleration. - # - # TODO: Is there a better formulation to discretize the Laplace operator? - # Because when using this formulation for the pressure acceleration, it is not - # energy conserving. - # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394 - visc = (volume_a^2 + volume_b^2) * dot(grad_kernel, pos_diff) * tmp / m_a - - return visc .* v_diff + return adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel, + m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon) end function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length, @@ -286,20 +274,56 @@ end return current_velocity(v, system, particle) end +@doc raw""" + ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.01) + +Viscosity model that extends the standard Adami formulation by incorporating a subgrid-scale (SGS) +eddy viscosity via a Smagorinsky-type closure. The effective kinematic viscosity is defined as + +``` +\nu_{\mathrm{eff}} = \nu_{\\mathrm{std}} + \nu_{\\mathrm{SGS}}, +``` + +with + +``` +\nu_{\mathrm{SGS}} = (C_S * h)^2 * |S|, +``` + +and an approximation for the strain rate magnitude given by + +``` +|S| \approx \frac{\|v_a - v_b\|}{\|r_a - r_b\| + \epsilon}, +``` + +where: +- $C_S$ is the Smagorinsky constant (typically 0.1 to 0.2), +- $h$ is the local smoothing length, and + +The effective dynamic viscosities are then computed as +``` +\eta_{a,\mathrm{eff}} = \rho_a\, \nu_{\mathrm{eff}}, +``` +and averaged as +``` +\bar{\eta}_{ab} = \frac{2 \eta_{a,\mathrm{eff}} \eta_{b,\mathrm{eff}}}{\eta_{a,\mathrm{eff}}+\eta_{b,\mathrm{eff}}}. +``` + +This model is appropriate for turbulent flows where unresolved scales contribute additional dissipation. + +# Keywords +- `nu`: Standard kinematic viscosity. +- `C_S`: Smagorinsky constant. +- `epsilon`: Epsilon for singularity prevention [e.g., 0.001] +""" struct ViscosityAdamiSGS{ELTYPE} nu::ELTYPE # Standard (molecular) kinematic viscosity [e.g., 1e-6 m²/s] C_S::ELTYPE # Smagorinsky constant [e.g., 0.1-0.2] - epsilon::ELTYPE # Epsilon for singularity prevention [e.g., 0.01] + epsilon::ELTYPE # Epsilon for singularity prevention [e.g., 0.001] end # Convenient constructor with default values for C_S and epsilon. -ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.01) = ViscosityAdamiSGS(nu, C_S, epsilon) - -function compute_SGS_nu(h, v_diff, distance; C_S=0.1, δ=1e-6) - # Estimate the strain rate magnitude |S| - S_mag = norm(v_diff) / (distance + δ) - return (C_S * h)^2 * S_mag -end +ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, epsilon) @propagate_inbounds function (viscosity::ViscosityAdamiSGS)(particle_system, neighbor_system, @@ -308,10 +332,19 @@ end particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - # Retrieve the local smoothing length for the reference particle. - h = smoothing_length(particle_system, particle) + epsilon = viscosity.epsilon + + smoothing_length_particle = smoothing_length(particle_system, particle) + smoothing_length_neighbor = smoothing_length(particle_system, neighbor) + smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2 + + nu_a = kinematic_viscosity(particle_system, + viscosity_model(neighbor_system, particle_system), + smoothing_length_particle, sound_speed) + nu_b = kinematic_viscosity(neighbor_system, + viscosity_model(particle_system, neighbor_system), + smoothing_length_neighbor, sound_speed) - # Retrieve the velocities v_a = viscous_velocity(v_particle_system, particle_system, particle) v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) v_diff = v_a - v_b @@ -319,24 +352,101 @@ end # ------------------------------------------------------------------------------ # SGS part: Compute the subgrid-scale eddy viscosity. # ------------------------------------------------------------------------------ - nu_SGS = compute_SGS_nu(h, v_diff, distance; C_S=viscosity.C_S, δ=1e-6) + # Estimate the strain rate magnitude |S| (rough approximation) + S_mag = norm(v_diff) / (distance + epsilon) + nu_SGS = (C_S * smoothing_length_average)^2 * S_mag # Effective kinematic viscosity is the sum of the standard and SGS parts. - nu_eff = viscosity.nu + nu_SGS + nu_a = nu_a + nu_SGS + nu_b = nu_b + nu_SGS + + return adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel, + m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon) +end + +function kinematic_viscosity(system, viscosity::ViscosityAdamiSGS, smoothing_length, + sound_speed) + return viscosity.nu +end + +@doc raw""" + ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) + +Subgrid-scale (SGS) viscosity model based on the formulation by [Morris (1997)](@cite Morris1997), +extended with a Smagorinsky-type eddy viscosity term for modeling turbulent flows. + +The acceleration on particle `a` due to viscosity interaction with particle `b` is calculated as: +```math +\frac{d v_a}{dt} = \sum_b m_b \frac{\mu_{a,\mathrm{eff}} + \mu_{b,\mathrm{eff}}}{\rho_a \rho_b} \frac{r_{ab} \cdot \nabla_a W_{ab}}{r_{ab}^2 + \epsilon h_{ab}^2} v_{ab} +where $v_{ab} = v_a - v_b$, $r_{ab} = r_a - r_b$, $r_{ab} = \|r_{ab}\|$, $h_{ab} = (h_a + h_b)/2$, $W_{ab}$ is the smoothing kernel, $\rho$ is density, $m$ is mass, and $\epsilon$ is a regularization parameter. + +The effective dynamic viscosity $\mu_{i,\mathrm{eff}}$ for each particle `i` (a or b) includes both the standard molecular viscosity and an SGS eddy viscosity contribution: +```math +\mu_{i,\mathrm{eff}} = \rho_i \nu_{i,\mathrm{eff}} = \rho_i (\nu_{\mathrm{std}} + \nu_{i,\mathrm{SGS}}) +``` +The standard kinematic viscosity $\nu_{\mathrm{std}}$ is provided by the `nu` parameter. The SGS kinematic viscosity is calculated using the Smagorinsky model: +```math +\nu_{i,\mathrm{SGS}} = (C_S h_{ab})^2 |\bar{S}_{ab}| +``` +where $C_S$ is the Smagorinsky constant. This implementation uses a simplified pairwise approximation for the strain rate tensor magnitude $|\bar{S}_{ab}|$: +```math +|\bar{S}_{ab}| \approx \frac{\|v_a - v_b\|}{\|r_a - r_b\| + \epsilon} +``` + +This model is appropriate for turbulent flows where unresolved scales contribute additional dissipation. + +# Keywords +- `nu`: Standard kinematic viscosity. +- `C_S`: Smagorinsky constant. +- `epsilon`: Epsilon for singularity prevention [e.g., 0.001] +""" +ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, epsilon) + +@propagate_inbounds function (viscosity::ViscosityMorrisSGS)(particle_system, + neighbor_system, + v_particle_system, + v_neighbor_system, + particle, neighbor, pos_diff, + distance, sound_speed, m_a, + m_b, + rho_a, rho_b, grad_kernel) + epsilon_val = viscosity.epsilon + smoothing_length_particle = smoothing_length(particle_system, particle) + smoothing_length_neighbor = smoothing_length(particle_system, neighbor) + smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2 + + nu_a = kinematic_viscosity(particle_system, + viscosity_model(neighbor_system, particle_system), + smoothing_length_particle, sound_speed) + nu_b = kinematic_viscosity(neighbor_system, + viscosity_model(particle_system, neighbor_system), + smoothing_length_neighbor, sound_speed) + + v_a = viscous_velocity(v_particle_system, particle_system, particle) + v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) + v_diff = v_a - v_b # ------------------------------------------------------------------------------ - # Common Adami formulation: Compute the viscous acceleration. + # SGS part: Compute the subgrid-scale eddy viscosity. # ------------------------------------------------------------------------------ - return adami_viscosity_force(h, pos_diff, distance, grad_kernel, m_a, m_b, - rho_a, rho_b, v_diff, nu_eff; ε=viscosity.epsilon) + # Estimate the strain rate magnitude |S| (rough approximation) + S_mag = norm(v_diff) / (distance + epsilon) + nu_SGS = (C_S * smoothing_length_average)^2 * S_mag + + # Effective viscosities include the SGS term. + nu_a_eff = nu_a + nu_SGS + nu_b_eff = nu_b + nu_SGS + + # For the Morris model, dynamic viscosities are: + mu_a = nu_a_eff * rho_a + mu_b = nu_b_eff * rho_b + + force_Morris = (mu_a + mu_b) / (rho_a * rho_b) * (dot(pos_diff, grad_kernel)) / + (distance^2 + epsilon_val * smoothing_length_average^2) * v_diff + return m_b * force_Morris end -# ------------------------------------------------------------------------------ -# A helper to compute the effective kinematic viscosity for diagnostics, if needed. -# ------------------------------------------------------------------------------ -function kinematic_viscosity(particle_system, viscosity::ViscosityAdamiSGS) - # In this diagnostic, we return the baseline (molecular) viscosity. - # The effective viscosity may vary locally. - (; smoothing_length) = particle_system +function kinematic_viscosity(system, viscosity::ViscosityMorrisSGS, smoothing_length, + sound_speed) return viscosity.nu end From 365c1601ac5447854ecda2607657f3b05e10c190 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Apr 2025 13:32:28 +0200 Subject: [PATCH 03/26] fix --- src/schemes/fluid/viscosity.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index b385d5cfd1..1dda8d102b 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -400,6 +400,12 @@ This model is appropriate for turbulent flows where unresolved scales contribute - `C_S`: Smagorinsky constant. - `epsilon`: Epsilon for singularity prevention [e.g., 0.001] """ +struct ViscosityMorrisSGS{ELTYPE} + nu::ELTYPE # Standard (molecular) kinematic viscosity [e.g., 1e-6 m²/s] + C_S::ELTYPE # Smagorinsky constant [e.g., 0.1-0.2] + epsilon::ELTYPE # Epsilon for singularity prevention [e.g., 0.001] +end + ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, epsilon) @propagate_inbounds function (viscosity::ViscosityMorrisSGS)(particle_system, From 1c9e306b90163639eaf909ed15a016a82ead2998 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Apr 2025 13:38:42 +0200 Subject: [PATCH 04/26] remove --- src/schemes/fluid/viscosity.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 1dda8d102b..1ce2c55393 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -270,10 +270,6 @@ function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length return viscosity.nu end -@propagate_inbounds function viscous_velocity(v, system, particle) - return current_velocity(v, system, particle) -end - @doc raw""" ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.01) From 19451d80d55a82ad18508807e075c0ab339870cb Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Apr 2025 13:39:59 +0200 Subject: [PATCH 05/26] fix doc --- src/schemes/fluid/viscosity.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 1ce2c55393..29ab8c9572 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -276,19 +276,19 @@ end Viscosity model that extends the standard Adami formulation by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type closure. The effective kinematic viscosity is defined as -``` +```math \nu_{\mathrm{eff}} = \nu_{\\mathrm{std}} + \nu_{\\mathrm{SGS}}, ``` with -``` +```math \nu_{\mathrm{SGS}} = (C_S * h)^2 * |S|, ``` and an approximation for the strain rate magnitude given by -``` +```math |S| \approx \frac{\|v_a - v_b\|}{\|r_a - r_b\| + \epsilon}, ``` @@ -297,11 +297,11 @@ where: - $h$ is the local smoothing length, and The effective dynamic viscosities are then computed as -``` +```math \eta_{a,\mathrm{eff}} = \rho_a\, \nu_{\mathrm{eff}}, ``` and averaged as -``` +```math \bar{\eta}_{ab} = \frac{2 \eta_{a,\mathrm{eff}} \eta_{b,\mathrm{eff}}}{\eta_{a,\mathrm{eff}}+\eta_{b,\mathrm{eff}}}. ``` From 22d62eaff16a63300664276ed9bc0a51dc3f352d Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Apr 2025 15:12:54 +0200 Subject: [PATCH 06/26] fix --- src/schemes/fluid/viscosity.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 29ab8c9572..1fa5a7eaea 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -270,6 +270,10 @@ function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length return viscosity.nu end +@propagate_inbounds function viscous_velocity(v, system, particle) + return current_velocity(v, system, particle) +end + @doc raw""" ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.01) From 6d7bec4bc725d2ba893cb2e83d4cdc316a89540f Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Apr 2025 15:41:39 +0200 Subject: [PATCH 07/26] fixes --- examples/fluid/dam_break_2d.jl | 11 +++++++++-- src/schemes/fluid/viscosity.jl | 4 ++-- src/visualization/write2vtk.jl | 2 +- 3 files changed, 12 insertions(+), 5 deletions(-) diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index c3fe33566d..e4684da591 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -14,7 +14,7 @@ W = 2 * H # ========================================================================================== # ==== Resolution -fluid_particle_spacing = H / 40 +fluid_particle_spacing = H / 60 # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model boundary_layers = 4 @@ -50,7 +50,9 @@ fluid_density_calculator = ContinuityDensity() viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) # nu = 0.02 * smoothing_length * sound_speed/8 # viscosity = ViscosityMorris(nu=nu) +# viscosity = ViscosityMorrisSGS(nu=nu) # viscosity = ViscosityAdami(nu=nu) +# viscosity = ViscosityAdamiSGS(nu=nu) # Alternatively the density diffusion model by Molteni & Colagrossi can be used, # which will run faster. # density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) @@ -67,12 +69,17 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, # ========================================================================================== # ==== Boundary boundary_density_calculator = AdamiPressureExtrapolation() +viscosity_wall = nothing +# For a no-slip condition the corresponding wall viscosity withouth SGS can be set +#viscosity_wall = ViscosityAdami(nu=nu) +#viscosity_wall = ViscosityMorris(nu=nu) boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, smoothing_kernel, smoothing_length, correction=nothing, - reference_particle_spacing=0) + reference_particle_spacing=0, + viscosity=viscosity_wall) boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, adhesion_coefficient=0.0) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 1fa5a7eaea..b5261b41df 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -354,7 +354,7 @@ ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, eps # ------------------------------------------------------------------------------ # Estimate the strain rate magnitude |S| (rough approximation) S_mag = norm(v_diff) / (distance + epsilon) - nu_SGS = (C_S * smoothing_length_average)^2 * S_mag + nu_SGS = (viscosity.C_S * smoothing_length_average)^2 * S_mag # Effective kinematic viscosity is the sum of the standard and SGS parts. nu_a = nu_a + nu_SGS @@ -437,7 +437,7 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e # ------------------------------------------------------------------------------ # Estimate the strain rate magnitude |S| (rough approximation) S_mag = norm(v_diff) / (distance + epsilon) - nu_SGS = (C_S * smoothing_length_average)^2 * S_mag + nu_SGS = (viscosity.C_S * smoothing_length_average)^2 * S_mag # Effective viscosities include the SGS term. nu_a_eff = nu_a + nu_SGS diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 8e49874a2c..700116ac8d 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -338,7 +338,7 @@ end write2vtk!(vtk, viscosity::Nothing) = vtk -function write2vtk!(vtk, viscosity::Union{ViscosityAdami, ViscosityMorris}) +function write2vtk!(vtk, viscosity::Union{ViscosityAdami, ViscosityMorris, ViscosityAdamiSGS, ViscosityMorrisSGS}) vtk["viscosity_nu"] = viscosity.nu vtk["viscosity_epsilon"] = viscosity.epsilon end From 5f90bf7f8fa3691d63e5fd160376a006482cce88 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Apr 2025 15:42:46 +0200 Subject: [PATCH 08/26] typo --- examples/fluid/dam_break_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index e4684da591..557997b1a1 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -70,7 +70,7 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, # ==== Boundary boundary_density_calculator = AdamiPressureExtrapolation() viscosity_wall = nothing -# For a no-slip condition the corresponding wall viscosity withouth SGS can be set +# For a no-slip condition the corresponding wall viscosity without SGS can be set #viscosity_wall = ViscosityAdami(nu=nu) #viscosity_wall = ViscosityMorris(nu=nu) boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, From 296320276ab825da9f0969cdce8323ae32b57776 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Apr 2025 15:43:05 +0200 Subject: [PATCH 09/26] tpo --- examples/fluid/dam_break_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 557997b1a1..7cf9ebc3ce 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -14,7 +14,7 @@ W = 2 * H # ========================================================================================== # ==== Resolution -fluid_particle_spacing = H / 60 +fluid_particle_spacing = H / 40 # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model boundary_layers = 4 From 82c9e9bca04b5590096ea1718733e50e7eb8f4f8 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 9 Apr 2025 15:48:03 +0200 Subject: [PATCH 10/26] format --- src/visualization/write2vtk.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 700116ac8d..729e021189 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -338,7 +338,9 @@ end write2vtk!(vtk, viscosity::Nothing) = vtk -function write2vtk!(vtk, viscosity::Union{ViscosityAdami, ViscosityMorris, ViscosityAdamiSGS, ViscosityMorrisSGS}) +function write2vtk!(vtk, + viscosity::Union{ViscosityAdami, ViscosityMorris, ViscosityAdamiSGS, + ViscosityMorrisSGS}) vtk["viscosity_nu"] = viscosity.nu vtk["viscosity_epsilon"] = viscosity.epsilon end From 5ff0c7437da3f4a05f05524f1e4a33ad1bb6c0d5 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 6 May 2025 11:57:32 +0200 Subject: [PATCH 11/26] add explanation --- src/schemes/fluid/viscosity.jl | 39 +++++++++++++++++++++++++++++++--- 1 file changed, 36 insertions(+), 3 deletions(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 1416f4d477..1ceab53b94 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -316,7 +316,24 @@ ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, eps # ------------------------------------------------------------------------------ # SGS part: Compute the subgrid-scale eddy viscosity. # ------------------------------------------------------------------------------ - # Estimate the strain rate magnitude |S| (rough approximation) + # In classical LES the Smagorinsky model defines: + # ν_SGS = (C_S Δ)^2 |S|, + # where |S| is the norm of the strain-rate tensor Sᵢⱼ = ½(∂ᵢvⱼ+∂ⱼvᵢ). + # + # In SPH, one could compute ∂ᵢvⱼ via kernel gradients, but this is costly. + # A common low-order surrogate is to approximate the strain‐rate magnitude by a + # finite-difference along each particle pair: + # + # |S| ≈ ‖vₐ–v_b‖ / (‖rₐ–r_b‖ + δ), + # + # where δ regularizes the denominator to avoid singularities when particles are very close. + # + # This yields: + # S_mag = norm(v_diff) / (distance + ε) + # + # and then the Smagorinsky eddy viscosity: + # ν_SGS = (C_S * h̄)^2 * S_mag. + # S_mag = norm(v_diff) / (distance + epsilon) nu_SGS = (viscosity.C_S * smoothing_length_average)^2 * S_mag @@ -399,8 +416,24 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e # ------------------------------------------------------------------------------ # SGS part: Compute the subgrid-scale eddy viscosity. # ------------------------------------------------------------------------------ - # Estimate the strain rate magnitude |S| (rough approximation) - S_mag = norm(v_diff) / (distance + epsilon) + # In classical LES the Smagorinsky model defines: + # ν_SGS = (C_S Δ)^2 |S|, + # where |S| is the norm of the strain-rate tensor Sᵢⱼ = ½(∂ᵢvⱼ+∂ⱼvᵢ). + # + # In SPH, one could compute ∂ᵢvⱼ via kernel gradients, but this is costly. + # A common low-order surrogate is to approximate the strain‐rate magnitude by a + # finite-difference along each particle pair: + # + # |S| ≈ ‖vₐ–v_b‖ / (‖rₐ–r_b‖ + δ), + # + # where δ regularizes the denominator to avoid singularities when particles are very close. + # + # This yields: + # S_mag = norm(v_diff) / (distance + ε) + # + # and then the Smagorinsky eddy viscosity: + # ν_SGS = (C_S * h̄)^2 * S_mag. + # S_mag = norm(v_diff) / (distance + epsilon) nu_SGS = (viscosity.C_S * smoothing_length_average)^2 * S_mag # Effective viscosities include the SGS term. From 77bf5ef8931e673ee5ee30edb68324794e7943f2 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 6 May 2025 13:04:29 +0200 Subject: [PATCH 12/26] adjust tests --- src/schemes/fluid/viscosity.jl | 10 ++-- test/schemes/fluid/viscosity.jl | 84 ++++++++++++++++++++++++++++----- 2 files changed, 78 insertions(+), 16 deletions(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 1ceab53b94..30c194c011 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -395,9 +395,9 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e v_neighbor_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, - m_b, - rho_a, rho_b, grad_kernel) - epsilon_val = viscosity.epsilon + m_b, rho_a, rho_b, grad_kernel) + epsilon = viscosity.epsilon + smoothing_length_particle = smoothing_length(particle_system, particle) smoothing_length_neighbor = smoothing_length(particle_system, neighbor) smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2 @@ -433,7 +433,7 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e # # and then the Smagorinsky eddy viscosity: # ν_SGS = (C_S * h̄)^2 * S_mag. - # S_mag = norm(v_diff) / (distance + epsilon) + S_mag = norm(v_diff) / (distance + epsilon) nu_SGS = (viscosity.C_S * smoothing_length_average)^2 * S_mag # Effective viscosities include the SGS term. @@ -445,7 +445,7 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e mu_b = nu_b_eff * rho_b force_Morris = (mu_a + mu_b) / (rho_a * rho_b) * (dot(pos_diff, grad_kernel)) / - (distance^2 + epsilon_val * smoothing_length_average^2) * v_diff + (distance^2 + epsilon * smoothing_length_average^2) * v_diff return m_b * force_Morris end diff --git a/test/schemes/fluid/viscosity.jl b/test/schemes/fluid/viscosity.jl index f11f4b05b7..66b770763f 100644 --- a/test/schemes/fluid/viscosity.jl +++ b/test/schemes/fluid/viscosity.jl @@ -1,3 +1,4 @@ +include("../../test_util.jl") @testset verbose=true "Viscosity" begin particle_spacing = 0.2 smoothing_length = 1.2 * particle_spacing @@ -9,8 +10,8 @@ fluid = rectangular_patch(particle_spacing, (3, 3), seed=1) - v_diff = [0.1, -0.75] - pos_diff = [-0.5 * smoothing_length, 0.75 * smoothing_length] + v_diff = [0.3, -1.0] + pos_diff = [-0.25 * smoothing_length, 0.375 * smoothing_length] distance = norm(pos_diff) rho_a = rho_b = rho_mean = 1000.0 @@ -29,8 +30,8 @@ rho_mean, rho_a, rho_b, smoothing_length, grad_kernel, 0.0, 0.0) - @test isapprox(dv[1], -0.007073849138494646, atol=6e-15) - @test isapprox(dv[2], 0.01061077370774197, atol=6e-15) + @test isapprox(dv[1], -0.02049217623299368, atol=6e-15) + @test isapprox(dv[2], 0.03073826434949052, atol=6e-15) end @testset verbose=true "`ViscosityMorris`" begin nu = 7e-3 @@ -41,13 +42,22 @@ grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff, distance, 1) + v = fluid.velocity - dv = viscosity(sound_speed, v_diff, pos_diff, distance, - rho_mean, rho_a, rho_b, smoothing_length, - grad_kernel, nu, nu) + m_a = 0.01 + m_b = 0.01 - @test isapprox(dv[1], -0.00018421886647594437, atol=6e-15) - @test isapprox(dv[2], 0.0013816414985695826, atol=6e-15) + v[1, 1] = v_diff[1] + v[2, 1] = v_diff[2] + v[1, 2] = 0.0 + v[2, 2] = 0.0 + + dv = viscosity(system_wcsph, system_wcsph, + v, v, 1, 2, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + + @test isapprox(dv[1], -1.0895602048035404e-5, atol=6e-15) + @test isapprox(dv[2], 3.631867349345135e-5, atol=6e-15) end @testset verbose=true "`ViscosityAdami`" begin viscosity = ViscosityAdami(nu=7e-3) @@ -71,7 +81,59 @@ v, v, 1, 2, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - @test isapprox(dv[1], -1.8421886647594435e-6, atol=6e-15) - @test isapprox(dv[2], 1.3816414985695826e-5, atol=6e-15) + @test isapprox(dv[1], -1.089560204803541e-5, atol=6e-15) + @test isapprox(dv[2], 3.6318673493451364e-5, atol=6e-15) + end + @testset verbose=true "`ViscosityMorrisSGS`" begin + nu = 7e-3 + viscosity = ViscosityMorrisSGS(nu=nu) + system_wcsph = WeaklyCompressibleSPHSystem(fluid, ContinuityDensity(), + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity) + + grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff, + distance, 1) + + v = fluid.velocity + + m_a = 0.01 + m_b = 0.01 + + v[1, 1] = v_diff[1] + v[2, 1] = v_diff[2] + v[1, 2] = 0.0 + v[2, 2] = 0.0 + + dv = viscosity(system_wcsph, system_wcsph, + v, v, 1, 2, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + + @test isapprox(dv[1], -2.032835697804103e-5, atol=6e-15) + @test isapprox(dv[2], 6.776118992680343e-5, atol=6e-15) + end + @testset verbose=true "`ViscosityAdamiSGS`" begin + viscosity = ViscosityAdamiSGS(nu=7e-3) + system_wcsph = WeaklyCompressibleSPHSystem(fluid, ContinuityDensity(), + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity) + + grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff, + distance, 1) + v = fluid.velocity + + m_a = 0.01 + m_b = 0.01 + + v[1, 1] = v_diff[1] + v[2, 1] = v_diff[2] + v[1, 2] = 0.0 + v[2, 2] = 0.0 + + dv = viscosity(system_wcsph, system_wcsph, + v, v, 1, 2, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + + @test isapprox(dv[1], -2.0328356978041036e-5, atol=6e-15) + @test isapprox(dv[2], 6.776118992680346e-5, atol=6e-15) end end From bad1954f23e22119bd4013b30fb983d0121baf52 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 13 May 2025 16:50:12 +0200 Subject: [PATCH 13/26] add to validation setup for wcsph --- examples/fluid/taylor_green_vortex_2d.jl | 4 ++-- .../validation_taylor_green_vortex_2d.jl | 19 +++++++++++++------ 2 files changed, 15 insertions(+), 8 deletions(-) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index c1f0f72fd4..4fb430a277 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -79,7 +79,7 @@ if wcsph state_equation, smoothing_kernel, pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure, smoothing_length, - viscosity=ViscosityAdami(; nu), + viscosity=ViscosityAdamiSGS(; nu, C_S=0.2), transport_velocity=TransportVelocityAdami(background_pressure)) else density_calculator = SummationDensity() @@ -100,7 +100,7 @@ ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) -saving_callback = SolutionSavingCallback(dt=0.02) +saving_callback = SolutionSavingCallback(dt=0.02, prefix="SGS") pp_callback = nothing diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl index a30f26f98c..70b47d4577 100644 --- a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -2,7 +2,8 @@ using TrixiParticles # ========================================================================================== # ==== Resolution -particle_spacings = [0.02, 0.01, 0.005] +# particle_spacings = [0.02, 0.01, 0.005] +particle_spacings = [0.02] # ========================================================================================== # ==== Experiment Setup @@ -12,10 +13,13 @@ reynolds_number = 100.0 density_calculators = [ContinuityDensity(), SummationDensity()] perturb_coordinates = [false, true] -function compute_l1v_error(v, u, t, system) +function compute_l1v_error(system, v_ode, u_ode, semi, t) v_analytical_avg = 0.0 v_avg = 0.0 + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + for particle in TrixiParticles.eachparticle(system) position = TrixiParticles.current_coords(u, system, particle) @@ -32,11 +36,14 @@ function compute_l1v_error(v, u, t, system) return v_avg /= v_analytical_avg end -function compute_l1p_error(v, u, t, system) +function compute_l1p_error(system, v_ode, u_ode, semi, t) p_max_exact = 0.0 L1p = 0.0 + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + for particle in TrixiParticles.eachparticle(system) position = TrixiParticles.current_coords(u, system, particle) @@ -57,8 +64,8 @@ end # The pressure plotted in the paper is the difference of the local pressure minus # the average of the pressure of all particles. -function diff_p_loc_p_avg(v, u, t, system) - p_avg_tot = avg_pressure(v, u, t, system) +function diff_p_loc_p_avg(system, v, u, semi, t) + p_avg_tot = avg_pressure(system, v, u, semi, t) return v[end, :] .- p_avg_tot end @@ -74,7 +81,7 @@ for density_calculator in density_calculators, perturbation in perturb_coordinat output_directory = joinpath("out_tgv", name_density_calculator * name_perturbation, wcsph ? "wcsph" : "edac", "validation_run_taylor_green_vortex_2d_nparticles_$(n_particles_xy)x$(n_particles_xy)") - saving_callback = SolutionSavingCallback(dt=0.02, + saving_callback = SolutionSavingCallback(dt=0.1, output_directory=output_directory, p_avg=diff_p_loc_p_avg) From bea7af621bace532e3dd523f32174bff3705cdff Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 14 May 2025 01:29:28 +0200 Subject: [PATCH 14/26] review fixes --- examples/fluid/dam_break_2d.jl | 17 +++++++++++++---- examples/fluid/taylor_green_vortex_2d.jl | 2 +- src/schemes/fluid/viscosity.jl | 16 ++++++++-------- .../validation_taylor_green_vortex_2d.jl | 7 +++++++ 4 files changed, 29 insertions(+), 13 deletions(-) diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 159eddaa6d..b1298afea9 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -47,12 +47,21 @@ smoothing_length = 1.75 * fluid_particle_spacing smoothing_kernel = WendlandC2Kernel{2}() fluid_density_calculator = ContinuityDensity() -viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) -# nu = 0.02 * smoothing_length * sound_speed/8 -# viscosity = ViscosityMorris(nu=nu) -# viscosity = ViscosityMorrisSGS(nu=nu) +alpha = 0.02 +viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) +# A typical formula to convert Artificial viscosity to a +# kinematic viscosity is provided by Monaghan as +# nu = alpha * smoothing_length * sound_speed/8 + +# Alternatively a kinematic viscosity for water can be set +# nu = 1.0e-6 + +# This allows the use of a physical viscosity model like: # viscosity = ViscosityAdami(nu=nu) +# or with additional diffusion through the Smagorinsky model # viscosity = ViscosityAdamiSGS(nu=nu) +# For more details see the documentation [Viscosity model overview](@ref viscosity_sph). + # Alternatively the density diffusion model by Molteni & Colagrossi can be used, # which will run faster. # density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 4fb430a277..0f36054de6 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -79,7 +79,7 @@ if wcsph state_equation, smoothing_kernel, pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure, smoothing_length, - viscosity=ViscosityAdamiSGS(; nu, C_S=0.2), + viscosity=ViscosityAdami(; nu), transport_velocity=TransportVelocityAdami(background_pressure)) else density_calculator = SummationDensity() diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 30c194c011..15590de45e 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -177,8 +177,7 @@ struct ViscosityAdami{ELTYPE} end function adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel, - m_a, - m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon) + m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon) eta_a = nu_a * rho_a eta_b = nu_b * rho_b @@ -241,8 +240,9 @@ end @doc raw""" ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.01) -Viscosity model that extends the standard Adami formulation by incorporating a subgrid-scale (SGS) -eddy viscosity via a Smagorinsky-type closure. The effective kinematic viscosity is defined as +Viscosity model that extends the standard [Adami formulation](@ref ViscosityAdami) +by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type closure. +The effective kinematic viscosity is defined as ```math \nu_{\mathrm{eff}} = \nu_{\\mathrm{std}} + \nu_{\\mathrm{SGS}}, @@ -278,10 +278,10 @@ This model is appropriate for turbulent flows where unresolved scales contribute # Keywords - `nu`: Standard kinematic viscosity. - `C_S`: Smagorinsky constant. -- `epsilon`: Epsilon for singularity prevention [e.g., 0.001] +- `epsilon=0.01`: Parameter to prevent singularities """ struct ViscosityAdamiSGS{ELTYPE} - nu::ELTYPE # Standard (molecular) kinematic viscosity [e.g., 1e-6 m²/s] + nu::ELTYPE # kinematic viscosity [e.g., 1e-6 m²/s] C_S::ELTYPE # Smagorinsky constant [e.g., 0.1-0.2] epsilon::ELTYPE # Epsilon for singularity prevention [e.g., 0.001] end @@ -379,10 +379,10 @@ This model is appropriate for turbulent flows where unresolved scales contribute # Keywords - `nu`: Standard kinematic viscosity. - `C_S`: Smagorinsky constant. -- `epsilon`: Epsilon for singularity prevention [e.g., 0.001] +- `epsilon=0.01`: Parameter to prevent singularities """ struct ViscosityMorrisSGS{ELTYPE} - nu::ELTYPE # Standard (molecular) kinematic viscosity [e.g., 1e-6 m²/s] + nu::ELTYPE # kinematic viscosity [e.g., 1e-6 m²/s] C_S::ELTYPE # Smagorinsky constant [e.g., 0.1-0.2] epsilon::ELTYPE # Epsilon for singularity prevention [e.g., 0.001] end diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl index 70b47d4577..a6a73f8477 100644 --- a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -2,7 +2,14 @@ using TrixiParticles # ========================================================================================== # ==== Resolution +# P. Ramachandran, K. Puri +# "Entropically damped artificial compressibility for SPH". +# In: Computers and Fluids, Volume 179 (2019), pages 579-594. +# https://doi.org/10.1016/j.compfluid.2018.11.023 +# The following resolutions are used in the paper: # particle_spacings = [0.02, 0.01, 0.005] + +# To save time only the configuration with 50x50 particles is used here. particle_spacings = [0.02] # ========================================================================================== From fa6df699c233b9a6d50f0a512c068eafa0d13bb1 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 14 May 2025 09:22:56 +0200 Subject: [PATCH 15/26] Update taylor_green_vortex_2d.jl --- examples/fluid/taylor_green_vortex_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index 0f36054de6..c1f0f72fd4 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -100,7 +100,7 @@ ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) -saving_callback = SolutionSavingCallback(dt=0.02, prefix="SGS") +saving_callback = SolutionSavingCallback(dt=0.02) pp_callback = nothing From 0458976eb3462ea27ee1b98a73d459daffa05698 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 14 May 2025 09:24:15 +0200 Subject: [PATCH 16/26] Update validation_taylor_green_vortex_2d.jl --- .../taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl index a6a73f8477..56b7ee5646 100644 --- a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -9,7 +9,7 @@ using TrixiParticles # The following resolutions are used in the paper: # particle_spacings = [0.02, 0.01, 0.005] -# To save time only the configuration with 50x50 particles is used here. +# To save time only the configuration with 50x50 particles are used here. particle_spacings = [0.02] # ========================================================================================== From de606e126153c844a4697c3d722afc53e5c50ae2 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 14 May 2025 09:26:42 +0200 Subject: [PATCH 17/26] Update dam_break_2d.jl --- examples/fluid/dam_break_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index b1298afea9..1f198b691d 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -58,7 +58,7 @@ viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) # This allows the use of a physical viscosity model like: # viscosity = ViscosityAdami(nu=nu) -# or with additional diffusion through the Smagorinsky model +# or with additional dissipation through a Smagorinsky model # viscosity = ViscosityAdamiSGS(nu=nu) # For more details see the documentation [Viscosity model overview](@ref viscosity_sph). From c45fdb6a98e0b2711b855845265829f0397a124d Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 14 May 2025 10:54:55 +0000 Subject: [PATCH 18/26] fix and add citations --- docs/src/refs.bib | 24 ++++++++++++++++++- examples/fluid/taylor_green_vortex_2d.jl | 2 +- src/schemes/fluid/viscosity.jl | 8 +++---- .../validation_taylor_green_vortex_2d.jl | 2 +- 4 files changed, 29 insertions(+), 7 deletions(-) diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 67f75b6925..9ed1f0bd96 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -704,7 +704,7 @@ @article{MeloEspinosa2014 journal = {Energy Conversion and Management}, volume = {84}, pages = {50--60}, - doi = {https://doi.org/10.1016/j.enconman.2014.08.004}, + doi = {10.1016/j.enconman.2014.08.004}, year = {2014} } @@ -723,4 +723,26 @@ @book{Poling2001 year = {2001}, publisher = {McGraw-Hill}, address = {New York} +} + +@article{Smagorinsky1963, + author = {Smagorinsky, Joseph}, + title = {General Circulation Experiments with the Primitive Equations. I. The Basic Experiment}, + journal = {Monthly Weather Review}, + volume = {91}, + number = {3}, + pages = {99--164}, + year = {1963}, + doi = {10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2} +} + +@inproceedings{Lilly1967, + author = {Lilly, Douglas K.}, + title = {The Representation of Small‑Scale Turbulence in Numerical Simulation Experiments}, + booktitle = {Proceedings of the {IBM} Scientific Computing Symposium on Environmental Sciences}, + editor = {Goldstine, Herman H.}, + address = {Yorktown Heights, New York}, + pages = {195--210}, + year = {1967}, + doi = {10.5065/D62R3PMM} } \ No newline at end of file diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index c1f0f72fd4..75e6387184 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -3,7 +3,7 @@ # P. Ramachandran, K. Puri # "Entropically damped artificial compressibility for SPH". # In: Computers and Fluids, Volume 179 (2019), pages 579-594. -# https://doi.org/10.1016/j.compfluid.2018.11.023 +# https://doi.org/10.1016/j.compfluid.2018.11.023 using TrixiParticles using OrdinaryDiffEq diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 15590de45e..d9be50dc09 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -241,7 +241,7 @@ end ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.01) Viscosity model that extends the standard [Adami formulation](@ref ViscosityAdami) -by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type closure. +by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type [Smagorinsky (1963)](@cite Smagorinsky1963) closure. The effective kinematic viscosity is defined as ```math @@ -316,7 +316,7 @@ ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, eps # ------------------------------------------------------------------------------ # SGS part: Compute the subgrid-scale eddy viscosity. # ------------------------------------------------------------------------------ - # In classical LES the Smagorinsky model defines: + # In classical LES [Lilly (1967)](@cite Lilly1967) the Smagorinsky model defines: # ν_SGS = (C_S Δ)^2 |S|, # where |S| is the norm of the strain-rate tensor Sᵢⱼ = ½(∂ᵢvⱼ+∂ⱼvᵢ). # @@ -354,7 +354,7 @@ end ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) Subgrid-scale (SGS) viscosity model based on the formulation by [Morris (1997)](@cite Morris1997), -extended with a Smagorinsky-type eddy viscosity term for modeling turbulent flows. +extended with a Smagorinsky-type eddy viscosity term [Smagorinsky (1963)](@cite Smagorinsky1963) for modeling turbulent flows. The acceleration on particle `a` due to viscosity interaction with particle `b` is calculated as: ```math @@ -416,7 +416,7 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e # ------------------------------------------------------------------------------ # SGS part: Compute the subgrid-scale eddy viscosity. # ------------------------------------------------------------------------------ - # In classical LES the Smagorinsky model defines: + # In classical LES [Lilly (1967)](@cite Lilly1967) the Smagorinsky model defines: # ν_SGS = (C_S Δ)^2 |S|, # where |S| is the norm of the strain-rate tensor Sᵢⱼ = ½(∂ᵢvⱼ+∂ⱼvᵢ). # diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl index 56b7ee5646..bb7d2d3773 100644 --- a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -5,7 +5,7 @@ using TrixiParticles # P. Ramachandran, K. Puri # "Entropically damped artificial compressibility for SPH". # In: Computers and Fluids, Volume 179 (2019), pages 579-594. -# https://doi.org/10.1016/j.compfluid.2018.11.023 +# https://doi.org/10.1016/j.compfluid.2018.11.023 # The following resolutions are used in the paper: # particle_spacings = [0.02, 0.01, 0.005] From 41938afd9fc0e8de95d5eb1707079eb227ececf5 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 14 May 2025 10:58:28 +0000 Subject: [PATCH 19/26] add news --- NEWS.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/NEWS.md b/NEWS.md index cb91788fef..a688a03ac8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,12 @@ TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1) used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Version 0.3.1 + +### Features + +- New Viscosity model 'ViscosityMorrisSGS' and 'ViscosityAdamiSGS' were added which use simplified Smagorinsky type SGS (#753). + ## Version 0.3 ### API Changes From af656b60f8bee4f290b71ed1f7b1282d6a165832 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 21 May 2025 11:33:44 +0200 Subject: [PATCH 20/26] implement suggestions --- NEWS.md | 2 +- examples/fluid/dam_break_2d.jl | 14 ++--- examples/fluid/dam_break_2phase_2d.jl | 2 +- examples/fluid/dam_break_oil_film_2d.jl | 2 +- src/schemes/fluid/viscosity.jl | 77 +++++++++++-------------- test/examples/examples_fluid.jl | 16 ++--- test/schemes/fluid/viscosity.jl | 1 - 7 files changed, 53 insertions(+), 61 deletions(-) diff --git a/NEWS.md b/NEWS.md index abfe525434..fc9b832e30 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,7 +8,7 @@ used in the Julia ecosystem. Notable changes will be documented in this file for ### Features -- New Viscosity model 'ViscosityMorrisSGS' and 'ViscosityAdamiSGS' were added which use simplified Smagorinsky type SGS (#753). +- New viscosity models `ViscosityMorrisSGS` and `ViscosityAdamiSGS` were added, which use a simplified Smagorinsky-type SGS (#753). - With all CPU backends, a new array type is used for the integration array, which defines broadcasting to be multithreaded, leading to speedups of up to 5x with large thread counts diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 1f198b691d..950f8c8060 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -48,7 +48,7 @@ smoothing_kernel = WendlandC2Kernel{2}() fluid_density_calculator = ContinuityDensity() alpha = 0.02 -viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) +viscosity_fluid = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) # A typical formula to convert Artificial viscosity to a # kinematic viscosity is provided by Monaghan as # nu = alpha * smoothing_length * sound_speed/8 @@ -57,10 +57,10 @@ viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) # nu = 1.0e-6 # This allows the use of a physical viscosity model like: -# viscosity = ViscosityAdami(nu=nu) +# viscosity_fluid = ViscosityAdami(nu=nu) # or with additional dissipation through a Smagorinsky model -# viscosity = ViscosityAdamiSGS(nu=nu) -# For more details see the documentation [Viscosity model overview](@ref viscosity_sph). +# viscosity_fluid = ViscosityAdamiSGS(nu=nu) +# For more details see the documentation "Viscosity model overview". # Alternatively the density diffusion model by Molteni & Colagrossi can be used, # which will run faster. @@ -69,7 +69,7 @@ density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, - smoothing_length, viscosity=viscosity, + smoothing_length, viscosity=viscosity_fluid, density_diffusion=density_diffusion, acceleration=(0.0, -gravity), correction=nothing, surface_tension=nothing, @@ -80,8 +80,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, boundary_density_calculator = AdamiPressureExtrapolation() viscosity_wall = nothing # For a no-slip condition the corresponding wall viscosity without SGS can be set -#viscosity_wall = ViscosityAdami(nu=nu) -#viscosity_wall = ViscosityMorris(nu=nu) +# viscosity_wall = ViscosityAdami(nu=nu) +# viscosity_wall = ViscosityMorris(nu=nu) boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, diff --git a/examples/fluid/dam_break_2phase_2d.jl b/examples/fluid/dam_break_2phase_2d.jl index de99a5f397..247462b52b 100644 --- a/examples/fluid/dam_break_2phase_2d.jl +++ b/examples/fluid/dam_break_2phase_2d.jl @@ -37,7 +37,7 @@ water_viscosity = ViscosityMorris(nu=nu_sim_water) trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), sol=nothing, fluid_particle_spacing=fluid_particle_spacing, - viscosity=water_viscosity, smoothing_length=smoothing_length, + viscosity_fluid=water_viscosity, smoothing_length=smoothing_length, gravity=gravity, tspan=tspan, density_diffusion=nothing, sound_speed=sound_speed, exponent=7, tank_size=(floor(5.366 * H / fluid_particle_spacing) * fluid_particle_spacing, diff --git a/examples/fluid/dam_break_oil_film_2d.jl b/examples/fluid/dam_break_oil_film_2d.jl index a9143b1128..5d71c40146 100644 --- a/examples/fluid/dam_break_oil_film_2d.jl +++ b/examples/fluid/dam_break_oil_film_2d.jl @@ -30,7 +30,7 @@ oil_viscosity = ViscosityMorris(nu=nu_sim_oil) # TODO: broken if both systems use surface tension trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), sol=nothing, fluid_particle_spacing=fluid_particle_spacing, - viscosity=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length, + viscosity_fluid=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length, gravity=gravity, density_diffusion=nothing, sound_speed=sound_speed, prefix="", reference_particle_spacing=fluid_particle_spacing) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index d9be50dc09..1d8dbc3c01 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -245,32 +245,32 @@ by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type [Sm The effective kinematic viscosity is defined as ```math -\nu_{\mathrm{eff}} = \nu_{\\mathrm{std}} + \nu_{\\mathrm{SGS}}, +\nu_{\text{eff}} = \nu_{\text{std}} + \nu_{\text{SGS}}, ``` with ```math -\nu_{\mathrm{SGS}} = (C_S * h)^2 * |S|, +\nu_{\text{SGS}} = (C_S h)^2 |S|, ``` and an approximation for the strain rate magnitude given by ```math -|S| \approx \frac{\|v_a - v_b\|}{\|r_a - r_b\| + \epsilon}, +|S| \approx \frac{\|v_ab\|}{\|r_ab\| + \epsilon}, ``` where: -- $C_S$ is the Smagorinsky constant (typically 0.1 to 0.2), -- $h$ is the local smoothing length, and +- ``C_S`` is the Smagorinsky constant (typically 0.1 to 0.2), +- ``h`` is the local smoothing length. The effective dynamic viscosities are then computed as ```math -\eta_{a,\mathrm{eff}} = \rho_a\, \nu_{\mathrm{eff}}, +\eta_{a,\text{eff}} = \rho_a\, \nu_{\text{eff}}, ``` and averaged as ```math -\bar{\eta}_{ab} = \frac{2 \eta_{a,\mathrm{eff}} \eta_{b,\mathrm{eff}}}{\eta_{a,\mathrm{eff}}+\eta_{b,\mathrm{eff}}}. +\bar{\eta}_{ab} = \frac{2 \eta_{a,\text{eff}} \eta_{b,\text{eff}}}{\eta_{a,\text{eff}}+\eta_{b,\text{eff}}}. ``` This model is appropriate for turbulent flows where unresolved scales contribute additional dissipation. @@ -281,12 +281,11 @@ This model is appropriate for turbulent flows where unresolved scales contribute - `epsilon=0.01`: Parameter to prevent singularities """ struct ViscosityAdamiSGS{ELTYPE} - nu::ELTYPE # kinematic viscosity [e.g., 1e-6 m²/s] - C_S::ELTYPE # Smagorinsky constant [e.g., 0.1-0.2] - epsilon::ELTYPE # Epsilon for singularity prevention [e.g., 0.001] + nu :: ELTYPE # kinematic viscosity [e.g., 1e-6 m²/s] + C_S :: ELTYPE # Smagorinsky constant [e.g., 0.1-0.2] + epsilon :: ELTYPE # Epsilon for singularity prevention [e.g., 0.001] end -# Convenient constructor with default values for C_S and epsilon. ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, epsilon) @propagate_inbounds function (viscosity::ViscosityAdamiSGS)(particle_system, @@ -322,9 +321,9 @@ ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, eps # # In SPH, one could compute ∂ᵢvⱼ via kernel gradients, but this is costly. # A common low-order surrogate is to approximate the strain‐rate magnitude by a - # finite-difference along each particle pair: + # finite difference along each particle pair: # - # |S| ≈ ‖vₐ–v_b‖ / (‖rₐ–r_b‖ + δ), + # |S| ≈ ‖v_ab‖ / (‖r_ab‖ + δ), # # where δ regularizes the denominator to avoid singularities when particles are very close. # @@ -354,24 +353,36 @@ end ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) Subgrid-scale (SGS) viscosity model based on the formulation by [Morris (1997)](@cite Morris1997), -extended with a Smagorinsky-type eddy viscosity term [Smagorinsky (1963)](@cite Smagorinsky1963) for modeling turbulent flows. +by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type [Smagorinsky (1963)](@cite Smagorinsky1963) closure. +The effective kinematic viscosity is defined as + +```math +\nu_{\text{eff}} = \nu_{\text{std}} + \nu_{\text{SGS}}, +``` + +with -The acceleration on particle `a` due to viscosity interaction with particle `b` is calculated as: ```math -\frac{d v_a}{dt} = \sum_b m_b \frac{\mu_{a,\mathrm{eff}} + \mu_{b,\mathrm{eff}}}{\rho_a \rho_b} \frac{r_{ab} \cdot \nabla_a W_{ab}}{r_{ab}^2 + \epsilon h_{ab}^2} v_{ab} -where $v_{ab} = v_a - v_b$, $r_{ab} = r_a - r_b$, $r_{ab} = \|r_{ab}\|$, $h_{ab} = (h_a + h_b)/2$, $W_{ab}$ is the smoothing kernel, $\rho$ is density, $m$ is mass, and $\epsilon$ is a regularization parameter. +\nu_{\text{SGS}} = (C_S h)^2 |S|, +``` + +and an approximation for the strain rate magnitude given by -The effective dynamic viscosity $\mu_{i,\mathrm{eff}}$ for each particle `i` (a or b) includes both the standard molecular viscosity and an SGS eddy viscosity contribution: ```math -\mu_{i,\mathrm{eff}} = \rho_i \nu_{i,\mathrm{eff}} = \rho_i (\nu_{\mathrm{std}} + \nu_{i,\mathrm{SGS}}) +|S| \approx \frac{\|v_ab\|}{\|r_ab\| + \epsilon}, ``` -The standard kinematic viscosity $\nu_{\mathrm{std}}$ is provided by the `nu` parameter. The SGS kinematic viscosity is calculated using the Smagorinsky model: + +where: +- ``C_S`` is the Smagorinsky constant (typically 0.1 to 0.2), +- ``h`` is the local smoothing length. + +The effective dynamic viscosities are then computed as ```math -\nu_{i,\mathrm{SGS}} = (C_S h_{ab})^2 |\bar{S}_{ab}| +\eta_{a,\text{eff}} = \rho_a\, \nu_{\text{eff}}, ``` -where $C_S$ is the Smagorinsky constant. This implementation uses a simplified pairwise approximation for the strain rate tensor magnitude $|\bar{S}_{ab}|$: +and averaged as ```math -|\bar{S}_{ab}| \approx \frac{\|v_a - v_b\|}{\|r_a - r_b\| + \epsilon} +\bar{\eta}_{ab} = \frac{2 \eta_{a,\text{eff}} \eta_{b,\text{eff}}}{\eta_{a,\text{eff}}+\eta_{b,\text{eff}}}. ``` This model is appropriate for turbulent flows where unresolved scales contribute additional dissipation. @@ -413,26 +424,8 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) v_diff = v_a - v_b - # ------------------------------------------------------------------------------ # SGS part: Compute the subgrid-scale eddy viscosity. - # ------------------------------------------------------------------------------ - # In classical LES [Lilly (1967)](@cite Lilly1967) the Smagorinsky model defines: - # ν_SGS = (C_S Δ)^2 |S|, - # where |S| is the norm of the strain-rate tensor Sᵢⱼ = ½(∂ᵢvⱼ+∂ⱼvᵢ). - # - # In SPH, one could compute ∂ᵢvⱼ via kernel gradients, but this is costly. - # A common low-order surrogate is to approximate the strain‐rate magnitude by a - # finite-difference along each particle pair: - # - # |S| ≈ ‖vₐ–v_b‖ / (‖rₐ–r_b‖ + δ), - # - # where δ regularizes the denominator to avoid singularities when particles are very close. - # - # This yields: - # S_mag = norm(v_diff) / (distance + ε) - # - # and then the Smagorinsky eddy viscosity: - # ν_SGS = (C_S * h̄)^2 * S_mag. + # See comments above for `ViscosityAdamiSGS`. S_mag = norm(v_diff) / (distance + epsilon) nu_SGS = (viscosity.C_S * smoothing_length_average)^2 * S_mag diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index 32b2ccbde1..73338a82e4 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -25,18 +25,18 @@ clip_negative_pressure=true), "WCSPH with ViscosityAdami" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityAdami(nu=0.0015),), + viscosity_fluid=ViscosityAdami(nu=0.0015),), "WCSPH with ViscosityMorris" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityMorris(nu=0.0015),), + viscosity_fluid=ViscosityMorris(nu=0.0015),), "WCSPH with ViscosityAdami and SummationDensity" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityAdami(nu=0.0015), + viscosity_fluid=ViscosityAdami(nu=0.0015), fluid_density_calculator=SummationDensity(), clip_negative_pressure=true), "WCSPH with ViscosityMorris and SummationDensity" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityMorris(nu=0.0015), + viscosity_fluid=ViscosityMorris(nu=0.0015), fluid_density_calculator=SummationDensity(), clip_negative_pressure=true), "WCSPH with smoothing_length=1.3" => (smoothing_length=1.3,), @@ -55,7 +55,7 @@ smoothing_kernel, smoothing_length, sound_speed, - viscosity=viscosity, + viscosity_fluid=viscosity, density_calculator=ContinuityDensity(), acceleration=(0.0, -gravity))), @@ -63,7 +63,7 @@ smoothing_kernel, smoothing_length, sound_speed, - viscosity=viscosity, + viscosity_fluid=viscosity, density_calculator=SummationDensity(), acceleration=(0.0, -gravity)),), @@ -71,7 +71,7 @@ smoothing_kernel, smoothing_length, sound_speed, - viscosity=ViscosityAdami(nu=0.0015), + viscosity_fluid=ViscosityAdami(nu=0.0015), density_calculator=ContinuityDensity(), acceleration=(0.0, -gravity)),), @@ -79,7 +79,7 @@ smoothing_kernel, smoothing_length, sound_speed, - viscosity=ViscosityMorris(nu=0.0015), + viscosity_fluid=ViscosityMorris(nu=0.0015), density_calculator=ContinuityDensity(), acceleration=(0.0, -gravity)),) diff --git a/test/schemes/fluid/viscosity.jl b/test/schemes/fluid/viscosity.jl index 66b770763f..852ebfbc34 100644 --- a/test/schemes/fluid/viscosity.jl +++ b/test/schemes/fluid/viscosity.jl @@ -1,4 +1,3 @@ -include("../../test_util.jl") @testset verbose=true "Viscosity" begin particle_spacing = 0.2 smoothing_length = 1.2 * particle_spacing From 6a1444341c65b573c4886cbc4a9f38a8c41dbaa1 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 21 May 2025 11:34:12 +0200 Subject: [PATCH 21/26] format --- examples/fluid/dam_break_oil_film_2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/fluid/dam_break_oil_film_2d.jl b/examples/fluid/dam_break_oil_film_2d.jl index 5d71c40146..6effdf218a 100644 --- a/examples/fluid/dam_break_oil_film_2d.jl +++ b/examples/fluid/dam_break_oil_film_2d.jl @@ -30,7 +30,8 @@ oil_viscosity = ViscosityMorris(nu=nu_sim_oil) # TODO: broken if both systems use surface tension trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), sol=nothing, fluid_particle_spacing=fluid_particle_spacing, - viscosity_fluid=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length, + viscosity_fluid=ViscosityMorris(nu=nu_sim_water), + smoothing_length=smoothing_length, gravity=gravity, density_diffusion=nothing, sound_speed=sound_speed, prefix="", reference_particle_spacing=fluid_particle_spacing) From 0ee856ef8510820c3e072ba901268187a900af8f Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 22 May 2025 09:16:19 +0200 Subject: [PATCH 22/26] rename viscosity --- examples/fluid/hydrostatic_water_column_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/fluid/hydrostatic_water_column_2d.jl b/examples/fluid/hydrostatic_water_column_2d.jl index b492290756..f60a76f965 100644 --- a/examples/fluid/hydrostatic_water_column_2d.jl +++ b/examples/fluid/hydrostatic_water_column_2d.jl @@ -32,7 +32,7 @@ smoothing_length = 1.2 * fluid_particle_spacing smoothing_kernel = SchoenbergCubicSplineKernel{2}() alpha = 0.02 -viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) +viscosity_fluid = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) fluid_density_calculator = ContinuityDensity() @@ -40,7 +40,7 @@ fluid_density_calculator = ContinuityDensity() system_acceleration = (0.0, -gravity) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, - smoothing_length, viscosity=viscosity, + smoothing_length, viscosity=viscosity_fluid, acceleration=system_acceleration, source_terms=nothing) From f038002171fdd4c41fcd7797d2bb56ba08b141f9 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 22 May 2025 11:47:24 +0200 Subject: [PATCH 23/26] too much renaming --- test/examples/examples_fluid.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index 21eac4eac3..9a26ed1230 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -55,7 +55,7 @@ smoothing_kernel, smoothing_length, sound_speed, - viscosity_fluid=viscosity, + viscosity=viscosity_fluid, density_calculator=ContinuityDensity(), acceleration=(0.0, -gravity))), @@ -63,7 +63,7 @@ smoothing_kernel, smoothing_length, sound_speed, - viscosity_fluid=viscosity, + viscosity=viscosity_fluid, density_calculator=SummationDensity(), acceleration=(0.0, -gravity)),), @@ -71,7 +71,7 @@ smoothing_kernel, smoothing_length, sound_speed, - viscosity_fluid=ViscosityAdami(nu=0.0015), + viscosity=ViscosityAdami(nu=0.0015), density_calculator=ContinuityDensity(), acceleration=(0.0, -gravity)),), @@ -79,7 +79,7 @@ smoothing_kernel, smoothing_length, sound_speed, - viscosity_fluid=ViscosityMorris(nu=0.0015), + viscosity=ViscosityMorris(nu=0.0015), density_calculator=ContinuityDensity(), acceleration=(0.0, -gravity)),) From 9c9186fb2bc8ac431b03ae4de30f2205f54bfdf3 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 23 May 2025 16:57:34 +0200 Subject: [PATCH 24/26] implement suggestions --- src/schemes/fluid/viscosity.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 1d8dbc3c01..ddd74a7481 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -257,7 +257,7 @@ with and an approximation for the strain rate magnitude given by ```math -|S| \approx \frac{\|v_ab\|}{\|r_ab\| + \epsilon}, +|S| \approx \frac{\|v_{ab}\|}{\|r_{ab}\| + \epsilon}, ``` where: @@ -369,7 +369,7 @@ with and an approximation for the strain rate magnitude given by ```math -|S| \approx \frac{\|v_ab\|}{\|r_ab\| + \epsilon}, +|S| \approx \frac{\|v_{ab}\|}{\|r_{ab}\| + \epsilon}, ``` where: From 2f3bbbb067310aecedf645704a3515f4aff0f31d Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 28 May 2025 08:54:37 +0200 Subject: [PATCH 25/26] Update validation_taylor_green_vortex_2d.jl --- .../taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl index db07b67a09..4543412306 100644 --- a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -85,7 +85,7 @@ for density_calculator in density_calculators, perturbation in perturb_coordinat output_directory = joinpath("out_tgv", name_density_calculator * name_perturbation, wcsph ? "wcsph" : "edac", "validation_run_taylor_green_vortex_2d_nparticles_$(n_particles_xy)x$(n_particles_xy)") - saving_callback = SolutionSavingCallback(dt=0.1, + saving_callback = SolutionSavingCallback(dt=0.02, output_directory=output_directory, p_avg=diff_p_loc_p_avg) From 6a83b83ff03887124a13fc5d842697d2fde55cd4 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 13 Jun 2025 14:09:23 +0200 Subject: [PATCH 26/26] fix tests again --- test/examples/gpu.jl | 10 +++++----- test/validation/validation.jl | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl index 9d551a5273..85d657b678 100644 --- a/test/examples/gpu.jl +++ b/test/examples/gpu.jl @@ -178,13 +178,13 @@ end clip_negative_pressure=true), "WCSPH with ViscosityAdami" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityAdami(nu=0.0015f0),), + viscosity_fluid=ViscosityAdami(nu=0.0015f0),), "WCSPH with ViscosityMorris" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityMorris(nu=0.0015f0),), + viscosity_fluid=ViscosityMorris(nu=0.0015f0),), "WCSPH with ViscosityAdami and SummationDensity" => ( # from 0.02*10.0*1.2*0.05/8 - viscosity=ViscosityAdami(nu=0.0015f0), + viscosity_fluid=ViscosityAdami(nu=0.0015f0), fluid_density_calculator=SummationDensity(), maxiters=38, # 38 time steps on CPU clip_negative_pressure=true), @@ -205,7 +205,7 @@ end smoothing_kernel, smoothing_length, sound_speed, - viscosity=viscosity, + viscosity=viscosity_fluid, density_calculator=ContinuityDensity(), acceleration=(0.0, -gravity))), @@ -213,7 +213,7 @@ end smoothing_kernel, smoothing_length, sound_speed, - viscosity=viscosity, + viscosity=viscosity_fluid, density_calculator=SummationDensity(), acceleration=(0.0, -gravity)),) diff --git a/test/validation/validation.jl b/test/validation/validation.jl index d1c8420842..70aceb4bf4 100644 --- a/test/validation/validation.jl +++ b/test/validation/validation.jl @@ -50,7 +50,7 @@ if Sys.ARCH === :aarch64 # MacOS ARM produces slightly different pressure values than x86. # Note that pressure values are in the order of 1e5. - @test isapprox(error_edac_P1, 0, atol=4e-6) + @test isapprox(error_edac_P1, 0, atol=5e-6) @test isapprox(error_edac_P2, 0, atol=4e-11) @test isapprox(error_wcsph_P1, 0, atol=400.0) @test isapprox(error_wcsph_P2, 0, atol=0.03)