Add Carreau-Yasuda Poiseuille validation

Author: Sherin Joseph

examples/fluid/poiseuille_carreau_2d.jl

@@ -0,0 +1,285 @@
+using TrixiParticles
+using OrdinaryDiffEqLowStorageRK
+
+# ==========================================================================================
+# 2D Periodic Poiseuille Flow with Carreau-Yasuda Viscosity
+#
+# This example simulates pressure-driven channel flow with the
+# `ViscosityCarreauYasuda` non-Newtonian viscosity model. The driving pressure
+# gradient is represented by an equivalent body acceleration.
+# ==========================================================================================
+
+# ==========================================================================================
+# ==== Resolution
+ny = 50
+
+# ==========================================================================================
+# ==== Experiment Setup
+t_end_factor = 0.1
+eps_factor = 1.0
+sound_speed_factor = 60.0
+initial_condition_mode = "analytical"
+power_law_index = 1.0
+
+channel_height = 1.0
+channel_length = 6.0 * channel_height
+particle_spacing = channel_height / ny
+boundary_layers = 5
+
+fluid_density = 1000.0
+nu0 = 1.0e-3
+nu_inf = 0.0
+lambda_exponent = 2.0
+reynolds_number = 200.0
+
+reference_velocity = reynolds_number * nu0 / channel_height
+pressure_gradient = 8.0 * fluid_density * reference_velocity^2 /
+                    (reynolds_number * channel_height)
+acceleration_x = pressure_gradient / fluid_density
+carreau_time_constant = channel_height / max(reference_velocity, eps())
+
+t_end = t_end_factor * channel_height / max(reference_velocity, eps())
+tspan = (0.0, t_end)
+
+if !(initial_condition_mode in ("newtonian", "analytical", "zero"))
+    throw(ArgumentError("initial condition mode must be \"newtonian\", " *
+                        "\"analytical\", or \"zero\""))
+end
+
+# ==========================================================================================
+# ==== Analytical Solution
+function linear_interpolation_clamped(x, y, interpolation_point)
+    interpolation_point <= first(x) && return first(y)
+    interpolation_point >= last(x) && return last(y)
+
+    i = searchsortedlast(x, interpolation_point)
+    x0, x1 = x[i], x[i + 1]
+    y0, y1 = y[i], y[i + 1]
+    return y0 + (y1 - y0) * (interpolation_point - x0) / (x1 - x0)
+end
+
+function carreau_yasuda_kinematic_viscosity(shear_rate, nu0, nu_inf,
+                                            time_constant, lambda_exponent,
+                                            power_law_index)
+    return nu_inf + (nu0 - nu_inf) *
+           (1.0 + (time_constant * shear_rate)^lambda_exponent)^((power_law_index -
+                                                                    1.0) /
+                                                                   lambda_exponent)
+end
+
+function solve_shear_rate_from_stress(shear_stress, density, nu0, nu_inf,
+                                      time_constant, lambda_exponent,
+                                      power_law_index)
+    shear_stress <= 0 && return 0.0
+
+    residual(shear_rate) = density *
+                           carreau_yasuda_kinematic_viscosity(shear_rate, nu0,
+                                                              nu_inf,
+                                                              time_constant,
+                                                              lambda_exponent,
+                                                              power_law_index) *
+                           shear_rate - shear_stress
+
+    lower = 0.0
+    upper = 1.0
+    while residual(upper) < 0.0
+        upper *= 2.0
+        upper > 1.0e12 &&
+            error("failed to bracket shear-rate root for shear stress $shear_stress")
+    end
+
+    for _ in 1:120
+        middle = 0.5 * (lower + upper)
+        residual_middle = residual(middle)
+
+        if abs(residual_middle) <= 1.0e-12 * max(shear_stress, 1.0)
+            return middle
+        elseif residual_middle > 0
+            upper = middle
+        else
+            lower = middle
+        end
+    end
+
+    return 0.5 * (lower + upper)
+end
+
+function analytical_ux_profile(y_positions, power_law_index, channel_height,
+                               density, nu0, nu_inf, time_constant,
+                               lambda_exponent, pressure_gradient)
+    distances_to_centerline = sort(unique(abs.(y_positions .- 0.5 * channel_height)))
+    shear_rates = similar(distances_to_centerline)
+
+    for i in eachindex(distances_to_centerline)
+        shear_stress = pressure_gradient * distances_to_centerline[i]
+        shear_rates[i] = solve_shear_rate_from_stress(shear_stress, density, nu0,
+                                                      nu_inf, time_constant,
+                                                      lambda_exponent,
+                                                      power_law_index)
+    end
+
+    velocity_at_distance = zeros(length(distances_to_centerline))
+    for i in (lastindex(distances_to_centerline) - 1):-1:firstindex(distances_to_centerline)
+        ds = distances_to_centerline[i + 1] - distances_to_centerline[i]
+        velocity_at_distance[i] = velocity_at_distance[i + 1] +
+                                  0.5 * (shear_rates[i + 1] + shear_rates[i]) * ds
+    end
+
+    velocity = Vector{Float64}(undef, length(y_positions))
+    for (i, y) in pairs(y_positions)
+        distance_to_centerline = abs(y - 0.5 * channel_height)
+        velocity[i] = linear_interpolation_clamped(distances_to_centerline,
+                                                   velocity_at_distance,
+                                                   distance_to_centerline)
+    end
+
+    return velocity
+end
+
+function newtonian_ux(y, channel_height, density, nu0, pressure_gradient)
+    return pressure_gradient / (2.0 * density * nu0) * y * (channel_height - y)
+end
+
+function l2_velocity_error(system::TrixiParticles.AbstractFluidSystem,
+                           dv_ode, du_ode, v_ode, u_ode, semi, t)
+    v = TrixiParticles.wrap_v(v_ode, system, semi)
+    u = TrixiParticles.wrap_u(u_ode, system, semi)
+
+    y_positions = [TrixiParticles.current_coords(u, system, particle)[2]
+                   for particle in TrixiParticles.eachparticle(system)]
+    analytical_velocity = analytical_ux_profile(y_positions, power_law_index,
+                                                channel_height, fluid_density, nu0,
+                                                nu_inf, carreau_time_constant,
+                                                lambda_exponent, pressure_gradient)
+
+    squared_error = 0.0
+    squared_reference = 0.0
+    for (i, particle) in enumerate(TrixiParticles.eachparticle(system))
+        ux = TrixiParticles.current_velocity(v, system, particle)[1]
+        squared_error += (ux - analytical_velocity[i])^2
+        squared_reference += analytical_velocity[i]^2
+    end
+
+    return sqrt(squared_error / nparticles(system)) /
+           (sqrt(squared_reference / nparticles(system)) + eps())
+end
+l2_velocity_error(system, dv_ode, du_ode, v_ode, u_ode, semi, t) = nothing
+
+function max_velocity_error(system::TrixiParticles.AbstractFluidSystem,
+                            dv_ode, du_ode, v_ode, u_ode, semi, t)
+    v = TrixiParticles.wrap_v(v_ode, system, semi)
+    u = TrixiParticles.wrap_u(u_ode, system, semi)
+
+    y_positions = [TrixiParticles.current_coords(u, system, particle)[2]
+                   for particle in TrixiParticles.eachparticle(system)]
+    analytical_velocity = analytical_ux_profile(y_positions, power_law_index,
+                                                channel_height, fluid_density, nu0,
+                                                nu_inf, carreau_time_constant,
+                                                lambda_exponent, pressure_gradient)
+
+    error = 0.0
+    for (i, particle) in enumerate(TrixiParticles.eachparticle(system))
+        ux = TrixiParticles.current_velocity(v, system, particle)[1]
+        error = max(error, abs(ux - analytical_velocity[i]))
+    end
+
+    return error
+end
+max_velocity_error(system, dv_ode, du_ode, v_ode, u_ode, semi, t) = nothing
+
+# ==========================================================================================
+# ==== Initial Condition
+initial_velocity = if initial_condition_mode == "zero"
+    (0.0, 0.0)
+elseif initial_condition_mode == "analytical"
+    y_reference = collect(range(0.0, channel_height; length=4 * ny + 1))
+    ux_reference = analytical_ux_profile(y_reference, power_law_index,
+                                         channel_height, fluid_density, nu0,
+                                         nu_inf, carreau_time_constant,
+                                         lambda_exponent, pressure_gradient)
+    x -> (linear_interpolation_clamped(y_reference, ux_reference, x[2]), 0.0)
+else
+    x -> (newtonian_ux(x[2], channel_height, fluid_density, nu0,
+                       pressure_gradient), 0.0)
+end
+
+# ==========================================================================================
+# ==== Fluid
+tank = RectangularTank(particle_spacing, (channel_length, channel_height),
+                       (channel_length, channel_height), fluid_density;
+                       n_layers=boundary_layers,
+                       faces=(false, false, true, true),
+                       velocity=initial_velocity,
+                       coordinates_eltype=Float64)
+
+smoothing_length = 1.2 * particle_spacing
+smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+
+sound_speed = sound_speed_factor * reference_velocity
+state_equation = StateEquationCole(; sound_speed,
+                                    reference_density=fluid_density,
+                                    exponent=7)
+
+viscosity = ViscosityCarreauYasuda(; nu0, nu_inf,
+                                   lambda=carreau_time_constant,
+                                   a=lambda_exponent,
+                                   n=power_law_index,
+                                   epsilon=max(0.5, eps_factor) * particle_spacing)
+
+fluid_system = WeaklyCompressibleSPHSystem(tank.fluid;
+                                           density_calculator=ContinuityDensity(),
+                                           state_equation,
+                                           smoothing_kernel,
+                                           smoothing_length,
+                                           acceleration=(acceleration_x, 0.0),
+                                           viscosity,
+                                           shifting_technique=nothing)
+
+# ==========================================================================================
+# ==== Boundary
+boundary_model = BoundaryModelDummyParticles(tank.boundary.density,
+                                             tank.boundary.mass,
+                                             AdamiPressureExtrapolation(),
+                                             smoothing_kernel,
+                                             smoothing_length;
+                                             state_equation,
+                                             viscosity)
+
+boundary_system = WallBoundarySystem(tank.boundary, boundary_model)
+
+# ==========================================================================================
+# ==== Simulation
+periodic_box = PeriodicBox(min_corner=[0.0, -10.0 * channel_height],
+                           max_corner=[channel_length, 10.0 * channel_height])
+neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box)
+
+semi = Semidiscretization(fluid_system, boundary_system;
+                          neighborhood_search,
+                          parallelization_backend=PolyesterBackend())
+
+ode = semidiscretize(semi, tspan)
+
+n_label = replace(string(power_law_index), "." => "p")
+output_directory = joinpath("out_poiseuille_carreau", "n_$power_law_index")
+result_filename = "validation_run_poiseuille_carreau_2d_n_$(n_label)_ny_$ny"
+
+info_callback = InfoCallback(interval=200)
+saving_callback = SolutionSavingCallback(; dt=t_end / 20,
+                                         prefix="",
+                                         output_directory)
+pp_callback = PostprocessCallback(; dt=t_end / 20,
+                                  output_directory,
+                                  filename=result_filename,
+                                  l2_velocity_error,
+                                  max_velocity_error,
+                                  write_csv=true)
+cfl_callback = StepsizeCallback(cfl=0.2)
+callbacks = CallbackSet(info_callback, saving_callback, pp_callback,
+                        cfl_callback, UpdateCallback())
+
+sol = solve(ode, RDPK3SpFSAL35();
+            abstol=1.0e-7,
+            reltol=1.0e-4,
+            save_everystep=false,
+            callback=callbacks,
+            maxiters=2_000_000);

validation/poiseuille_carreau_2d/README.md

@@ -0,0 +1,112 @@
+# 2D Poiseuille Flow with Carreau-Yasuda Viscosity
+
+This folder contains a periodic 2D Poiseuille validation case for the
+`ViscosityCarreauYasuda` non-Newtonian viscosity model in TrixiParticles.jl.
+
+The case checks whether the numerical solution approaches the steady
+one-dimensional Carreau-Yasuda velocity profile in a pressure-driven channel.
+It also includes `n = 1`, which recovers the Newtonian parabolic Poiseuille
+profile. The setup is inspired by the Poiseuille validation in section 3.1 of
+Coclite et al. (2020), but uses WCSPH with an equivalent body acceleration
+instead of the MRT-LBM pressure-gradient implementation used there.
+
+## Files
+
+- `../../examples/fluid/poiseuille_carreau_2d.jl`: reusable example setup for a
+  single Carreau-Yasuda power-law index.
+- `validation_poiseuille_carreau_2d.jl`: validation runner. It runs one or more
+  Carreau-Yasuda power-law indices, writes VTU output plus error histories, and
+  checks the final relative L2 errors against validation bounds.
+- `plot_carreau_comparison.jl`: helper script for manual inspection of the
+  saved VTU profiles and error trends.
+
+## Setup
+
+- Channel height: `H = 1.0`
+- Channel length: `L = 6H`
+- Periodic direction: `x`
+- Solid walls: lower and upper `y` boundaries
+- Reference density: `rho0 = 1000.0`
+- Zero-shear kinematic viscosity: `nu0 = 1.0e-3`
+- Infinite-shear kinematic viscosity: `nu_inf = 0.0`
+- Yasuda transition parameter: `a = 2.0`
+- Default power-law indices: `n = (1.0, 1.5, 0.5, 0.25)`
+
+## Analytical Profile
+
+For each `n`, the script computes the steady profile by solving the implicit
+Carreau-Yasuda stress relation
+
+```text
+tau(y) = rho0 * nu(gammadot) * gammadot
+```
+
+where `tau(y) = dpdx * abs(y - H / 2)`. The resulting shear-rate profile is
+integrated from the wall to the centerline to obtain `u_x(y)`.
+
+The validation records two quantities with `PostprocessCallback`:
+
+- relative L2 velocity error against the analytical profile
+- maximum absolute velocity error against the analytical profile
+
+## Running
+
+Default run:
+
+```bash
+julia --project=. validation/poiseuille_carreau_2d/validation_poiseuille_carreau_2d.jl
+```
+
+The default uses `ny = 50`, `t_end_factor = 0.1`, and analytical initial
+conditions so it can be run quickly while still checking the numerical profile
+against the analytical solution.
+
+High-resolution run matching the paper-inspired channel resolution:
+
+```bash
+julia --project=. validation/poiseuille_carreau_2d/validation_poiseuille_carreau_2d.jl 200 0.02 1.0 60.0 analytical 0.25,0.5,1.0,1.5
+```
+
+Optional command line arguments are:
+
+```text
+ny t_end_factor eps_factor sound_speed_factor initial_condition_mode n_values
+```
+
+where `initial_condition_mode` is one of `newtonian`, `analytical`, or `zero`,
+and `n_values` is a comma-separated list such as `1.0,0.5,0.25`.
+
+Results are written to:
+
+```text
+out_poiseuille_carreau/n_<n>/
+```
+
+Each case contains `fluid_1.pvd`/VTU output and
+`validation_run_poiseuille_carreau_2d_*.csv/.json` error histories.
+
+The validation runner checks the final relative L2 error with the following
+bounds:
+
+```text
+n = 0.25: relative L2 <= 0.06
+n = 0.5:  relative L2 <= 0.06
+n = 1.0:  relative L2 <= 0.06
+n = 1.5:  relative L2 <= 0.06
+```
+
+The relative L2 error is the main quantitative comparison metric. The maximum
+absolute error is also written to the CSV/JSON files, but it should not be
+compared directly across different `n` values without normalization because the
+velocity scale changes strongly for shear-thinning cases such as `n = 0.25`.
+
+## Plotting
+
+After a run, create comparison plots with:
+
+```bash
+julia --project=. -e 'push!(LOAD_PATH, "test"); append!(ARGS, ["out_poiseuille_carreau"]); include("validation/poiseuille_carreau_2d/plot_carreau_comparison.jl")'
+```
+
+The plotting helper uses optional plotting packages that are available through
+the test environment. They are not runtime dependencies of TrixiParticles.jl.