RMC01_PS.cppGo to the documentation of this file.00001 //================================================================================ 00002 //# COPY LEFT and RIGHT: # 00003 //# Commercial use of this program without express permission of the # 00004 //# University of California, is strictly encouraged. Copyright and Copyleft # 00005 //# are covered by the following clause: # 00006 //# # 00007 //# Woody's license: # 00008 //# ``This source code is Copyrighted in U.S., by the The Regents of the # 00009 //# University of California, for an indefinite period, and anybody caught # 00010 //# using it without our permission, will be mighty good friends of ourn, # 00011 //# cause we don't give a darn. Hack it. Compile it. Debug it. Run it. Yodel # 00012 //# it. Enjoy it. We wrote it, that's all we wanted to do.'' bj # 00013 //# # 00014 //# # 00015 //# # 00016 //# PROJECT: Object Oriented Finite Element Program # 00017 //# PURPOSE: Rounded Mohr Coulomb Potential Surface # 00018 //# CLASS: # 00019 //# # 00020 //# VERSION: # 00021 //# LANGUAGE: C++ # 00022 //# TARGET OS: DOS || UNIX || . . . # 00023 //# DESIGNER(S): Boris Jeremic jeremic@ucdavis.edu # 00024 //# Zhao Cheng, # 00025 //# PROGRAMMER(S): Zhao Cheng, Boris Jeremic # 00026 //# # 00027 //# # 00028 //# DATE: 12 Feb. 2003 # 00029 //# UPDATE HISTORY: Feb 25 2003 # 00030 //# # 00031 //# # 00032 //# # 00033 //# # 00034 //# # 00035 //################################################################################ 00036 //*/ 00037 00038 #ifndef RMC01_PS_CPP 00039 #define RMC01_PS_CPP 00040 00041 #include "RMC01_PS.h" 00042 00043 //================================================================================ 00044 //create a clone of itself 00045 //================================================================================ 00046 00047 PotentialSurface * RMC01PotentialSurface::newObj() 00048 { 00049 PotentialSurface *new_YS = new RMC01PotentialSurface(); 00050 return new_YS; 00051 } 00052 00053 00054 00055 //================================================================================ 00056 // tensor dQ/dsigma_ij 00057 //================================================================================ 00058 00059 tensor RMC01PotentialSurface::dQods(const EPState *EPS) const 00060 { 00061 tensor dQoverds( 2, def_dim_2, 0.0); 00062 // double p = EPS->getStress().p_hydrostatic(); // p 00063 double q = EPS->getStress().q_deviatoric(); // q 00064 double theta = EPS->getStress().theta(); // theta 00065 // double temp_phi = EPS->getScalarVar(1)*3.14159265358979/180.0; // frictional angle 00066 // double temp_cohesive = EPS->getScalarVar(2); // cohesion 00067 tensor DpoDs = EPS->getStress().dpoverds(); // dp/ds 00068 tensor DqoDs = EPS->getStress().dqoverds(); // dq/ds 00069 tensor DthetaoDs = EPS->getStress().dthetaoverds(); // d(theta)/ds 00070 double alfa = EPS->getScalarVar(1); // Take alfa & k as internal variables 00071 // double k = EPS->getScalarVar(2); // instead of phi & conhesive 00072 double a1 = (3.0*1.7320508076*alfa) / (2.0+1.7320508076*alfa); 00073 // double a1 = -6*sin(temp_phi)/(3.0-sin(temp_phi)); 00074 // double a2 = -6*temp_cohesive*cos(temp_phi)/(3.0-sin(temp_phi)); 00075 double e = (3.0-a1)/(3.0+a1); 00076 double Frou = g_0(theta, e); // r(theta) 00077 double Frou_prime = g_prime(theta, e); // r'(theta) 00078 // double dQoverdp = a1; // dQ/dp 00079 // double dQoverdq = Frou; // dQ/dq 00080 // double dQoverdtheta = q*Frou_prime; // dQ/d(theta) 00081 double dQoverdp = alfa*(-3.0); // dQ/dp 00082 double dQoverdq = Frou/1.7320508076; // dQ/dq 00083 double dQoverdtheta = q*Frou_prime/1.7320508076; // dQ/d(theta) 00084 00085 dQoverds = DpoDs * dQoverdp + 00086 DqoDs * dQoverdq + 00087 DthetaoDs * dQoverdtheta; 00088 00089 return dQoverds; 00090 00091 } 00092 00093 00094 //================================================================================ 00095 // tensor d2Q/dsigma_ij_2 00096 //================================================================================ 00097 00098 tensor RMC01PotentialSurface::d2Qods2( const EPState *EPS ) const 00099 { 00100 tensor d2Qoverds2( 4, def_dim_4, 0.0); // d2Q/ds2(pqmn) 00101 00102 // double p = EPS->getStress().p_hydrostatic(); // p 00103 double q = EPS->getStress().q_deviatoric(); // q 00104 double theta = EPS->getStress().theta(); // theta 00105 // double temp_phi = EPS->getScalarVar(1)*3.14159265358979/180.0; // frictional angle 00106 // double temp_cohesive = EPS->getScalarVar(2); // conhesion 00107 double alfa = EPS->getScalarVar(1); 00108 // double k = EPS->getScalarVar(2); 00109 double a1 = (3.0*1.7320508076*alfa) / (2.0+1.7320508076*alfa); 00110 tensor DpoDs = EPS->getStress().dpoverds(); // dp/ds 00111 tensor DqoDs = EPS->getStress().dqoverds(); // dq/ds 00112 tensor DthetaoDs = EPS->getStress().dthetaoverds(); // d(theta)/ds 00113 tensor D2poDs2 = EPS->getStress().d2poverds2(); // d2p/ds2 00114 tensor D2qoDs2 = EPS->getStress().d2qoverds2(); // d2q/ds2 00115 tensor D2thetaoDs2 = EPS->getStress().d2thetaoverds2(); // d2(theta)/ds2 00116 // double a1 = -6*sin(temp_phi)/(3.0-sin(temp_phi)); 00117 // double a2 = -6*temp_cohesive*cos(temp_phi)/(3.0-sin(temp_phi)); 00118 double e = (3.0-a1)/(3.0+a1); 00119 double Frou = g_0(theta, e); // r(theta) 00120 double Frou_prime = g_prime(theta, e); // r'(theta) 00121 double Frou_second = g_second(theta, e); // r"(theta) 00122 // double dQoverdp = a1; // dQ/dp 00123 // double dQoverdq = Frou; // dQ/dq 00124 // double dQoverdtheta = q*Frou_prime; // dQ/d(theta) 00125 double dQoverdp = alfa*(-3.0); // dQ/dp 00126 double dQoverdq = Frou/1.7320508076; // dQ/dq 00127 double dQoverdtheta = q*Frou_prime/1.7320508076; // dQ/d(theta) 00128 // double a23 = Frou_prime; // d2Q/dqd(theta) 00129 // double a32 = a23; // d2Q/d(theta)dq 00130 // double a33 = q * Frou_second; // d2Q/d(theta)2 00131 double a23 = Frou_prime/1.7320508076; // d2Q/dqd(theta) 00132 double a32 = a23/1.7320508076; // d2Q/d(theta)dq 00133 double a33 = q * Frou_second/1.7320508076; // d2Q/d(theta)2 00134 00135 // d2Qoverds2 = DthetaoDs("mn") * DqoDs("pq") * a23 + 00136 // DqoDs("mn") * DthetaoDs("pq") * a32 + 00137 // DthetaoDs("mn") * DthetaoDs("pq") * a33 + 00138 // D2poDs2("pqmn") * dQoverdp + 00139 // D2qoDs2("pqmn") * dQoverdq + 00140 // D2thetaoDs2("pqmn") * dQoverdtheta; 00141 d2Qoverds2 = DqoDs("pq") * DthetaoDs("mn") * a23 + 00142 DthetaoDs("pq") * DqoDs("mn") * a32 + 00143 DthetaoDs("pq") * DthetaoDs("mn") * a33 + 00144 D2poDs2("pqmn") * dQoverdp + 00145 D2qoDs2("pqmn") * dQoverdq + 00146 D2thetaoDs2("pqmn") * dQoverdtheta; 00147 return d2Qoverds2; 00148 00149 } 00150 00151 // For Consistent Algorithm, Z Cheng, Jan 2004 00152 tensor RMC01PotentialSurface::d2Qodsds1(const EPState *EPS) const 00153 { 00154 tensor I("I", 2, def_dim_2); 00155 tensor d2Qoverdsds1 = I; 00156 return d2Qoverdsds1; 00157 } 00158 00159 //================================================================================ 00160 // friend ostream functions for output 00161 //================================================================================ 00162 00163 OPS_Stream & operator<< (OPS_Stream & os, const RMC01PotentialSurface & PS) 00164 { 00165 os << "ROUNDED MC Potential Surface Parameters: " << endln; 00166 return os; 00167 } 00168 00169 00170 #endif 00171 |