correction of description of the linear hardening implementation and solving of an...

git-svn-id: https://svn.dealii.org/trunk@28456 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc
index 57edd6d89518ae99ab7eff35fe7fbf00dee61b49..530f2a57f3e710cb252ee436944d3f12b6886eef 100644 (file)
--- 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<dim> &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<< " " <<std::endl;
         pcout<< "   Newton iteration " << j <<std::endl;
         pcout<< "      Updating active set..." <<std::endl;
@@ -1288,6 +1308,7 @@ namespace Step42
 
         int damped = 0;
         tmp_vector = old_solution;
+//        constraints.distribute (old_solution);
         double a = 0;
         for (unsigned int i=0; (i<10)&&(!damped); i++)
           {