Update NorSand for Jefferies and Been (2019)

[jgiven100]

Update NorSand for Jefferies and Been (2019)

Summary

Current NorSand formulation is based primarily on Jefferies (1993) with some additions from Jefferies and Shuttle (2012). The 2nd edition of Soil Liquefaction (Jefferies and Been, 2019) provides updates for NorSand. These updates include the introduction of soil state-dilatancy property \chi, simplified critical friction ratio M, and changes to volumetric coupling property N.

Additional internal parameters include \psi_i, \chi_i, and M_i or the state paramter image, state-dilatancy image, and critical state ratio image. The resulting new yield surface is "bullet" shaped, similar to original Cam Clay.

Motivation

NorSand was initially published with Nova's rule:

D^p = ( M_{tc} - \eta ) / ( 1 - N )

but Dafalias has shown that M_{tc} is strain dependent. Thus stress-dilatancy should be updated to:

D^p = M_i - \eta 

NorSand should be updated prior to exploring anisotropic stresses within the model as briefly discussed here -> https://github.com/cb-geo/mpm/discussions/713

Design Detail

Generally, the existing NorSand material class will be updated based on the improved yield surface and hardening rule. This will require:

New Private Variable

chi_image_ will be added to the material construtor (beneath Mtc_). This will be a constant since chi_, lambda_, and Mtc_ are all currently implemented as constants.

chi_image_ = chi_ / (1. - ((chi_ * lambda_) / Mtc_))

New State Parameters

Both M_imgae and psi_image will be added to the state_vars dense_map. They will be initialized as

// M_image
{"M_image", 0.},
 // State parameter image
{"psi_image", void_ratio_initial_ - (gamma_ - lambda_ * log(p_image_initial_ / reference_pressure_))}

New Functions

M_image and psi_image will be updated via the following new function.

//! Compute image parameters
template <unsigned Tdim>
void mpm::NorSand<Tdim>::compute_image_parameters(mpm::dense_map* state_vars,
                                                  const double void_ratio,
                                                  const double e_image,
                                                  const double M_theta) {
  // Compute state parameter image
  const double psi_image_ = void_ratio - e_image;
  (*state_vars).at("psi_image") = psi_image_;

  // Compute critical state coefficient image
  (*state_vars).at("M_image") =
      M_theta * (1. - ((chi_image_ * N_ * std::fabs(psi_image_)) / Mtc_));
}

Updated Functions

• compute_yield_state(...) will most importantly have updates to the yield function,

  // Compute yield functions
  (*yield_function) = deviatoric_q / (mean_p + p_cohesion) - M_image +
                      M_image * std::log((mean_p + p_cohesion) /
                                         (p_image + p_cohesion + p_dilation));

• compute_state_variables(...) will have updates to how p_image is computed, and

    // Compute and update pressure image
    p_image =
        std::pow(
            std::exp(1 - (deviatoric_q / ((mean_p + p_cohesion) * M_image))),
            -1) *
            (mean_p + p_cohesion) -
        (p_cohesion + p_dilation);
    (*state_vars).at("p_image") = p_image;

• compute_plastic_tensor(...) will have the bulk of updates with the derivatives and plastic tensor construction changing.

Drawbacks

• These improvements should have limited impact on users. Additional parameters are strictly for internal use (i.e., no additional input paramters).
• These imporvements should have limited impact on other processes since the global NorSand behavior will be the same.

Rationale and Alternatives

• Do nothing and continue to use current NorSand implementation.

Prior Art

"Original" NorSand is already implemented within CB-Geo. Thus, most of the hard work has been completed. These changes simply look to update NorSand to the latest recommendation of Jefferies and Been (2019).

Unresolved questions

First

Previous implementation has assumed:

e_image = (gamma_ - lambda_ * log(p_image_initial_ / reference_pressure_)

I'm not sure where this is coming from and if it is still an acceptable way to find the void ratio image.

Second

I am still working through how the stiffness tenosr will be updated due to the new yield surface and hardening relationship. This is the remaining key element to be completed.

Third

What is a good way to showcase differences between existing and updated model prior to merging?

Changelog

Edit 1: add Table 3.1 from Jefferies and Been (2019)

[jgiven100]

One thing I've come across while working on the updates is that function for dp_dsigma in the material utilities (code here) requires passing a 6x1 stress tensor:

//! Compute derivative of p in terms of stress sigma
inline const Eigen::Matrix<double, 6, 1> mpm::materials::dp_dsigma(
    const Eigen::Matrix<double, 6, 1>& stress) {

  Eigen::Matrix<double, 6, 1> dp_dsigma = Eigen::Matrix<double, 6, 1>::Zero();
  dp_dsigma(0) = 1. / 3.;
  dp_dsigma(1) = 1. / 3.;
  dp_dsigma(2) = 1. / 3.;

  return dp_dsigma;
}

Since this function returns a constant tensor, this can be simplified to

//! Compute derivative of p in terms of stress sigma
inline const Eigen::Matrix<double, 6, 1> mpm::materials::dp_dsigma() {

  Eigen::Matrix<double, 6, 1> dp_dsigma = Eigen::Matrix<double, 6, 1>::Zero();
  dp_dsigma(0) = 1. / 3.;
  dp_dsigma(1) = 1. / 3.;
  dp_dsigma(2) = 1. / 3.;

  return dp_dsigma;
}

Are there any reasons not to make this simplification?

[kks32]

@jgiven100 Considering we have merged #717 can this be closed?

[jgiven100]

@kks32 Yes.