diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml index 5e093c7b7..4919ff9ef 100644 --- a/.github/workflows/Documenter.yml +++ b/.github/workflows/Documenter.yml @@ -17,8 +17,6 @@ on: workflow_dispatch: concurrency: - # Skip intermediate builds: always. - # Cancel intermediate builds: only if it is a pull request build. group: ${{ github.workflow }}-${{ github.ref }} cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} diff --git a/NEWS.md b/NEWS.md index e60c63e9f..2fb07858b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,24 +3,51 @@ 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.2.4 + +### Features + +- Support for surface tension was added to EDAC (#539) +- A method to prevent penetration of fast moving particles with solids was added (#498) +- Added the callback `SteadyStateReachedCallback` to detect convergence of static simulations (#601) +- Added Ideal Gas State Equation (#607) + +### Documentation + +- GPU Support Documentation was added (#660) +- A new user tutorial was added (#514) + +### Fixes + +- Diverse Doc fixes (#663, #659, #637, #658, #664) +- Simulations can be run with `Float32` (#662) + +### Refactored + +- Surface normal calculation was moved from surface_tension.jl to surface_normal_sph.jl (#539) + ## Version 0.2.3 ### Highlights + Transport Velocity Formulation (TVF) based on the work of Ramachandran et al. "Entropically damped artificial compressibility for SPH" (2019) was added. ## Version 0.2.2 ### Highlights + Hotfix for threaded sampling of complex geometries. ## Version 0.2.1 ### Highlights + Particle sampling of complex geometries from `.stl` and `.asc` files. ## Version 0.2.0 ### Removed + Use of the internal neighborhood search has been removed and replaced with PointNeighbors.jl. ## Development Cycle 0.1 @@ -28,26 +55,33 @@ Use of the internal neighborhood search has been removed and replaced with Point ### Highlights #### Discrete Element Method + A basic implementation of the discrete element method was added. #### Surface Tension and Adhesion Model + A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids" (2013) was added to WCSPH. #### Support for Open Boundaries + Open boundaries using the method of characteristics based on the work of Lastiwka et al., "Permeable and non-reflecting boundary conditions in SPH" (2009) were added for WCSPH and EDAC. ## Pre Initial Release (v0.1.0) + This section summarizes the initial features that TrixiParticles.jl was released with. ### Highlights #### EDAC + An implementation of EDAC (Entropically Damped Artificial Compressibility) was added, which allows for more stable simulations compared to basic WCSPH and reduces spurious pressure oscillations. #### WCSPH + An implementation of WCSPH (Weakly Compressible Smoothed Particle Hydrodynamics), which is the classical SPH approach. Features: + - Correction schemes (Shepard (0. Order) ... MixedKernelGradient (1. Order)) - Density reinitialization - Kernel summation and Continuity equation density formulations @@ -57,4 +91,5 @@ Features: #### TLSPH + An implementation of TLSPH (Total Lagrangian Smoothed Particle Hydrodynamics) for solid bodies enabling FSI (Fluid Structure Interactions). diff --git a/Project.toml b/Project.toml index 398268160..b757a53e5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TrixiParticles" uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a" authors = ["erik.faulhaber <44124897+efaulhaber@users.noreply.github.com>"] -version = "0.2.4-dev" +version = "0.2.4" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/docs/make.jl b/docs/make.jl index 1fb16e935..a2592d2bb 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -139,8 +139,8 @@ makedocs(sitename="TrixiParticles.jl", "Util" => joinpath("general", "util.md") ], "Systems" => [ - "Discrete Element Method (Solid)" => joinpath("systems", - "dem.md"), + "Fluid Models" => joinpath("systems", "fluid.md"), + "Discrete Element Method (Solid)" => joinpath("systems", "dem.md"), "Weakly Compressible SPH (Fluid)" => joinpath("systems", "weakly_compressible_sph.md"), "Entropically Damped Artificial Compressibility for SPH (Fluid)" => joinpath("systems", diff --git a/docs/src/refs.bib b/docs/src/refs.bib index e1b287b07..8927830a9 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -329,6 +329,19 @@ @Article{Hormann2001 publisher = {Elsevier BV}, } +@Article{Huber2016, + title = {On the physically based modeling of surface tension and moving contact lines with dynamic contact angles on the continuum scale}, + journal = {Journal of Computational Physics}, + volume = {310}, + pages = {459-477}, + year = {2016}, + issn = {0021-9991}, + doi = {https://doi.org/10.1016/j.jcp.2016.01.030}, + url = {https://www.sciencedirect.com/science/article/pii/S0021999116000383}, + author = {M. Huber and F. Keller and W. Säckel and M. Hirschler and P. Kunz and S.M. Hassanizadeh and U. Nieken}, + keywords = {SPH, Two-phase flow, Surface tension, Moving contact line, Dynamic contact angle, CSF, CLF} +} + @Article{Jacobson2013, author = {Jacobson, Alec and Kavan, Ladislav and Sorkine-Hornung, Olga}, @@ -491,6 +504,19 @@ @Article{Morris1997 publisher = {Elsevier BV}, } +@article{Morris2000, + author = {Morris, Joseph P.}, + title = {Simulating surface tension with smoothed particle hydrodynamics}, + journal = {International Journal for Numerical Methods in Fluids}, + volume = {33}, + number = {3}, + pages = {333-353}, + keywords = {interfacial flow, meshless methods, surface tension}, + doi = {https://doi.org/10.1002/1097-0363(20000615)33:3<333::AID-FLD11>3.0.CO;2-7}, + year = {2000} +} + + @InProceedings{Mueller2003, author = {M{\"u}ller, Matthias and Charypar, David and Gross, Markus}, booktitle = {Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation}, diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md new file mode 100644 index 000000000..d559790aa --- /dev/null +++ b/docs/src/systems/fluid.md @@ -0,0 +1,212 @@ + +# [Fluid Models](@id fluid_models) +Currently available fluid methods are the [weakly compressible SPH method](@ref wcsph) and the [entropically damped artificial compressibility for SPH](@ref edac). +This page lists models and techniques that apply to both of these methods. + +## [Viscosity](@id viscosity_wcsph) + +TODO: Explain viscosity. + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "fluid", "viscosity.jl")] +``` + +## [Corrections](@id corrections) + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("general", "corrections.jl")] +``` + + +--- + +## [Surface Normals](@id surface_normal) + +### Overview of Surface Normal Calculation in SPH + +Surface normals are essential for modeling surface tension as they provide the directionality of forces acting at the fluid interface. They are calculated based on the particle properties and their spatial distribution within the smoothed particle hydrodynamics (SPH) framework. + +#### Color Field and Gradient-Based Surface Normals + +The surface normal at a particle is derived from the color field, a scalar field assigned to particles to distinguish between different fluid phases or between fluid and air. The color field gradients point towards the interface, and the normalized gradient defines the surface normal direction. + +The simplest SPH formulation for surface normal, \( n_a \), is given as: +```math +n_a = \sum_b m_b \frac{c_b}{\rho_b} \nabla_a W_{ab}, +``` +where: +- \( c_b \) is the color field value for particle \( b \), +- \( m_b \) is the mass of particle \( b \), +- \( \rho_b \) is the density of particle \( b \), +- \( \nabla_a W_{ab} \) is the gradient of the smoothing kernel \( W_{ab} \) with respect to particle \( a \). + +#### Normalization of Surface Normals + +The calculated normals are normalized to unit vectors: +```math +\hat{n}_a = \frac{n_a}{\Vert n_a \Vert}. +``` +Normalization ensures that the magnitude of the normals does not bias the curvature calculations or the resulting surface tension forces. + +#### Handling Noise and Errors in Normal Calculation + +In regions distant from the interface, the calculated normals may be small or inaccurate due to the smoothing kernel's support radius. To mitigate this: +1. Normals below a threshold are excluded from further calculations. +2. Curvature calculations use a corrected formulation to reduce errors near interface fringes. + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "fluid", "surface_normal_sph.jl")] +``` + +--- + +## [Surface Tension](@id surface_tension) + +### Introduction to Surface Tension in SPH + +Surface tension is a key phenomenon in fluid dynamics, influencing the behavior of droplets, bubbles, and fluid interfaces. In SPH, surface tension is modeled as forces arising due to surface curvature and particle interactions, ensuring realistic simulation of capillary effects, droplet coalescence, and fragmentation. + +### Akinci-Based Intra-Particle Force Surface Tension and Wall Adhesion Model + +The Akinci model divides surface tension into distinct force components: + +#### Cohesion Force + +The cohesion force captures the attraction between particles at the fluid interface, creating the effect of surface tension. It is defined by the distance between particles and the support radius \( h_c \), using a kernel-based formulation. + +**Key Features:** +- Particles within half the support radius experience a repulsive force to prevent clustering. +- Particles beyond half the radius but within the support radius experience an attractive force to simulate cohesion. + +Mathematically: +```math +F_{\text{cohesion}} = -\sigma m_b C(r) \frac{r}{\Vert r \Vert}, +``` +where \( C(r) \), the cohesion kernel, is defined as: +```math +C(r)=\frac{32}{\pi h_c^9} +\begin{cases} +(h_c-r)^3 r^3, & \text{if } 2r > h_c, \\ +2(h_c-r)^3 r^3 - \frac{h^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c, \\ +0, & \text{otherwise.} +\end{cases} +``` + +#### Surface Area Minimization Force + +The surface area minimization force models the curvature reduction effects, aligning particle motion to reduce the interface's total area. It acts based on the difference in surface normals: +```math +F_{\text{curvature}} = -\sigma (n_a - n_b), +``` +where \( n_a \) and \( n_b \) are the surface normals of the interacting particles. + +#### Wall Adhesion Force + +This force models the interaction between fluid and solid boundaries, simulating adhesion effects at walls. It uses a custom kernel with a peak at 0.75 times the support radius: +```math +F_{\text{adhesion}} = -\beta m_b A(r) \frac{r}{\Vert r \Vert}, +``` +where \( A(r) \) is the adhesion kernel: +```math +A(r) = \frac{0.007}{h_c^{3.25}} +\begin{cases} +\sqrt[4]{-\frac{4r^2}{h_c} + 6r - 2h_c}, & \text{if } 2r > h_c \text{ and } r \leq h_c, \\ +0, & \text{otherwise.} +\end{cases} +``` + +--- + +### Morris-Based Momentum-Conserving Surface Tension Model + +In addition to the Akinci model, Morris (2000) introduced a momentum-conserving approach to surface tension. This model uses stress tensors to ensure exact conservation of linear momentum, providing a robust method for high-resolution simulations. + +#### Stress Tensor Formulation + +The force is calculated as: +```math +F_{\text{surface tension}} = \nabla \cdot S, +``` +where \( S \) is the stress tensor: +```math +S = \sigma \delta_s (I - \hat{n} \otimes \hat{n}), +``` +with: +- \( \delta_s \): Surface delta function, +- \( \hat{n} \): Unit normal vector, +- \( I \): Identity matrix. + +#### Advantages and Limitations + +While momentum conservation makes this model attractive, it requires additional computational effort and stabilization techniques to address instabilities in high-density regions. + +# [Huber Model](@id huber_model) + +## Introduction to the Huber Contact Force Model + +The **Huber Model**, introduced in [Huber2016](@cite), provides a physically-based approach for simulating surface tension and contact line dynamics in Smoothed Particle Hydrodynamics (SPH). It is specifically designed to address challenges in modeling wetting phenomena, including dynamic contact angles and their influence on fluid behavior. + +The Huber Model introduces a **Contact Line Force (CLF)** that complements the **Continuum Surface Force (CSF)** model, allowing for accurate representation of fluid-fluid and fluid-solid interactions. The dynamic evolution of contact angles emerges naturally from the force balance without requiring artificial adjustments or fitting parameters. + +--- + +## Key Features of the Huber Model + +1. **Dynamic Contact Angles**: + - Captures the transition between static and dynamic contact angles based on interfacial forces and contact line velocities. + - Removes the need for predefined constitutive equations, as the contact angle is derived from the system's dynamics. + +2. **Volume Reformulation**: + - Transforms line-based forces (e.g., those acting at the contact line) into volume-based forces for compatibility with SPH formulations. + - Ensures smooth force distribution near the contact line, reducing numerical artifacts. + +3. **Momentum Balance**: + - Extends the Navier-Stokes equations to include contributions from the contact line, ensuring accurate modeling of wetting and spreading dynamics. + +4. **No Fitting Parameters**: + - Fully physics-driven, requiring only measurable inputs like surface tension coefficients and static contact angles. + +--- + +## Mathematical Formulation + +### Contact Line Force (CLF) + +The force acting along the contact line is derived from the unbalanced Young Force: +```math +f_{\text{CLF}} = \sigma_{\text{wn}} [\cos(\alpha_s) - \cos(\alpha_d)] \hat{\nu}, +``` +where: +- \( \sigma_{\text{wn}} \): Surface tension coefficient of the fluid-fluid interface, +- \( \alpha_s \): Static contact angle, +- \( \alpha_d \): Dynamic contact angle, +- \( \hat{\nu} \): Tangential unit vector along the fluid-solid interface. + +### Volume Reformulation + +To incorporate the CLF into SPH, it is reformulated as a volume force: +```math +F_{\text{CLF}} = f_{\text{CLF}} \delta_{\text{CL}}, +``` +where \( \delta_{\text{CL}} \) is a Dirac delta function approximated by SPH kernels, ensuring the force is applied locally near the contact line. + +--- + +## Applications + +1. **Droplet Dynamics**: + - Simulates droplet spreading, recoiling, and merging with accurate contact line evolution. + +2. **Capillary Action**: + - Models fluid behavior in porous media and confined geometries, where contact lines play a critical role. + +3. **Wetting Phenomena**: + - Predicts equilibrium shapes and transient states of droplets and films on solid surfaces. + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "fluid", "surface_tension.jl")] +``` diff --git a/docs/src/systems/weakly_compressible_sph.md b/docs/src/systems/weakly_compressible_sph.md index e561d2cb9..8848dfeff 100644 --- a/docs/src/systems/weakly_compressible_sph.md +++ b/docs/src/systems/weakly_compressible_sph.md @@ -44,15 +44,6 @@ Modules = [TrixiParticles] Pages = [joinpath("schemes", "fluid", "weakly_compressible_sph", "state_equations.jl")] ``` -## [Viscosity](@id viscosity_wcsph) - -TODO: Explain viscosity. - -```@autodocs -Modules = [TrixiParticles] -Pages = [joinpath("schemes", "fluid", "viscosity.jl")] -``` - ## [Density Diffusion](@id density_diffusion) Density diffusion can be used with [`ContinuityDensity`](@ref) to remove the noise in the @@ -118,72 +109,3 @@ term is very expensive and adds about 40--50% of computational cost. Modules = [TrixiParticles] Pages = [joinpath("schemes", "fluid", "weakly_compressible_sph", "density_diffusion.jl")] ``` - -## [Corrections](@id corrections) - -```@autodocs -Modules = [TrixiParticles] -Pages = [joinpath("general", "corrections.jl")] -``` - -## [Surface Tension](@id surface_tension) - -### Akinci-based intra-particle force surface tension and wall adhesion model -The work by Akinci proposes three forces: -- a cohesion force -- a surface area minimization force -- a wall adhesion force - -The classical model is composed of the curvature minimization and cohesion force. - -#### Cohesion force -The model calculates the cohesion force based on the distance between particles and the support radius ``h_c``. -This force is determined using two distinct regimes within the support radius: -- For particles closer than half the support radius, - a repulsive force is calculated to prevent particle clustering too tightly, - enhancing the simulation's stability and realism. -- Beyond half the support radius and within the full support radius, - an attractive force is computed, simulating the effects of surface tension that draw particles together. -The cohesion force, ``F_{\text{cohesion}}``, for a pair of particles is given by: -```math -F_{\text{cohesion}} = -\sigma m_b C(r) \frac{r}{\Vert r \Vert}, -``` -where: -- ``\sigma`` represents the surface tension coefficient, adjusting the overall strength of the cohesion effect. -- ``C`` is a scalar function of the distance between particles. - -The cohesion kernel ``C`` is defined as -```math -C(r)=\frac{32}{\pi h_c^9} -\begin{cases} -(h_c-r)^3 r^3, & \text{if } 2r > h_c \\ -2(h_c-r)^3 r^3 - \frac{h^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c \\ -0, & \text{otherwise} -\end{cases} -``` - -#### Surface area minimization force -To model the minimization of the surface area and curvature of the fluid, a curvature force is used, which is calculated as -```math -F_{\text{curvature}} = -\sigma (n_a - n_b) -``` - -#### Wall adhesion force -The wall adhesion model proposed by Akinci et al. is based on a kernel function which is 0 from 0.0 to 0.5 support radiia with a maximum at 0.75. -With the force calculated with an adhesion coefficient ``\beta`` as -```math -F_{\text{adhesion}} = -\beta m_b A(r) \frac{r}{\Vert r \Vert}, -``` -with ``A`` being the adhesion kernel defined as -```math -A(r)= \frac{0.007}{h_c^{3.25}} -\begin{cases} -\sqrt[4]{- \frac{4r^2}{h_c} + 6r - 2h_c}, & \text{if } 2r > h_c \text{ and } r \leq h_c \\ -0, & \text{otherwise.} -\end{cases} -``` - -```@autodocs -Modules = [TrixiParticles] -Pages = [joinpath("schemes", "fluid", "surface_tension.jl")] -``` diff --git a/examples/fluid/calibration_water_drop_shape.jl b/examples/fluid/calibration_water_drop_shape.jl new file mode 100644 index 000000000..916033223 --- /dev/null +++ b/examples/fluid/calibration_water_drop_shape.jl @@ -0,0 +1,248 @@ +# In this example we try to approach the static shape of a water droplet on a horizontal plane. +# The shape of a static droplet can be calculated from the Young-Laplace equation. +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.0001 + +boundary_layers = 4 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 5.0) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.0, 0.0) +tank_size = (0.05, 0.01) + +fluid_density = 1000.0 +sound_speed = 50 + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation=state_equation) + +sphere_radius = 0.0015 + +sphere1_center = (0.5, sphere_radius) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=VoxelSphere()) + +# ========================================================================================== +# ==== Fluid +# fluid_smoothing_length = 1.2 * fluid_particle_spacing +# fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}() + +fluid_smoothing_length = 2.5 * fluid_particle_spacing +fluid_smoothing_kernel = WendlandC2Kernel{2}() + +fluid_smoothing_length_2 = 2.5 * fluid_particle_spacing +fluid_smoothing_kernel_2 = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() + +# water at 20C +#nu=0.00089 +nu = 0.00089 +# too much 0.00089 -> 0.00045 -> 0.0002 -> 0.0001 +# no impact 0.00005 + +# the following term only holds for 2d sims +# alpha = 8 * nu / (fluid_smoothing_length * sound_speed) +viscosity = ViscosityAdami(nu=nu) +# density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1) + +# with increased smoothing length surface_tension is too small +sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator, + state_equation, fluid_smoothing_kernel, + fluid_smoothing_length, + viscosity=viscosity, + acceleration=(0.0, -gravity), + surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.008), + correction=AkinciFreeSurfaceCorrection(fluid_density)) + +# 0.001 +# 0.2 > 90 +# 0.15 > 90 +# 0.125 > 90 +# 0.11 60-90 +# 0.1 ~ 45 + +# 0.0008 +# 0.11 ~45-50 +# 0.1125 ~90 +# 0.115 ~90 +# 0.12 > 90 + +# 0.0005 +#0.11 ~100 + +# 0.00025 +# 0.115 ~60 grad aber zu tief in der mitte 0.006 +# 0.115 and increase nu to 0.0002 + +# 0.0001 deg 90 (x-axis: 2mm, y-axis: 1.8mm) +# start x-axis ~1.7mm y-axis: 2.2mm ~110deg distance to boundary to large at adhesion coefficient 1.0 +# increase adhesion coefficient to 1.1 x-axis ~1.8mm y-axis: 2mm ~100deg distance still to high +# increase adhesion coefficient to 1.2 x-axis ~2.8 y-axis: 1.7 distance still to high +# decrease adhesion coefficient to 1.15 x-axis 2.4mm y-axis: 1.9 distance high +# decrease adhesion coefficient to 1.125 x-axis 2-2.2mm y-axis: 1.9 distance high ~90+ +# decrease surface tension coefficient from 0.115 to 0.11 x-axis: 2.6mm y-axis: 1.7 -> invalid too much adhesion +# increase surface tension coefficient from 0.11 to 0.1125 x-axis: 2-2.2mm y-axis: 1.9 +# increase surface tension coefficient from 0.1125 to 0.12 x-axis: 1.9-2mm y-axis: 2.0 +# increase viscosity from 0.0001 to 0.0002 x-axis: 1.8mm y-axis:2.2 +# increase adhesion coefficient from 1.125 to 1.15 x-axis: 1.8mm y-axis:2.1 +# increase adhesion coefficient from 1.15 to 1.2 x-axis: 1.8mm y-axis:2.1 +# increase adhesion coefficient from 1.2 to 1.25 x-axis: 2.6mm y-axis:1.7 +# increase viscosity from 0.0002 to 0.0003 x-axis: 2.1mm y-axis:1.9 <=== (adh=1.25, surft=0.12, nu=0.0003) +# increase adhesion coefficient from 1.25 to 1.275 x-axis: 2.6mm y-axis:1.7 +# decrease adhesion coefficient from 1.275 to 1.26 x-axis: 2.2mm y-axis:1.8 +# decrease adhesion coefficient from 1.26 to 1.255 x-axis: 1.8-2.6mm y-axis:1.8 +# increase viscosity from 0.0003 to 0.00031 x-axis: 2.2-2.8mm y-axis:1.7 +# increase viscosity from 0.00031 to 0.00032 x-axis: 2.2-2.8mm y-axis:1.7 + +# 0.0001 deg 75 (x-axis: 2.2mm, y-axis: 1.6mm) +#(adh=1.25, surft=0.12, nu=0.0003) +# decrease surft=0.11 x-axis: 2.4-2.6mm y-axis:1.6 +# decrease adh=1.2 x-axis: 2.2-2.4mm y-axis:1.8 +# increase viscosity = 0.00035 x-axis: 2.4-3.2mm y-axis:1.6 +# increase viscosity = 0.0004 x-axis: 2.6-2.8mm y-axis:1.4 +# increase viscosity = 0.00045 x-axis: 2.4-3.2mm y-axis:1.4 +# decrease adh=1.15 x-axis: 1.8mm y-axis:2.1 +# increase adh=1.175 x-axis: 1.8mm y-axis:2.2 +# increase adh=1.19 x-axis: 1.8mm y-axis:2.1 +#(adh=1.25, surft=0.12, nu=0.0003) +# increase adh=1.3 and x-axis: 2.4mm y-axis:1.8 +# decrease adh=1.275 and x-axis: 2.4mm y-axis:1.8 +# decrease adh=1.26 and x-axis: 2.4mm y-axis:1.8 +# decrease adh=1.255 and x-axis: 2.4mm y-axis:1.8 +# decrease smoothing length from 4 -> 3.75 x-axis: 2mm y-axis:1.9 + +# 60deg x=2.4, y=1.3 +# adh = 1.15, surft=0.115, nu=0.0003, h=3.75, x=1.7-1.8, y=2.2 +# adh = 1.2, surft=0.115, nu=0.0003, h=3.75, x=2-2.2, y=1.8 +# adh = 1.25, surft=0.115, nu=0.0003, h=3.75, x=2.2, y=1.8 +# adh = 1.3, surft=0.115, nu=0.0003, h=3.75, x=2.6, y=1.6 +# adh = 1.3, surft=0.12, nu=0.0003, h=3.75, x=2.2-2.6, y=1.6 +# adh = 1.3, surft=0.125, nu=0.0003, h=3.75, x=2, y=1.8 +# adh = 1.35, surft=0.125, nu=0.0003, h=3.75, x=2.5, y=1.6 +# adh = 1.35, surft=0.13, nu=0.0003, h=3.75, x=2-2.2, y=1.8 +# adh = 1.4, surft=0.13, nu=0.0003, h=3.75, x=2-2.2, y=1.8 +# adh = 1.45, surft=0.13, nu=0.0003, h=3.75, x=2.6, y=1.6 +# not possible with smoothing_length -> lower to 3 and switch from C4 to C2 +# adh = 1.45, surft=0.13, nu=0.0003, h=3.0, x=1.4, y=2.6 +# adh = 1.45, surft=0.12, nu=0.0003, h=3.0, x=1.4, y=2.4 +# adh = 1.5, surft=0.12, nu=0.0003, h=3.0, x=1.6, y=2.6 +# adh = 1.5, surft=0.11, nu=0.0003, h=3.0, x=1.4, y=2.4 +# adh = 1.5, surft=0.10, nu=0.0003, h=3.0, x=1.4, y=2.4 +# adh = 1.5, surft=0.05, nu=0.0003, h=3.0, x=3.5, y=2.0 +# adh = 1.4, surft=0.075, nu=0.0003, h=3.0, x=3, y=2.0 +# adh = 1.4, surft=0.075, nu=0.0003, h=2.0, x=3, y=2.0 crap +# adh = 1.4, surft=0.05, nu=0.0003, h=2.0, x=3, y=2.0 crap +# adh = 1.4, surft=0.01, nu=0.0003, h=2.0, x=3, y=2.0 crap +# adh = 1.4, surft=0.01, nu=0.0001, h=2.0, x=3, y=2.0 crap +# adh = 1.4, surft=0.001, nu=0.00089, h=2.0, x=3, y=2.0 crap +# adh = 1.2, surft=0.002, nu=0.00089, h=2.0, x=3, y=2.0 crap +# adh = 1.2, surft=0.002, nu=0.00089, h=2.5, x=2.5, y=1.4 close but too much adh +# adh = 0.6, surft=0.002, nu=0.00089, h=2.5, x=2.7, y=1.3 close but too much adh +# adh = 0.3, surft=0.002, nu=0.00089, h=2.5, x=2.8, y=1.2 +# adh = 0.15, surft=0.002, nu=0.00089, h=2.5, x=2.8, y=1.2 +# adh = 0.075, surft=0.002, nu=0.00089, h=2.5, x=3.4, y=1.2 +# adh = 0.075, surft=0.004, nu=0.00089, h=2.5, x=2.8, y=1.4 +# adh = 0.075, surft=0.01, nu=0.00089, h=2.5, x=2.2, y=1.7 +# adh = 0.075, surft=0.007, nu=0.00089, h=2.5, x=2.6, y=1.4 <- okish but form needs more adh +# adh = 0.075, surft=0.008, nu=0.00089, h=2.5, x=2.6, y=1.4 +# adh = 0.1, surft=0.01, nu=0.00089, h=2.5, x=2.4, y=1.6 <- okish but form needs more adh +# adh = 0.15, surft=0.01, nu=0.00089, wall_nu=0.00089, h=2.5, x=2.4, y=1.8 +# adh = 0.15, surft=0.01, nu=0.00089, h=2.5, wall_nu=0.5*0.00089, x=2.2, y=1.7 <-lower wall viscosity +# adh = 0.15, surft=0.01, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.6, y=1.6 <-lower wall viscosity +# adh = 0.15, surft=0.015, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.2, y=1.8 +# adh = 0.2, surft=0.015, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.2, y=1.8 +# adh = 0.2, surft=0.0125, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.6, y=1.7 +# adh = 0.25, surft=0.014, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.4, y=1.6 +# adh = 0.25, surft=0.014, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.8, y=1.5 diverged for up to 3s +# adh = 0.25, surft=0.015, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.5, y=1.6 diverged at 2.722 <---- +# adh = 0.25, surft=0.016, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.5, y=1.6 +# adh = 0.1, surft=0.0125, nu=0.00089, h=2.5, wall_nu=20*0.00089, x=2.2, y=1.7 <- 75d +# adh = 0.15, surft=0.0125, nu=0.00089, h=2.5, wall_nu=20*0.00089, x=2.2, y=1.7 +# adh = 0.15, surft=0.0125, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.3, y=1.7 +# adh = 0.2, surft=0.0125, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.6, y=1.6 +# adh = 0.15, surft=0.0125, nu=0.00089, h=2.5, wall_nu=0.00089, x=2.6, y=1.5 +# adh = 0.15, surft=0.0125, nu=0.00089, h=2.5, wall_nu=5 * 0.00089, x=2.4, y=1.6 +# adh = 0.13, surft=0.011, nu=0.00089, h=2.5, wall_nu=5*0.00089, x=2.6, y=1.5 +# adh = 0.11, surft=0.011, nu=0.00089, h=2.5, wall_nu=5*0.00089, x=2.4, y=1.6 + +# 75deg (x-axis: 2.2mm, y-axis: 1.6mm) +# adh = 0.20, surft=0.015, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.5, y=1.6 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=20*0.00089, x=2, y=1.7 <- 90d +# adh = 0.1, surft=0.0125, nu=0.00089, h=2.5, wall_nu=20*0.00089, x=2.2, y=1.7 +# adh = 0.1, surft=0.0125, nu=0.00089, h=2.5, wall_nu=7*0.00089, x=2.2, y=1.6 +# adh = 0.08, surft=0.011, nu=0.00089, h=2.5, wall_nu=7*0.00089, x=2.3, y=1.6 +# adh = 0.07, surft=0.011, nu=0.00089, h=2.5, wall_nu=7*0.00089, x=2.2, y=1.6 +# adh = 0.07, surft=0.011, nu=0.00089, h=2.5, wall_nu=5*0.00089, x=2.3, y=1.6 +# adh = 0.05, surft=0.008, nu=0.00089, h=2.5, wall_nu=2*0.00089, x=2.2, y=1.6 + +# 90deg (x-axis: 2mm, y-axis: 1.8mm) +# adh = 0.15, surft=0.015, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.5, y=1.6 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=0.25*0.00089, x=2.2, y=1.7 after capallariy test +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=0.4*0.00089, x=2.2, y=1.6 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=0.00089, x=2.2, y=1.6 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=2*0.00089, x=2.2, y=1.7 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=4*0.00089, x=2.2, y=1.7 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.2, y=1.7 +# adh = 0.05, surft=0.0125, nu=0.00089, h=2.5, wall_nu=20*0.00089, x=2, y=1.7 <- does not work in capallariy test +# adh = 0.04, surft=0.0125, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.1, y=1.7 <- further testing showed surface tension is about 10-20% too high +# adh = 0.04, surft=0.011, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.2, y=1.7 +# adh = 0.03, surft=0.011, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.2, y=1.7 **** <- seems to be still too high +# adh = 0.02, surft=0.011, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.2, y=1.8 +# adh = 0.02, surft=0.011, nu=0.00089, h=2.5, wall_nu=15*0.00089, x=2.2, y=1.7 +# adh = 0.01, surft=0.011, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.2, y=1.7 +# adh = 0.01, surft=0.008, nu=0.00089, h=2.5, wall_nu=5*0.00089, x=2.2, y=1.6 +# adh = 0.01, surft=0.008, nu=0.00089, h=2.5, wall_nu=10*0.00089, x=2.2, y=1.6 + +# sphere = WeaklyCompressibleSPHSystem(sphere2, fluid_density_calculator, +# state_equation, fluid_smoothing_kernel_2, +# fluid_smoothing_length_2, viscosity=viscosity, +# density_diffusion=density_diffusion, +# acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel_2, + fluid_smoothing_length_2, + viscosity=ViscosityAdami(nu=10 * nu)) + +# adhesion_coefficient = 1.0 and surface_tension_coefficient=0.01 for perfect wetting +# adhesion_coefficient = 0.001 and surface_tension_coefficient=2.0 for no wetting +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, + adhesion_coefficient=0.01) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(boundary_system, sphere_surface_tension) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="", + write_meta_data=true) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-5, # Default abstol is 1e-6 + reltol=1e-4, # Default reltol is 1e-3 + dt=1e-4, maxiters=1E12, + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 117dde772..829ab3fa5 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -61,7 +61,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, smoothing_length, viscosity=viscosity, density_diffusion=density_diffusion, acceleration=(0.0, -gravity), correction=nothing, - surface_tension=nothing) + surface_tension=nothing, + reference_particle_spacing=fluid_particle_spacing) # ========================================================================================== # ==== Boundary diff --git a/examples/fluid/dam_break_oil_film_2d.jl b/examples/fluid/dam_break_oil_film_2d.jl index b221d8cd6..222085375 100644 --- a/examples/fluid/dam_break_oil_film_2d.jl +++ b/examples/fluid/dam_break_oil_film_2d.jl @@ -25,18 +25,31 @@ nu_ratio = nu_water / nu_oil nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil) nu_sim_water = nu_ratio * nu_sim_oil +oil_viscosity = ViscosityMorris(nu=nu_sim_oil) +# Physical values +nu_water = 8.9E-7 +nu_oil = 6E-5 +nu_ratio = nu_water / nu_oil + +nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil) +nu_sim_water = nu_ratio * nu_sim_oil + oil_viscosity = ViscosityMorris(nu=nu_sim_oil) +# TODO: broken if both are set to 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, - gravity=gravity, density_diffusion=nothing, sound_speed=sound_speed) + gravity=gravity, tspan=tspan, density_diffusion=nothing, + sound_speed=sound_speed, prefix="") # ========================================================================================== # ==== Setup oil layer oil_size = (W, 0.1 * H) oil_density = 700.0 +oil_eos = StateEquationCole(; sound_speed, reference_density=oil_density, exponent=1, + clip_negative_pressure=false) oil = RectangularShape(fluid_particle_spacing, round.(Int, oil_size ./ fluid_particle_spacing), @@ -50,12 +63,19 @@ end oil_state_equation = StateEquationCole(; sound_speed, reference_density=oil_density, exponent=1, clip_negative_pressure=false) oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator, - oil_state_equation, smoothing_kernel, + oil_eos, smoothing_kernel, smoothing_length, viscosity=oil_viscosity, - density_diffusion=density_diffusion, acceleration=(0.0, -gravity), - surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.02), - correction=AkinciFreeSurfaceCorrection(oil_density)) + surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.01), + correction=AkinciFreeSurfaceCorrection(oil_density), + reference_particle_spacing=fluid_particle_spacing) + +# oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator, +# oil_eos, smoothing_kernel, +# smoothing_length, viscosity=oil_viscosity, +# acceleration=(0.0, -gravity), +# surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.03), +# reference_particle_spacing=fluid_particle_spacing) # ========================================================================================== # ==== Simulation diff --git a/examples/fluid/drops_on_box.jl b/examples/fluid/drops_on_box.jl new file mode 100644 index 000000000..54974e4dd --- /dev/null +++ b/examples/fluid/drops_on_box.jl @@ -0,0 +1,143 @@ +# In this example two circles of water drop to the floor demonstrating the difference +# between the behavior with and without surface tension modelling. +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.0005 + +boundary_layers = 3 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 1.0) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.0, 0.0) +tank_size = (1.0, 0.5) + +fluid_density = 1000.0 +sound_speed = 100 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation=state_equation) + +box = RectangularShape(fluid_particle_spacing, (300, 125), (0.485, 0.0), + density=fluid_density) + +# box = SphereShape(fluid_particle_spacing, 1.0, (0.5, 0.0), +# fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0)) + +sphere_radius = 0.0025 + +sphere1_center = (0.5, 0.05) +sphere2_center = (0.5, 0.1) +sphere3_center = (0.5, 0.15) +sphere4_center = (0.5, 0.2) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0)) +sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0)) +sphere3 = SphereShape(fluid_particle_spacing, sphere_radius, sphere3_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0)) +sphere4 = SphereShape(fluid_particle_spacing, sphere_radius, sphere4_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0)) + +water = union(sphere1, sphere2, sphere3, sphere4) + +# ========================================================================================== +# ==== Fluid +fluid_smoothing_length = 1.0 * fluid_particle_spacing +fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}() + +fluid_density_calculator = ContinuityDensity() + +nu = 0.00089 +alpha = 0.75 * 8 * nu / (fluid_smoothing_length * sound_speed) +# viscosity = ViscosityAdami(nu=nu) +viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) + +# density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1) + +sphere_surface_tension = WeaklyCompressibleSPHSystem(water, fluid_density_calculator, + state_equation, fluid_smoothing_kernel, + fluid_smoothing_length, + viscosity=viscosity, + acceleration=(0.0, -gravity), + surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.1 * + 0.011), + correction=AkinciFreeSurfaceCorrection(fluid_density), + reference_particle_spacing=fluid_particle_spacing) + +# sphere_surface_tension = WeaklyCompressibleSPHSystem(water, fluid_density_calculator, +# state_equation, fluid_smoothing_kernel, +# fluid_smoothing_length, +# viscosity=viscosity, +# acceleration=(0.0, -gravity)) + +# sphere_surface_tension2 = WeaklyCompressibleSPHSystem(sphere2, fluid_density_calculator, +# state_equation, fluid_smoothing_kernel, +# fluid_smoothing_length, +# viscosity=viscosity, +# acceleration=(0.0, -gravity), +# surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.011), +# correction=AkinciFreeSurfaceCorrection(fluid_density)) + +# sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator, +# state_equation, fluid_smoothing_kernel, +# fluid_smoothing_length, +# viscosity=viscosity, +# acceleration=(0.0, -gravity), +# surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.01), +# correction=AkinciFreeSurfaceCorrection(fluid_density)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, fluid_smoothing_length, + viscosity=ViscosityAdami(nu=0.5 * nu)) + +boundary_model_box = BoundaryModelDummyParticles(box.density, box.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, + fluid_smoothing_length, + viscosity=ViscosityAdami(nu=0.5 * nu)) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, + adhesion_coefficient=0.01) + +boundary_system2 = BoundarySPHSystem(box, boundary_model_box, + adhesion_coefficient=0.01) + +# boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +# boundary_system2 = BoundarySPHSystem(box, boundary_model_box) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(boundary_system, boundary_system2, sphere_surface_tension) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.0025, output_directory="out", + prefix="four_drop_terminal_m20_90d_01surft_05wallnu_075viscmon_001adh", + write_meta_data=true) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-5, # Default abstol is 1e-6 + reltol=1e-3, # Default reltol is 1e-3 + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl index cd41d6da8..0de249587 100644 --- a/examples/fluid/falling_water_spheres_2d.jl +++ b/examples/fluid/falling_water_spheres_2d.jl @@ -40,7 +40,7 @@ sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center, # ========================================================================================== # ==== Fluid -fluid_smoothing_length = 1.0 * fluid_particle_spacing +fluid_smoothing_length = 1.0 * fluid_particle_spacing - eps() fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}() fluid_density_calculator = ContinuityDensity() @@ -50,13 +50,13 @@ alpha = 8 * nu / (fluid_smoothing_length * sound_speed) viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1) -sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator, - state_equation, fluid_smoothing_kernel, +sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel, fluid_smoothing_length, - viscosity=viscosity, + sound_speed, viscosity=viscosity, + density_calculator=ContinuityDensity(), acceleration=(0.0, -gravity), surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05), - correction=AkinciFreeSurfaceCorrection(fluid_density)) + reference_particle_spacing=fluid_particle_spacing) sphere = WeaklyCompressibleSPHSystem(sphere2, fluid_density_calculator, state_equation, fluid_smoothing_kernel, @@ -75,16 +75,17 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundar viscosity=ViscosityAdami(nu=wall_viscosity)) boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, - adhesion_coefficient=0.001) + adhesion_coefficient=1.0, + surface_normal_method=StaticNormals((0.0, 1.0))) # ========================================================================================== # ==== Simulation -semi = Semidiscretization(boundary_system, sphere_surface_tension, sphere) +semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=50) -saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="", - write_meta_data=true) +saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", + prefix="", write_meta_data=true) callbacks = CallbackSet(info_callback, saving_callback) diff --git a/examples/fluid/falling_water_spheres_3d.jl b/examples/fluid/falling_water_spheres_3d.jl index a4ebb799e..078defee7 100644 --- a/examples/fluid/falling_water_spheres_3d.jl +++ b/examples/fluid/falling_water_spheres_3d.jl @@ -3,30 +3,30 @@ using OrdinaryDiffEq # ========================================================================================== # ==== Resolution -fluid_particle_spacing = 0.008 +fluid_particle_spacing = 0.005 # ========================================================================================== # ==== Experiment Setup gravity = 9.81 -nu = 0.01 +nu = 0.001 fluid_density = 1000.0 sound_speed = 50 sphere1_radius = 0.05 -sphere1_center = (0.5, 0.5, 0.2) -sphere2_center = (1.5, 0.5, 0.2) +sphere1_center = (0.5, 0.5, 0.075) +sphere2_center = (1.5, 0.5, 0.075) sphere1 = SphereShape(fluid_particle_spacing, sphere1_radius, sphere1_center, - fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -2.0)) + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -1.0)) sphere2 = SphereShape(fluid_particle_spacing, sphere1_radius, sphere2_center, - fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -2.0)) + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -1.0)) # `compact_support` needs to be `2.0 * particle_spacing` to be correct fluid_smoothing_length = 1.0 * fluid_particle_spacing trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "falling_water_spheres_2d.jl"), - fluid_particle_spacing=fluid_particle_spacing, tspan=(0.0, 0.2), + fluid_particle_spacing=fluid_particle_spacing, tspan=(0.0, 0.1), initial_fluid_size=(0.0, 0.0, 0.0), tank_size=(2.0, 1.0, 0.1), sound_speed=sound_speed, faces=(true, true, true, true, true, false), @@ -34,4 +34,4 @@ trixi_include(@__MODULE__, fluid_smoothing_length=fluid_smoothing_length, fluid_smoothing_kernel=SchoenbergCubicSplineKernel{3}(), nu=nu, alpha=10 * nu / (fluid_smoothing_length * sound_speed), - surface_tension_coefficient=1.5, adhesion_coefficient=0.5) + surface_tension_coefficient=10, adhesion_coefficient=0.1) diff --git a/examples/fluid/falling_water_spheres_morris_2d.jl b/examples/fluid/falling_water_spheres_morris_2d.jl new file mode 100644 index 000000000..b328aa9b2 --- /dev/null +++ b/examples/fluid/falling_water_spheres_morris_2d.jl @@ -0,0 +1,129 @@ +# In this example two circles of water drop to the floor demonstrating the difference +# between the behavior with and without surface tension modelling. +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.002 + +boundary_layers = 4 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 0.2) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.0, 0.0) +tank_size = (2.0, 0.5) + +fluid_density = 1000.0 +sound_speed = 100 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation=state_equation) + +sphere_radius = 0.05 + +sphere1_center = (0.5, 0.2) +sphere2_center = (1.5, 0.2) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -1.0)) +sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -1.0)) + +# ========================================================================================== +# ==== Fluid +fluid_smoothing_length = 3.5 * fluid_particle_spacing +fluid_smoothing_kernel = WendlandC2Kernel{2}() + +nu = 0.01 +viscosity = ViscosityMorris(nu=nu) +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +# TODO: sinks into wall with EDAC and MixedKernelGradientCorrection +# sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, state_equation, +# fluid_smoothing_kernel, +# fluid_smoothing_length, +# viscosity=viscosity, +# density_calculator=SummationDensity(), +# acceleration=(0.0, -gravity), +# reference_particle_spacing=fluid_particle_spacing, +# surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.0728), +# correction=MixedKernelGradientCorrection(), +# surface_normal_method=ColorfieldSurfaceNormal(ideal_density_threshold=0.0, +# interface_threshold=0.025)) + +# sphere = WeaklyCompressibleSPHSystem(sphere2, state_equation, fluid_smoothing_kernel, +# fluid_smoothing_length, +# viscosity=viscosity, +# density_calculator=SummationDensity(), +# acceleration=(0.0, -gravity), +# reference_particle_spacing=fluid_particle_spacing, +# surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.0728), +# surface_normal_method=ColorfieldSurfaceNormal(ideal_density_threshold=0.0, +# interface_threshold=0.025)) + +sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel, + fluid_smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=ContinuityDensity(), + acceleration=(0.0, -gravity), + reference_particle_spacing=fluid_particle_spacing, + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.0728), + surface_normal_method=ColorfieldSurfaceNormal(ideal_density_threshold=0.9, + interface_threshold=0.001, boundary_contact_threshold=0.1)) + +sphere = EntropicallyDampedSPHSystem(sphere2, fluid_smoothing_kernel, + fluid_smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=SummationDensity(), + acceleration=(0.0, -gravity), + reference_particle_spacing=fluid_particle_spacing, + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.0728), + surface_normal_method=ColorfieldSurfaceNormal(ideal_density_threshold=0.0, + interface_threshold=0.025)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +wall_viscosity = nu +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, fluid_smoothing_length, + viscosity=ViscosityAdami(nu=wall_viscosity)) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, + surface_normal_method=StaticNormals((0.0, 1.0))) + +# ========================================================================================== +# ==== Simulation + +# min_corner = (-1.0, -1.0) +# max_corner = (2.5, 1.0) +# cell_list = TrixiParticles.PointNeighbors.FullGridCellList(; min_corner, max_corner) +# semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system, +# neighborhood_search=TrivialNeighborhoodSearch{2}()) +semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", + prefix="", write_meta_data=true) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-6, # Default abstol is 1e-6 + reltol=1e-4, # Default reltol is 1e-3 + dt=1e-6, + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/falling_water_spheres_wetting_2d.jl b/examples/fluid/falling_water_spheres_wetting_2d.jl new file mode 100644 index 000000000..24a59d36c --- /dev/null +++ b/examples/fluid/falling_water_spheres_wetting_2d.jl @@ -0,0 +1,130 @@ +# In this example two circles of water drop to the floor demonstrating the difference +# between the behavior with and without surface tension modelling. +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +# fluid_particle_spacing = 0.001 +fluid_particle_spacing = 0.0005 +boundary_layers = 4 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 1.0) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.0, 0.0) +sphere_radius = 0.01 + +tank_size = (40 * sphere_radius, 0.5) + +fluid_density = 1000.0 +sound_speed = 100 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation=state_equation) + +sphere1_center = (tank_size[1] / 4, sphere_radius) +sphere2_center = (3 * tank_size[1] / 4, sphere_radius) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -0.0)) +sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -0.0)) + +# ========================================================================================== +# ==== Fluid +fluid_smoothing_length = 3.5 * fluid_particle_spacing +fluid_smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() + +# nu = 0.005 +# alpha = 8 * nu / (fluid_smoothing_length * sound_speed) +# viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) +nu = 0.005 +viscosity = ViscosityMorris(nu=nu) +# density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1) + +# sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, state_equation, +# fluid_smoothing_kernel, +# fluid_smoothing_length, +# viscosity=viscosity, +# density_calculator=SummationDensity(), +# acceleration=(0.0, -gravity), +# reference_particle_spacing=fluid_particle_spacing, +# correction=MixedKernelGradientCorrection(), +# surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.02, +# contact_model=HuberContactModel()), +# surface_normal_method=ColorfieldSurfaceNormal(ideal_density_threshold=0.0, +# interface_threshold=0.025)) + +# sphere = WeaklyCompressibleSPHSystem(sphere2, state_equation, fluid_smoothing_kernel, +# fluid_smoothing_length, +# viscosity=viscosity, +# density_calculator=SummationDensity(), +# acceleration=(0.0, -gravity), +# reference_particle_spacing=fluid_particle_spacing, +# correction=MixedKernelGradientCorrection(), +# surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.02), +# surface_normal_method=ColorfieldSurfaceNormal(ideal_density_threshold=0.0, +# interface_threshold=0.025)) + +sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel, + fluid_smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=ContinuityDensity(), + acceleration=(0.0, -gravity), + reference_particle_spacing=fluid_particle_spacing, + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.0728, + contact_model=HuberContactModel()), + surface_normal_method=ColorfieldSurfaceNormal(ideal_density_threshold=0.0, + interface_threshold=0.025)) + +sphere = EntropicallyDampedSPHSystem(sphere2, fluid_smoothing_kernel, + fluid_smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=ContinuityDensity(), + acceleration=(0.0, -gravity), + reference_particle_spacing=fluid_particle_spacing, + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.0728), + surface_normal_method=ColorfieldSurfaceNormal(ideal_density_threshold=0.0, + interface_threshold=0.025)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +wall_viscosity = nu +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, fluid_smoothing_length, + viscosity=ViscosityAdami(nu=wall_viscosity)) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, + surface_normal_method=StaticNormals((0.0, 1.0)), + static_contact_angle=90) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.0001, output_directory="out", + prefix="", write_meta_data=true) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-6, # Default abstol is 1e-6 + reltol=1e-4, # Default reltol is 1e-3 + dt=1e-6, + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/falling_water_spheres_zf.jl b/examples/fluid/falling_water_spheres_zf.jl new file mode 100644 index 000000000..a9112df81 --- /dev/null +++ b/examples/fluid/falling_water_spheres_zf.jl @@ -0,0 +1,163 @@ +using TrixiParticles +using OrdinaryDiffEq +using Random # Import Random module for generating random positions and radii + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.00125 + +boundary_layers = 4 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 0.5) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.0, 0.0) +tank_size = (2.0, 0.5) + +fluid_density = 1000.0 +sound_speed = 100 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation=state_equation) + +box = RectangularShape(fluid_particle_spacing, + (Int(0.25 / fluid_particle_spacing), + Int(0.1 / fluid_particle_spacing)), (0.5 * tank_size[1], 0.0), + density=fluid_density) + +# Sphere radius range +sphere_radius_min = 0.005 # Minimum radius +sphere_radius_max = 0.025 # Maximum radius +sphere_radius_range = (sphere_radius_min, sphere_radius_max) + +# Set random seed for reproducibility (optional) +Random.seed!(1234) + +# Define the tank boundaries to avoid placing spheres too close to the walls +x_range = (0.0, tank_size[1]) +y_range = (0.1, tank_size[2] + 0.1) + +# Function to generate non-overlapping spheres with random radii +function generate_non_overlapping_spheres(num_spheres, radius_range, x_range, y_range) + spheres = [] + positions_radii = [] + attempt = 0 + max_attempts = 1000 # Prevent infinite loops + while length(spheres) < num_spheres && attempt < max_attempts + attempt += 1 + # Generate random radius within the specified range + radius = rand() * (radius_range[2] - radius_range[1]) + radius_range[1] + # Adjust x and y ranges to account for sphere radius and avoid walls + x_min = x_range[1] + radius + 0.01 # Add small buffer to avoid wall overlap + x_max = x_range[2] - radius - 0.01 + y_min = y_range[1] + radius + 0.01 + y_max = y_range[2] - radius - 0.01 + if x_min >= x_max || y_min >= y_max + continue # Skip if adjusted ranges are invalid + end + # Generate random position within adjusted ranges + position = (rand() * (x_max - x_min) + x_min, + rand() * (y_max - y_min) + y_min) + # Check for overlap with existing spheres + no_overlap = true + for (pos_existing, radius_existing) in positions_radii + dist = sqrt((position[1] - pos_existing[1])^2 + + (position[2] - pos_existing[2])^2) + if dist < (radius + radius_existing) + 0.01 # Add small buffer + no_overlap = false + break + end + end + if no_overlap + # Create the sphere + sphere = SphereShape(fluid_particle_spacing, radius, position, + fluid_density, sphere_type=VoxelSphere(), + velocity=(0.0, -1.0)) + push!(spheres, sphere) + push!(positions_radii, (position, radius)) + end + end + if length(spheres) < num_spheres + error("Could not generate the required number of non-overlapping spheres within the maximum attempts.") + end + return spheres +end + +num_spheres = 20 +spheres = generate_non_overlapping_spheres(num_spheres, sphere_radius_range, x_range, + y_range) + +# Combine all spheres using the union function +combined_spheres = reduce(union, spheres) + +# ========================================================================================== +# ==== Fluid +fluid_smoothing_length = 3.5 * fluid_particle_spacing +fluid_smoothing_kernel = WendlandC2Kernel{2}() + +nu = 0.0025 +viscosity = ViscosityMorris(nu=nu) +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +sphere_surface_tension = EntropicallyDampedSPHSystem(combined_spheres, + fluid_smoothing_kernel, + fluid_smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=ContinuityDensity(), + acceleration=(0.0, -gravity), + reference_particle_spacing=fluid_particle_spacing, + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.0728), + surface_normal_method=ColorfieldSurfaceNormal(ideal_density_threshold=0.9, + interface_threshold=0.001, + boundary_contact_threshold=0.0)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +wall_viscosity = nu +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, fluid_smoothing_length, + viscosity=ViscosityAdami(nu=wall_viscosity), + correction=MixedKernelGradientCorrection()) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, + surface_normal_method=StaticNormals((0.0, 1.0))) + +box_model = BoundaryModelDummyParticles(box.density, box.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, fluid_smoothing_length, + viscosity=ViscosityAdami(nu=wall_viscosity), + correction=MixedKernelGradientCorrection()) + +box_system = BoundarySPHSystem(box, box_model, + surface_normal_method=StaticNormals((0.0, 1.0))) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(sphere_surface_tension, boundary_system, box_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", + prefix="", write_meta_data=true) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error-based) time step size control. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-6, # Default abstol is 1e-6 + reltol=1e-4, # Default reltol is 1e-3 + dt=1e-6, + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/periodic_field_2d.jl b/examples/fluid/periodic_field_2d.jl new file mode 100644 index 000000000..692696f5a --- /dev/null +++ b/examples/fluid/periodic_field_2d.jl @@ -0,0 +1,213 @@ +using TrixiParticles +using OrdinaryDiffEq + +# Domain size +domain_width = 0.5 +domain_height = 1.0 +no_particles = 50 + +# Particle spacing +particle_spacing = domain_width / no_particles # Adjust for resolution + +# Numerical settings +smoothing_length = 1.2 * particle_spacing +sound_speed = 20.0 + +# Fluid properties +fluid_density = 1.0 +# No gravity +gravity = (0.0, 0.0) + +# Time span +tspan = (0.0, 1.0) + +rect_size = (domain_width, domain_width) + +color0 = RectangularShape(particle_spacing, + round.(Int, rect_size ./ particle_spacing), + zeros(length(rect_size)), density=fluid_density) + +color1 = RectangularShape(particle_spacing, + round.(Int, rect_size ./ particle_spacing), + (0.0, domain_width), density=fluid_density) + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1, + clip_negative_pressure=true) + +wcsph_color0 = WeaklyCompressibleSPHSystem(color0, SummationDensity(), + state_equation, SchoenbergCubicSplineKernel{2}(), + smoothing_length, + reference_particle_spacing=particle_spacing, + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=1.0), + color_value=0) + +wcsph_color1 = WeaklyCompressibleSPHSystem(color1, SummationDensity(), + state_equation, SchoenbergCubicSplineKernel{2}(), + smoothing_length, + reference_particle_spacing=particle_spacing, + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=1.0), + color_value=1) + +periodic_box = PeriodicBox(min_corner=[0.0, 0.0], max_corner=[domain_width, domain_height]) +semi = Semidiscretization(wcsph_color0, wcsph_color1, + neighborhood_search=GridNeighborhoodSearch{2}(update_strategy=nothing, + periodic_box=periodic_box)) +ode = semidiscretize(semi, tspan, data_type=nothing) + +info_callback = InfoCallback(interval=100) + +solution_prefix = "" +saving_callback = SolutionSavingCallback(dt=0.02, prefix=solution_prefix) + +# This can be overwritten with `trixi_include` +extra_callback = nothing + +use_reinit = false +stepsize_callback = StepsizeCallback(cfl=0.9) + +callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback, extra_callback) + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # This is overwritten by the stepsize callback + save_everystep=false, callback=callbacks); + +# # Create particle coordinates +# nx = Int(domain_width / particle_spacing) + 1 +# ny = Int(domain_height / particle_spacing) + 1 +# x_coords = range(0.0, stop=domain_width, length=nx) +# y_coords = range(0.0, stop=domain_height, length=ny) + +# coordinates = [ [x, y] for y in y_coords, x in x_coords ] +# coordinates = reduce(vcat, coordinates)' + +# # Assign colors based on y-coordinate +# colors = [ coord[2] > 0.5 ? 1 : 0 for coord in eachcol(coordinates) ] + +# # Create particle properties +# n_particles = size(coordinates, 2) +# particle_mass = fluid_density * particle_spacing^2 + +# particles = Particles( +# coordinates = coordinates, +# velocities = zeros(2, n_particles), +# masses = fill(particle_mass, n_particles), +# densities = fill(fluid_density, n_particles), +# colors = colors +# ) + +# # Initialize random velocities +# # Internal energy per particle +# internal_energy = 0.5 * sound_speed^2 + +# # Desired kinetic energy per particle +# desired_ke_per_particle = internal_energy * 1e-6 + +# # Generate random velocities +# using Random +# Random.seed!(1234) # For reproducibility + +# velocities = zeros(2, n_particles) +# for i in 1:n_particles +# # Random velocity direction +# theta = 2 * π * rand() +# # Random velocity magnitude +# v_mag = sqrt(2 * desired_ke_per_particle) +# velocities[:, i] = v_mag * [cos(theta), sin(theta)] +# end + +# # Assign velocities to particles +# particles.velocities = velocities + +# # Use appropriate density calculator +# fluid_density_calculator = ContinuityDensity() + +# # Exclude viscosity +# viscosity = nothing + +# # Create the fluid system +# fluid_system = WeaklyCompressibleSPHSystem( +# particles, +# fluid_density_calculator, +# StateEquationCole( +# sound_speed = sound_speed, +# reference_density = fluid_density, +# exponent = 7.0, +# clip_negative_pressure = false +# ), +# SmoothingKernelCubicSpline(), +# smoothing_length, +# viscosity = viscosity, +# acceleration = gravity, +# surface_tension = SurfaceTensionMomentumMorris() +# ) + +# # Define the periodic box matching the simulation domain +# periodic_box = PeriodicBox( +# min_corner = [0.0, 0.0], +# max_corner = [domain_width, domain_height] +# ) + +# # Configure the neighborhood search with the periodic box +# neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box) + +# # Set up the semidiscretization and ODE problem +# semi = Semidiscretization( +# fluid_system, +# neighborhood_search = neighborhood_search +# ) + +# ode = semidiscretize(semi, tspan) + +# # Define callbacks to record kinetic and internal energy +# kinetic_energies = Float64[] +# internal_energies = Float64[] +# times = Float64[] + +# function record_energies!(integrator) +# u = integrator.u +# v = integrator.cache + +# # Compute kinetic energy +# velocities = v.velocities +# masses = particles.masses +# ke = 0.5 * sum(masses .* sum(velocities .^ 2, dims=1)) +# push!(kinetic_energies, ke) + +# # Compute internal energy +# densities = v.densities +# pressures = v.pressures +# ie = sum(pressures .* masses / (fluid_density * sound_speed^2)) +# push!(internal_energies, ie) + +# push!(times, integrator.t) +# return +# end + +# cb = CallbackSet(SavingCallback(record_energies!)) + +# # Time integration settings +# dtmax = 0.25 * particle_spacing / (sound_speed + maximum(abs.(velocities))) + +# # Solve the ODE problem +# sol = solve( +# ode, +# RDPK3SpFSAL35(), +# callback = cb, +# save_everystep = false, +# abstol = 1e-6, +# reltol = 1e-6, +# dtmax = dtmax +# ) + +# # Compute the ratio of kinetic energy to internal energy +# energy_ratios = kinetic_energies ./ internal_energies + +# # Plot the ratio over time +# plot( +# times, energy_ratios, +# xlabel = "Time", +# ylabel = "KE / Internal Energy", +# title = "Parasitic Currents Test", +# legend = false +# ) diff --git a/examples/fluid/rain_flat.jl b/examples/fluid/rain_flat.jl new file mode 100644 index 000000000..e3786efb4 --- /dev/null +++ b/examples/fluid/rain_flat.jl @@ -0,0 +1,130 @@ +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.00004 # x75 visc +# fluid_particle_spacing = 0.0001 # no improvement +boundary_particle_spacing = 0.00005 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 0.02) +fluid_density = 1000.0 +sound_speed = 100 + +# Boundary geometry and initial fluid particle positions +plate_length = 0.01 +plate_width = 0.01 + +sphere_radius = 0.0005 +#terminal_velocity = 75 * 5.8 # match Re +terminal_velocity = 5.8 # real +#terminal_velocity = 1.0 + +boundary_thickness = 4 * boundary_particle_spacing + +# ground floor +plate_size = (plate_length, plate_width, boundary_thickness) + +plate = RectangularShape(boundary_particle_spacing, + round.(Int, plate_size ./ boundary_particle_spacing), + (0.0, 0.0, 0.0), density=fluid_density) + +sphere1_center = (0.005, 0.005, sphere_radius + boundary_thickness) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=VoxelSphere(), + velocity=(0.0, 0.0, -terminal_velocity)) + +# box = RectangularTank(fluid_particle_spacing, (0.3, 0.125, 0.0), (0.35, 0.075, 0.2), fluid_density, +# n_layers=8, spacing_ratio=spacing_ratio, acceleration=(0.0, -gravity, 0.0), state_equation=state_equation, +# faces=(false, false, true, false, true, true), +# min_coordinates=(tank_size[1], 0.0, 0.0)) + +# # move to the end +# for i in axes(end_tank.boundary.coordinates, 2) +# end_tank.boundary.coordinates[:, i] .+= [tank_size[1], 0.0] +# end + +#tank = union(tank, tank_end, tank_side_mz, tank_side_pz, tank_side_mz_end, tank_side_pz_end) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 3.5 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{3}() + +fluid_density_calculator = ContinuityDensity() +# viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +# kinematic viscosity of water at 20C +nu_water = 8.9E-7 + +# Morris has a higher velocity with same viscosity. +# viscosity = ViscosityMorris(nu=75*nu_water) +viscosity = ViscosityAdami(nu=20 * nu_water) + +density_diffusion = DensityDiffusionAntuono(sphere1, delta=0.1) + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, clip_negative_pressure=false) + +fluid_system = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity, + density_diffusion=density_diffusion, + acceleration=(0.0, 0.0, -gravity), + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.0728, + free_surface_threshold=0.6), # 0.55 too many # 0.8 even more + reference_particle_spacing=fluid_particle_spacing, + surface_normal_method=ColorfieldSurfaceNormal(smoothing_kernel, + smoothing_length, + boundary_contact_threshold=1.0)) + +# fluid_system = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator, +# state_equation, smoothing_kernel, +# smoothing_length, viscosity=viscosity, +# density_diffusion=density_diffusion, +# acceleration=(0.0, 0.0, -gravity)) + +# fluid_system = EntropicallyDampedSPHSystem(sphere1, smoothing_kernel, +# smoothing_length, +# sound_speed, +# viscosity=viscosity, +# density_calculator=fluid_density_calculator, +# reference_particle_spacing=fluid_particle_spacing, +# acceleration=(0.0, 0.0, -gravity)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(plate.density, plate.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity) + +boundary_system = BoundarySPHSystem(plate, boundary_model) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.00005, + prefix="rainfall_morris_alpha001_h40micro_nu20_lowerdt") +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +# Limiting of the maximum stepsize is necessary to prevent crashing. +# When particles are approaching a wall in a uniform way, they can be advanced +# with large time steps. Close to the wall, the stepsize has to be reduced drastically. +# Sometimes, the method fails to do so because forces become extremely large when +# fluid particles are very close to boundary particles, and the time integration method +# interprets this as an instability. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + dtmax=1e-4, # Limit stepsize to prevent crashing + dt=1e-8, + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/sphere_surface_tension_2d.jl b/examples/fluid/sphere_surface_tension_2d.jl index 9006ebf0e..1b0b68bae 100644 --- a/examples/fluid/sphere_surface_tension_2d.jl +++ b/examples/fluid/sphere_surface_tension_2d.jl @@ -5,7 +5,9 @@ using OrdinaryDiffEq fluid_density = 1000.0 -particle_spacing = 0.1 +particle_spacing = 0.025 +# Use a higher resolution for a better result +# particle_spacing = 0.025 # Note: Only square shapes will result in a sphere. # Furthermore, changes of the coefficients might be necessary for higher resolutions or larger squares. @@ -20,35 +22,58 @@ state_equation = StateEquationCole(; sound_speed, reference_density=fluid_densit # smoothing_kernel = WendlandC2Kernel{2}() # nu = 0.01 -smoothing_length = 1.0 * particle_spacing -smoothing_kernel = SchoenbergCubicSplineKernel{2}() -nu = 0.025 +smoothing_length = 3.5 * particle_spacing +fluid_smoothing_kernel = WendlandC2Kernel{2}() +# nu = 0.001 # SurfaceTensionMomentumMorris +nu = 0.05 # SurfaceTensionMorris fluid = RectangularShape(particle_spacing, round.(Int, fluid_size ./ particle_spacing), zeros(length(fluid_size)), density=fluid_density) alpha = 8 * nu / (smoothing_length * sound_speed) source_terms = SourceTermDamping(; damping_coefficient=0.5) -fluid_system = WeaklyCompressibleSPHSystem(fluid, SummationDensity(), - state_equation, smoothing_kernel, +# fluid_system = WeaklyCompressibleSPHSystem(fluid, SummationDensity(), +# state_equation, fluid_smoothing_kernel, +# smoothing_length, +# reference_particle_spacing=particle_spacing, +# viscosity=ArtificialViscosityMonaghan(alpha=alpha, +# beta=0.0), +# surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.02), +# correction=AkinciFreeSurfaceCorrection(fluid_density), +# source_terms=source_terms) + +fluid_system = EntropicallyDampedSPHSystem(fluid, fluid_smoothing_kernel, smoothing_length, - viscosity=ArtificialViscosityMonaghan(alpha=alpha, - beta=0.0), - surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.02), - correction=AkinciFreeSurfaceCorrection(fluid_density), - source_terms=source_terms) + sound_speed, + viscosity=ViscosityMorris(nu=nu), + density_calculator=ContinuityDensity(), + reference_particle_spacing=particle_spacing, + acceleration=zeros(length(fluid_size)), + surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1), + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=5 * + 0.0728)) + +# fluid_system = EntropicallyDampedSPHSystem(fluid, fluid_smoothing_kernel, +# smoothing_length, +# sound_speed, +# viscosity=ViscosityMorris(nu=nu), +# density_calculator=ContinuityDensity(), +# reference_particle_spacing=particle_spacing, +# acceleration=zeros(length(fluid_size)), +# surface_normal_method=ColorfieldSurfaceNormal(), +# surface_tension=SurfaceTensionMomentumMorris(surface_tension_coefficient=1.0)) # ========================================================================================== # ==== Simulation semi = Semidiscretization(fluid_system) -tspan = (0.0, 3.0) +tspan = (0.0, 50.0) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) # For overwriting via `trixi_include` -saving_callback = SolutionSavingCallback(dt=0.02) +saving_callback = SolutionSavingCallback(dt=1.0) stepsize_callback = StepsizeCallback(cfl=1.0) diff --git a/examples/fluid/sphere_surface_tension_3d.jl b/examples/fluid/sphere_surface_tension_3d.jl index 58a506294..aff151fc0 100644 --- a/examples/fluid/sphere_surface_tension_3d.jl +++ b/examples/fluid/sphere_surface_tension_3d.jl @@ -5,7 +5,7 @@ using OrdinaryDiffEq fluid_density = 1000.0 -particle_spacing = 0.1 +particle_spacing = 0.15 fluid_size = (0.9, 0.9, 0.9) sound_speed = 20.0 @@ -13,13 +13,12 @@ sound_speed = 20.0 # For all surface tension simulations, we need a compact support of `2 * particle_spacing` smoothing_length = 1.0 * particle_spacing -nu = 0.01 +nu = 0.04 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "sphere_surface_tension_2d.jl"), - dt=0.1, cfl=1.2, surface_tension_coefficient=0.1, - tspan=(0.0, 10.0), nu=nu, - alpha=10 * nu / (smoothing_length * sound_speed), - smoothing_kernel=SchoenbergCubicSplineKernel{3}(), + surface_tension_coefficient=0.5, dt=0.25, + tspan=(0.0, 100.0), nu=nu, smoothing_length=3.0 * particle_spacing, + fluid_smoothing_kernel=WendlandC2Kernel{3}(), particle_spacing=particle_spacing, sound_speed=sound_speed, fluid_density=fluid_density, fluid_size=fluid_size) diff --git a/examples/fluid/sphere_surface_tension_wall_2d.jl b/examples/fluid/sphere_surface_tension_wall_2d.jl index 3e1cbdb80..4a4eef531 100644 --- a/examples/fluid/sphere_surface_tension_wall_2d.jl +++ b/examples/fluid/sphere_surface_tension_wall_2d.jl @@ -1,5 +1,3 @@ -# In this example we try to approach the static shape of a water droplet on a horizontal plane. -# The shape of a static droplet can be calculated from the Young-Laplace equation. using TrixiParticles using OrdinaryDiffEq @@ -9,6 +7,7 @@ fluid_particle_spacing = 0.0025 # ========================================================================================== # ==== Experiment Setup +gravity = 9.81 tspan = (0.0, 0.5) # Boundary geometry and initial fluid particle positions @@ -17,6 +16,9 @@ tank_size = (0.5, 0.1) fluid_density = 1000.0 sound_speed = 120.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + sphere_radius = 0.05 sphere1_center = (0.25, sphere_radius) @@ -25,7 +27,8 @@ sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, # ========================================================================================== # ==== Fluid -fluid_smoothing_length = 1.0 * fluid_particle_spacing +fluid_smoothing_length = 1.0 * fluid_particle_spacing - eps() +fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}() # For perfect wetting # nu = 0.0005 @@ -35,10 +38,22 @@ nu = 0.001 alpha = 8 * nu / (fluid_smoothing_length * sound_speed) # `adhesion_coefficient = 1.0` and `surface_tension_coefficient = 0.01` for perfect wetting # `adhesion_coefficient = 0.001` and `surface_tension_coefficient = 2.0` for no wetting + +viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) +sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, ContinuityDensity(), + state_equation, fluid_smoothing_kernel, + fluid_smoothing_length, + viscosity=viscosity, + acceleration=(0.0, -gravity), + surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=2.0), + correction=AkinciFreeSurfaceCorrection(fluid_density), + reference_particle_spacing=fluid_particle_spacing) + trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "falling_water_spheres_2d.jl"), sphere=nothing, sphere1=sphere1, adhesion_coefficient=0.001, - wall_viscosity=4.0 * nu, surface_tension_coefficient=2.0, alpha=alpha, + wall_viscosity=4.0 * nu, surface_tension_coefficient=0.9, alpha=alpha, sound_speed=sound_speed, fluid_density=fluid_density, nu=nu, fluid_particle_spacing=fluid_particle_spacing, tspan=tspan, - tank_size=tank_size) + tank_size=tank_size, fluid_smoothing_length=fluid_smoothing_length, + sphere_surface_tension=sphere_surface_tension) diff --git a/examples/fluid/static_spheres.jl b/examples/fluid/static_spheres.jl new file mode 100644 index 000000000..e5fb8683b --- /dev/null +++ b/examples/fluid/static_spheres.jl @@ -0,0 +1,98 @@ +# In this example two circles of water drop to the floor demonstrating the difference +# between the behavior with and without surface tension modelling. +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.005 + +boundary_layers = 3 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 0.0 +tspan = (0.0, 0.3) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.0, 0.0) +tank_size = (2.5, 0.5) + +fluid_density = 1000.0 +sound_speed = 100 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation=state_equation) + +sphere_radius = 0.05 + +sphere1_center = (0.5, 0.5) +sphere2_center = (1.5, 0.5) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=RoundSphere(), velocity=(0.0, 0.0)) +sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center, + fluid_density, sphere_type=RoundSphere(), velocity=(0.0, 0.0)) + +# ========================================================================================== +# ==== Fluid +fluid_smoothing_length = 3.5 * fluid_particle_spacing +fluid_smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() + +nu = 0.005 +alpha = 8 * nu / (fluid_smoothing_length * sound_speed) +viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) +# density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1) + +sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel, + fluid_smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=ContinuityDensity(), + reference_particle_spacing=fluid_particle_spacing, + acceleration=(0.0, -gravity), + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.0728)) + +sphere = EntropicallyDampedSPHSystem(sphere2, fluid_smoothing_kernel, + fluid_smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=ContinuityDensity(), + reference_particle_spacing=fluid_particle_spacing, + acceleration=(0.0, -gravity), + surface_normal_method=ColorfieldSurfaceNormal(fluid_smoothing_kernel, + fluid_smoothing_length)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +wall_viscosity = nu +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, fluid_smoothing_length, + viscosity=ViscosityAdami(nu=wall_viscosity)) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", + prefix="static", write_meta_data=true) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-7, # Default abstol is 1e-6 + reltol=1e-4, # Default reltol is 1e-3 + dt=1e-6, + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/waterfall.jl b/examples/fluid/waterfall.jl new file mode 100644 index 000000000..516ffc6c7 --- /dev/null +++ b/examples/fluid/waterfall.jl @@ -0,0 +1,127 @@ +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.008 + +# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +boundary_layers = 7 +spacing_ratio = 1 + +boundary_particle_spacing = fluid_particle_spacing / spacing_ratio + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 5.0) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.2, 0.1, 0.2) +tank_size = (0.3, 0.15, 0.20) +end_size = (0.3, 0.05, 0.20) + +fluid_density = 1000.0 +sound_speed = 40 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, background_pressure=1000) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=4, spacing_ratio=spacing_ratio, + acceleration=(0.0, -gravity, 0.0), state_equation=state_equation, + faces=(true, false, true, false, true, true)) + +end_tank = RectangularTank(fluid_particle_spacing, (0.0, 0.0, 0.0), end_size, fluid_density, + n_layers=8, spacing_ratio=spacing_ratio, + acceleration=(0.0, -gravity, 0.0), state_equation=state_equation, + faces=(false, false, true, false, true, true), + min_coordinates=(tank_size[1], 0.0, 0.0)) + +# # move to the end +# for i in axes(end_tank.boundary.coordinates, 2) +# end_tank.boundary.coordinates[:, i] .+= [tank_size[1], 0.0] +# end + +# tank = union(tank, end_tank) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 3.5 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{3}() + +fluid_density_calculator = ContinuityDensity() +# viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +nu_water = 8.9E-7 + +# Morris has a higher velocity with same viscosity. +# viscosity = ViscosityMorris(nu=100*nu_water) +viscosity = ViscosityAdami(nu=100 * nu_water) + +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) + +# fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, +# state_equation, smoothing_kernel, +# smoothing_length, viscosity=viscosity, +# density_diffusion=density_diffusion, +# acceleration=(0.0, -gravity, 0.0)) + +fluid_system = EntropicallyDampedSPHSystem(tank.fluid, smoothing_kernel, + smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=ContinuityDensity(), + acceleration=(0.0, -gravity, 0.0), + # surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05), + reference_particle_spacing=fluid_particle_spacing, + buffer_size=1) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +boundary_model2 = BoundaryModelDummyParticles(end_tank.boundary.density, + end_tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity) + +boundary_system2 = BoundarySPHSystem(end_tank.boundary, boundary_model2) + +outflow = BoundaryZone(; plane=([0.4, -0.2, -0.1], [0.4, -0.2, 0.3], [0.8, -0.2, 0.3]), + plane_normal=[0.0, -1.0, 0.0], open_boundary_layers=1, + density=2 * eps(), particle_spacing=fluid_particle_spacing, + boundary_type=:outflow) + +open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, + boundary_model=BasicOutlet()) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, boundary_system, boundary_system2, + open_boundary_out) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.01, prefix="lower_visc") +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +# Limiting of the maximum stepsize is necessary to prevent crashing. +# When particles are approaching a wall in a uniform way, they can be advanced +# with large time steps. Close to the wall, the stepsize has to be reduced drastically. +# Sometimes, the method fails to do so because forces become extremely large when +# fluid particles are very close to boundary particles, and the time integration method +# interprets this as an instability. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + dtmax=1e-2, # Limit stepsize to prevent crashing + dt=1e-6, + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/waterfallv2.jl b/examples/fluid/waterfallv2.jl new file mode 100644 index 000000000..021392510 --- /dev/null +++ b/examples/fluid/waterfallv2.jl @@ -0,0 +1,171 @@ +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.0025 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 5.0) +fluid_density = 1000.0 +sound_speed = 80 + +# Boundary geometry and initial fluid particle positions +canal_width = 0.1 +canal_height = 0.1 +canal_length = 0.1 +outflow_length = 0.3 + +initial_fluid_size = (canal_length, canal_height, canal_width) +boundary_thickness = 4 * fluid_particle_spacing + +# ground floor +plate = (canal_length + outflow_length, boundary_thickness, canal_width) +# end to close the canal +plate_end = (boundary_thickness, canal_height + boundary_thickness, canal_width) + +plate_side = (1.5 * canal_length, canal_height + boundary_thickness, boundary_thickness) +# goes from the end of plate_side to the end of the plate +plate_side_end = (plate[1] - 1.5 * canal_length + boundary_thickness, + 0.5 * canal_height + boundary_thickness, boundary_thickness) + +fluid = RectangularShape(fluid_particle_spacing, + round.(Int, initial_fluid_size ./ fluid_particle_spacing), + (0.0, plate[2], 0.0), density=fluid_density) + +tank = RectangularShape(fluid_particle_spacing, + round.(Int, plate ./ fluid_particle_spacing), + (0.0, 0.0, 0.0), density=fluid_density) + +tank_end = RectangularShape(fluid_particle_spacing, + round.(Int, plate_end ./ fluid_particle_spacing), + (-plate_end[1], 0.0, 0.0), density=fluid_density) + +tank_side_pz = RectangularShape(fluid_particle_spacing, + round.(Int, plate_side ./ fluid_particle_spacing), + (-plate_end[1], 0.0, initial_fluid_size[1]), + density=fluid_density) + +tank_side_mz = RectangularShape(fluid_particle_spacing, + round.(Int, plate_side ./ fluid_particle_spacing), + (-plate_end[1], 0.0, -plate_end[1]), density=fluid_density) + +tank_side_pz_end = RectangularShape(fluid_particle_spacing, + round.(Int, plate_side_end ./ fluid_particle_spacing), + (plate_side[1] - boundary_thickness, 0.0, + initial_fluid_size[1]), density=fluid_density) + +tank_side_mz_end = RectangularShape(fluid_particle_spacing, + round.(Int, plate_side_end ./ fluid_particle_spacing), + (plate_side[1] - boundary_thickness, 0.0, + -plate_end[1]), density=fluid_density) + +# box = RectangularTank(fluid_particle_spacing, (0.3, 0.125, 0.0), (0.35, 0.075, 0.2), fluid_density, +# n_layers=8, spacing_ratio=spacing_ratio, acceleration=(0.0, -gravity, 0.0), state_equation=state_equation, +# faces=(false, false, true, false, true, true), +# min_coordinates=(tank_size[1], 0.0, 0.0)) + +# # move to the end +# for i in axes(end_tank.boundary.coordinates, 2) +# end_tank.boundary.coordinates[:, i] .+= [tank_size[1], 0.0] +# end + +tank = union(tank, tank_end, tank_side_mz, tank_side_pz, tank_side_mz_end, tank_side_pz_end) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 3.25 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{3}() + +fluid_density_calculator = ContinuityDensity() +# viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +nu_water = 8.9E-7 + +# Morris has a higher velocity with same viscosity. +# viscosity = ViscosityMorris(nu=100*nu_water) +viscosity = ViscosityAdami(nu=50 * nu_water) + +# density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, clip_negative_pressure=false) + +# fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, +# state_equation, smoothing_kernel, +# smoothing_length, viscosity=viscosity, +# density_diffusion=density_diffusion, +# acceleration=(0.0, -gravity, 0.0)) + +fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity, + density_diffusion=density_diffusion, + acceleration=(0.0, -gravity, 0.0), + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=50 * + 0.0728, + free_surface_threshold=0.7), + reference_particle_spacing=fluid_particle_spacing, + surface_normal_method=ColorfieldSurfaceNormal(smoothing_kernel, + smoothing_length, + boundary_contact_threshold=0.1)) + +# fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, +# smoothing_length, +# sound_speed, viscosity=viscosity, +# density_calculator=ContinuityDensity(), +# acceleration=(0.0, -gravity, 0.0), +# # surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05), +# reference_particle_spacing=fluid_particle_spacing) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.density, tank.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity) + +boundary_system = BoundarySPHSystem(tank, boundary_model) + +# boundary_model2 = BoundaryModelDummyParticles(end_tank.boundary.density, +# end_tank.boundary.mass, +# state_equation=state_equation, +# boundary_density_calculator, +# smoothing_kernel, smoothing_length, +# viscosity=viscosity) + +# boundary_system2 = BoundarySPHSystem(end_tank.boundary, boundary_model2) + +# outflow = BoundaryZone(; plane=([0.4, -0.2, -0.1], [0.4, -0.2, 0.3], [0.8, -0.2, 0.3]), +# plane_normal=[0.0, -1.0, 0.0], open_boundary_layers=1, +# density=2 * eps(), particle_spacing=fluid_particle_spacing, +# boundary_type=:outflow) + +# open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, +# boundary_model=BasicOutlet()) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.01, prefix="waterfall_v2_less_height") +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +# Limiting of the maximum stepsize is necessary to prevent crashing. +# When particles are approaching a wall in a uniform way, they can be advanced +# with large time steps. Close to the wall, the stepsize has to be reduced drastically. +# Sometimes, the method fails to do so because forces become extremely large when +# fluid particles are very close to boundary particles, and the time integration method +# interprets this as an instability. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + dtmax=1e-2, # Limit stepsize to prevent crashing + dt=1e-6, + save_everystep=false, callback=callbacks); diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index d85096860..8aeffa5ee 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -87,6 +87,8 @@ export kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure, max_density, min_density, avg_density export interpolate_line, interpolate_point, interpolate_plane_3d, interpolate_plane_2d, interpolate_plane_2d_vtk -export SurfaceTensionAkinci, CohesionForceAkinci +export SurfaceTensionAkinci, CohesionForceAkinci, SurfaceTensionMorris, + SurfaceTensionMomentumMorris, HuberContactModel +export ColorfieldSurfaceNormal, StaticNormals end # module diff --git a/src/general/corrections.jl b/src/general/corrections.jl index b1c66a06a..f583de7ea 100644 --- a/src/general/corrections.jl +++ b/src/general/corrections.jl @@ -413,3 +413,92 @@ function correction_matrix_inversion_step!(corr_matrix, system) return corr_matrix end + +create_cache_correction(correction, density, NDIMS, nparticles) = (;) + +function create_cache_correction(::ShepardKernelCorrection, density, NDIMS, n_particles) + return (; kernel_correction_coefficient=similar(density)) +end + +function create_cache_correction(::KernelCorrection, density, NDIMS, n_particles) + dw_gamma = Array{Float64}(undef, NDIMS, n_particles) + return (; kernel_correction_coefficient=similar(density), dw_gamma) +end + +function create_cache_correction(::Union{GradientCorrection, BlendedGradientCorrection}, + density, + NDIMS, n_particles) + correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) + return (; correction_matrix) +end + +function create_cache_correction(::MixedKernelGradientCorrection, density, NDIMS, + n_particles) + dw_gamma = Array{Float64}(undef, NDIMS, n_particles) + correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) + + return (; kernel_correction_coefficient=similar(density), dw_gamma, correction_matrix) +end + +function kernel_correct_density!(system::FluidSystem, v, u, v_ode, u_ode, + semi, correction, density_calculator) + return system +end + +function kernel_correct_density!(system::FluidSystem, v, u, v_ode, u_ode, + semi, corr::ShepardKernelCorrection, ::SummationDensity) + system.cache.density ./= system.cache.kernel_correction_coefficient +end + +function compute_gradient_correction_matrix!(correction, system::FluidSystem, u, + v_ode, u_ode, semi) + return system +end + +function compute_gradient_correction_matrix!(corr::Union{GradientCorrection, + BlendedGradientCorrection, + MixedKernelGradientCorrection}, + system::FluidSystem, u, + v_ode, u_ode, semi) + (; cache, correction, smoothing_kernel, smoothing_length) = system + (; correction_matrix) = cache + + system_coords = current_coordinates(u, system) + + compute_gradient_correction_matrix!(correction_matrix, system, system_coords, + v_ode, u_ode, semi, correction, smoothing_length, + smoothing_kernel) +end + +function reinit_density!(vu_ode, semi) + v_ode, u_ode = vu_ode.x + + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + reinit_density!(system, v, u, v_ode, u_ode, semi) + end + + return vu_ode +end + +function reinit_density!(system::FluidSystem, v, u, v_ode, u_ode, semi) + # Compute density with `SummationDensity` and store the result in `v`, + # overwriting the previous integrated density. + summation_density!(system, semi, u, u_ode, v[end, :]) + + # Apply `ShepardKernelCorrection` + kernel_correction_coefficient = zeros(size(v[end, :])) + compute_shepard_coeff!(system, current_coordinates(u, system), v_ode, u_ode, semi, + kernel_correction_coefficient) + v[end, :] ./= kernel_correction_coefficient + + compute_pressure!(system, v) + + return system +end + +function reinit_density!(system, v, u, v_ode, u_ode, semi) + return system +end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 600d573d5..032bdab60 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -407,7 +407,7 @@ end function calculate_dt(v_ode, u_ode, cfl_number, semi::Semidiscretization) (; systems) = semi - return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system), systems) + return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system, semi), systems) end function drift!(du_ode, v_ode, u_ode, semi, t) @@ -818,13 +818,49 @@ function update!(neighborhood_search, system::GPUSystem, x, y; points_moving=(tr end function check_configuration(systems) + min_color = 0 + max_color = 0 + no_fluid_and_bnd_systems = 0 + uses_surface_tension_model = false + foreach_system(systems) do system check_configuration(system, systems) + + if system isa FluidSystem && !isnothing(system.surface_tension) + uses_surface_tension_model = true + end + + if system isa FluidSystem || system isa BoundarySPHSystem + no_fluid_and_bnd_systems += 1 + if system.color < min_color + min_color = system.color + end + + if system.color > max_color + max_color = system.color + end + end + end + + if max_color == 0 && min_color == 0 && no_fluid_and_bnd_systems > 1 && + uses_surface_tension_model + throw(ArgumentError("If a surface tension model is used the values of at least one system needs to have a color different than 0.")) end end check_configuration(system, systems) = nothing +function check_configuration(fluid_system::FluidSystem, systems) + if !isnothing(fluid_system.surface_tension) + foreach_system(systems) do neighbor + if neighbor isa FluidSystem && isnothing(fluid_system.surface_tension) && + isnothing(fluid_system.surface_normal_method) + throw(ArgumentError("All `FluidSystem` need to use a surface tension model or a surface normal method.")) + end + end + end +end + function check_configuration(boundary_system::BoundarySPHSystem, systems) (; boundary_model) = boundary_system diff --git a/src/general/system.jl b/src/general/system.jl index 48d556eba..b87441742 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -151,3 +151,17 @@ end # Only for systems requiring a mandatory callback reset_callback_flag!(system) = system + +# Assuming a constant particle spacing one can calculate the number of neighbors within the +# compact support for an undisturbed particle distribution. +function ideal_neighbor_count(::Val{D}, particle_spacing, compact_support) where {D} + throw(ArgumentError("Unsupported dimension: $D")) +end + +@inline function ideal_neighbor_count(::Val{2}, particle_spacing, compact_support) + return floor(Int, pi * compact_support^2 / particle_spacing^2) +end + +@inline @fastpow function ideal_neighbor_count(::Val{3}, particle_spacing, compact_support) + return floor(Int, 4.0 / 3.0 * pi * compact_support^3 / particle_spacing^3) +end diff --git a/src/schemes/boundary/boundary.jl b/src/schemes/boundary/boundary.jl index 6d50eec9b..70ef62bb3 100644 --- a/src/schemes/boundary/boundary.jl +++ b/src/schemes/boundary/boundary.jl @@ -5,3 +5,5 @@ include("open_boundary/method_of_characteristics.jl") include("open_boundary/system.jl") # Monaghan-Kajtar repulsive boundary particles require the `BoundarySPHSystem` # and the `TotalLagrangianSPHSystem` and are therefore included later. + +@inline Base.ndims(boundary_model::BoundaryModelDummyParticles) = ndims(boundary_model.smoothing_kernel) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 8b7e507ed..eadd37b91 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -2,7 +2,8 @@ BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, density_calculator, smoothing_kernel, smoothing_length; viscosity=nothing, - state_equation=nothing, correction=nothing) + state_equation=nothing, correction=nothing, + reference_particle_spacing=0.0) Boundary model for `BoundarySPHSystem`. @@ -16,12 +17,12 @@ Boundary model for `BoundarySPHSystem`. - `smoothing_length`: Smoothing length should be the same as for the adjacent fluid system. # Keywords -- `state_equation`: This should be the same as for the adjacent fluid system - (see e.g. [`StateEquationCole`](@ref)). -- `correction`: Correction method of the adjacent fluid system (see [Corrections](@ref corrections)). -- `viscosity`: Slip (default) or no-slip condition. See description below for further - information. - +- `state_equation`: This should be the same as for the adjacent fluid system + (see e.g. [`StateEquationCole`](@ref)). +- `correction`: Correction method of the adjacent fluid system (see [Corrections](@ref corrections)). +- `viscosity`: Slip (default) or no-slip condition. See description below for further + information. +- `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, which currently is only needed when using surface tension. # Examples ```jldoctest; output = false, setup = :(densities = [1.0, 2.0, 3.0]; masses = [0.1, 0.2, 0.3]; smoothing_kernel = SchoenbergCubicSplineKernel{2}(); smoothing_length = 0.1) # Free-slip condition @@ -38,15 +39,16 @@ BoundaryModelDummyParticles(AdamiPressureExtrapolation, ViscosityAdami) ``` """ struct BoundaryModelDummyParticles{DC, ELTYPE <: Real, VECTOR, SE, K, V, COR, C} - pressure :: VECTOR # Vector{ELTYPE} - hydrodynamic_mass :: VECTOR # Vector{ELTYPE} - state_equation :: SE - density_calculator :: DC - smoothing_kernel :: K - smoothing_length :: ELTYPE - viscosity :: V - correction :: COR - cache :: C + pressure :: VECTOR # Vector{ELTYPE} + hydrodynamic_mass :: VECTOR # Vector{ELTYPE} + state_equation :: SE + density_calculator :: DC + smoothing_kernel :: K + smoothing_length :: ELTYPE + ideal_neighbor_count :: Int + viscosity :: V + correction :: COR + cache :: C end # The default constructor needs to be accessible for Adapt.jl to work with this struct. @@ -54,19 +56,34 @@ end function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, density_calculator, smoothing_kernel, smoothing_length; viscosity=nothing, - state_equation=nothing, correction=nothing) + state_equation=nothing, correction=nothing, + reference_particle_spacing=0.0) pressure = initial_boundary_pressure(initial_density, density_calculator, state_equation) NDIMS = ndims(smoothing_kernel) + ELTYPE = eltype(smoothing_length) n_particles = length(initial_density) cache = (; create_cache_model(viscosity, n_particles, NDIMS)..., create_cache_model(initial_density, density_calculator)..., - create_cache_model(correction, initial_density, NDIMS, n_particles)...) + create_cache_model(correction, initial_density, NDIMS, n_particles)..., + (; colorfield_bnd=zeros(ELTYPE, n_particles), + colorfield=zeros(ELTYPE, n_particles), + neighbor_count=zeros(ELTYPE, n_particles))...) + + # If the `reference_density_spacing` is set calculate the `ideal_neighbor_count` + ideal_neighbor_count_ = 0 + if reference_particle_spacing > 0.0 + ideal_neighbor_count_ = ideal_neighbor_count(Val(ndims(boundary_model)), + reference_particle_spacing, + compact_support(smoothing_kernel, + smoothing_length)) + end return BoundaryModelDummyParticles(pressure, hydrodynamic_mass, state_equation, density_calculator, smoothing_kernel, - smoothing_length, viscosity, correction, cache) + smoothing_length, ideal_neighbor_count_, viscosity, + correction, cache) end @doc raw""" diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index a3f47f611..aa01b5505 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -9,36 +9,53 @@ The interaction between fluid and boundary particles is specified by the boundar - `boundary_model`: Boundary model (see [Boundary Models](@ref boundary_models)) # Keyword Arguments -- `movement`: For moving boundaries, a [`BoundaryMovement`](@ref) can be passed. -- `adhesion_coefficient`: Coefficient specifying the adhesion of a fluid to the surface. - Note: currently it is assumed that all fluids have the same adhesion coefficient. +- `movement` : For moving boundaries, a [`BoundaryMovement`](@ref) can be passed. +- `adhesion_coefficient` : Coefficient specifying the adhesion of a fluid to the surface. + Note: currently it is assumed that all fluids have the same adhesion coefficient. +- `static_contact_angle` : The static contact or also sometimes equilibrium contact angle is + used by the wetting model in connection with a surface tension model + to model the correct contact mechanics between a fluid and a solid. + The contact angle should be provided in degrees. + Note: currently it is assumed that all fluids have the same adhesion coefficient. +- `surface_normal_method`: The surface normal method is used to calculate the surface normals for the boundary. + Currently available are the following models: `StaticNormals` +- `color` : The color is used to differentiate different phases/materials/fluids. + Typically all boundaries should have the same value. (By default boundaries::color=0 fluids::color=1) + """ -struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, IC, CO, M, IM, +struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, IC, CO, M, IM, SRFN, CA} <: BoundarySystem{NDIMS, IC} - initial_condition :: IC - coordinates :: CO # Array{ELTYPE, 2} - boundary_model :: BM - movement :: M - ismoving :: IM # Ref{Bool} (to make a mutable field compatible with GPUs) - adhesion_coefficient :: ELTYPE - cache :: CA - buffer :: Nothing + initial_condition :: IC + coordinates :: CO # Array{ELTYPE, 2} + boundary_model :: BM + movement :: M + ismoving :: IM # Ref{Bool} (to make a mutable field compatible with GPUs) + adhesion_coefficient :: ELTYPE + static_contact_angle :: ELTYPE # sometimes named equilibrium contact angle + surface_normal_method :: SRFN + color :: Int + cache :: CA + buffer :: Nothing # This constructor is necessary for Adapt.jl to work with this struct. # See the comments in general/gpu.jl for more details. function BoundarySPHSystem(initial_condition, coordinates, boundary_model, movement, - ismoving, adhesion_coefficient, cache, buffer) + ismoving, adhesion_coefficient, static_contact_angle, + surface_normal_method, color_value, cache, buffer) ELTYPE = eltype(coordinates) new{typeof(boundary_model), size(coordinates, 1), ELTYPE, typeof(initial_condition), typeof(coordinates), typeof(movement), typeof(ismoving), + typeof(surface_normal_method), typeof(cache)}(initial_condition, coordinates, boundary_model, movement, - ismoving, adhesion_coefficient, cache, buffer) + ismoving, adhesion_coefficient, static_contact_angle, + surface_normal_method, color_value, cache, buffer) end end function BoundarySPHSystem(initial_condition, model; movement=nothing, - adhesion_coefficient=0.0) + adhesion_coefficient=0.0, static_contact_angle=0.0, + surface_normal_method=nothing, color_value=0) coordinates = copy(initial_condition.coordinates) ismoving = Ref(!isnothing(movement)) @@ -52,9 +69,46 @@ function BoundarySPHSystem(initial_condition, model; movement=nothing, movement.moving_particles .= collect(1:nparticles(initial_condition)) end + if static_contact_angle < 0 || static_contact_angle > 180 + throw(ArgumentError("The `static_contact_angle` must be between 0 and 180.")) + end + + # convert degrees to radians + static_contact_rad = static_contact_angle * pi / 180 + # Because of dispatches boundary model needs to be first! return BoundarySPHSystem(initial_condition, coordinates, model, movement, - ismoving, adhesion_coefficient, cache, nothing) + ismoving, adhesion_coefficient, static_contact_rad, + surface_normal_method, color_value, cache, nothing) +end + +function Base.show(io::IO, system::BoundarySPHSystem) + @nospecialize system # reduce precompilation time + + print(io, "BoundarySPHSystem{", ndims(system), "}(") + print(io, system.boundary_model) + print(io, ", ", system.movement) + print(io, ", ", system.adhesion_coefficient) + print(io, ", ", system.color) + print(io, ") with ", nparticles(system), " particles") +end + +function Base.show(io::IO, ::MIME"text/plain", system::BoundarySPHSystem) + @nospecialize system # reduce precompilation time + + if get(io, :compact, false) + show(io, system) + else + summary_header(io, "BoundarySPHSystem{$(ndims(system))}") + summary_line(io, "#particles", nparticles(system)) + summary_line(io, "boundary model", system.boundary_model) + summary_line(io, "movement function", + isnothing(system.movement) ? "nothing" : + string(system.movement.movement_function)) + summary_line(io, "adhesion coefficient", system.adhesion_coefficient) + summary_line(io, "color", system.color) + summary_footer(io) + end end """ @@ -162,33 +216,6 @@ function create_cache_boundary(::BoundaryMovement, initial_condition) return (; velocity, acceleration, initial_coordinates) end -function Base.show(io::IO, system::BoundarySPHSystem) - @nospecialize system # reduce precompilation time - - print(io, "BoundarySPHSystem{", ndims(system), "}(") - print(io, system.boundary_model) - print(io, ", ", system.movement) - print(io, ", ", system.adhesion_coefficient) - print(io, ") with ", nparticles(system), " particles") -end - -function Base.show(io::IO, ::MIME"text/plain", system::BoundarySPHSystem) - @nospecialize system # reduce precompilation time - - if get(io, :compact, false) - show(io, system) - else - summary_header(io, "BoundarySPHSystem{$(ndims(system))}") - summary_line(io, "#particles", nparticles(system)) - summary_line(io, "boundary model", system.boundary_model) - summary_line(io, "movement function", - isnothing(system.movement) ? "nothing" : - string(system.movement.movement_function)) - summary_line(io, "adhesion coefficient", system.adhesion_coefficient) - summary_footer(io) - end -end - timer_name(::Union{BoundarySPHSystem, BoundaryDEMSystem}) = "boundary" @inline function Base.eltype(system::Union{BoundarySPHSystem, BoundaryDEMSystem}) @@ -336,6 +363,11 @@ end function update_quantities!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi, t) (; boundary_model) = system + (; colorfield, colorfield_bnd) = boundary_model.cache + + # Reset to the constant boundary interpolated color values + colorfield .= colorfield_bnd + update_density!(boundary_model, system, v, u, v_ode, u_ode, semi) @@ -395,6 +427,35 @@ end # viscosity model has to be used. @inline viscosity_model(system::BoundarySPHSystem, neighbor_system::FluidSystem) = neighbor_system.viscosity -function calculate_dt(v_ode, u_ode, cfl_number, system::BoundarySystem) +function calculate_dt(v_ode, u_ode, cfl_number, system::BoundarySystem, semi) return Inf end + +function initialize!(system::BoundarySPHSystem, neighborhood_search) + initialize_colorfield!(system, system.boundary_model, neighborhood_search) + return system +end + +function initialize_colorfield!(system, boundary_model, neighborhood_search) + return system +end + +function initialize_colorfield!(system, ::BoundaryModelDummyParticles, neighborhood_search) + system_coords = system.coordinates + (; smoothing_kernel, smoothing_length) = system.boundary_model + + foreach_point_neighbor(system, system, + system_coords, system_coords, + neighborhood_search, + points=eachparticle(system)) do particle, neighbor, pos_diff, + distance + system.boundary_model.cache.colorfield_bnd[particle] += system.initial_condition.mass[particle] / + system.initial_condition.density[particle] * + system.color * + kernel(smoothing_kernel, + distance, + smoothing_length) + system.boundary_model.cache.neighbor_count[particle] += 1 + end + return system +end diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 3644c3630..c5514fe97 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -8,6 +8,9 @@ function interact!(dv, v_particle_system, u_particle_system, system_coords = current_coordinates(u_particle_system, particle_system) neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + surface_tension_a = surface_tension_model(particle_system) + surface_tension_b = surface_tension_model(neighbor_system) + # Loop over all pairs of particles and neighbors within the kernel cutoff. foreach_point_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords, @@ -15,8 +18,8 @@ function interact!(dv, v_particle_system, u_particle_system, # Only consider particles with a distance > 0. distance < sqrt(eps()) && return - rho_a = particle_density(v_particle_system, particle_system, particle) - rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + rho_a = @inbounds particle_density(v_particle_system, particle_system, particle) + rho_b = @inbounds particle_density(v_neighbor_system, neighbor_system, neighbor) p_a = particle_pressure(v_particle_system, particle_system, particle) p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) @@ -27,8 +30,8 @@ function interact!(dv, v_particle_system, u_particle_system, # Note that the return value is zero when not using EDAC with TVF. p_avg = average_pressure(particle_system, particle) - m_a = hydrodynamic_mass(particle_system, particle) - m_b = hydrodynamic_mass(neighbor_system, neighbor) + m_a = @inbounds hydrodynamic_mass(particle_system, particle) + m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance) @@ -48,8 +51,19 @@ function interact!(dv, v_particle_system, u_particle_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) + # TODO: only applies to fluids with the same color + dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b, + particle_system, neighbor_system, + particle, neighbor, pos_diff, distance, + rho_a, rho_b, grad_kernel) + + dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system, + particle, neighbor, pos_diff, distance) + for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_convection[i] + @inbounds dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + + dv_convection[i] + dv_surface_tension[i] + + dv_adhesion[i] end v_diff = current_velocity(v_particle_system, particle_system, particle) - @@ -66,12 +80,21 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system, grad_kernel) end + for particle in each_moving_particle(particle_system) + F = contact_force(particle_system.surface_tension.contact_model, + particle_system, particle) + for i in 1:ndims(particle_system) + @inbounds dv[i, particle] += F[i] + end + end + return dv end +@inline -@inline function pressure_evolution!(dv, particle_system, v_diff, grad_kernel, particle, - pos_diff, distance, sound_speed, m_a, m_b, - p_a, p_b, rho_a, rho_b) +function pressure_evolution!(dv, particle_system, v_diff, grad_kernel, particle, + pos_diff, distance, sound_speed, m_a, m_b, + p_a, p_b, rho_a, rho_b) (; smoothing_length) = particle_system volume_a = m_a / rho_a @@ -109,19 +132,20 @@ end return dv end -# We need a separate method for EDAC since the density is stored in `v[end-1,:]`. -@inline function continuity_equation!(dv, density_calculator::ContinuityDensity, - vdiff, particle, m_b, rho_a, rho_b, - particle_system::EntropicallyDampedSPHSystem, - grad_kernel) +@inline# We need a separate method for EDAC since the density is stored in `v[end-1,:]`. +function continuity_equation!(dv, density_calculator::ContinuityDensity, + vdiff, particle, m_b, rho_a, rho_b, + particle_system::EntropicallyDampedSPHSystem, + grad_kernel) dv[end - 1, particle] += rho_a / rho_b * m_b * dot(vdiff, grad_kernel) return dv end +@inline -@inline function continuity_equation!(dv, density_calculator, - vdiff, particle, m_b, rho_a, rho_b, - particle_system::EntropicallyDampedSPHSystem, - grad_kernel) +function continuity_equation!(dv, density_calculator, + vdiff, particle, m_b, rho_a, rho_b, + particle_system::EntropicallyDampedSPHSystem, + grad_kernel) return dv end diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index eb2aa1c1d..9116593a6 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -42,34 +42,42 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more [`BoundaryModelDummyParticles`](@ref) and [`AdamiPressureExtrapolation`](@ref). The keyword argument `acceleration` should be used instead for gravity-like source terms. +- `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) """ -struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, - PF, ST, B, C} <: FluidSystem{NDIMS, IC} +struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, COR, + PF, ST, SRFT, SRFN, B, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: M # Vector{ELTYPE}: [particle] density_calculator :: DC smoothing_kernel :: K smoothing_length :: ELTYPE + ideal_neighbor_count :: Int + color :: Int sound_speed :: ELTYPE viscosity :: V nu_edac :: ELTYPE acceleration :: SVector{NDIMS, ELTYPE} - correction :: Nothing + correction :: COR pressure_acceleration_formulation :: PF transport_velocity :: TV source_terms :: ST + surface_tension :: SRFT + surface_normal_method :: SRFN buffer :: B cache :: C function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, sound_speed; pressure_acceleration=inter_particle_averaged_pressure, - density_calculator=SummationDensity(), + density_calculator=ContinuityDensity(), transport_velocity=nothing, alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), - source_terms=nothing, buffer_size=nothing) + correction=nothing, + source_terms=nothing, surface_tension=nothing, + surface_normal_method=nothing, buffer_size=nothing, + reference_particle_spacing=0.0, color_value=1) buffer = isnothing(buffer_size) ? nothing : SystemBuffer(nparticles(initial_condition), buffer_size) @@ -79,6 +87,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, ELTYPE = eltype(initial_condition) mass = copy(initial_condition.mass) + n_particles = length(initial_condition.mass) if ndims(smoothing_kernel) != NDIMS throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) @@ -89,24 +98,59 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end + if surface_tension !== nothing && surface_normal_method === nothing + surface_normal_method = ColorfieldSurfaceNormal() + end + + if surface_normal_method !== nothing && reference_particle_spacing < eps() + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) + end + + if correction isa ShepardKernelCorrection && + density_calculator isa ContinuityDensity + throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`")) + end + + ideal_neighbor_count_ = 0 + if reference_particle_spacing > 0.0 + ideal_neighbor_count_ = ideal_neighbor_count(Val(NDIMS), + reference_particle_spacing, + compact_support(smoothing_kernel, + smoothing_length)) + end + pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, density_calculator, NDIMS, ELTYPE, - nothing) + correction) nu_edac = (alpha * smoothing_length * sound_speed) / 8 cache = create_cache_density(initial_condition, density_calculator) - cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...) + cache = (; + create_cache_correction(correction, initial_condition.density, NDIMS, + n_particles)..., + create_cache_edac(initial_condition, transport_velocity)..., + create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, + n_particles)..., + create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, + n_particles)..., + cache...) new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), - typeof(transport_velocity), typeof(pressure_acceleration), typeof(source_terms), - typeof(buffer), - typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, - smoothing_length, sound_speed, viscosity, nu_edac, acceleration_, - nothing, pressure_acceleration, transport_velocity, source_terms, - buffer, cache) + typeof(transport_velocity), typeof(correction), typeof(pressure_acceleration), + typeof(source_terms), + typeof(surface_tension), typeof(surface_normal_method), + typeof(buffer), typeof(cache)}(initial_condition, mass, density_calculator, + smoothing_kernel, smoothing_length, + ideal_neighbor_count_, + color_value, sound_speed, viscosity, nu_edac, + acceleration_, + correction, pressure_acceleration, + transport_velocity, source_terms, + surface_tension, surface_normal_method, buffer, + cache) end end @@ -115,9 +159,12 @@ function Base.show(io::IO, system::EntropicallyDampedSPHSystem) print(io, "EntropicallyDampedSPHSystem{", ndims(system), "}(") print(io, system.density_calculator) + print(io, ", ", system.correction) print(io, ", ", system.viscosity) print(io, ", ", system.smoothing_kernel) print(io, ", ", system.acceleration) + print(io, ", ", system.surface_tension) + print(io, ", ", system.surface_normal_method) print(io, ") with ", nparticles(system), " particles") end @@ -136,12 +183,16 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end summary_line(io, "density calculator", system.density_calculator |> typeof |> nameof) + summary_line(io, "correction method", + system.correction |> typeof |> nameof) summary_line(io, "viscosity", system.viscosity |> typeof |> nameof) summary_line(io, "ν₍EDAC₎", "≈ $(round(system.nu_edac; digits=3))") summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) summary_line(io, "tansport velocity formulation", system.transport_velocity |> typeof |> nameof) summary_line(io, "acceleration", system.acceleration) + summary_line(io, "surface tension", system.surface_tension) + summary_line(io, "surface normal method", system.surface_normal_method) summary_footer(io) end end @@ -206,6 +257,32 @@ end function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) compute_density!(system, u, u_ode, semi, system.density_calculator) +end + +function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) + (; density_calculator, correction) = system + + compute_correction_values!(system, correction, u, v_ode, u_ode, semi) + + compute_gradient_correction_matrix!(correction, system, u, v_ode, u_ode, semi) + + # `kernel_correct_density!` only performed for `SummationDensity` + kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction, + density_calculator) + compute_surface_normal!(system, system.surface_normal_method, v, u, v_ode, u_ode, semi, + t) + compute_surface_delta_function!(system, system.surface_tension) + compute_wall_contact_values!(system, system.surface_tension.contact_model, v, u, v_ode, + u_ode, semi, t) +end + +function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) + (; surface_tension) = system + + # Surface normal of neighbor and boundary needs to have been calculated already + compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) + compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) update_average_pressure!(system, system.transport_velocity, v_ode, u_ode, semi) end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index fd0a9c3a5..80515dc7e 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -55,16 +55,22 @@ function compute_density!(system, u, u_ode, semi, ::SummationDensity) summation_density!(system, semi, u, u_ode, density) end -function calculate_dt(v_ode, u_ode, cfl_number, system::FluidSystem) - (; smoothing_length, viscosity, acceleration) = system +function calculate_dt(v_ode, u_ode, cfl_number, system::FluidSystem, semi) + (; smoothing_length, viscosity, acceleration, surface_tension) = system - dt_viscosity = 0.125 * smoothing_length^2 / kinematic_viscosity(system, viscosity) + if !(system.viscosity isa Nothing) + dt_viscosity = 0.125 * smoothing_length^2 / kinematic_viscosity(system, viscosity) + else + dt_viscosity = 0.125 * smoothing_length^2 + end # TODO Adami et al. (2012) just use the gravity here, but Antuono et al. (2012) # are using a per-particle acceleration. Is that supposed to be the previous RHS? + # Morris (2000) also uses the acceleration and cites Monaghan (1992) dt_acceleration = 0.25 * sqrt(smoothing_length / norm(acceleration)) # TODO Everyone seems to be doing this differently. + # Morris (2000) uses the additional condition CFL < 0.25. # Sun et al. (2017) only use h / c (because c depends on v_max as c >= 10 v_max). # Adami et al. (2012) use h / (c + v_max) with a fixed CFL of 0.25. # Antuono et al. (2012) use h / (c + v_max + h * pi_max), where pi is the viscosity coefficient. @@ -73,13 +79,48 @@ function calculate_dt(v_ode, u_ode, cfl_number, system::FluidSystem) # See docstring of the callback for the references. dt_sound_speed = cfl_number * smoothing_length / system_sound_speed(system) - return min(dt_viscosity, dt_acceleration, dt_sound_speed) + # Eq. 28 in Morris (2000) + dt_surface_tension = max(dt_viscosity, dt_acceleration, dt_sound_speed) + if surface_tension isa SurfaceTensionMorris || + surface_tension isa SurfaceTensionMomentumMorris + v = wrap_v(v_ode, system, semi) + dt_surface_tension = sqrt(particle_density(v, system, 1) * smoothing_length^3 / + (2 * pi * surface_tension.surface_tension_coefficient)) + end + + return min(dt_viscosity, dt_acceleration, dt_sound_speed, dt_surface_tension) +end + +@inline function surface_tension_model(system::FluidSystem) + return system.surface_tension +end + +@inline function surface_tension_model(system) + return nothing +end + +@inline function surface_normal_method(system::FluidSystem) + return system.surface_normal_method +end + +@inline function surface_normal_method(system) + return nothing +end + +@inline function correction_matrix(system::FluidSystem, particle) + extract_smatrix(system.cache.correction_matrix, system, particle) +end + +@inline function curvature(particle_system::FluidSystem, particle) + (; cache) = particle_system + return cache.curvature[particle] end include("pressure_acceleration.jl") include("viscosity.jl") include("transport_velocity.jl") include("surface_tension.jl") +include("surface_normal_sph.jl") include("weakly_compressible_sph/weakly_compressible_sph.jl") include("entropically_damped_sph/entropically_damped_sph.jl") diff --git a/src/schemes/fluid/pressure_acceleration.jl b/src/schemes/fluid/pressure_acceleration.jl index 5b0b82e11..2316ab762 100644 --- a/src/schemes/fluid/pressure_acceleration.jl +++ b/src/schemes/fluid/pressure_acceleration.jl @@ -58,6 +58,20 @@ end return -volume_term * pressure_tilde * W_a end +# TODO: some one please check this that this correct... +@inline function inter_particle_averaged_pressure(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a, + W_b) + volume_a = m_a / rho_a + volume_b = m_b / rho_b + volume_term = (volume_a^2 + volume_b^2) / m_a + + # Inter-particle averaged pressure + pressure_tilde = (rho_b * p_a * W_a + rho_a * p_b * W_b) / + (rho_a * W_a + rho_b * W_b) + + return -volume_term * pressure_tilde +end + function choose_pressure_acceleration_formulation(pressure_acceleration, density_calculator, NDIMS, ELTYPE, correction) @@ -90,8 +104,7 @@ end function choose_pressure_acceleration_formulation(pressure_acceleration::Nothing, density_calculator::SummationDensity, - NDIMS, ELTYPE, - correction) + NDIMS, ELTYPE, correction) # Choose the pressure acceleration formulation corresponding to the density calculator. return pressure_acceleration_summation_density @@ -99,8 +112,7 @@ end function choose_pressure_acceleration_formulation(pressure_acceleration::Nothing, density_calculator::ContinuityDensity, - NDIMS, ELTYPE, - correction) + NDIMS, ELTYPE, correction) # Choose the pressure acceleration formulation corresponding to the density calculator. return pressure_acceleration_continuity_density diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl new file mode 100644 index 000000000..32bd72367 --- /dev/null +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -0,0 +1,603 @@ +# TODO: handle corners with StaticNormals +# TODO: add method that determines local static normals +# TODO: add method for TLSPH +@doc raw""" + StaticNormals(normal_vectors::Tuple{Vararg{ELTYPE, NDIMS}}) where {NDIMS, ELTYPE <: Real} + +Represents a unit normal vector and its corresponding unit tangential vector in +2D or 3D space. The input normal vector is normalized to ensure it has a length +of 1, and the tangential vector is computed as a vector perpendicular to the +normal vector. + +# Keywords +- `normal_vectors::Tuple{Vararg{ELTYPE, NDIMS}}`: A tuple representing the normal + vector in NDIMS-dimensional space. It will be normalized internally. + +# Tangential Vector Calculation +- In 2D: The tangential vector is calculated as `[-n[2], n[1]]`, ensuring + it is perpendicular to the normal vector. +- In 3D: The tangential vector is computed using a cross product with a + reference vector that is not parallel to the normal vector. The result is + normalized to ensure unit length. + +# Errors +- Throws `ArgumentError` if the provided normal vector has a length of 0. + +# Example +```julia +sn2d = StaticNormals((3.0, 4.0)) # Normalizes to (0.6, 0.8) and computes tangential (-0.8, 0.6) +sn3d = StaticNormals((0.0, 1.0, 0.0)) # Computes normal and tangential vectors in 3D +``` +""" +struct StaticNormals{NDIMS, ELTYPE <: Real} + normal_vectors::SVector{NDIMS, ELTYPE} + tangential_vectors::SVector{NDIMS, ELTYPE} +end + +function StaticNormals(normal_vectors::Tuple{Vararg{ELTYPE, NDIMS}}) where {NDIMS, + ELTYPE <: Real} + norm_value = sqrt(dot(normal_vectors, normal_vectors)) + if norm_value == 0 + throw(ArgumentError("Normal vector cannot be zero-length.")) + end + normalized_vector = SVector{NDIMS, ELTYPE}(normal_vectors) / norm_value + tangential_vector = calculate_tangential_vector(normalized_vector) + return StaticNormals(normalized_vector, tangential_vector) +end + +# Helper function to calculate tangential vector +function calculate_tangential_vector(n::SVector{2, ELTYPE}) where {ELTYPE} + # Perpendicular vector in 2D + return SVector(-n[2], n[1]) +end + +function calculate_tangential_vector(n::SVector{3, ELTYPE}) where {ELTYPE} + # Cross product with a reference vector to get a perpendicular vector in 3D + ref = abs(n[1]) < abs(n[2]) ? SVector(1.0, 0.0, 0.0) : SVector(0.0, 1.0, 0.0) + t = cross(ref, n) # Perpendicular to 'n' + t_norm = norm(t) + if t_norm == 0 + throw(ArgumentError("Cannot compute tangential vector; invalid input normal vector.")) + end + return t / t_norm # Normalize +end + +function wall_tangential(particle_system::BoundarySystem, particle) + wall_tangential(particle_system, particle, particle_system.surface_normal_method) +end + +function wall_tangential(particle_system::BoundarySystem, particle, surface_normal_method) + return zero(SVector{ndims(particle_system), eltype(particle_system)}) +end + +function wall_tangential(::BoundarySystem, particle, surface_normal_method::StaticNormals) + return surface_normal_method.tangential_vectors +end + +@doc raw""" + ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01, + ideal_density_threshold=0.0) + +Implements a model for computing surface normals based on a color field representation. +This approach is commonly used in fluid simulations to determine interface normals +for multiphase flows or free surfaces. + +# Keywords +- `boundary_contact_threshold=0.0`: The threshold value used to determine contact with boundaries. + Adjust this to refine the detection of surface interfaces near boundaries. +- `interface_threshold=0.01`: The threshold value that defines the presence of an interface. + Lower values can improve sensitivity but may introduce noise. +- `ideal_density_threshold=0.0`: The ideal density threshold used for interface calculations. + This value can be tuned based on the density variations in the simulation. +""" +struct ColorfieldSurfaceNormal{ELTYPE} + boundary_contact_threshold::ELTYPE + interface_threshold::ELTYPE + ideal_density_threshold::ELTYPE +end + +function ColorfieldSurfaceNormal(; boundary_contact_threshold=0.0, interface_threshold=0.01, + ideal_density_threshold=0.0) + return ColorfieldSurfaceNormal(boundary_contact_threshold, interface_threshold, + ideal_density_threshold) +end + +function create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, nparticles) + return (;) +end + +function create_cache_surface_normal(::ColorfieldSurfaceNormal, ELTYPE, NDIMS, nparticles) + surface_normal = Array{ELTYPE, 2}(undef, NDIMS, nparticles) + neighbor_count = Array{ELTYPE, 1}(undef, nparticles) + colorfield = Array{ELTYPE, 1}(undef, nparticles) + return (; surface_normal, neighbor_count, colorfield) +end + +@inline function surface_normal(particle_system::FluidSystem, particle) + (; cache) = particle_system + return extract_svector(cache.surface_normal, particle_system, particle) +end + +@inline function surface_normal(particle_system::BoundarySystem, particle) + (; surface_normal_method) = particle_system + return surface_normal(particle_system, particle, surface_normal_method) +end + +@inline function surface_normal(particle_system::BoundarySystem, particle, + surface_normal_method) + return zero(SVector{ndims(particle_system), eltype(particle_system)}) +end + +@inline function surface_normal(::BoundarySystem, particle, + surface_normal_method::StaticNormals) + return surface_normal_method.normal_vectors +end + +function calc_normal!(system, neighbor_system, u_system, v, v_neighbor_system, + u_neighbor_system, semi, surfn, nsurfn) + # Normal not needed + return system +end + +# Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" +# and Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function calc_normal!(system::FluidSystem, neighbor_system::FluidSystem, v, u, + v_neighbor_system, u_neighbor_system, semi, surfn, + ::ColorfieldSurfaceNormal) + (; cache) = system + + # TODO: actually calculate if there is more than one colored fluid + cache.colorfield .= system.color + system_coords = current_coordinates(u, system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + nhs = get_neighborhood_search(system, neighbor_system, semi) + + foreach_point_neighbor(system, neighbor_system, + system_coords, neighbor_system_coords, + nhs) do particle, neighbor, pos_diff, distance + m_b = hydrodynamic_mass(neighbor_system, neighbor) + density_neighbor = particle_density(v_neighbor_system, + neighbor_system, neighbor) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance) + for i in 1:ndims(system) + cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] + end + + cache.neighbor_count[particle] += 1 + end + + return system +end + +# Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" +# and Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function calc_normal!(system::FluidSystem, neighbor_system::BoundarySystem, v, u, + v_neighbor_system, u_neighbor_system, semi, surfn, nsurfn) + (; cache) = system + (; colorfield, colorfield_bnd) = neighbor_system.boundary_model.cache + (; boundary_contact_threshold) = surfn + + system_coords = current_coordinates(u, system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + nhs = get_neighborhood_search(neighbor_system, system, semi) + + # First we need to calculate the smoothed colorfield values of the boundary + # TODO: move colorfield to extra step + # TODO: this is only correct for a single fluid + + # # Reset to the constant boundary interpolated color values + # colorfield .= colorfield_bnd + + # Accumulate fluid neighbors + foreach_point_neighbor(neighbor_system, system, neighbor_system_coords, system_coords, + nhs, points=eachparticle(neighbor_system)) do particle, neighbor, pos_diff, distance + colorfield[particle] += hydrodynamic_mass(system, neighbor) / + particle_density(v, system, neighbor) * system.color * + smoothing_kernel(system, distance) + end + + if boundary_contact_threshold < eps() + return system + end + + foreach_point_neighbor(system, neighbor_system, + system_coords, neighbor_system_coords, + nhs) do particle, neighbor, pos_diff, distance + # we assume that we are in contact with the boundary if the color of the boundary particle is larger than the threshold + if colorfield[neighbor] > boundary_contact_threshold + m_b = hydrodynamic_mass(system, particle) + density_neighbor = particle_density(v, system, particle) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance) + for i in 1:ndims(system) + cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] + end + cache.neighbor_count[particle] += 1 + end + end + + return system +end + +function remove_invalid_normals!(system::FluidSystem, surface_tension, surfn) + (; cache) = system + + # We remove invalid normals (too few neighbors) to reduce the impact of underdefined normals + for particle in each_moving_particle(system) + # A corner has that many neighbors assuming a regular 2 * r distribution and a compact_support of 4r + if cache.neighbor_count[particle] < 2^ndims(system) + 1 + cache.surface_normal[1:ndims(system), particle] .= 0 + end + end + + return system +end + +# see Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +# Note: this also normalizes the normals +function remove_invalid_normals!(system::FluidSystem, + surface_tension::Union{SurfaceTensionMorris, + SurfaceTensionMomentumMorris}, + surfn::ColorfieldSurfaceNormal) + (; cache, smoothing_length, smoothing_kernel, ideal_neighbor_count) = system + (; ideal_density_threshold, interface_threshold) = surfn + + # We remove invalid normals i.e. they have a small norm (eq. 20) + normal_condition2 = (interface_threshold / + compact_support(smoothing_kernel, smoothing_length))^2 + + for particle in each_moving_particle(system) + + # heuristic condition if there is no gas phase to find the free surface + # We remove normals for particles which have alot of support e.g. they are in the inside + if ideal_density_threshold > 0 && + ideal_density_threshold * ideal_neighbor_count < cache.neighbor_count[particle] + if surface_tension.contact_model isa HuberContactModel + cache.normal_v[1:ndims(system), particle] .= 0 + end + cache.surface_normal[1:ndims(system), particle] .= 0 + continue + end + + particle_surface_normal = cache.surface_normal[1:ndims(system), particle] + norm2 = dot(particle_surface_normal, particle_surface_normal) + + # see eq. 21 + if norm2 > normal_condition2 + if surface_tension.contact_model isa HuberContactModel + cache.normal_v[1:ndims(system), particle] = particle_surface_normal + end + + cache.surface_normal[1:ndims(system), particle] = particle_surface_normal / + sqrt(norm2) + else + if surface_tension.contact_model isa HuberContactModel + cache.normal_v[1:ndims(system), particle] .= 0 + end + cache.surface_normal[1:ndims(system), particle] .= 0 + end + end + + return system +end + +function compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) + return system +end + +function compute_surface_normal!(system::FluidSystem, + surface_normal_method_::ColorfieldSurfaceNormal, + v, u, v_ode, u_ode, semi, t) + (; cache, surface_tension) = system + + # Reset surface normal + set_zero!(cache.surface_normal) + set_zero!(cache.neighbor_count) + + # TODO: if color values are set only different systems need to be called + @trixi_timeit timer() "compute surface normal" foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + + calc_normal!(system, neighbor_system, v, u, v_neighbor_system, + u_neighbor_system, semi, surface_normal_method_, + surface_normal_method(neighbor_system)) + end + remove_invalid_normals!(system, surface_tension, surface_normal_method_) + + return system +end + +function calc_curvature!(system, neighbor_system, u_system, v, + v_neighbor_system, u_neighbor_system, semi, surfn, nsurfn) +end + +# Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function calc_curvature!(system::FluidSystem, neighbor_system::FluidSystem, u_system, v, + v_neighbor_system, u_neighbor_system, semi, + surfn::ColorfieldSurfaceNormal, nsurfn::ColorfieldSurfaceNormal) + (; cache) = system + (; curvature) = cache + + system_coords = current_coordinates(u_system, system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + nhs = get_neighborhood_search(system, neighbor_system, semi) + correction_factor = fill(eps(eltype(system)), n_moving_particles(system)) + + no_valid_neighbors = 0 + + foreach_point_neighbor(system, neighbor_system, + system_coords, neighbor_system_coords, + nhs) do particle, neighbor, pos_diff, distance + m_b = hydrodynamic_mass(neighbor_system, neighbor) + rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + n_a = surface_normal(system, particle) + n_b = surface_normal(neighbor_system, neighbor) + v_b = m_b / rho_b + + # eq. 22 we can test against eps() here since the surface normals that are invalid have been reset + if dot(n_a, n_a) > eps() && dot(n_b, n_b) > eps() + w = smoothing_kernel(system, distance) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + + for i in 1:ndims(system) + curvature[particle] += v_b * (n_b[i] - n_a[i]) * grad_kernel[i] + end + # eq. 24 + correction_factor[particle] += v_b * w + # prevent NaNs from systems that are entirely skipped + no_valid_neighbors += 1 + end + end + + # eq. 23 + if no_valid_neighbors > 0 + for i in 1:n_moving_particles(system) + curvature[i] /= correction_factor[i] + end + end + + return system +end + +function compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) + return system +end + +function compute_curvature!(system::FluidSystem, surface_tension::SurfaceTensionMorris, v, + u, v_ode, u_ode, semi, t) + (; cache, surface_tension) = system + + # Reset surface curvature + set_zero!(cache.curvature) + + @trixi_timeit timer() "compute surface curvature" foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + + calc_curvature!(system, neighbor_system, u, v, v_neighbor_system, + u_neighbor_system, semi, surface_normal_method(system), + surface_normal_method(neighbor_system)) + end + return system +end + +function calc_stress_tensors!(system::FluidSystem, neighbor_system::FluidSystem, u_system, + v, + v_neighbor_system, u_neighbor_system, semi, + surfn::ColorfieldSurfaceNormal, + nsurfn::ColorfieldSurfaceNormal) + (; cache) = system + (; smoothing_kernel, smoothing_length) = surfn + (; stress_tensor, delta_s) = cache + + neighbor_cache = neighbor_system.cache + neighbor_delta_s = neighbor_cache.delta_s + + NDIMS = ndims(system) + max_delta_s = max(maximum(delta_s), maximum(neighbor_delta_s)) + + for particle in each_moving_particle(system) + normal = surface_normal(system, particle) + delta_s_particle = delta_s[particle] + if delta_s_particle > eps() + for i in 1:NDIMS + for j in 1:NDIMS + delta_ij = (i == j) ? 1.0 : 0.0 + stress_tensor[i, j, particle] = delta_s_particle * + (delta_ij - normal[i] * normal[j]) - + delta_ij * max_delta_s + end + end + end + end + + return system +end + +function compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) + return system +end + +# Section 6 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function compute_stress_tensors!(system::FluidSystem, ::SurfaceTensionMomentumMorris, + v, u, v_ode, u_ode, semi, t) + (; cache) = system + (; delta_s, stress_tensor) = cache + + # Reset surface stress_tensor + set_zero!(stress_tensor) + + max_delta_s = maximum(delta_s) + NDIMS = ndims(system) + + @trixi_timeit timer() "compute surface stress tensor" for particle in each_moving_particle(system) + normal = surface_normal(system, particle) + delta_s_particle = delta_s[particle] + if delta_s_particle > eps() + for i in 1:NDIMS + for j in 1:NDIMS + delta_ij = (i == j) ? 1.0 : 0.0 + stress_tensor[i, j, particle] = delta_s_particle * + (delta_ij - normal[i] * normal[j]) - + delta_ij * max_delta_s + end + end + end + end + + return system +end + +function compute_surface_delta_function!(system, surface_tension) + return system +end + +# eq. 6 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function compute_surface_delta_function!(system, ::SurfaceTensionMomentumMorris) + (; cache) = system + (; delta_s) = cache + + set_zero!(delta_s) + + for particle in each_moving_particle(system) + delta_s[particle] = norm(surface_normal(system, particle)) + end + return system +end + +function calc_wall_distance_vector!(system, neighbor_system, + v, u, v_neighbor_system, u_neighbor_system, + semi, contact_model) +end + +# Computes the wall distance vector \( \mathbf{d}_a \) and its normalized version \( \hat{\mathbf{d}}_a \) for each particle in the fluid system. +# The distance vector \( \mathbf{d}_a \) is used to calculate dynamic contact angles and the delta function in wall-contact models. +function calc_wall_distance_vector!(system::FluidSystem, neighbor_system::BoundarySystem, + v, u, v_neighbor_system, u_neighbor_system, + semi, contact_model::HuberContactModel) + cache = system.cache + (; d_hat, d_vec) = cache + NDIMS = ndims(system) + + system_coords = current_coordinates(u, system) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + nhs = get_neighborhood_search(system, neighbor_system, semi) + + foreach_point_neighbor(system, neighbor_system, + system_coords, neighbor_coords, + nhs) do particle, neighbor, pos_diff, distance + V_b = hydrodynamic_mass(neighbor_system, neighbor) / + particle_density(v_neighbor_system, neighbor_system, neighbor) + W_ab = smoothing_kernel(system, distance) + + # equation 51 + for i in 1:NDIMS + d_vec[i, particle] += V_b * pos_diff[i] * W_ab + end + end + + # d_hat is the unit version of d_vec + for particle in each_moving_particle(system) + norm_d = sqrt(sum(d_vec[i, particle]^2 for i in 1:NDIMS)) + if norm_d > eps() + for i in 1:NDIMS + d_hat[i, particle] = d_vec[i, particle] / norm_d + end + else + for i in 1:NDIMS + d_hat[i, particle] = 0.0 + end + end + end + + return system +end + +function calc_wall_contact_values!(system, neighbor_system, + v, u, v_neighbor_system, u_neighbor_system, + semi, contact_model) +end + +function calc_wall_contact_values!(system::FluidSystem, neighbor_system::BoundarySystem, + v, u, v_neighbor_system, u_neighbor_system, + semi, contact_model::HuberContactModel) + cache = system.cache + (; d_hat, delta_wns, normal_v, nu_hat, d_vec) = cache + NDIMS = ndims(system) + + system_coords = current_coordinates(u, system) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + nhs = get_neighborhood_search(system, neighbor_system, semi) + + foreach_point_neighbor(system, neighbor_system, + system_coords, neighbor_coords, + nhs) do particle, neighbor, pos_diff, distance + V_b = hydrodynamic_mass(neighbor_system, neighbor) / + particle_density(v_neighbor_system, neighbor_system, neighbor) + + grad_W_ab = smoothing_kernel_grad(system, pos_diff, distance) + + # equation 53 computing the tangential direction vector + for i in 1:NDIMS + nu_hat[i, particle] = dot(d_vec[:, particle], d_vec[:, particle]) * + normal_v[i, particle] - + (d_vec[i, particle] * normal_v[i, particle]) * + d_vec[i, particle] + end + + nu_norm = norm(nu_hat[:, particle]) + if nu_norm > eps() + for i in 1:NDIMS + nu_hat[i, particle] /= nu_norm + end + else + for i in 1:NDIMS + nu_hat[i, particle] = 0.0 + end + end + + # equation 54 delta function + dot_nu_n = sum(nu_hat[i, particle] * normal_v[i, particle] for i in 1:NDIMS) + dot_d_gradW = sum(d_hat[i, particle] * grad_W_ab[i] for i in 1:NDIMS) + + delta_wns[particle] += 2 * V_b * dot_d_gradW * dot_nu_n + end + + return system +end + +function compute_wall_contact_values!(system, contact_model, v, u, v_ode, u_ode, semi, t) +end + +function compute_wall_contact_values!(system::FluidSystem, contact_model::HuberContactModel, + v, u, v_ode, u_ode, semi, t) + cache = system.cache + (; d_hat, delta_wns) = cache + + # Reset d_hat and delta_wns_partial + set_zero!(d_hat) + set_zero!(delta_wns) + + # Loop over boundary neighbor systems + @trixi_timeit timer() "compute wall contact values" foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + + calc_wall_distance_vector!(system, neighbor_system, v, u, + v_neighbor_system, u_neighbor_system, semi, + contact_model) + end + + # Loop over boundary neighbor systems + @trixi_timeit timer() "compute wall contact values" foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + + # Call calc_wall_distance_vector! + calc_wall_contact_values!(system, neighbor_system, v, u, + v_neighbor_system, u_neighbor_system, semi, + contact_model) + end + + return system +end diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 66140775e..6a43414c5 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -1,4 +1,5 @@ -abstract type AkinciTypeSurfaceTension end +abstract type SurfaceTension end +abstract type AkinciTypeSurfaceTension <: SurfaceTension end @doc raw""" CohesionForceAkinci(surface_tension_coefficient=1.0) @@ -38,6 +39,94 @@ struct SurfaceTensionAkinci{ELTYPE} <: AkinciTypeSurfaceTension end end +@doc raw""" + SurfaceTensionMorris(surface_tension_coefficient=1.0) + +This model implements the surface tension approach described by [Morris2000](@cite). +It calculates surface tension forces based on the curvature of the fluid interface +using particle normals and their divergence, making it suitable for simulating +phenomena like droplet formation and capillary wave dynamics. + +# Details +The method estimates curvature by combining particle color gradients and smoothing +functions to derive surface normals. The curvature is then used to compute forces +acting perpendicular to the interface. While this method provides accurate +surface tension forces, it does not conserve momentum explicitly. + +# Keywords +- `surface_tension_coefficient=1.0`: Adjusts the magnitude of the surface tension + forces, enabling tuning of fluid surface behaviors in simulations. +""" +struct SurfaceTensionMorris{ELTYPE, CM} <: SurfaceTension + surface_tension_coefficient :: ELTYPE + contact_model :: CM + + function SurfaceTensionMorris(; surface_tension_coefficient=1.0, contact_model=nothing) + new{typeof(surface_tension_coefficient), typeof(contact_model)}(surface_tension_coefficient, + contact_model) + end +end + +function create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, nparticles) + return (;) +end + +function create_cache_surface_tension(st::SurfaceTensionMorris, ELTYPE, NDIMS, nparticles) + curvature = Array{ELTYPE, 1}(undef, nparticles) + if st.contact_model isa Nothing + return (; curvature) + else + return (; + create_cache_contact_model(st.contact_model, ELTYPE, NDIMS, nparticles)..., + curvature) + end +end + +@doc raw""" + SurfaceTensionMomentumMorris(surface_tension_coefficient=1.0) + +This model implements the momentum-conserving surface tension approach outlined by +[Morris2000](@cite). It calculates surface tension forces using the gradient of a stress +tensor, ensuring exact conservation of linear momentum. This method is particularly +useful for simulations where momentum conservation is critical, though it may require +numerical adjustments at higher resolutions. + +# Details +The stress tensor approach replaces explicit curvature calculations, avoiding the +singularities associated with resolution increases. However, the method is computationally +intensive and may require stabilization techniques to handle tensile instability at high +particle densities. + +# Keywords +- `surface_tension_coefficient=1.0`: A parameter to adjust the strength of surface tension + forces, allowing fine-tuning to replicate physical behavior. +""" +struct SurfaceTensionMomentumMorris{ELTYPE, CM} <: SurfaceTension + surface_tension_coefficient::ELTYPE + contact_model::CM + + function SurfaceTensionMomentumMorris(; surface_tension_coefficient=1.0, + contact_model=nothing) + new{typeof(surface_tension_coefficient), typeof(contact_model)}(surface_tension_coefficient, + contact_model) + end +end + +function create_cache_surface_tension(::SurfaceTensionMomentumMorris, ELTYPE, NDIMS, + nparticles) + # Allocate stress tensor for each particle: NDIMS x NDIMS x nparticles + delta_s = Array{ELTYPE, 1}(undef, nparticles) + stress_tensor = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, nparticles) + + if st.contact_model isa Nothing + return (; stress_tensor, delta_s) + else + return (; + create_cache_contact_model(st.contact_model, ELTYPE, NDIMS, nparticles)..., + delta_s, stress_tensor) + end +end + # Note that `floating_point_number^integer_literal` is lowered to `Base.literal_pow`. # Currently, specializations reducing this to simple multiplications exist only up # to a power of three, see @@ -89,47 +178,18 @@ end return adhesion_force end -# Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" -# Note: most of the time this only leads to an approximation of the surface normal - -function calc_normal_akinci!(system, neighbor_system::FluidSystem, - surface_tension::SurfaceTensionAkinci, u_system, - v_neighbor_system, u_neighbor_system, - neighborhood_search) - (; smoothing_length, cache) = system - - system_coords = current_coordinates(u_system, system) - neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) - - foreach_point_neighbor(system, neighbor_system, - system_coords, neighbor_system_coords, - neighborhood_search) do particle, neighbor, pos_diff, distance - m_b = hydrodynamic_mass(neighbor_system, neighbor) - density_neighbor = particle_density(v_neighbor_system, - neighbor_system, neighbor) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, - particle) - for i in 1:ndims(system) - cache.surface_normal[i, particle] += m_b / density_neighbor * - grad_kernel[i] * smoothing_length - end - end - - return system -end - -function calc_normal_akinci!(system, neighbor_system, surface_tension, u_system, - v_neighbor_system, u_neighbor_system, - neighborhood_search) - # Normal not needed - return system +# Skip +@inline function surface_tension_force(surface_tension_a, surface_tension_b, + particle_system, neighbor_system, particle, neighbor, + pos_diff, distance, rho_a, rho_b, grad_kernel) + return zero(pos_diff) end @inline function surface_tension_force(surface_tension_a::CohesionForceAkinci, surface_tension_b::CohesionForceAkinci, particle_system::FluidSystem, neighbor_system::FluidSystem, particle, neighbor, - pos_diff, distance) + pos_diff, distance, rho_a, rho_b, grad_kernel) (; smoothing_length) = particle_system # No cohesion with oneself distance < sqrt(eps()) && return zero(pos_diff) @@ -144,7 +204,7 @@ end surface_tension_b::SurfaceTensionAkinci, particle_system::FluidSystem, neighbor_system::FluidSystem, particle, neighbor, - pos_diff, distance) + pos_diff, distance, rho_a, rho_b, grad_kernel) (; smoothing_length, smoothing_kernel) = particle_system (; surface_tension_coefficient) = surface_tension_a @@ -152,20 +212,47 @@ end distance < sqrt(eps()) && return zero(pos_diff) m_b = hydrodynamic_mass(neighbor_system, neighbor) - n_a = surface_normal(surface_tension_a, particle_system, particle) - n_b = surface_normal(surface_tension_b, neighbor_system, neighbor) + n_a = surface_normal(particle_system, particle) + n_b = surface_normal(neighbor_system, neighbor) support_radius = compact_support(smoothing_kernel, smoothing_length) return cohesion_force_akinci(surface_tension_a, support_radius, m_b, pos_diff, distance) .- - (surface_tension_coefficient * (n_a - n_b)) + (surface_tension_coefficient * (n_a - n_b) * smoothing_length) end -# Skip -@inline function surface_tension_force(surface_tension_a, surface_tension_b, - particle_system, neighbor_system, particle, neighbor, - pos_diff, distance) - return zero(pos_diff) +@inline function surface_tension_force(surface_tension_a::SurfaceTensionMorris, + surface_tension_b::SurfaceTensionMorris, + particle_system::FluidSystem, + neighbor_system::FluidSystem, particle, neighbor, + pos_diff, distance, rho_a, rho_b, grad_kernel) + (; surface_tension_coefficient) = surface_tension_a + + # This force only applies to itself + !(particle == neighbor) && return zero(pos_diff) + + n_a = surface_normal(particle_system, particle) + curvature_a = curvature(particle_system, particle) + + return -surface_tension_coefficient / rho_a * curvature_a * n_a +end + +@inline function surface_tension_force(surface_tension_a::SurfaceTensionMomentumMorris, + surface_tension_b::SurfaceTensionMomentumMorris, + particle_system::FluidSystem, + neighbor_system::FluidSystem, particle, neighbor, + pos_diff, distance, rho_a, rho_b, grad_kernel) + (; surface_tension_coefficient) = surface_tension_a + + # No surface tension with oneself + distance < sqrt(eps()) && return zero(pos_diff) + + S_a = particle_system.cache.stress_tensor[:, :, particle] + S_b = neighbor_system.cache.stress_tensor[:, :, neighbor] + + m_b = hydrodynamic_mass(neighbor_system, neighbor) + + return surface_tension_coefficient * m_b * (S_a + S_b) / (rho_a * rho_b) * grad_kernel end @inline function adhesion_force(surface_tension::AkinciTypeSurfaceTension, @@ -192,3 +279,63 @@ end neighbor, pos_diff, distance) return zero(pos_diff) end + +@doc raw""" + HuberContactModel{ELTYPE} + +Represents the contact force model inspired by [Huber2016](@cite). This model is +designed to calculate physically accurate surface tension and contact line +forces in fluid simulations, specifically those involving Smoothed Particle +Hydrodynamics (SPH). + +The model enables the simulation of wetting phenomena, including dynamic contact +angles and their influence on the system's dynamics, using physically-based +parameters without introducing fitting coefficients. +""" +struct HuberContactModel end + +function create_cache_contact_model(contact_model, ELTYPE, NDIMS, nparticles) + return (;) +end + +function create_cache_contact_model(::HuberContactModel, ELTYPE, NDIMS, nparticles) + d_vec = Array{ELTYPE}(undef, NDIMS, nparticles) + d_hat = Array{ELTYPE}(undef, NDIMS, nparticles) + nu_hat = Array{ELTYPE}(undef, NDIMS, nparticles) + nu_hat .= 0.0 + delta_wns = Array{ELTYPE}(undef, nparticles) + # non-normal vector + normal_v = Array{ELTYPE}(undef, NDIMS, nparticles) + return (; d_hat, delta_wns, normal_v, nu_hat, d_vec) +end + +@inline function contact_force(contact_model, particle_system, particle) + return zero(ndims(particle_system)) +end + +@inline function contact_force(::HuberContactModel, particle_system::FluidSystem, particle) + cache = particle_system.cache + sigma = particle_system.surface_tension.surface_tension_coefficient + + # Retrieve the normalized dynamic direction vector d̂ₐ for particle 'a' (Equation 51, normalized) + d_hat_particle = cache.d_hat[:, particle] + + if norm(d_hat_particle) < eps() + return zero(ndims(particle_system)) + end + + nu_hat_particle = cache.nu_hat[:, particle] + delta_wns_particle = cache.delta_wns[particle] + n_hat_particle = surface_normal(particle_system, particle) + + dot_d_n = dot(d_hat_particle, n_hat_particle) + + # Equation 52 + static_contact_angle = 1.5708 #neighbor_system.static_contact_angle #TODO: FIX + f_wns_a = sigma * (cos(static_contact_angle) + dot_d_n) .* nu_hat_particle + + # Compute the force contribution (Equation 57) (needs to be in acceleration form) + F_a = f_wns_a * delta_wns_particle / hydrodynamic_mass(particle_system, particle) + + return F_a +end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 0fe58ab21..d289fb261 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -6,7 +6,7 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, particle_system::WeaklyCompressibleSPHSystem, neighbor_system) - (; density_calculator, state_equation, correction, surface_tension) = particle_system + (; density_calculator, state_equation, correction) = particle_system (; sound_speed) = state_equation surface_tension_a = surface_tension_model(particle_system) @@ -65,17 +65,24 @@ function interact!(dv, v_particle_system, u_particle_system, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + # TODO: only applies to fluids with the same color dv_surface_tension = surface_tension_correction * surface_tension_force(surface_tension_a, surface_tension_b, particle_system, neighbor_system, - particle, neighbor, pos_diff, distance) + particle, neighbor, pos_diff, distance, + rho_a, rho_b, grad_kernel) - dv_adhesion = adhesion_force(surface_tension, particle_system, neighbor_system, + dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system, particle, neighbor, pos_diff, distance) + # dv_contact_force = contact_force(surface_tension_a.contact_model, + # particle_system, neighbor_system, particle, + # neighbor, pos_diff, distance, rho_a, rho_b) + for i in 1:ndims(particle_system) @inbounds dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_surface_tension[i] + dv_adhesion[i] + # + dv_contact_force[i] # Debug example # debug_array[i, particle] += dv_pressure[i] end @@ -94,6 +101,14 @@ function interact!(dv, v_particle_system, u_particle_system, # TODO: This call should use public API. This requires some additional changes to simplify the calls. # trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug", iter=iter += 1) + for particle in each_moving_particle(particle_system) + F = contact_force(particle_system.surface_tension.contact_model, + particle_system, particle) + for i in 1:ndims(particle_system) + @inbounds dv[i, particle] += F[i] + end + end + return dv end diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 6d66ba06a..2fecc66a5 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -45,7 +45,8 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. """ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, - V, DD, COR, PF, ST, B, SRFT, C} <: FluidSystem{NDIMS, IC} + V, DD, COR, PF, ST, B, SRFT, SRFN, C} <: + FluidSystem{NDIMS, IC} initial_condition :: IC mass :: MA # Array{ELTYPE, 1} pressure :: P # Array{ELTYPE, 1} @@ -53,6 +54,8 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, state_equation :: SE smoothing_kernel :: K smoothing_length :: ELTYPE + ideal_neighbor_count :: Int + color :: Int acceleration :: SVector{NDIMS, ELTYPE} viscosity :: V density_diffusion :: DD @@ -61,22 +64,24 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, transport_velocity :: Nothing # TODO source_terms :: ST surface_tension :: SRFT + surface_normal_method :: SRFN buffer :: B cache :: C end # The default constructor needs to be accessible for Adapt.jl to work with this struct. # See the comments in general/gpu.jl for more details. -function WeaklyCompressibleSPHSystem(initial_condition, - density_calculator, state_equation, +function WeaklyCompressibleSPHSystem(initial_condition, state_equation, smoothing_kernel, smoothing_length; + density_calculator=ContinuityDensity(), pressure_acceleration=nothing, buffer_size=nothing, viscosity=nothing, density_diffusion=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), correction=nothing, source_terms=nothing, - surface_tension=nothing) + surface_tension=nothing, surface_normal_method=nothing, + reference_particle_spacing=0.0, color_value=1) buffer = isnothing(buffer_size) ? nothing : SystemBuffer(nparticles(initial_condition), buffer_size) @@ -104,6 +109,21 @@ function WeaklyCompressibleSPHSystem(initial_condition, throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`")) end + if surface_tension !== nothing && surface_normal_method === nothing + surface_normal_method = ColorfieldSurfaceNormal() + end + + if surface_normal_method !== nothing && reference_particle_spacing < eps() + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) + end + + ideal_neighbor_count_ = 0 + if reference_particle_spacing > 0.0 + ideal_neighbor_count_ = ideal_neighbor_count(Val(NDIMS), reference_particle_spacing, + compact_support(smoothing_kernel, + smoothing_length)) + end + pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, density_calculator, NDIMS, ELTYPE, @@ -111,48 +131,23 @@ function WeaklyCompressibleSPHSystem(initial_condition, cache = create_cache_density(initial_condition, density_calculator) cache = (; - create_cache_wcsph(correction, initial_condition.density, NDIMS, - n_particles)..., cache...) - cache = (; - create_cache_wcsph(surface_tension, ELTYPE, NDIMS, n_particles)..., + create_cache_correction(correction, initial_condition.density, NDIMS, + n_particles)..., + create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, + n_particles)..., + create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, + n_particles)..., cache...) return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, smoothing_kernel, smoothing_length, + ideal_neighbor_count_, color_value, acceleration_, viscosity, density_diffusion, correction, pressure_acceleration, nothing, - source_terms, surface_tension, buffer, cache) -end - -create_cache_wcsph(correction, density, NDIMS, nparticles) = (;) - -function create_cache_wcsph(::ShepardKernelCorrection, density, NDIMS, n_particles) - return (; kernel_correction_coefficient=similar(density)) -end - -function create_cache_wcsph(::KernelCorrection, density, NDIMS, n_particles) - dw_gamma = Array{Float64}(undef, NDIMS, n_particles) - return (; kernel_correction_coefficient=similar(density), dw_gamma) -end - -function create_cache_wcsph(::Union{GradientCorrection, BlendedGradientCorrection}, density, - NDIMS, n_particles) - correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) - return (; correction_matrix) -end - -function create_cache_wcsph(::MixedKernelGradientCorrection, density, NDIMS, n_particles) - dw_gamma = Array{Float64}(undef, NDIMS, n_particles) - correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) - - return (; kernel_correction_coefficient=similar(density), dw_gamma, correction_matrix) -end - -function create_cache_wcsph(::SurfaceTensionAkinci, ELTYPE, NDIMS, nparticles) - surface_normal = Array{ELTYPE, 2}(undef, NDIMS, nparticles) - return (; surface_normal) + source_terms, surface_tension, surface_normal_method, + buffer, cache) end function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) @@ -166,6 +161,10 @@ function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) print(io, ", ", system.viscosity) print(io, ", ", system.density_diffusion) print(io, ", ", system.surface_tension) + print(io, ", ", system.surface_normal_method) + if system.surface_normal_method isa ColorfieldSurfaceNormal + print(io, ", ", system.color) + end print(io, ", ", system.acceleration) print(io, ", ", system.source_terms) print(io, ") with ", nparticles(system), " particles") @@ -193,6 +192,10 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst summary_line(io, "viscosity", system.viscosity) summary_line(io, "density diffusion", system.density_diffusion) summary_line(io, "surface tension", system.surface_tension) + summary_line(io, "surface normal method", system.surface_normal_method) + if system.surface_normal_method isa ColorfieldSurfaceNormal + summary_line(io, "color", system.color) + end summary_line(io, "acceleration", system.acceleration) summary_line(io, "source terms", system.source_terms |> typeof |> nameof) summary_footer(io) @@ -220,7 +223,7 @@ end function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) - (; density_calculator, density_diffusion, correction) = system + (; density_calculator, density_diffusion) = system compute_density!(system, u, u_ode, semi, density_calculator) @@ -232,7 +235,7 @@ function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, end function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) - (; density_calculator, correction, surface_tension) = system + (; density_calculator, correction, surface_normal_method, surface_tension) = system compute_correction_values!(system, correction, u, v_ode, u_ode, semi) @@ -242,74 +245,20 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction, density_calculator) compute_pressure!(system, v) - compute_surface_normal!(system, surface_tension, v, u, v_ode, u_ode, semi, t) - - return system -end - -function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, - semi, correction, density_calculator) + compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) + compute_surface_delta_function!(system, surface_tension) + compute_wall_contact_values!(system, system.surface_tension.contact_model, v, + u, v_ode, u_ode, semi, t) return system end -function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, - semi, corr::ShepardKernelCorrection, ::SummationDensity) - system.cache.density ./= system.cache.kernel_correction_coefficient -end - -function compute_gradient_correction_matrix!(correction, - system::WeaklyCompressibleSPHSystem, u, - v_ode, u_ode, semi) - return system -end +function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) + (; surface_tension) = system -function compute_gradient_correction_matrix!(corr::Union{GradientCorrection, - BlendedGradientCorrection, - MixedKernelGradientCorrection}, - system::WeaklyCompressibleSPHSystem, u, - v_ode, u_ode, semi) - (; cache, correction, smoothing_kernel, smoothing_length) = system - (; correction_matrix) = cache - - system_coords = current_coordinates(u, system) - - compute_gradient_correction_matrix!(correction_matrix, system, system_coords, - v_ode, u_ode, semi, correction, smoothing_length, - smoothing_kernel) -end - -function reinit_density!(vu_ode, semi) - v_ode, u_ode = vu_ode.x - - foreach_system(semi) do system - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - - reinit_density!(system, v, u, v_ode, u_ode, semi) - end - - return vu_ode -end - -function reinit_density!(system::WeaklyCompressibleSPHSystem, v, u, - v_ode, u_ode, semi) - # Compute density with `SummationDensity` and store the result in `v`, - # overwriting the previous integrated density. - summation_density!(system, semi, u, u_ode, v[end, :]) - - # Apply `ShepardKernelCorrection` - kernel_correction_coefficient = zeros(size(v[end, :])) - compute_shepard_coeff!(system, current_coordinates(u, system), v_ode, u_ode, semi, - kernel_correction_coefficient) - v[end, :] ./= kernel_correction_coefficient - - compute_pressure!(system, v) - - return system -end - -function reinit_density!(system, v, u, v_ode, u_ode, semi) - return system + # Surface normal of neighbor and boundary needs to have been calculated already + compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) + compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) end function compute_pressure!(system, v) @@ -353,43 +302,3 @@ function restart_with!(system, ::ContinuityDensity, v, u) return system end - -@inline function correction_matrix(system::WeaklyCompressibleSPHSystem, particle) - extract_smatrix(system.cache.correction_matrix, system, particle) -end - -function compute_surface_normal!(system, surface_tension, v, u, v_ode, u_ode, semi, t) - return system -end - -function compute_surface_normal!(system, surface_tension::SurfaceTensionAkinci, v, u, v_ode, - u_ode, semi, t) - (; cache) = system - - # Reset surface normal - set_zero!(cache.surface_normal) - - @trixi_timeit timer() "compute surface normal" foreach_system(semi) do neighbor_system - u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) - v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) - nhs = get_neighborhood_search(system, semi) - - calc_normal_akinci!(system, neighbor_system, surface_tension, u, v_neighbor_system, - u_neighbor_system, nhs) - end - return system -end - -@inline function surface_normal(::SurfaceTensionAkinci, particle_system::FluidSystem, - particle) - (; cache) = particle_system - return extract_svector(cache.surface_normal, particle_system, particle) -end - -@inline function surface_tension_model(system::WeaklyCompressibleSPHSystem) - return system.surface_tension -end - -@inline function surface_tension_model(system) - return nothing -end diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index dfd71a2d6..6bda711cf 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -201,6 +201,7 @@ end return system.boundary_model.hydrodynamic_mass[particle] end +# TODO: move TLSPH correction matrix into cache so that there is no conflict with the other schemes @inline function correction_matrix(system, particle) extract_smatrix(system.correction_matrix, system, particle) end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 9f657f733..f3e3e96c5 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -247,6 +247,67 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) vtk["pressure"] = [particle_pressure(v, system, particle) for particle in active_particles(system)] + if system.surface_normal_method !== nothing + vtk["surf_normal"] = [surface_normal(system, particle) + for particle in eachparticle(system)] + vtk["neighbor_count"] = system.cache.neighbor_count + vtk["color"] = system.color + end + + if system.surface_tension isa SurfaceTensionMorris || + system.surface_tension isa SurfaceTensionMomentumMorris + surft = zeros((ndims(system), n_moving_particles(system))) + system_coords = current_coordinates(u, system) + + surface_tension_a = surface_tension_model(system) + surface_tension_b = surface_tension_model(system) + nhs = create_neighborhood_search(nothing, system, system) + + foreach_point_neighbor(system, system, + system_coords, system_coords, + nhs) do particle, neighbor, pos_diff, + distance + rho_a = particle_density(v, system, particle) + rho_b = particle_density(v, system, neighbor) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + + surft[1:ndims(system), particle] .+= surface_tension_force(surface_tension_a, + surface_tension_b, + system, system, + particle, neighbor, + pos_diff, distance, + rho_a, rho_b, + grad_kernel) + end + vtk["surface_tension"] = surft + vtk["colorfield"] = system.cache.colorfield + + if system.surface_tension isa SurfaceTensionMorris + vtk["curvature"] = system.cache.curvature + end + if system.surface_tension isa SurfaceTensionMomentumMorris + vtk["surface_stress_tensor"] = system.cache.stress_tensor + end + if system.surface_tension isa SurfaceTensionMorris || + system.surface_tension isa SurfaceTensionMomentumMorris + if !(system.surface_tension.contact_model isa Nothing) + clf = zeros((ndims(system), n_moving_particles(system))) + for particle in each_moving_particle(system) + clf[:, particle] .= contact_force(system.surface_tension.contact_model, + system, particle) + end + + vtk["clf"] = clf + + vtk["d_hat"] = system.cache.d_hat + vtk["delta_wns"] = system.cache.delta_wns + vtk["nu_hat"] = system.cache.nu_hat + vtk["normal_v"] = system.cache.normal_v + vtk["d_vec"] = system.cache.d_vec + end + end + end + if write_meta_data vtk["acceleration"] = system.acceleration vtk["viscosity"] = type2string(system.viscosity) @@ -256,6 +317,8 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) vtk["density_calculator"] = type2string(system.density_calculator) if system isa WeaklyCompressibleSPHSystem + vtk["solver"] = "WCSPH" + vtk["correction_method"] = type2string(system.correction) if system.correction isa AkinciFreeSurfaceCorrection vtk["correction_rho0"] = system.correction.rho0 @@ -371,20 +434,30 @@ function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, system; write_meta_data=true) if write_meta_data vtk["boundary_model"] = "BoundaryModelDummyParticles" - vtk["smoothing_kernel"] = type2string(system.boundary_model.smoothing_kernel) - vtk["smoothing_length"] = system.boundary_model.smoothing_length - vtk["density_calculator"] = type2string(system.boundary_model.density_calculator) - vtk["state_equation"] = type2string(system.boundary_model.state_equation) + vtk["smoothing_kernel"] = type2string(model.smoothing_kernel) + vtk["smoothing_length"] = model.smoothing_length + vtk["density_calculator"] = type2string(model.density_calculator) + vtk["state_equation"] = type2string(model.state_equation) vtk["viscosity_model"] = type2string(model.viscosity) end vtk["hydrodynamic_density"] = [particle_density(v, system, particle) for particle in eachparticle(system)] vtk["pressure"] = model.pressure + vtk["colorfield_bnd"] = model.cache.colorfield_bnd + vtk["colorfield"] = model.cache.colorfield + vtk["neighbor_count"] = model.cache.neighbor_count if model.viscosity isa ViscosityAdami vtk["wall_velocity"] = view(model.cache.wall_velocity, 1:ndims(system), :) end + + if !(system.surface_normal_method isa Nothing) + vtk["normal"] = [surface_normal(system, particle) + for particle in eachparticle(system)] + vtk["tangential"] = [wall_tangential(system, particle) + for particle in eachparticle(system)] + end end function write2vtk!(vtk, v, u, t, system::BoundaryDEMSystem; write_meta_data=true) diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index 11a2d7f56..232b6b772 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -49,7 +49,11 @@ struct BoundaryModelMock end # Mock fluid system - struct FluidSystemMock <: TrixiParticles.FluidSystem{2, Nothing} end + struct FluidSystemMock <: TrixiParticles.FluidSystem{2, Nothing} + surface_tension::Nothing + color::Int64 + FluidSystemMock() = new(nothing, 0) + end kernel = Val(:smoothing_kernel) Base.ndims(::Val{:smoothing_kernel}) = 2 diff --git a/test/schemes/fluid/fluid.jl b/test/schemes/fluid/fluid.jl index cfb75a5c7..570c3bc86 100644 --- a/test/schemes/fluid/fluid.jl +++ b/test/schemes/fluid/fluid.jl @@ -1,6 +1,7 @@ include("weakly_compressible_sph/weakly_compressible_sph.jl") include("rhs.jl") include("pressure_acceleration.jl") +include("surface_normal_sph.jl") include("transport_velocity.jl") include("surface_tension.jl") include("viscosity.jl") diff --git a/test/schemes/fluid/surface_normal_sph.jl b/test/schemes/fluid/surface_normal_sph.jl new file mode 100644 index 000000000..fdf861ddb --- /dev/null +++ b/test/schemes/fluid/surface_normal_sph.jl @@ -0,0 +1,218 @@ +function create_boundary_system(coordinates, particle_spacing, state_equation, kernel, + smoothing_length, NDIMS, walldistance) + # Compute bounding box of fluid particles + xmin = minimum(coordinates[1, :]) + xmax = maximum(coordinates[1, :]) + ymin = minimum(coordinates[2, :]) + ymax = maximum(coordinates[2, :]) + + wall_thickness = 4 * particle_spacing + + if NDIMS == 2 + wall_width = xmax - xmin + wall_size = (wall_width, wall_thickness) + wall_coord = (xmin, ymin - walldistance) + elseif NDIMS == 3 + zmin = minimum(coordinates[3, :]) + wall_width_x = xmax - xmin + wall_width_y = ymax - ymin + wall_size = (wall_width_x, wall_width_y, wall_thickness) + wall_coord = (xmin, ymin, zmin - walldistance) + end + + # Create the wall shape + wall = RectangularShape(particle_spacing, + round.(Int, wall_size ./ particle_spacing), + wall_coord, + density=1000.0) + + boundary_model = BoundaryModelDummyParticles(wall.density, + wall.mass, + state_equation=state_equation, + AdamiPressureExtrapolation(), + kernel, + smoothing_length, + correction=nothing) + + boundary_system = BoundarySPHSystem(wall, boundary_model, adhesion_coefficient=0.0) + return boundary_system +end + +function create_fluid_system(coordinates, velocity, mass, density, particle_spacing, + surface_tension; + surface_normal_method=ColorfieldSurfaceNormal(), + NDIMS=2, smoothing_length=1.0, wall=false, walldistance=0.0, + smoothing_kernel=SchoenbergCubicSplineKernel{NDIMS}()) + tspan = (0.0, 0.01) + + fluid = InitialCondition(coordinates=coordinates, + velocity=velocity, + mass=mass, + density=density, + particle_spacing=particle_spacing) + + state_equation = StateEquationCole(sound_speed=10.0, + reference_density=1000.0, + exponent=1) + + system = WeaklyCompressibleSPHSystem(fluid, SummationDensity(), state_equation, + smoothing_kernel, smoothing_length; + surface_normal_method=surface_normal_method, + reference_particle_spacing=particle_spacing, + surface_tension=surface_tension) + + if wall + boundary_system = create_boundary_system(coordinates, particle_spacing, + state_equation, smoothing_kernel, + smoothing_length, + NDIMS, walldistance) + semi = Semidiscretization(system, boundary_system) + else + semi = Semidiscretization(system) + boundary_system = nothing + end + + ode = semidiscretize(semi, tspan) + TrixiParticles.update_systems_and_nhs(ode.u0.x..., semi, 0.0) + + return system, boundary_system, semi, ode +end + +# TODO: change to rect +function compute_and_test_surface_normals(system, semi, ode; NDIMS=2) + v0_ode, u0_ode = ode.u0.x + v = TrixiParticles.wrap_v(v0_ode, system, semi) + u = TrixiParticles.wrap_u(u0_ode, system, semi) + + # Compute the surface normals + TrixiParticles.compute_surface_normal!(system, system.surface_normal_method, v, u, + v0_ode, u0_ode, semi, 0.0) + + TrixiParticles.remove_invalid_normals!(system, SurfaceTensionMorris(), + system.surface_normal_method) + + # After computation, check that surface normals have been computed and are not NaN or Inf + @test all(isfinite.(system.cache.surface_normal)) + @test all(isfinite.(system.cache.neighbor_count)) + @test size(system.cache.surface_normal, 1) == NDIMS + + nparticles = size(u, 2) + + # Check that the threshold has been applied correctly + threshold = 2^ndims(system) + 1 + + # Test the surface normals based on validity conditions + + for i in 1:nparticles + if system.cache.neighbor_count[i] < threshold + @test all(system.cache.surface_normal[:, i] .== 0.0) + else + # For the linear arrangement, surface normals may still be zero + @test true + end + end +end + +@testset "Surface Normal Computation" begin + # Test case 1: Simple linear arrangement of particles + nparticles = 5 + particle_spacing = 1.0 + NDIMS = 2 + + coordinates = zeros(NDIMS, nparticles) + for i in 1:nparticles + coordinates[:, i] = [i * particle_spacing, 0.0] + end + velocity = zeros(NDIMS, nparticles) + mass = fill(1.0, nparticles) + fluid_density = 1000.0 + density = fill(fluid_density, nparticles) + + system, bnd_system, semi, ode = create_fluid_system(coordinates, velocity, mass, + density, + particle_spacing, nothing; + NDIMS=NDIMS) + + compute_and_test_surface_normals(system, semi, ode; NDIMS=NDIMS) +end + +@testset "Sphere Surface Normals" begin + # Define each variation as a tuple of parameters: + # (NDIMS, smoothing_kernel, particle_spacing, smoothing_length_multiplier, radius, center) + variations = [ + (2, SchoenbergCubicSplineKernel{2}(), 0.25, 3.0, 1.0, (0.0, 0.0)), + (2, SchoenbergCubicSplineKernel{2}(), 0.1, 3.5, 1.0, (0.0, 0.0)), + (3, SchoenbergCubicSplineKernel{3}(), 0.25, 3.0, 1.0, (0.0, 0.0, 0.0)), + (2, WendlandC2Kernel{2}(), 0.3, 3.0, 1.0, (0.0, 0.0)), + (3, WendlandC2Kernel{3}(), 0.3, 3.0, 1.0, (0.0, 0.0, 0.0)) + ] + + for (NDIMS, smoothing_kernel, particle_spacing, smoothing_length_mult, radius, center) in variations + @testset "NDIMS: $(NDIMS), Kernel: $(typeof(smoothing_kernel)), spacing: $(particle_spacing)" begin + smoothing_length = smoothing_length_mult * particle_spacing + + # Create a `SphereShape`, which is a disk in 2D + sphere_ic = SphereShape(particle_spacing, radius, center, 1000.0) + + coordinates = sphere_ic.coordinates + velocity = zeros(NDIMS, size(coordinates, 2)) + mass = sphere_ic.mass + density = sphere_ic.density + + system, bnd_system, semi, ode = create_fluid_system(coordinates, velocity, mass, + density, + particle_spacing, nothing; + NDIMS=NDIMS, + smoothing_length=smoothing_length, + smoothing_kernel=smoothing_kernel, + surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, + ideal_density_threshold=0.9), + wall=true, walldistance=2.0) + + compute_and_test_surface_normals(system, semi, ode; NDIMS=NDIMS) + + nparticles = size(coordinates, 2) + expected_normals = zeros(NDIMS, nparticles) + surface_particles = Int[] + + # Compute expected normals and identify surface particles + for i in 1:nparticles + pos = coordinates[:, i] + r = pos .- center + norm_r = norm(r) + + # If particle is on the circumference of the circle + if abs(norm_r - radius) < particle_spacing + expected_normals[:, i] = -r / norm_r + push!(surface_particles, i) + else + expected_normals[:, i] .= 0.0 + end + end + + # Normalize computed normals + computed_normals = copy(system.cache.surface_normal) + for i in surface_particles + norm_computed = norm(computed_normals[:, i]) + if norm_computed > 0 + computed_normals[:, i] /= norm_computed + end + end + + # Boundary system + bnd_color = bnd_system.boundary_model.cache.colorfield_bnd + # this is only true since it assumed that the color is 1 + @test all(bnd_color .>= 0.0) + + # Test that computed normals match expected normals + for i in surface_particles + @test isapprox(computed_normals[:, i], expected_normals[:, i], atol=0.04) + end + + # Optionally, test that interior particles have near-zero normals + # for i in setdiff(1:nparticles, surface_particles) + # @test isapprox(norm(system.cache.surface_normal[:, i]), 0.0, atol=1e-4) + # end + end + end +end diff --git a/test/systems/boundary_system.jl b/test/systems/boundary_system.jl index e483ec3f5..89f526dd5 100644 --- a/test/systems/boundary_system.jl +++ b/test/systems/boundary_system.jl @@ -111,7 +111,7 @@ system = BoundarySPHSystem(initial_condition, model) - show_compact = "BoundarySPHSystem{2}((hydrodynamic_mass = 3,), nothing, 0.0) with 2 particles" + show_compact = "BoundarySPHSystem{2}((hydrodynamic_mass = 3,), nothing, 0.0, 0) with 2 particles" @test repr(system) == show_compact show_box = """ @@ -122,6 +122,7 @@ │ boundary model: ……………………………………… (hydrodynamic_mass = 3,) │ │ movement function: ……………………………… nothing │ │ adhesion coefficient: ……………………… 0.0 │ + │ color: ……………………………………………………………… 0 │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box end diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index 9be9feeb9..887e77444 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -128,7 +128,7 @@ system = EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, sound_speed) - show_compact = "EntropicallyDampedSPHSystem{2}(SummationDensity(), nothing, Val{:smoothing_kernel}(), [0.0, 0.0]) with 2 particles" + show_compact = "EntropicallyDampedSPHSystem{2}(SummationDensity(), nothing, Val{:smoothing_kernel}(), [0.0, 0.0], nothing, nothing) with 2 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -141,6 +141,8 @@ │ smoothing kernel: ………………………………… Val │ │ tansport velocity formulation: Nothing │ │ acceleration: …………………………………………… [0.0, 0.0] │ + │ surface tension: …………………………………… nothing │ + │ surface normal method: …………………… nothing │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box end diff --git a/test/systems/wcsph_system.jl b/test/systems/wcsph_system.jl index 46b4900cd..533af89f4 100644 --- a/test/systems/wcsph_system.jl +++ b/test/systems/wcsph_system.jl @@ -193,7 +193,7 @@ smoothing_length, density_diffusion=density_diffusion) - show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), nothing, [0.0, 0.0], nothing) with 2 particles" + show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), nothing, nothing, [0.0, 0.0], nothing) with 2 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -207,6 +207,7 @@ │ viscosity: …………………………………………………… nothing │ │ density diffusion: ……………………………… Val{:density_diffusion}() │ │ surface tension: …………………………………… nothing │ + │ surface normal method: …………………… nothing │ │ acceleration: …………………………………………… [0.0, 0.0] │ │ source terms: …………………………………………… Nothing │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""