diff --git a/.gitignore b/.gitignore index be4f7441f..f7151b63f 100644 --- a/.gitignore +++ b/.gitignore @@ -18,3 +18,4 @@ run *.vtu *.pvd *.mp4 +*.csv diff --git a/docs/src/install.md b/docs/src/install.md index c633634f3..110aef029 100644 --- a/docs/src/install.md +++ b/docs/src/install.md @@ -30,11 +30,28 @@ The advantage of using a separate `run` directory is that you can also add other related packages (e.g., OrdinaryDiffEq.jl, see above) to the project in the `run` folder and always have a reproducible environment at hand to share with others. -## Optional Software/Packages -- [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) -- a Julia package of ordinary differential equation solvers that is used in the examples +## Optional software/packages +- [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) -- A Julia package of ordinary differential equation solvers that is used in the examples - [Plots.jl](https://github.com/JuliaPlots/Plots.jl) -- Julia Plotting library that is used in some examples - [PythonPlot.jl](https://github.com/JuliaPy/PythonPlot.jl) -- Plotting library that can be used instead of Plots.jl - [ParaView](https://www.paraview.org/) -- Software that can be used for visualization of results ## Common issues +If you followed the [installation instructions for developers](@ref for-developers) and you +run into any problems with packages when pulling the latest version of TrixiParticles.jl, +start Julia with the project in the `run` folder, +```bash + julia --project=run +``` +update all packages in that project, resolve all conflicts in the project, and install all +new dependencies: +```julia +julia> using Pkg + +julia> Pkg.update() + +julia> Pkg.resolve() + +julia> Pkg.instantiate() +``` diff --git a/examples/fluid/oscillating_drop_2d.jl b/examples/fluid/oscillating_drop_2d.jl new file mode 100644 index 000000000..da4e855ac --- /dev/null +++ b/examples/fluid/oscillating_drop_2d.jl @@ -0,0 +1,93 @@ +# Oscillating drop in a central force field, based on +# +# J. J. Monaghan, Ashkan Rafiee. +# "A Simple SPH Algorithm for Multi-Fluid Flow with High Density Ratios." +# In: International Journal for Numerical Methods in Fluids 71, no. 5 (2013), pages 537-61. +# https://doi.org/10.1002/fld.3671. + +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.05 + +# ========================================================================================== +# ==== Experiment Setup +radius = 1.0 +sigma = 0.5 +# Make this a constant because global variables in the source terms are slow +const OMEGA = 1.0 + +source_terms = (coords, velocity, density, pressure) -> -OMEGA^2 * coords + +# 1 period in the exact solution as computed below (but integrated with a small timestep) +period = 4.567375 +tspan = (0.0, 1period) + +fluid_density = 1000.0 +sound_speed = 10.0 +state_equation = StateEquationCole(; sound_speed, exponent=7, + reference_density=fluid_density) + +# Equation A.19 in the paper rearranged. +# sigma^2 = Q / rho * (a^2 + b^2) / (a^2 b^2) - OMEGA^2. +Q = (sigma^2 + OMEGA^2) * fluid_density / 2 +pressure = coords -> Q * (1 - coords[1]^2 - coords[2]^2) +density = coords -> TrixiParticles.inverse_state_equation(state_equation, pressure(coords)) + +fluid = SphereShape(fluid_particle_spacing, radius, (0.0, 0.0), + density, pressure=pressure, + sphere_type=RoundSphere(), + velocity=coords -> sigma .* (coords[1], -coords[2])) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 3.0 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() +viscosity = ArtificialViscosityMonaghan(alpha=0.01, beta=0.0) + +density_diffusion = DensityDiffusionAntuono(fluid, delta=0.1) +fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity, + density_diffusion=density_diffusion, + source_terms=source_terms) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.04, prefix="") + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +sol = solve(ode, RDPK3SpFSAL49(), + abstol=1e-7, # Default abstol is 1e-6 (may need to be tuned to prevent intabilities) + reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent intabilities) + save_everystep=false, callback=callbacks); + +@inline function exact_solution_rhs(u, p, t) + sigma, A, B = u + + dsigma = (sigma^2 + OMEGA^2) * ((B^2 - A^2) / (B^2 + A^2)) + dA = sigma * A + dB = -sigma * B + + return SVector(dsigma, dA, dB) +end + +exact_u0 = SVector(sigma, radius, radius) +exact_solution_ode = ODEProblem(exact_solution_rhs, exact_u0, tspan) + +# Use the same time integrator to avoid compilation of another integrator in CI +sol_exact = solve(exact_solution_ode, RDPK3SpFSAL49(), save_everystep=false) + +# Error in the semi-major axis of the elliptical drop +error_A = maximum(sol.u[end].x[2]) + 0.5fluid_particle_spacing - + maximum(sol_exact.u[end][2:3]) diff --git a/examples/fsi/dam_break_2d.jl b/examples/fsi/dam_break_2d.jl index 0bb668fb6..8126e9759 100644 --- a/examples/fsi/dam_break_2d.jl +++ b/examples/fsi/dam_break_2d.jl @@ -88,8 +88,8 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) # ========================================================================================== # ==== Solid -solid_smoothing_length = sqrt(2) * solid_particle_spacing -solid_smoothing_kernel = SchoenbergCubicSplineKernel{2}() +solid_smoothing_length = 2 * sqrt(2) * solid_particle_spacing +solid_smoothing_kernel = WendlandC2Kernel{2}() # For the FSI we need the hydrodynamic masses and densities in the solid boundary model hydrodynamic_densites = fluid_density * ones(size(solid.density)) diff --git a/examples/fsi/dam_break_gate_2d.jl b/examples/fsi/dam_break_gate_2d.jl index 6ba326e50..f82d780eb 100644 --- a/examples/fsi/dam_break_gate_2d.jl +++ b/examples/fsi/dam_break_gate_2d.jl @@ -117,8 +117,8 @@ boundary_system_gate = BoundarySPHSystem(gate, boundary_model_gate, movement=gat # ========================================================================================== # ==== Solid -solid_smoothing_length = sqrt(2) * solid_particle_spacing -solid_smoothing_kernel = SchoenbergCubicSplineKernel{2}() +solid_smoothing_length = 2 * sqrt(2) * solid_particle_spacing +solid_smoothing_kernel = WendlandC2Kernel{2}() # For the FSI we need the hydrodynamic masses and densities in the solid boundary model hydrodynamic_densites = fluid_density * ones(size(solid.density)) diff --git a/examples/solid/oscillating_beam_2d.jl b/examples/solid/oscillating_beam_2d.jl index c396eec65..b542d5ed7 100644 --- a/examples/solid/oscillating_beam_2d.jl +++ b/examples/solid/oscillating_beam_2d.jl @@ -45,12 +45,12 @@ solid = union(beam, fixed_particles) # ========================================================================================== # ==== Solid +# The kernel in the reference uses a differently scaled smoothing length, +# so this is equivalent to the smoothing length of `sqrt(2) * particle_spacing` used in the paper. +smoothing_length = 2 * sqrt(2) * particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() -smoothing_length = sqrt(2) * particle_spacing -smoothing_kernel = SchoenbergCubicSplineKernel{2}() - -solid_system = TotalLagrangianSPHSystem(solid, - smoothing_kernel, smoothing_length, +solid_system = TotalLagrangianSPHSystem(solid, smoothing_kernel, smoothing_length, E, nu, n_fixed_particles=nparticles(fixed_particles), acceleration=(0.0, -gravity)) @@ -61,7 +61,28 @@ semi = Semidiscretization(solid_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) -saving_callback = SolutionSavingCallback(dt=0.02, prefix="") + +# Track the position of the particle in the middle of the tip of the beam. +particle_id = Int(n_particles_per_dimension[1] * (n_particles_per_dimension[2] + 1) / 2) + +shift_x = beam.coordinates[1, particle_id] +shift_y = beam.coordinates[2, particle_id] + +function x_deflection(v, u, t, system) + particle_position = TrixiParticles.current_coords(u, system, particle_id) + + return particle_position[1] - shift_x +end + +function y_deflection(v, u, t, system) + particle_position = TrixiParticles.current_coords(u, system, particle_id) + + return particle_position[2] - shift_y +end + +saving_callback = SolutionSavingCallback(dt=0.02, prefix="", + x_deflection=x_deflection, + y_deflection=y_deflection) callbacks = CallbackSet(info_callback, saving_callback) diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index 53496603a..996c31433 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -100,9 +100,6 @@ struct InitialCondition{ELTYPE} if velocity isa AbstractMatrix velocities = velocity - if size(coordinates) != size(velocities) - throw(ArgumentError("`coordinates` and `velocities` must be of the same size")) - end else # Assuming `velocity` is a scalar or a function velocity_fun = wrap_function(velocity, Val(NDIMS)) @@ -111,7 +108,10 @@ struct InitialCondition{ELTYPE} "for $NDIMS-dimensional `coordinates`")) end velocities_svector = velocity_fun.(coordinates_svector) - velocities = reinterpret(reshape, ELTYPE, velocities_svector) + velocities = stack(velocities_svector) + end + if size(coordinates) != size(velocities) + throw(ArgumentError("`coordinates` and `velocities` must be of the same size")) end if density isa AbstractVector diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index bce3fba90..43f87215f 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -33,15 +33,15 @@ The zero subscript on quantities denotes that the quantity is to be measured in The difference in the initial coordinates is denoted by $\bm{X}_{ab} = \bm{X}_a - \bm{X}_b$, the difference in the current coordinates is denoted by $\bm{x}_{ab} = \bm{x}_a - \bm{x}_b$. -For the computation of the PK1 stress tensor, the deformation gradient $\bm{J}$ is computed per particle as +For the computation of the PK1 stress tensor, the deformation gradient $\bm{F}$ is computed per particle as ```math -\bm{J}_a = \sum_b \frac{m_{0b}}{\rho_{0b}} \bm{x}_{ba} (\bm{L}_{0a}\nabla_{0a} W(\bm{X}_{ab}))^T \\ +\bm{F}_a = \sum_b \frac{m_{0b}}{\rho_{0b}} \bm{x}_{ba} (\bm{L}_{0a}\nabla_{0a} W(\bm{X}_{ab}))^T \\ \qquad = -\left(\sum_b \frac{m_{0b}}{\rho_{0b}} \bm{x}_{ab} (\nabla_{0a} W(\bm{X}_{ab}))^T \right) \bm{L}_{0a}^T ``` with $1 \leq i,j \leq d$. From the deformation gradient, the Green-Lagrange strain ```math -\bm{E} = \frac{1}{2}(\bm{J}^T\bm{J} - \bm{I}) +\bm{E} = \frac{1}{2}(\bm{F}^T\bm{F} - \bm{I}) ``` and the second Piola-Kirchhoff stress tensor ```math @@ -49,7 +49,7 @@ and the second Piola-Kirchhoff stress tensor ``` are computed to obtain the PK1 stress tensor as ```math -\bm{P} = \bm{J}\bm{S}. +\bm{P} = \bm{F}\bm{S}. ``` Here, @@ -285,8 +285,8 @@ end calc_deformation_grad!(deformation_grad, neighborhood_search, system) @threaded for particle in eachparticle(system) - J_particle = deformation_gradient(system, particle) - pk1_particle = pk1_stress_tensor(J_particle, system) + F_particle = deformation_gradient(system, particle) + pk1_particle = pk1_stress_tensor(F_particle, system) pk1_particle_corrected = pk1_particle * correction_matrix(system, particle) @inbounds for j in 1:ndims(system), i in 1:ndims(system) @@ -332,18 +332,18 @@ end end # First Piola-Kirchhoff stress tensor -@inline function pk1_stress_tensor(J, system) - S = pk2_stress_tensor(J, system) +@inline function pk1_stress_tensor(F, system) + S = pk2_stress_tensor(F, system) - return J * S + return F * S end # Second Piola-Kirchhoff stress tensor -@inline function pk2_stress_tensor(J, system) +@inline function pk2_stress_tensor(F, system) (; lame_lambda, lame_mu) = system # Compute the Green-Lagrange strain - E = 0.5 * (transpose(J) * J - I) + E = 0.5 * (transpose(F) * F - I) return lame_lambda * tr(E) * I + 2 * lame_mu * E end @@ -411,3 +411,47 @@ end function viscosity_model(system::TotalLagrangianSPHSystem) return system.boundary_model.viscosity end + +# An explanation of these equation can be found in +# J. Lubliner, 2008. Plasticity theory. +# See here below Equation 5.3.21 for the equation for the equivalent stress. +# The von-Mises stress is one form of equivalent stress, where sigma is the deviatoric stress. +# See pages 32 and 123. +function von_mises_stress(system::TotalLagrangianSPHSystem) + von_mises_stress = zeros(eltype(system.pk1_corrected), nparticles(system)) + + @threaded for particle in each_moving_particle(system) + F = deformation_gradient(system, particle) + J = det(F) + P = pk1_corrected(system, particle) + sigma = (1.0 / J) * P * F' + + # Calculate deviatoric stress tensor + s = sigma - (1.0 / 3.0) * tr(sigma) * I + + von_mises_stress[particle] = sqrt(3.0 / 2.0 * sum(s .^ 2)) + end + + return von_mises_stress +end + +# An explanation of these equation can be found in +# J. Lubliner, 2008. Plasticity theory. +# See here page 473 for the relation between the `pk1`, the first Piola-Kirchhoff tensor, +# and the Cauchy stress. +function cauchy_stress(system::TotalLagrangianSPHSystem) + NDIMS = ndims(system) + + cauchy_stress_tensors = zeros(eltype(system.pk1_corrected), NDIMS, NDIMS, + nparticles(system)) + + @threaded for particle in each_moving_particle(system) + F = deformation_gradient(system, particle) + J = det(F) + P = pk1_corrected(system, particle) + sigma = (1.0 / J) * P * F' + cauchy_stress_tensors[:, :, particle] = sigma + end + + return cauchy_stress_tensors +end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 11b88524a..bad071243 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -216,11 +216,25 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d vtk["velocity"] = hcat(view(v, 1:ndims(system), :), zeros(ndims(system), n_fixed_particles)) + vtk["jacobian"] = [det(deformation_gradient(system, particle)) + for particle in eachparticle(system)] + + vtk["von_mises_stress"] = von_mises_stress(system) + + sigma = cauchy_stress(system) + vtk["sigma_11"] = sigma[1, 1, :] + vtk["sigma_22"] = sigma[2, 2, :] + if ndims(system) == 3 + vtk["sigma_33"] = sigma[3, 3, :] + end + vtk["material_density"] = system.material_density vtk["young_modulus"] = system.young_modulus vtk["poisson_ratio"] = system.poisson_ratio vtk["lame_lambda"] = system.lame_lambda vtk["lame_mu"] = system.lame_mu + vtk["smoothing_kernel"] = type2string(system.smoothing_kernel) + vtk["smoothing_length"] = system.smoothing_length write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data) end diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 3df9b7d40..ac4f6218f 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -2,6 +2,15 @@ # but without checking the correctness of the solution. @testset verbose=true "Examples" begin @testset verbose=true "Fluid" begin + @trixi_testset "fluid/oscillating_drop_2d.jl" begin + @test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "oscillating_drop_2d.jl")) + @test sol.retcode == ReturnCode.Success + # This is the error on an Apple M2 Pro. We need this tolerance to make CI pass. + @test isapprox(error_A, 0.0001717690010767381, atol=1e-8) + end + @trixi_testset "fluid/rectangular_tank_2d.jl" begin @test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", diff --git a/test/systems/solid_system.jl b/test/systems/solid_system.jl index a58720991..0d512c205 100644 --- a/test/systems/solid_system.jl +++ b/test/systems/solid_system.jl @@ -338,4 +338,38 @@ @test v0 == velocity end + + @testset verbose=true "compute_von_mises_stress" begin + # System setup + coordinates = [1.0 2.0; 1.0 2.0] + velocity = zero(coordinates) + mass = [1.25, 1.5] + material_densities = [990.0, 1000.0] + smoothing_kernel = Val(:smoothing_kernel) + smoothing_length = 0.362 + nu = 0.25 # Poisson's ratio + E = 2.5 # Young's modulus + boundary_model = Val(:boundary_model) + + initial_condition = InitialCondition(; coordinates, velocity, mass, + density=material_densities) + system = TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, + smoothing_length, E, nu, boundary_model) + + # Initialize deformation_grad and pk1_corrected with arbitrary values + for particle in TrixiParticles.eachparticle(system) + system.deformation_grad[:, :, particle] = [1.0 0.2; 0.2 1.0] + system.pk1_corrected[:, :, particle] = [1.0 0.5; 0.5 1.0] + end + + von_mises_stress = TrixiParticles.von_mises_stress(system) + cauchy_stress = TrixiParticles.cauchy_stress(system) + + reference_stress_tensor = [1.145833 0.729167; 0.729167 1.145833;;; + 1.145833 0.729167; 0.729167 1.145833] + + # Verify against calculation by hand + @test isapprox(von_mises_stress[1], 1.4257267477533202, atol=1e-14) + @test isapprox(reference_stress_tensor, cauchy_stress, atol=1e-6) + end end