elastic_solid.h

Go to the documentation of this file.

/* -----------------------------------------------------------------------------*
 *                               SPHinXsys                                      *
 * -----------------------------------------------------------------------------*
 * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle    *
 * Hydrodynamics for industrial compleX systems. It provides C++ APIs for       *
 * physical accurate simulation and aims to model coupled industrial dynamic    *
 * systems including fluid, solid, multi-body dynamics and beyond with SPH      *
 * (smoothed particle hydrodynamics), a meshless computational method using     *
 * particle discretization.                                                     *
 *                                                                              *
 * SPHinXsys is partially funded by German Research Foundation                  *
 * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1,               *
 * HU1527/12-1 and HU1527/12-4.                                                 *
 *                                                                              *
 * Portions copyright (c) 2017-2022 Technical University of Munich and          *
 * the authors' affiliations.                                                   *
 *                                                                              *
 * Licensed under the Apache License, Version 2.0 (the "License"); you may      *
 * not use this file except in compliance with the License. You may obtain a    *
 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.           *
 *                                                                              *
 * -----------------------------------------------------------------------------*/
#ifndef ELASTIC_SOLID_H
#define ELASTIC_SOLID_H

#include "base_material.h"
#include <fstream>

namespace SPH
{
    class ElasticSolid : public Solid
    {
    protected:
        Real E0_;  /*< Youngs or tensile modules  */
        Real G0_;  /*< shear modules  */
        Real K0_;  /*< bulk modules  */
        Real nu_;  /*< Poisson ratio  */
        Real c0_;  /*< sound wave speed */
        Real ct0_; /*< tensile wave speed */
        Real cs0_; /*< shear wave speed */

        void setSoundSpeeds();

    public:
        explicit ElasticSolid(Real rho0)
            : Solid(rho0), c0_(0.0), ct0_(0.0), cs0_(0.0),
              E0_(0.0), G0_(0.0), K0_(0.0), nu_(0.0)
        {
            material_type_name_ = "ElasticSolid";
        };
        virtual ~ElasticSolid(){};

        Real ReferenceSoundSpeed() { return c0_; };
        Real TensileWaveSpeed() { return ct0_; };
        Real ShearWaveSpeed() { return cs0_; };
        Real YoungsModulus() { return E0_; };
        Real ShearModulus() { return G0_; };
        Real BulkModulus() { return K0_; };
        Real PoissonRatio() { return nu_; };

        virtual Matd StressPK2(Matd &deformation, size_t particle_index_i) = 0;
        virtual Matd StressCauchy(Matd &almansi_strain, Matd &F, size_t particle_index_i) = 0;
        virtual Matd NumericalDampingRightCauchy(Matd &deformation, Matd &deformation_rate, Real smoothing_length, size_t particle_index_i);
        virtual Matd NumericalDampingLeftCauchy(Matd &deformation, Matd &deformation_rate, Real smoothing_length, size_t particle_index_i);
        virtual Real PairNumericalDamping(Real dE_dt_ij, Real smoothing_length);

        virtual Matd DeviatoricKirchhoff(const Matd &deviatoric_be);
        virtual Real VolumetricKirchhoff(Real J) = 0;
        virtual std::string getRelevantStressMeasureName() = 0;

        virtual ElasticSolid *ThisObjectPtr() override { return this; };
    };

    class LinearElasticSolid : public ElasticSolid
    {
    public:
        explicit LinearElasticSolid(Real rho0, Real youngs_modulus, Real poisson_ratio);
        virtual ~LinearElasticSolid(){};

        virtual Matd StressPK2(Matd &deformation, size_t particle_index_i) override;
        virtual Matd StressCauchy(Matd &almansi_strain, Matd &F, size_t particle_index_i) override;
        virtual Real VolumetricKirchhoff(Real J) override;
        virtual std::string getRelevantStressMeasureName() override { return "PK2"; };

        Real getYoungsModulus() { return E0_; };
        Real getPoissonRatio() { return nu_; };
        Real getDensity() { return rho0_; };

    protected:
        Real lambda0_; /*< first Lame parameter */
        Real getBulkModulus(Real youngs_modulus, Real poisson_ratio);
        Real getShearModulus(Real youngs_modulus, Real poisson_ratio);
        Real getLambda(Real youngs_modulus, Real poisson_ratio);
    };

    class NeoHookeanSolid : public LinearElasticSolid
    {
    public:
        explicit NeoHookeanSolid(Real rho0, Real youngs_modulus, Real poisson_ratio)
            : LinearElasticSolid(rho0, youngs_modulus, poisson_ratio)
        {
            material_type_name_ = "NeoHookeanSolid";
        };
        virtual ~NeoHookeanSolid(){};

        virtual Matd StressPK2(Matd &deformation, size_t particle_index_i) override;
        virtual Matd StressCauchy(Matd &almansi_strain, Matd &F, size_t particle_index_i) override;
        virtual Real VolumetricKirchhoff(Real J) override;
        virtual std::string getRelevantStressMeasureName() override { return "Cauchy"; };
    };

    class NeoHookeanSolidIncompressible : public LinearElasticSolid
    {
    public:
        NeoHookeanSolidIncompressible(Real rho_0, Real Youngs_modulus, Real poisson)
            : LinearElasticSolid(rho_0, Youngs_modulus, poisson)
        {
            material_type_name_ = "NeoHookeanSolidIncompressible";
        };
        virtual ~NeoHookeanSolidIncompressible() {};
    
        virtual Matd StressPK2(Matd &deformation, size_t particle_index_i) override;
        virtual Matd StressCauchy(Matd &almansi_strain, Matd &F, size_t particle_index_i) override;
        virtual Real VolumetricKirchhoff(Real J) override;
    };

    class OrthotropicSolid : public LinearElasticSolid
    {
    public:
        OrthotropicSolid(Real rho_0, std::array<Vecd, 3> a, std::array<Real, 3> E, std::array<Real, 3> G, std::array<Real, 3> poisson);

        virtual Matd StressPK2(Matd& deformation, size_t particle_index_i) override;
        virtual Real VolumetricKirchhoff(Real J) override;

    protected:
        // input data
        std::array<Vecd, 3> a_;
        std::array<Real, 3> E_;
        std::array<Real, 3> G_;
        std::array<Real, 3> poisson_;
        // calculated data
        Real Mu_[3];
        Matd Lambda_;
        Matd A_[3];

        virtual void CalculateAllMu();
        virtual void CalculateAllLambda();
        virtual void CalculateA0();
    };

    class FeneNeoHookeanSolid : public LinearElasticSolid
    {
    protected:
        Real j1_m_; 
    public:
        explicit FeneNeoHookeanSolid(Real rho0, Real youngs_modulus, Real poisson_ratio)
            : LinearElasticSolid(rho0, youngs_modulus, poisson_ratio), j1_m_(1.0)
        {
            material_type_name_ = "FeneNeoHookeanSolid";
        };
        virtual ~FeneNeoHookeanSolid() {};
        virtual Matd StressPK2(Matd& deformation, size_t particle_index_i) override;
        virtual std::string getRelevantStressMeasureName() override { return "Cauchy"; };
    };

    class Muscle : public NeoHookeanSolid
    {
    public:
        explicit Muscle(Real rho0, Real bulk_modulus,
                        const Vecd &f0, const Vecd &s0, const Real (&a0)[4], const Real (&b0)[4])
            : NeoHookeanSolid(rho0, this->getYoungsModulus(bulk_modulus, a0, b0), this->getPoissonRatio(bulk_modulus, a0, b0)),
              f0_(f0), s0_(s0), f0f0_(SimTK::outer(f0_, f0_)), s0s0_(SimTK::outer(s0_, s0_)),
              f0s0_(SimTK::outer(f0_, s0_))
        {
            material_type_name_ = "Muscle";
            std::copy(a0, a0 + 4, a0_);
            std::copy(b0, b0 + 4, b0_);
        };
        virtual ~Muscle(){};

        virtual Matd MuscleFiberDirection(size_t particle_index_i) { return f0f0_; };
        virtual Matd StressPK2(Matd &deformation, size_t particle_index_i) override;
        virtual Real VolumetricKirchhoff(Real J) override;
        virtual std::string getRelevantStressMeasureName() override { return "Cauchy"; };

        virtual Muscle *ThisObjectPtr() override { return this; };

    protected:
        Vecd f0_, s0_;            
        Matd f0f0_, s0s0_, f0s0_; 
        Real a0_[4], b0_[4];      
    private:
        Real getPoissonRatio(Real bulk_modulus, const Real (&a0)[4], const Real (&b0)[4]);
        Real getShearModulus(const Real (&a0)[4], const Real (&b0)[4]);
        Real getYoungsModulus(Real bulk_modulus, const Real (&a0)[4], const Real (&b0)[4]);
    };

    class LocallyOrthotropicMuscle : public Muscle
    {
    protected:
        StdLargeVec<Matd> local_f0f0_, local_s0s0_, local_f0s0_; 
        virtual void initializeFiberAndSheet();
        void initializeFiberAndSheetTensors();

    public:
        StdLargeVec<Vecd> local_f0_; 
        StdLargeVec<Vecd> local_s0_; 
        explicit LocallyOrthotropicMuscle(Real rho0, Real bulk_modulus,
                        const Vecd &f0, const Vecd &s0, const Real (&a0)[4], const Real (&b0)[4])
            : Muscle(rho0, bulk_modulus, f0, s0, a0, b0)
        {
            material_type_name_ = "LocallyOrthotropicMuscle";
        };
        virtual ~LocallyOrthotropicMuscle(){};

        virtual void assignBaseParticles(BaseParticles *base_particles) override;
        virtual Matd MuscleFiberDirection(size_t particle_index_i) override { return local_f0f0_[particle_index_i]; };
        virtual Matd StressPK2(Matd& deformation, size_t particle_index_i) override;
        virtual void readFromXmlForLocalParameters(const std::string &filefullpath) override;
        virtual std::string getRelevantStressMeasureName() override { return "Cauchy"; };
    };
}
#endif //ELASTIC_SOLID_H