1 /*! 2 * \file src/NUMODIS/FrankRead.cxx 3 * \brief 4 * \author Laurent Dupuy 5 * \date 9/06/2017 6 * \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All rights 7 * reserved. 8 * This project is publicly released under either the GNU GPL Licence 9 * or the CECILL-A licence. A copy of thoses licences are delivered 10 * with the sources of TFEL. CEA or EDF may also distribute this 11 * project under specific licensing conditions. 12 */ 13 14 #include <cmath> 15 #include"TFEL/Math/General/IEEE754.hxx" 16 #include"TFEL/Math/General/CubicRoots.hxx" 17 #include"NUMODIS/FrankRead.hxx" 18 19 namespace numodis{ 20 21 //====================================================== 22 // FrankRead::FrankRead 23 //------------------------------------------------------ 24 //! Constructor 25 //------------------------------------------------------ 26 /*! 27 \param mu shear modulus 28 \param burgers Burgers vector 29 \param nu Poisson coefficient 30 */ 31 //====================================================== FrankRead(const IsotropicLineTensionModel & linetension)32 FrankRead::FrankRead(const IsotropicLineTensionModel& linetension) 33 :_linetension(linetension) 34 {} 35 36 //================================================================= 37 // FrankRead::computeBetaF 38 //----------------------------------------------------------------- 39 //! Compute the final tangent angle knowing the first one 40 //----------------------------------------------------------------- 41 /*! 42 \param alpha angle between ... and ... 43 \param beta0 first tangent angle in range[0;pi] 44 \param tol tolerance / precision of the calculation 45 \param Nmax maximum number of iterations 46 \return betaF final tangent angle 47 */ 48 //================================================================= computeBetaF(double alpha,double beta0,double tol,int Nmax) const49 double FrankRead::computeBetaF(double alpha,double beta0,double tol,int Nmax) const 50 { 51 52 // pi = 3.14159... 53 double pi = atan(1.0)*4.0; 54 55 // range for betaF 56 double bFmin = beta0 - pi; 57 double bFmax = 0.0; 58 if ( bFmin > bFmax ) 59 throw -1; // throw: ill-defined beta0 value 60 61 double Eq4min = this->Equation4(alpha,beta0,bFmin); 62 double Eq4max = this->Equation4(alpha,beta0,bFmax); 63 if ( Eq4min*Eq4max > 0.0 ) 64 throw -2; // throw: both sides have positive values 65 66 // find betaF using a simple bisection method 67 int N = 0; 68 while ( N < Nmax ) { 69 double bFmid = ( bFmax + bFmin ) / 2.0; 70 double Eq4mid = this->Equation4(alpha,beta0,bFmid); 71 if ((tfel::math::ieee754::fpclassify(Eq4mid)==FP_ZERO) || 72 (( bFmax - bFmin ) / 2.0 < tol)){ 73 return bFmid; 74 } 75 ++N; 76 77 if ( std::signbit(Eq4mid) == std::signbit(Eq4min) ){ 78 bFmin = bFmid; 79 Eq4min = Eq4mid; 80 } else { 81 bFmax = bFmid; 82 Eq4max = Eq4mid; 83 } 84 } 85 throw -3; // throw: calculation did not converge 86 87 } 88 89 //================================================================= 90 // FrankRead::computeActivationAngle 91 //----------------------------------------------------------------- 92 //! Compute the tangent angle beta0 at activation stress 93 //----------------------------------------------------------------- 94 /*! 95 This method was directly taken out of the code used by 96 Dupuy & Fivel in Acta Materiala 2002 (jojo5). 97 98 \param phi0 angle between the initial source and the Burgers 99 \return activation angle betaF (in range[0:pi]) 100 101 */ 102 //================================================================= computeActivationAngle(double phi0) const103 double FrankRead::computeActivationAngle(double phi0) const 104 { 105 using tfel::math::CubicRoots; 106 double nu = _linetension.getNu(); 107 double alpha = 1. - nu *cos(phi0)*cos(phi0); 108 double beta = 1. - nu *sin(phi0)*sin(phi0); 109 double gamma = -nu*sin(2.*phi0); 110 double Q = (2.*alpha-beta)/(3.*alpha); 111 double R = -gamma/(2.*alpha); 112 double D = Q*Q*Q + R*R; 113 double SC = R+sqrt(D); 114 double TC = R-sqrt(D); 115 double x0 = CubicRoots::cbrt(SC) + CubicRoots::cbrt(TC); 116 return atan2( 1.,x0 ); 117 } 118 119 } // end of namespace numodis 120