diff --git a/diffmpm/materials/__init__.py b/diffmpm/materials/__init__.py index ce35083..028e715 100644 --- a/diffmpm/materials/__init__.py +++ b/diffmpm/materials/__init__.py @@ -1,3 +1,4 @@ from diffmpm.materials._base import _Material from diffmpm.materials.simple import SimpleMaterial from diffmpm.materials.linear_elastic import LinearElastic +from diffmpm.materials.newtonian import Newtonian diff --git a/diffmpm/materials/newtonian.py b/diffmpm/materials/newtonian.py index e5a0d9b..558832f 100644 --- a/diffmpm/materials/newtonian.py +++ b/diffmpm/materials/newtonian.py @@ -1 +1,114 @@ -#!/usr/bin/env python3 +import jax.numpy as jnp +from jax import Array, lax +from jax.typing import ArrayLike + +from ._base import _Material + + +class Newtonian(_Material): + """Newtonian fluid material model.""" + + _props = ("density", "bulk_modulus", "dynamic_viscosity") + state_vars = ("pressure",) + + def __init__(self, material_properties: dict): + """Create a Newtonian material. + + Parameters + ---------- + material_properties: dict + Dictionary with material properties. For newtonian + materials, `density`, `bulk_modulus` and `dynamic_viscosity` + are required keys. + """ + self.validate_props(material_properties) + compressibility = 1 + + if material_properties.get("incompressible", False): + compressibility = 0 + + self.properties = { + **material_properties, + "compressibility": compressibility, + } + + def __repr__(self): + return f"Newtonian(props={self.properties})" + + def initialize_state_variables(self, nparticles: int) -> dict: + """Return initial state variables dictionary. + + Parameters + ---------- + nparticles : int + Number of particles being simulated with this material. + + Returns + ------- + dict + Dictionary of state variables initialized with values + decided by material type. + """ + state_vars_dict = {var: jnp.zeros((nparticles, 1)) for var in self.state_vars} + return state_vars_dict + + def _thermodynamic_pressure(self, volumetric_strain: ArrayLike) -> Array: + return -self.properties["bulk_modulus"] * volumetric_strain + + def compute_stress(self, particles): + """Compute material stress.""" + ndim = particles.loc.shape[-1] + if ndim not in {2, 3}: + raise ValueError(f"Cannot compute stress for {ndim}-d Newotonian material.") + volumetric_strain_rate = ( + particles.strain_rate[:, 0] + particles.strain_rate[:, 1] + ) + particles.state_vars["pressure"] = ( + particles.state_vars["pressure"] + .at[:] + .add( + self.properties["compressibility"] + * self._thermodynamic_pressure(particles.dvolumetric_strain) + ) + ) + + volumetric_stress_component = self.properties["compressibility"] * ( + -particles.state_vars["pressure"] + - (2 * self.properties["dynamic_viscosity"] * volumetric_strain_rate / 3) + ) + + stress = jnp.zeros_like(particles.stress) + stress = stress.at[:, 0].set( + volumetric_stress_component + + 2 * self.properties["dynamic_viscosity"] * particles.strain_rate[:, 0] + ) + stress = stress.at[:, 1].set( + volumetric_stress_component + + 2 * self.properties["dynamic_viscosity"] * particles.strain_rate[:, 1] + ) + + extra_component_2 = lax.select( + ndim == 3, + 2 * self.properties["dynamic_viscosity"] * particles.strain_rate[:, 2], + jnp.zeros_like(particles.strain_rate[:, 2]), + ) + stress = stress.at[:, 2].set(volumetric_stress_component + extra_component_2) + + stress = stress.at[:, 3].set( + self.properties["dynamic_viscosity"] * particles.strain_rate[:, 3] + ) + + component_4 = lax.select( + ndim == 3, + self.properties["dynamic_viscosity"] * particles.strain_rate[:, 4], + jnp.zeros_like(particles.strain_rate[:, 4]), + ) + stress = stress.at[:, 4].set(component_4) + component_5 = lax.select( + ndim == 3, + self.properties["dynamic_viscosity"] * particles.strain_rate[:, 5], + jnp.zeros_like(particles.strain_rate[:, 5]), + ) + stress = stress.at[:, 5].set(component_5) + + return stress diff --git a/tests/test_newtonian.py b/tests/test_newtonian.py new file mode 100644 index 0000000..518a246 --- /dev/null +++ b/tests/test_newtonian.py @@ -0,0 +1,90 @@ +import jax.numpy as jnp +import pytest +from diffmpm.constraint import Constraint +from diffmpm.element import Quadrilateral4Node +from diffmpm.materials import Newtonian +from diffmpm.node import Nodes +from diffmpm.particle import Particles + +particles_element_targets = [ + ( + Particles( + jnp.array([[0.5, 0.5]]).reshape(1, 1, 2), + Newtonian( + { + "density": 1000, + "bulk_modulus": 8333333.333333333, + "dynamic_viscosity": 8.9e-4, + } + ), + jnp.array([0]), + ), + Quadrilateral4Node( + (1, 1), + 1, + (4.0, 4.0), + [(0, Constraint(0, 0.02)), (0, Constraint(1, 0.03))], + Nodes(4, jnp.array([-2, -2, 2, -2, -2, 2, 2, 2]).reshape((4, 1, 2))), + ), + jnp.array( + [ + -52083.3333338896, + -52083.3333355583, + -52083.3333305521, + -0.0000041719, + 0, + 0, + ] + ).reshape(1, 6, 1), + ), + ( + Particles( + jnp.array([[0.5, 0.5]]).reshape(1, 1, 2), + Newtonian( + { + "density": 1000, + "bulk_modulus": 8333333.333333333, + "dynamic_viscosity": 8.9e-4, + "incompressible": True, + } + ), + jnp.array([0]), + ), + Quadrilateral4Node( + (1, 1), + 1, + (4.0, 4.0), + [(0, Constraint(0, 0.02)), (0, Constraint(1, 0.03))], + Nodes(4, jnp.array([-2, -2, 2, -2, -2, 2, 2, 2]).reshape((4, 1, 2))), + ), + jnp.array( + [ + -0.0000033375, + -0.00000500625, + 0, + -0.0000041719, + 0, + 0, + ] + ).reshape(1, 6, 1), + ), +] + + +@pytest.mark.parametrize( + "particles, element, target", + particles_element_targets, +) +def test_compute_stress(particles, element, target): + dt = 1 + particles.update_natural_coords(element) + if element.constraints: + element.apply_boundary_constraints() + particles.compute_strain(element, dt) + stress = particles.material.compute_stress(particles) + assert jnp.allclose(stress, target) + + +def test_init(): + with pytest.raises(KeyError): + Newtonian({"dynamic_viscosity": 1, "density": 1})