From bd68a725b8b28b2cf70a70fcb8df155ce80dc2f2 Mon Sep 17 00:00:00 2001 From: frohne Date: Mon, 18 Feb 2013 21:16:37 +0000 Subject: [PATCH] correction of description of the linear hardening implementation and solving of an elastic problem for a better start solution git-svn-id: https://svn.dealii.org/trunk@28456 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/step-42.cc | 51 ++++++++++++++++++++--------- 1 file changed, 36 insertions(+), 15 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index 57edd6d895..530f2a57f3 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -350,6 +350,7 @@ namespace Step42 inline SymmetricTensor<2,dim> get_strain (const FEValues &fe_values, const unsigned int shape_func, const unsigned int q_point) const; + void set_sigma_0 (double sigma_hlp) {sigma_0 = sigma_hlp;} private: SymmetricTensor<4,dim> stress_strain_tensor_mu; @@ -404,7 +405,6 @@ namespace Step42 { SymmetricTensor<2,dim> stress_tensor; stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)*strain_tensor; - double tmp = E/((1+nu)*(1-2*nu)); SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor); @@ -413,10 +413,12 @@ namespace Step42 yield = 0; stress_strain_tensor = stress_strain_tensor_mu; double beta = 1.0; - if (deviator_stress_tensor_norm >= sigma_0) + if (deviator_stress_tensor_norm > + sigma_0) { - beta = (sigma_0 + gamma)/deviator_stress_tensor_norm; - stress_strain_tensor *= beta; +// beta = (sigma_0 + gamma/(2*mu + gamma)*(deviator_stress_tensor_norm - sigma_0))/deviator_stress_tensor_norm; + beta = sigma_0/deviator_stress_tensor_norm; + stress_strain_tensor *= (gamma + (1 - gamma)*beta); yield = 1; plast_points += 1; } @@ -436,23 +438,27 @@ namespace Step42 { SymmetricTensor<2,dim> stress_tensor; stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)*strain_tensor; - double tmp = E/((1+nu)*(1-2*nu)); - - stress_strain_tensor = stress_strain_tensor_mu; - stress_strain_tensor_linearized = stress_strain_tensor_mu; SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor); double deviator_stress_tensor_norm = deviator_stress_tensor.norm (); + stress_strain_tensor = stress_strain_tensor_mu; + stress_strain_tensor_linearized = stress_strain_tensor_mu; double beta = 1.0; - if (deviator_stress_tensor_norm >= sigma_0) + if (deviator_stress_tensor_norm > + sigma_0) { - beta = (sigma_0 + gamma)/deviator_stress_tensor_norm; - stress_strain_tensor *= beta; - stress_strain_tensor_linearized *= beta; +// beta = (sigma_0 + gamma/(2*mu + gamma)*(deviator_stress_tensor_norm - sigma_0))/deviator_stress_tensor_norm; + beta = sigma_0/deviator_stress_tensor_norm; + stress_strain_tensor *= (gamma + (1 - gamma)*beta); + stress_strain_tensor_linearized *= (gamma + (1 - gamma)*beta); deviator_stress_tensor /= deviator_stress_tensor_norm; - stress_strain_tensor_linearized -= beta*2*mu*outer_product(deviator_stress_tensor, deviator_stress_tensor); + stress_strain_tensor_linearized -= (1 - gamma)*beta*2*mu*outer_product(deviator_stress_tensor, deviator_stress_tensor); + +// std::cout << "Plastisch mit beta = " << beta +// << ", sigma_0 = " << sigma_0 +// << std::endl; } stress_strain_tensor += stress_strain_tensor_kappa; @@ -605,8 +611,8 @@ namespace Step42 pcout (std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), sigma_0 (400), - gamma (1.e-2), - e_modul (2.0e5), + gamma (0.01), + e_modul (2.0e+5), nu (0.3) { // double _E, double _nu, double _sigma_0, double _gamma @@ -1061,6 +1067,12 @@ namespace Step42 double solution_index_z = solution (index_z); double gap = obstacle_value - point (2); + +// std::cout << "lambda = " << lambda (index_z) +// << ", solution_index_z - gap = " << solution_index_z - gap +// << ", diag_mass_matrix_vector_relevant = " << diag_mass_matrix_vector_relevant (index_z) +// << std::endl; + if (lambda (index_z) + c * diag_mass_matrix_vector_relevant (index_z) * @@ -1235,6 +1247,8 @@ namespace Step42 ComponentMask(), constant_modes); + double sigma_hlp = sigma_0; + additional_data.elliptic = true; additional_data.n_cycles = 1; additional_data.w_cycle = false; @@ -1247,6 +1261,12 @@ namespace Step42 unsigned int number_assemble_system = 0; for (; j<=100; j++) { + // Solve an elastic problem to obtain a better start solution + if (j == 0) + plast_lin_hard->set_sigma_0 (1e+10); + else if (j == 1) + plast_lin_hard->set_sigma_0 (sigma_hlp); + pcout<< " " <