Dear all,
I'm trying to extend step-42, according to what the authors Frohne,
Heister, and Bangerth suggested in the related paper, ie: modify the
current one-step analysis to a sequence of loading steps (quasi-static
analysis).
Apart from a for-loop over the pseudo-time parameter,
i) I recast the procedures computing the constitutive laws in incremental
form (see the attached files "model.cpp" and "model.hpp"), by introducing
also the eta parameter which accounts for the load history as accumulated
plastic strain;
ii) I coherently modified the sphere obstacle function by increasing for
each time step the position of the sphere center;
ii) I added in the residuals the contribution of the stress stored at the
previous time-step (local_q_points_history[q_point].old_stress) as
*cell_rhs(i) -= (- local_q_points_history[q_point].old_stress +
strain_tensors[q_point] * stress_strain_tensor ) *
fe_values[displacement].symmetric_gradient(i, q_point) *
fe_values.JxW(q_point);*
The variables eta and old_stress are both stored as local_q_points_history.
I tried with different mesh refinements and different time-step sizes, but
the analyses provide unsatisfactory results. By summarizing, the size of
active set reduces to zero after the first time step, and then increases.
The convergence of the Newton solver slows down since the initial loops,
leading to loose convergence (for max loop reached). The evaluation of the
resulting contact force fails giving negative values.
At the moment, I'm really stuck. Is there anybody who can give me some
hints?
Thanks in advance,
Giovanni
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/d65d2b44-10ed-4d2f-ab73-3dbff7932680n%40googlegroups.com.
//
// model.cpp
//
#include "model.hpp"
////////////////////////////////////////////////////////////////////////////////
ConstitutiveLaw::ConstitutiveLaw(const double E,
const double nu,
const double sigma_0,
const double gamma,
const double eta_0)
: kappa(E / (3 * (1 - 2 * nu)))
, mu(E / (2 * (1 + nu)))
, sigma_0(sigma_0)
, gamma(gamma)
, eta(eta_0)
, stress_strain_tensor_kappa(kappa *
outer_product(unit_symmetric_tensor<dim>(),
unit_symmetric_tensor<dim>()))
, stress_strain_tensor_mu(2 * mu *
(identity_tensor<dim>() - outer_product(unit_symmetric_tensor<dim>(),
unit_symmetric_tensor<dim>()) / 3.0))
{}
////////////////////////////////////////////////////////////////////////////////
void ConstitutiveLaw::set_sigma_0(double sigma_zero)
{
sigma_0 = sigma_zero;
}
void ConstitutiveLaw::set_eta(double new_eta)
{
eta = new_eta;
}
////////////////////////////////////////////////////////////////////////////////
bool ConstitutiveLaw::get_stress_strain_tensor(const SymmetricTensor<2, dim> &stress_tensor,
const SymmetricTensor<2, dim> &strain_tensor,
SymmetricTensor<4, dim> &stress_strain_tensor) const
{
Assert(dim == 3, ExcNotImplemented());
SymmetricTensor<2, dim> current_stress_tensor;
current_stress_tensor = stress_tensor + (stress_strain_tensor_kappa + stress_strain_tensor_mu) * strain_tensor;
const SymmetricTensor<2, dim> deviator_stress_tensor = deviator(current_stress_tensor);
const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
const double gamma_iso = 2.0 * mu * gamma / (1.0 - gamma);
const double gamma_eta = 1.0 - 2.0 * mu * eta / sigma_0;
stress_strain_tensor = stress_strain_tensor_mu;
if (deviator_stress_tensor_norm > sigma_0 + gamma_iso * eta)
{
const double beta = sigma_0 / deviator_stress_tensor_norm;
stress_strain_tensor *= (gamma + (1 - gamma * gamma_eta) * beta);
}
stress_strain_tensor += stress_strain_tensor_kappa;
return (deviator_stress_tensor_norm > sigma_0);
}
////////////////////////////////////////////////////////////////////////////////
void ConstitutiveLaw::get_linearized_stress_strain_tensors(const SymmetricTensor<2, dim> &stress_tensor,
const SymmetricTensor<2, dim> &strain_tensor,
SymmetricTensor<4, dim> &stress_strain_tensor_linearized,
SymmetricTensor<4, dim> &stress_strain_tensor) const
{
Assert(dim == 3, ExcNotImplemented());
SymmetricTensor<2, dim> current_stress_tensor;
current_stress_tensor = stress_tensor + (stress_strain_tensor_kappa + stress_strain_tensor_mu) * strain_tensor;
stress_strain_tensor = stress_strain_tensor_mu;
stress_strain_tensor_linearized = stress_strain_tensor_mu;
SymmetricTensor<2, dim> deviator_stress_tensor = deviator(current_stress_tensor);
const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
const double gamma_iso = 2.0 * mu * gamma / (1.0 - gamma);
const double gamma_eta = 1.0 - 2.0 * mu * eta / sigma_0;
if (deviator_stress_tensor_norm > sigma_0 + gamma_iso * eta)
{
const double beta = sigma_0 / deviator_stress_tensor_norm;
stress_strain_tensor *= (gamma + (1 - gamma * gamma_eta) * beta);
stress_strain_tensor_linearized *= (gamma + (1 - gamma * gamma_eta) * beta);
deviator_stress_tensor /= deviator_stress_tensor_norm;
stress_strain_tensor_linearized -= (1 - gamma) * beta * 2 * mu * outer_product(deviator_stress_tensor, deviator_stress_tensor);
}
stress_strain_tensor += stress_strain_tensor_kappa;
stress_strain_tensor_linearized += stress_strain_tensor_kappa;
}
////////////////////////////////////////////////////////////////////////////////
bool ConstitutiveLaw::get_current_stress_tensor(const SymmetricTensor<2, dim> &stress_tensor,
const SymmetricTensor<2, dim> &incremental_strain_tensor,
SymmetricTensor<2, dim> ¤t_stress_tensor) const
{
Assert(dim == 3, ExcNotImplemented());
// predict incremental elastic stress
current_stress_tensor = stress_tensor;
current_stress_tensor += (stress_strain_tensor_kappa + stress_strain_tensor_mu) * incremental_strain_tensor;
// get current stress-strain tensor
const SymmetricTensor<2, dim> deviator_stress_tensor = deviator(current_stress_tensor);
const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
const double gamma_iso = 2.0 * mu * gamma / (1.0 - gamma);
const double gamma_eta = 1.0 - 2.0 * mu * eta / sigma_0;
SymmetricTensor<4, dim> stress_strain_tensor = stress_strain_tensor_mu;
if (deviator_stress_tensor_norm > sigma_0 + gamma_iso * eta)
{
const double beta = sigma_0 / deviator_stress_tensor_norm;
stress_strain_tensor *= (gamma + (1 - gamma * gamma_eta) * beta);
}
stress_strain_tensor += stress_strain_tensor_kappa;
// get current stress tensor
current_stress_tensor = stress_tensor;
current_stress_tensor += stress_strain_tensor * incremental_strain_tensor;
// return if in plasticity
return (deviator_stress_tensor_norm > sigma_0);
}
////////////////////////////////////////////////////////////////////////////////
double ConstitutiveLaw::get_current_eta(const SymmetricTensor<2, dim> &stress_tensor) const
{
// get current stress-strain tensor
const SymmetricTensor<2, dim> deviator_stress_tensor = deviator(stress_tensor);
const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
const double gamma_iso = 2.0 * mu * gamma / (1.0 - gamma);
return (deviator_stress_tensor_norm - sigma_0 + 2.0 * mu * eta) / ( 2.0 * mu + gamma_iso);
}
//
// model.hpp
//
#ifndef model_hpp
#define model_hpp
#include "general.hpp"
/// physical model
class ConstitutiveLaw
{
public:
/// constructor
ConstitutiveLaw(const double E,
const double nu,
const double sigma_0,
const double gamma,
const double eta_0);
/// set yield stress value
void set_sigma_0(double sigma_zero);
/// set the internal hardening variable
void set_eta(double new_eta);
/// get stress strain tensor
bool get_stress_strain_tensor(const SymmetricTensor<2, dim> &stress_tensor,
const SymmetricTensor<2, dim> &strain_tensor,
SymmetricTensor<4, dim> & stress_strain_tensor) const;
/// get linearized stress strain tensor
void get_linearized_stress_strain_tensors(const SymmetricTensor<2, dim> &stress_tensor,
const SymmetricTensor<2, dim> &strain_tensor,
SymmetricTensor<4, dim> & stress_strain_tensor_linearized,
SymmetricTensor<4, dim> & stress_strain_tensor) const;
/// get current strain tensor
bool get_current_stress_tensor(const SymmetricTensor<2, dim> &stress_tensor,
const SymmetricTensor<2, dim> &incremental_strain_tensor,
SymmetricTensor<2, dim> ¤t_stress_tensor) const;
/// get current internal hardening variable
double get_current_eta(const SymmetricTensor<2, dim> &stress_tensor) const;
private:
const double kappa;
const double mu;
double sigma_0;
const double gamma;
double eta;
public:
const SymmetricTensor<4, dim> stress_strain_tensor_kappa;
const SymmetricTensor<4, dim> stress_strain_tensor_mu;
};
#endif /* model_hpp */
