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