diff --git a/NEWS.md b/NEWS.md index f004f1410..fc9b832e3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,8 @@ used in the Julia ecosystem. Notable changes will be documented in this file for ### Features +- 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 when combined with thread pinning (#722). diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 10c1b2ef8..b857c7a12 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -715,7 +715,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} } @@ -736,6 +736,28 @@ @book{Poling2001 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} +} + @article{Zhu2021, author = {Zhu, Yujie and Zhang, Chi and Yu, Yongchuan and Hu, Xiangyu}, journal = {Journal of Hydrodynamics}, diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 35493a8c2..950f8c806 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -47,10 +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 = ViscosityAdami(nu=nu) +alpha = 0.02 +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 + +# Alternatively a kinematic viscosity for water can be set +# nu = 1.0e-6 + +# This allows the use of a physical viscosity model like: +# viscosity_fluid = ViscosityAdami(nu=nu) +# or with additional dissipation through a Smagorinsky model +# 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. # density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) @@ -58,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, @@ -67,12 +78,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 without 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/examples/fluid/dam_break_2phase_2d.jl b/examples/fluid/dam_break_2phase_2d.jl index de99a5f39..247462b52 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 320f18509..f301f9964 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, tspan=tspan, - 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/examples/fluid/hydrostatic_water_column_2d.jl b/examples/fluid/hydrostatic_water_column_2d.jl index b49229075..f60a76f96 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) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index c1f0f72fd..75e638718 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/TrixiParticles.jl b/src/TrixiParticles.jl index bab7e7a81..ff23f9490 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -73,7 +73,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 tensile_instability_control diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index 0e5aaf4c1..6a9663c68 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -357,7 +357,9 @@ 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 diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index e39e7205d..ddd74a748 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -176,6 +176,32 @@ struct ViscosityAdami{ELTYPE} end 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) + 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 + 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 +end + @inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, @@ -185,6 +211,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), @@ -197,36 +224,225 @@ 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 + 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 - eta_tilde = 2 * (eta_a * eta_b) / (eta_a + eta_b) +function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length, + sound_speed) + 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) + +Viscosity model that extends the standard [Adami formulation](@ref ViscosityAdami) +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 + +```math +\nu_{\text{SGS}} = (C_S h)^2 |S|, +``` + +and an approximation for the strain rate magnitude given by + +```math +|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. + +The effective dynamic viscosities are then computed as +```math +\eta_{a,\text{eff}} = \rho_a\, \nu_{\text{eff}}, +``` +and averaged as +```math +\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. + +# Keywords +- `nu`: Standard kinematic viscosity. +- `C_S`: Smagorinsky constant. +- `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] +end + +ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, epsilon) + +@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) + 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 - tmp = eta_tilde / (distance^2 + epsilon * smoothing_length_average^2) - volume_a = m_a / rho_a - volume_b = m_b / rho_b + 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) - # 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. + 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. + # ------------------------------------------------------------------------------ + # 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ᵢ). # - # 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 + # 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_ab‖ / (‖r_ab‖ + δ), + # + # 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 - return visc .* v_diff + # Effective kinematic viscosity is the sum of the standard and SGS parts. + 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::ViscosityAdami, smoothing_length, +function kinematic_viscosity(system, viscosity::ViscosityAdamiSGS, smoothing_length, sound_speed) return viscosity.nu end -@propagate_inbounds function viscous_velocity(v, system, particle) - return current_velocity(v, system, particle) +@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), +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 + +```math +\nu_{\text{SGS}} = (C_S h)^2 |S|, +``` + +and an approximation for the strain rate magnitude given by + +```math +|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. + +The effective dynamic viscosities are then computed as +```math +\eta_{a,\text{eff}} = \rho_a\, \nu_{\text{eff}}, +``` +and averaged as +```math +\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. + +# Keywords +- `nu`: Standard kinematic viscosity. +- `C_S`: Smagorinsky constant. +- `epsilon=0.01`: Parameter to prevent singularities +""" +struct ViscosityMorrisSGS{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] +end + +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 = 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 + + # SGS part: Compute the subgrid-scale eddy viscosity. + # See comments above for `ViscosityAdamiSGS`. + S_mag = norm(v_diff) / (distance + epsilon) + nu_SGS = (viscosity.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 * smoothing_length_average^2) * v_diff + return m_b * force_Morris +end + +function kinematic_viscosity(system, viscosity::ViscosityMorrisSGS, smoothing_length, + sound_speed) + return viscosity.nu end diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index 4c557ee46..141d60278 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=viscosity_fluid, density_calculator=ContinuityDensity(), acceleration=(0.0, -gravity))), @@ -63,7 +63,7 @@ smoothing_kernel, smoothing_length, sound_speed, - viscosity=viscosity, + viscosity=viscosity_fluid, density_calculator=SummationDensity(), acceleration=(0.0, -gravity)),), diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl index 9d551a527..85d657b67 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/schemes/fluid/viscosity.jl b/test/schemes/fluid/viscosity.jl index f11f4b05b..852ebfbc3 100644 --- a/test/schemes/fluid/viscosity.jl +++ b/test/schemes/fluid/viscosity.jl @@ -9,8 +9,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 +29,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 +41,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 +80,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 diff --git a/test/validation/validation.jl b/test/validation/validation.jl index d1c842084..70aceb4bf 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) 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 73884e38c..454341230 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,6 +2,10 @@ 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 paper provides reference data for particle spacings particle_spacings = [0.02, 0.01, 0.005] particle_spacing = 0.02