diff --git a/NEWS.md b/NEWS.md index f004f1410e..fc9b832e30 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 10c1b2ef8c..b857c7a128 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 35493a8c2b..950f8c8060 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 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 320f185092..f301f99644 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 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) 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/TrixiParticles.jl b/src/TrixiParticles.jl index bab7e7a813..ff23f94902 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 0e5aaf4c15..6a9663c68c 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 e39e7205d0..ddd74a7481 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 4c557ee465..141d602787 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 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/schemes/fluid/viscosity.jl b/test/schemes/fluid/viscosity.jl index f11f4b05b7..852ebfbc34 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 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) 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 73884e38c2..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 @@ -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