Albany: a Trilinos-based PDE code: J2FiberModel

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include <Intrepid_MiniTensor.h>
00008 #include "Teuchos_TestForException.hpp"
00009 #include "Phalanx_DataLayout.hpp"
00010 
00011 namespace LCM
00012 {
00013 
00014 //------------------------------------------------------------------------------
00015 template<typename EvalT, typename Traits>
00016 J2FiberModel<EvalT, Traits>::
00017 J2FiberModel(Teuchos::ParameterList* p,
00018     const Teuchos::RCP<Albany::Layouts>& dl) :
00019         LCM::ConstitutiveModel<EvalT, Traits>(p, dl),
00020         k_f1_(p->get<RealType>("Fiber 1 k", 1.0)),
00021         q_f1_(p->get<RealType>("Fiber 1 q", 1.0)),
00022         volume_fraction_f1_(p->get<RealType>("Fiber 1 volume fraction", 0.0)),
00023         max_damage_f1_(p->get<RealType>("Fiber 1 maximum damage", 1.0)),
00024         saturation_f1_(p->get<RealType>("Fiber 1 damage saturation", 0.0)),
00025         k_f2_(p->get<RealType>("Fiber 2 k", 1.0)),
00026         q_f2_(p->get<RealType>("Fiber 2 q", 1.0)),
00027         volume_fraction_f2_(p->get<RealType>("Fiber 2 volume fraction", 0.0)),
00028         max_damage_f2_(p->get<RealType>("Fiber 2 maximum damage", 1.0)),
00029         saturation_f2_(p->get<RealType>("Fiber 2 damage saturation", 0.0)),
00030         volume_fraction_m_(p->get<RealType>("Matrix volume fraction", 1.0)),
00031         max_damage_m_(p->get<RealType>("Matrix maximum damage", 1.0)),
00032         saturation_m_(p->get<RealType>("Matrix damage saturation", 0.0)),
00033         sat_mod_(p->get<RealType>("Saturation Modulus", 0.0)),
00034         sat_exp_(p->get<RealType>("Saturation Exponent", 1.0)),
00035         local_coord_flag_(p->get<bool>("Local Coordinate Flag", false)),
00036         direction_f1_(
00037             p->get<Teuchos::Array<RealType> >("Fiber 1 Orientation Vector")
00038                 .toVector()),
00039         direction_f2_(
00040             p->get<Teuchos::Array<RealType> >("Fiber 2 Orientation Vector")
00041                 .toVector()),
00042         ring_center_(
00043             p->get<Teuchos::Array<RealType> >("Ring Center Vector").toVector())
00044 {
00045   // define the dependent fields
00046   this->dep_field_map_.insert(std::make_pair("F", dl->qp_tensor));
00047   this->dep_field_map_.insert(std::make_pair("J", dl->qp_scalar));
00048   this->dep_field_map_.insert(std::make_pair("Poissons Ratio", dl->qp_scalar));
00049   this->dep_field_map_.insert(std::make_pair("Elastic Modulus", dl->qp_scalar));
00050   this->dep_field_map_.insert(std::make_pair("Yield Strength", dl->qp_scalar));
00051   this->dep_field_map_.insert(
00052       std::make_pair("Hardening Modulus", dl->qp_scalar));
00053 
00054   if (local_coord_flag_) {
00055     need_integration_pt_locations_ = true;
00056   }
00057 
00058   // retrieve appropriate field name strings
00059   std::string cauchy_string = (*field_name_map_)["Cauchy_Stress"];
00060   std::string Fp_string = (*field_name_map_)["Fp"];
00061   std::string eqps_string = (*field_name_map_)["eqps"];
00062   std::string matrix_energy_string = (*field_name_map_)["Matrix_Energy"];
00063   std::string f1_energy_string = (*field_name_map_)["F1_Energy"];
00064   std::string f2_energy_string = (*field_name_map_)["F2_Energy"];
00065   std::string matrix_damage_string = (*field_name_map_)["Matrix_Damage"];
00066   std::string f1_damage_string = (*field_name_map_)["F1_Damage"];
00067   std::string f2_damage_string = (*field_name_map_)["F2_Damage"];
00068 
00069   // define the evaluated fields
00070   this->eval_field_map_.insert(std::make_pair(cauchy_string, dl->qp_tensor));
00071   this->eval_field_map_.insert(
00072       std::make_pair(matrix_energy_string, dl->qp_scalar));
00073   this->eval_field_map_.insert(std::make_pair(f1_energy_string, dl->qp_scalar));
00074   this->eval_field_map_.insert(std::make_pair(f2_energy_string, dl->qp_scalar));
00075   this->eval_field_map_.insert(
00076       std::make_pair(matrix_damage_string, dl->qp_scalar));
00077   this->eval_field_map_.insert(std::make_pair(f1_damage_string, dl->qp_scalar));
00078   this->eval_field_map_.insert(std::make_pair(f2_damage_string, dl->qp_scalar));
00079   this->eval_field_map_.insert(std::make_pair(Fp_string, dl->qp_tensor));
00080   this->eval_field_map_.insert(std::make_pair(eqps_string, dl->qp_scalar));
00081 
00082   // define the state variables
00083   //
00084   // stress
00085   this->num_state_variables_++;
00086   this->state_var_names_.push_back(cauchy_string);
00087   this->state_var_layouts_.push_back(dl->qp_tensor);
00088   this->state_var_init_types_.push_back("scalar");
00089   this->state_var_init_values_.push_back(0.0);
00090   this->state_var_old_state_flags_.push_back(false);
00091   this->state_var_output_flags_.push_back(true);
00092   //
00093   // Fp
00094   this->num_state_variables_++;
00095   this->state_var_names_.push_back(Fp_string);
00096   this->state_var_layouts_.push_back(dl->qp_tensor);
00097   this->state_var_init_types_.push_back("identity");
00098   this->state_var_init_values_.push_back(1.0);
00099   this->state_var_old_state_flags_.push_back(true);
00100   this->state_var_output_flags_.push_back(false);
00101   //
00102   // eqps
00103   this->num_state_variables_++;
00104   this->state_var_names_.push_back(eqps_string);
00105   this->state_var_layouts_.push_back(dl->qp_scalar);
00106   this->state_var_init_types_.push_back("scalar");
00107   this->state_var_init_values_.push_back(0.0);
00108   this->state_var_old_state_flags_.push_back(true);
00109   this->state_var_output_flags_.push_back(true);
00110   //
00111   // matrix energy
00112   this->num_state_variables_++;
00113   this->state_var_names_.push_back(matrix_energy_string);
00114   this->state_var_layouts_.push_back(dl->qp_scalar);
00115   this->state_var_init_types_.push_back("scalar");
00116   this->state_var_init_values_.push_back(0.0);
00117   this->state_var_old_state_flags_.push_back(true);
00118   this->state_var_output_flags_.push_back(true);
00119   //
00120   // fiber 1 energy
00121   this->num_state_variables_++;
00122   this->state_var_names_.push_back(f1_energy_string);
00123   this->state_var_layouts_.push_back(dl->qp_scalar);
00124   this->state_var_init_types_.push_back("scalar");
00125   this->state_var_init_values_.push_back(0.0);
00126   this->state_var_old_state_flags_.push_back(true);
00127   this->state_var_output_flags_.push_back(true);
00128   //
00129   // fiber 2 energy
00130   this->num_state_variables_++;
00131   this->state_var_names_.push_back(f2_energy_string);
00132   this->state_var_layouts_.push_back(dl->qp_scalar);
00133   this->state_var_init_types_.push_back("scalar");
00134   this->state_var_init_values_.push_back(0.0);
00135   this->state_var_old_state_flags_.push_back(true);
00136   this->state_var_output_flags_.push_back(true);
00137   //
00138   // matrix damage
00139   this->num_state_variables_++;
00140   this->state_var_names_.push_back(matrix_damage_string);
00141   this->state_var_layouts_.push_back(dl->qp_scalar);
00142   this->state_var_init_types_.push_back("scalar");
00143   this->state_var_init_values_.push_back(0.0);
00144   this->state_var_old_state_flags_.push_back(true);
00145   this->state_var_output_flags_.push_back(true);
00146   //
00147   // fiber 1 damage
00148   this->num_state_variables_++;
00149   this->state_var_names_.push_back(f1_damage_string);
00150   this->state_var_layouts_.push_back(dl->qp_scalar);
00151   this->state_var_init_types_.push_back("scalar");
00152   this->state_var_init_values_.push_back(0.0);
00153   this->state_var_old_state_flags_.push_back(true);
00154   this->state_var_output_flags_.push_back(true);
00155   //
00156   // fiber 2 damage
00157   this->num_state_variables_++;
00158   this->state_var_names_.push_back(f2_damage_string);
00159   this->state_var_layouts_.push_back(dl->qp_scalar);
00160   this->state_var_init_types_.push_back("scalar");
00161   this->state_var_init_values_.push_back(0.0);
00162   this->state_var_old_state_flags_.push_back(true);
00163   this->state_var_output_flags_.push_back(true);
00164 
00165 }
00166 //------------------------------------------------------------------------------
00167 template<typename EvalT, typename Traits>
00168 void J2FiberModel<EvalT, Traits>::
00169 computeState(typename Traits::EvalData workset,
00170     std::map<std::string, Teuchos::RCP<PHX::MDField<ScalarT> > > dep_fields,
00171     std::map<std::string, Teuchos::RCP<PHX::MDField<ScalarT> > > eval_fields)
00172 {
00173   // extract dependent MDFields
00174   PHX::MDField<ScalarT> def_grad = *dep_fields["F"];
00175   PHX::MDField<ScalarT> J = *dep_fields["J"];
00176   PHX::MDField<ScalarT> poissons_ratio = *dep_fields["Poissons Ratio"];
00177   PHX::MDField<ScalarT> elastic_modulus = *dep_fields["Elastic Modulus"];
00178   PHX::MDField<ScalarT> yield_strength = *dep_fields["Yield Strength"];
00179   PHX::MDField<ScalarT> hardening_modulus = *dep_fields["Hardening Modulus"];
00180   PHX::MDField<MeshScalarT, Cell, QuadPoint, Dim> gpt_location;
00181 
00182   // for now, force using global fiber direction
00183   if (local_coord_flag_) {
00184     gpt_location = this->coord_vec_;
00185   }
00186 
00187   // retrive appropriate field name strings
00188   std::string cauchy_string = (*field_name_map_)["Cauchy_Stress"];
00189   std::string Fp_string = (*field_name_map_)["Fp"];
00190   std::string eqps_string = (*field_name_map_)["eqps"];
00191   std::string matrix_energy_string = (*field_name_map_)["Matrix_Energy"];
00192   std::string f1_energy_string = (*field_name_map_)["F1_Energy"];
00193   std::string f2_energy_string = (*field_name_map_)["F2_Energy"];
00194   std::string matrix_damage_string = (*field_name_map_)["Matrix_Damage"];
00195   std::string f1_damage_string = (*field_name_map_)["F1_Damage"];
00196   std::string f2_damage_string = (*field_name_map_)["F2_Damage"];
00197 
00198   // extract evaluated MDFields
00199   PHX::MDField<ScalarT> stress = *eval_fields[cauchy_string];
00200   PHX::MDField<ScalarT> Fp = *eval_fields[Fp_string];
00201   PHX::MDField<ScalarT> eqps = *eval_fields[eqps_string];
00202   PHX::MDField<ScalarT> energy_m = *eval_fields[matrix_energy_string];
00203   PHX::MDField<ScalarT> energy_f1 = *eval_fields[f1_energy_string];
00204   PHX::MDField<ScalarT> energy_f2 = *eval_fields[f2_energy_string];
00205   PHX::MDField<ScalarT> damage_m = *eval_fields[matrix_damage_string];
00206   PHX::MDField<ScalarT> damage_f1 = *eval_fields[f1_damage_string];
00207   PHX::MDField<ScalarT> damage_f2 = *eval_fields[f2_damage_string];
00208 
00209   // previous state
00210   Albany::MDArray Fp_old =
00211       (*workset.stateArrayPtr)[Fp_string + "_old"];
00212   Albany::MDArray eqps_old =
00213       (*workset.stateArrayPtr)[eqps_string + "_old"];
00214   Albany::MDArray energy_m_old =
00215       (*workset.stateArrayPtr)[matrix_energy_string + "_old"];
00216   Albany::MDArray energy_f1_old =
00217       (*workset.stateArrayPtr)[f1_energy_string + "_old"];
00218   Albany::MDArray energy_f2_old =
00219       (*workset.stateArrayPtr)[f2_energy_string + "_old"];
00220 
00221   ScalarT kappa, mu, mubar, Jm23, K, Y, smag, p;
00222   ScalarT Phi, dgam;
00223   ScalarT trbe, trbeby3;
00224   ScalarT I4_f1, I4_f2;
00225   ScalarT alpha_f1, alpha_f2, alpha_m;
00226   ScalarT sq23 = std::sqrt(2. / 3.);
00227 
00228   // Define some tensors for use
00229   Intrepid::Tensor<ScalarT> F(num_dims_), Fpn(num_dims_), Fpnew(num_dims_);
00230   Intrepid::Tensor<ScalarT> Cpinv(num_dims_), be(num_dims_);
00231   Intrepid::Tensor<ScalarT> s(num_dims_), N(num_dims_);
00232   Intrepid::Tensor<ScalarT> expA(num_dims_);
00233   Intrepid::Tensor<ScalarT> sigma_m(num_dims_);
00234   Intrepid::Tensor<ScalarT> C(num_dims_);
00235   Intrepid::Tensor<ScalarT> sigma_f1(num_dims_), sigma_f2(num_dims_);
00236   Intrepid::Tensor<ScalarT> M1dyadM1(num_dims_), M2dyadM2(num_dims_);
00237   Intrepid::Tensor<ScalarT> S0_f1(num_dims_), S0_f2(num_dims_);
00238   Intrepid::Tensor<ScalarT> I(Intrepid::eye<ScalarT>(num_dims_));
00239 
00240   Intrepid::Vector<ScalarT> M1(num_dims_), M2(num_dims_);
00241 
00242   volume_fraction_m_ = 1.0 - volume_fraction_f1_ - volume_fraction_f2_;
00243 
00244   for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00245     for (std::size_t pt = 0; pt < num_pts_; ++pt) {
00246       // local parameters
00247       kappa = elastic_modulus(cell, pt)
00248           / (3. * (1. - 2. * poissons_ratio(cell, pt)));
00249       mu = elastic_modulus(cell, pt) / (2. * (1. + poissons_ratio(cell, pt)));
00250       K = hardening_modulus(cell, pt);
00251       Y = yield_strength(cell, pt);
00252       Jm23 = std::pow(J(cell, pt), -2. / 3.);
00253 
00254       // fill local tensors
00255       F.fill(&def_grad(cell, pt, 0, 0));
00256       //Fpn.fill( &Fpold(cell,pt,std::size_t(0),std::size_t(0)) );
00257       for (std::size_t i(0); i < num_dims_; ++i) {
00258         for (std::size_t j(0); j < num_dims_; ++j) {
00259           Fpn(i, j) = static_cast<ScalarT>(Fp_old(cell, pt, i, j));
00260         }
00261       }
00262 
00263       // compute Cpinv = inv(Fp) * inv(Fp)^T
00264       Cpinv = Intrepid::inverse(Fpn)
00265           * Intrepid::transpose(Intrepid::inverse(Fpn));
00266 
00267       // compute trial state
00268       be = Jm23 * F * Cpinv * Intrepid::transpose(F);
00269       trbe = Intrepid::trace(be);
00270       trbeby3 = trbe / num_dims_;
00271       mubar = trbeby3 * mu;
00272 
00273       // compute trial deviatoric stress
00274       s = mu * Intrepid::dev(be);
00275 
00276       // check for yielding
00277       smag = Intrepid::norm(s);
00278       Phi = smag - sq23 * (Y + K * eqps_old(cell, pt)
00279           + sat_mod_ * (1.0 - std::exp(-sat_exp_ * eqps_old(cell, pt))));
00280 
00281       // if yielding, find plastic increment via return mapping alg.
00282       if (Phi > 1.0e-11) {
00283         //return mapping algorithm
00284         ScalarT H = 0.0;
00285         ScalarT dH = 0.0;
00286         ScalarT g = Phi;
00287         ScalarT dg = (-2.0 * mubar) * (1.0 + dH / (3.0 * mubar));
00288         ScalarT alpha = 0.0;
00289         ScalarT norm_r = 0.0;
00290         ScalarT relative_r = 0.0;
00291         int iter = 0.0;
00292         dgam = 0.0;
00293 
00294         while (true) {
00295           iter++;
00296 
00297           dgam = dgam - g / dg;
00298           alpha = eqps_old(cell, pt) + sq23 * dgam;
00299           H = K * alpha + sat_mod_ * (1.0 - std::exp(-sat_exp_ * alpha));
00300           dH = K + sat_exp_ * sat_mod_ * std::exp(-sat_exp_ * alpha);
00301 
00302           g = smag - (2.0 * mubar * dgam + sq23 * (Y + H));
00303           dg = -2.0 * mubar * (1.0 + dH / (3.0 * mubar));
00304 
00305           norm_r = std::abs(g);
00306           relative_r = norm_r / Phi;
00307 
00308           if (norm_r < 1.e-11 || relative_r < 1.0e-11)
00309             break;
00310 
00311           if (iter > 25)
00312             break;
00313         }
00314 
00315         // plastic flow direction
00316         N = (1.0 / smag) * s;
00317 
00318         // adjust deviatoric stress to account for plastic increment
00319         s = s - 2.0 * mubar * dgam * N;
00320 
00321         // update eqps
00322         eqps(cell, pt) = alpha;
00323 
00324         // exponential map to get Fp
00325         expA = Intrepid::exp(dgam * N);
00326         Fpnew = expA * Fpn;
00327 
00328         for (std::size_t i(0); i < num_dims_; ++i)
00329           for (std::size_t j(0); j < num_dims_; ++j)
00330             Fp(cell, pt, i, j) = Fpnew(i, j);
00331 
00332       } else {
00333         // elasticity, set state variables to old values
00334         eqps(cell, pt) = eqps_old(cell, pt);
00335         for (std::size_t i(0); i < num_dims_; ++i)
00336           for (std::size_t j(0); j < num_dims_; ++j)
00337             Fp(cell, pt, i, j) = Fpn(i, j);
00338 
00339       }  // end of return mapping
00340 
00341       // compute pressure
00342       p = 0.5 * kappa * (J(cell, pt) - 1. / (J(cell, pt)));
00343 
00344       // compute Cauchy stress for matrix
00345       sigma_m = s / J(cell, pt) + p * I;
00346 
00347       // compute energy for matrix
00348       energy_m(cell, pt) = volume_fraction_m_ * (0.5 * kappa
00349           * (0.5 * (J(cell, pt) * J(cell, pt) - 1.0) - std::log(J(cell, pt)))
00350           + 0.5 * mu * (trbe - 3.0));
00351 
00352       // damage term in matrix
00353       alpha_m = energy_m_old(cell, pt);
00354       if (energy_m(cell, pt) > alpha_m) alpha_m = energy_m(cell, pt);
00355 
00356       damage_m(cell, pt) = max_damage_m_
00357           * (1.0 - std::exp(-alpha_m / saturation_m_));
00358 
00359       //-----------compute stress in Fibers
00360 
00361       // Right Cauchy-Green Tensor C = F^{T} * F
00362       C = Intrepid::transpose(F) * F;
00363 
00364       // Fiber orientation vectors
00365       if (local_coord_flag_) {
00366         // compute fiber orientation based on local coordinates
00367         // special case of plane strain M1(3) = 0; M2(3) = 0;
00368         Intrepid::Vector<ScalarT> gpt(gpt_location(cell, pt, 0),
00369             gpt_location(cell, pt, 1), gpt_location(cell, pt, 2));
00370         Intrepid::Vector<ScalarT> OA(gpt(0) - ring_center_[0],
00371             gpt(1) - ring_center_[1], 0);
00372 
00373         M1 = OA / Intrepid::norm(OA);
00374         M2(0) = -M1(1);
00375         M2(1) = M1(0);
00376         M2(2) = M1(2);
00377       } else {
00378         for (std::size_t i = 0; i < num_dims_; ++i) {
00379           M1(i) = direction_f1_[i];
00380           M2(i) = direction_f2_[i];
00381         }
00382         M1 = M1 / Intrepid::norm(M1);
00383         M2 = M2 / Intrepid::norm(M2);
00384       }
00385 
00386       // Anisotropic invariants I4 = M_{i} * C * M_{i}
00387       I4_f1 = Intrepid::dot(M1, Intrepid::dot(C, M1));
00388       I4_f2 = Intrepid::dot(M2, Intrepid::dot(C, M2));
00389       M1dyadM1 = Intrepid::dyad(M1, M1);
00390       M2dyadM2 = Intrepid::dyad(M2, M2);
00391 
00392       // undamaged stress (2nd PK stress)
00393       S0_f1 = (4.0 * k_f1_ * (I4_f1 - 1.0)
00394           * std::exp(q_f1_ * (I4_f1 - 1) * (I4_f1 - 1))) * M1dyadM1;
00395       S0_f2 = (4.0 * k_f2_ * (I4_f2 - 1.0)
00396           * std::exp(q_f2_ * (I4_f2 - 1) * (I4_f2 - 1))) * M2dyadM2;
00397 
00398       // compute energy for fibers
00399       energy_f1(cell, pt) = volume_fraction_f1_ * (k_f1_
00400           * (std::exp(q_f1_ * (I4_f1 - 1) * (I4_f1 - 1)) - 1) / q_f1_);
00401       energy_f2(cell, pt) = volume_fraction_f2_ * (k_f2_
00402           * (std::exp(q_f2_ * (I4_f2 - 1) * (I4_f2 - 1)) - 1) / q_f2_);
00403 
00404       // Fiber Cauchy stress
00405       sigma_f1 = (1.0 / J(cell, pt))
00406           * Intrepid::dot(F, Intrepid::dot(S0_f1, Intrepid::transpose(F)));
00407       sigma_f2 = (1.0 / J(cell, pt))
00408           * Intrepid::dot(F, Intrepid::dot(S0_f2, Intrepid::transpose(F)));
00409 
00410       // maximum thermodynamic forces
00411       alpha_f1 = energy_f1_old(cell, pt);
00412       alpha_f2 = energy_f2_old(cell, pt);
00413 
00414       if (energy_f1(cell, pt) > alpha_f1) alpha_f1 = energy_f1(cell, pt);
00415 
00416       if (energy_f2(cell, pt) > alpha_f2) alpha_f2 = energy_f2(cell, pt);
00417 
00418       // damage term in fibers
00419       damage_f1(cell, pt) =
00420           max_damage_f1_ * (1 - std::exp(-alpha_f1 / saturation_f1_));
00421       damage_f2(cell, pt) =
00422           max_damage_f1_ * (1 - std::exp(-alpha_f2 / saturation_f2_));
00423 
00424       // total Cauchy stress (M, Fibers)
00425       for (std::size_t i(0); i < num_dims_; ++i) {
00426         for (std::size_t j(0); j < num_dims_; ++j) {
00427           stress(cell, pt, i, j) =
00428               volume_fraction_m_ * (1 - damage_m(cell, pt)) * sigma_m(i, j)
00429                   + volume_fraction_f1_ * (1 - damage_f1(cell, pt))
00430                       * sigma_f1(i, j)
00431                   + volume_fraction_f2_ * (1 - damage_f2(cell, pt))
00432                       * sigma_f2(i, j);
00433         }
00434       }
00435     }
00436   }
00437 }
00438 //------------------------------------------------------------------------------
00439 }
00440
J2FiberModel_Def.hpp
