Bingham plastic dynamics

Palabos General
#1

I have implemented a very simple Bingham plastic dynamics.

I use it as:


Bingham::omegaMin =  1. / 5.;
Bingham::plasticViscosity = 0.02;
Bingham::yieldStress = 0.02;
new Bingham::BinghamDynamics<T,D>(new BGKdynamics<T,D>(1.0)),

Here it is:


#ifndef BINGHAM_DYNAMICS_H
#define BINGHAM_DYNAMICS_H

namespace plb {
namespace Bingham {
    
double yieldStress;
double plasticViscosity;
double omegaMin;

template<typename T, template<typename U> class D>
class BinghamDynamics : public OmegaFromPiDynamics<T,D> {
public:
    BinghamDynamics(Dynamics<T,D>* baseDynamics, bool automaticPrepareCollision = true)
    : OmegaFromPiDynamics<T,D>(baseDynamics, automaticPrepareCollision){}
    
    BinghamDynamics<T,D>* clone() const { return new BinghamDynamics<T,D>(*this); }
    
    virtual int getId() const { return id; }
    
    virtual T getOmegaFromPiAndRhoBar(Array<T,SymmetricTensor<T,D>::n> const& PiNeq, T rhoBar) const {
        PLB_PRECONDITION(yieldStress > -1e-15);
        PLB_PRECONDITION(plasticViscosity > 0);
        PLB_PRECONDITION(omegaMin > 0);

        T ic2 = D<T>::invCs2;
        T omega;
        if (yieldStress < (T)1e-18) {
            omega = (T)1 / (plasticViscosity * ic2 + (T)0.5);
        } else {
            if (PiNeq[0] == 0) {
                omega = 1;
            } else {
                T piNeqNorm = sqrt(SymmetricTensor<T,D>::tensorNormSqr(PiNeq));
                T rho = D<T>::fullRho(rhoBar);
                T om = ((T)1 - sqrt((T)2) * yieldStress * rho / piNeqNorm) / (plasticViscosity * ic2 + (T)0.5);
                if (om != om || om < omegaMin) {
                    omega = omegaMin;
                } else {
                    omega = om;
                }
            }
        }
        return omega;
    }
private:
    static int id;
};

template<typename T, template<typename U> class D>    
int BinghamDynamics<T,D>::id = meta::registerCompositeDynamics<T,D, BinghamDynamics<T,D> >("BinghamDynamics");

}
}
#endif

EDIT:
Here I provide derivation of the equation, so that I can remember that as well :slight_smile: I followed thesis of Orestis Malaspinas chapter 5.1. The only exception is that I used sqrt(2) in Eq. 5.5 instead of 2.


veff = Sy / gdot + vpl
veff = c2 * (tau - 0.5)
gdot = 1 / sqrt(2) / c2 / rho / tau * sqrt(PP)
c2 * (tau - 0.5) = Sy / gdot + vpl
c2 * (tau - 0.5) = Sy / (sqrt(PP) / sqrt(2) / c2 / rho / tau) + vpl
tau = (Sy / (sqrt(PP) / sqrt(2) / c2 / rho / tau) + vpl) / c2 + 0.5
tau = Sy * sqrt(2) * rho * tau / sqrt(PP) + vpl / c2 + 0.5
tau * (1 - Sy * sqrt(2) * rho / sqrt(PP)) = vpl / c2 + 0.5
tau = (vpl / c2 + 0.5) / (1 - Sy * sqrt(2) * rho / sqrt(PP))
omega = (1 - Sy * sqrt(2) * rho / sqrt(PP)) / (vpl / c2 + 0.5)