Re: [deal.II] how to extend step-42 to quasi-static problem

[gif...@gmail.com]

Dear Wolfgang

thank you for your quick reply. All your suggestions are very useful.
I have already done some debugging you proposed, but I will go back on my 
steps, and try again to do it more carefully.

Best

Giovanni

On Tuesday, April 5, 2022 at 4:49:30 PM UTC+2 Wolfgang Bangerth wrote:

> On 4/5/22 07:21, gif...@gmail.com wrote:
> > 
> > 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?
>
> Giovanni,
> short of finding someone who is (i) an expert in time dependent 
> plasticity, 
> (ii) has plenty of spare time to look through your program, you are left 
> with 
> debugging the problem yourself. Here is how I would approach things:
>
> * Simplify. Instead of looking at a whole load history, just consider two 
> steps.
> * Pick a load history so that what step-42 produces corresponds to your 
> first 
> step. Make sure you get the same result. If it doesn't, you've just 
> figured 
> out that already your first load step has a problem.
> * Then double the load in the second step. You know that that should 
> result in 
> a larger indentation, and a larger active set.
> * If it doesn't, investigate why that is. One approach would be to 
> simplify 
> the set up so that instead of indenting a complex shape (or just a 
> sphere), 
> you impose a constant load across the top surface of the domain. You might 
> even be able to compute the solution for this problem by hand because of 
> the 
> symmetries of the deformation you expect.
>
>
> It is not uncommon for a complex program to be wrong in the beginning -- 
> in 
> fact, I'd say that's how nearly every program starts out. The question is 
> whether you can build the mental tools to break things down into smaller 
> pieces that can be debugged more easily. Simplifying the situation to 
> something for which it is easier to reason about the behavior is a key 
> first 
> step in this process.
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>

-- 
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 dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/5d3f2d6a-0c22-48b7-af1f-2bf9cc0f4998n%40googlegroups.com.

[deal.II] how to extend step-42 to quasi-static problem

[gif...@gmail.com]

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 dealii+unsubscr...@googlegroups.com.
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(),
   unit_symmetric_tensor()))
, stress_strain_tensor_mu(2 * mu *
  (identity_tensor() - outer_product(unit_symmetric_tensor(),
  unit_symmetric_tensor()) / 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> _tensor,
   const SymmetricTensor<2, dim> _tensor,
 SymmetricTensor<4, dim> _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> _tensor,
   const SymmetricTensor<2, dim> _tensor,
 SymmetricTensor<4, dim> _strain_tensor_linearized,