DM04_Elastic.cppGo to the documentation of this file.00001 00002 // COPYLEFT (C): Woody's viral GPL-like license (by BJ): 00003 // ``This source code is Copyrighted in 00004 // U.S., for an indefinite period, and anybody 00005 // caught using it without our permission, will be 00006 // mighty good friends of ourn, cause we don't give 00007 // a darn. Hack it. Compile it. Debug it. Run it. 00008 // Yodel it. Enjoy it. We wrote it, that's all we 00009 // wanted to do.'' 00010 // 00011 // 00012 // COPYRIGHT (C): :-)) 00013 // PROJECT: Object Oriented Finite Element Program 00014 // FILE: 00015 // CLASS: 00016 // MEMBER FUNCTIONS: 00017 // 00018 // MEMBER VARIABLES 00019 // 00020 // PURPOSE: 00021 // 00022 // RETURN: 00023 // VERSION: 00024 // LANGUAGE: C++ 00025 // TARGET OS: 00026 // DESIGNER: Zhao Cheng, Boris Jeremic 00027 // PROGRAMMER: Zhao Cheng, 00028 // DATE: Fall 2005 00029 // UPDATE HISTORY: 00030 // 00032 // 00033 00034 // This is based on G = G0*Pat*[(2.97-e)^2/(1+e)]*(p/Pat)^0.5 00035 // Richart et al 1970, Li & Dafalias 2000, Dafalias & Mazari 2004 00036 // Parameters: 00037 // 1: G0: reference shear mudulus factor (no unit); 00038 // 2: v: Poisson's ratio; 00039 // 3: Pat: Atmospheric pressure; 00040 // 4: k_c; cut-off factor, for p < k_c*Pat, let p = k_c*Pat to calculate G; 00041 // 5: e0; initial void ratio 00042 00043 #ifndef DM04_Elastic_CPP 00044 #define DM04_Elastic_CPP 00045 00046 #include "DM04_Elastic.h" 00047 00048 //BJtensor DM04_Elastic::StiffnessH(4, def_dim_4, 0.0); 00049 00050 DM04_Elastic::DM04_Elastic(int G0_in, 00051 int v_in, 00052 int Pat_in, 00053 int k_c_in, 00054 int e0_in, 00055 const stresstensor& initialStress, 00056 const straintensor& initialStrain) 00057 : ElasticState(initialStress, initialStrain), 00058 G0_index(G0_in), 00059 v_index(v_in), 00060 Pat_index(Pat_in), 00061 k_c_index(k_c_in), 00062 e0_index(e0_in) 00063 { 00064 00065 } 00066 00067 // Create a new 00068 ElasticState* DM04_Elastic::newObj() 00069 { 00070 ElasticState * Els = new DM04_Elastic(this->G0_index, 00071 this->v_index, 00072 this->Pat_index, 00073 this->k_c_index, 00074 this->e0_index, 00075 this->Stress, 00076 this->Strain); 00077 return Els; 00078 } 00079 00080 // Get Stiffness Tensor 00081 const BJtensor& DM04_Elastic::getElasticStiffness(const MaterialParameter &MaterialParameter_in) const 00082 { 00083 // Kronecker delta tensor 00084 BJtensor I2("I", 2, def_dim_2); 00085 00086 BJtensor I_ijkl = I2("ij")*I2("kl"); 00087 I_ijkl.null_indices(); 00088 BJtensor I_ikjl = I_ijkl.transpose0110(); 00089 BJtensor I_iljk = I_ijkl.transpose0111(); 00090 BJtensor I4s = (I_ikjl+I_iljk)*0.5; 00091 00092 double G0 = getG0(MaterialParameter_in); 00093 double v = getv(MaterialParameter_in); 00094 double Pat = getPat(MaterialParameter_in); 00095 double k_c = getk_c(MaterialParameter_in); 00096 double e0 = gete0(MaterialParameter_in); 00097 00098 double epsilon_v = this->getStrain().Iinvariant1(); 00099 double e = e0 + (1 + e0) *epsilon_v; 00100 double ef = (2.97-e)*(2.97-e)/(1.0+e); 00101 double p_cal = this->getStress().p_hydrostatic(); 00102 if (p_cal < 0.0) 00103 p_cal = 0.0; 00104 double p_cut = Pat *k_c; 00105 00106 double p = (p_cal > p_cut) ? p_cal : p_cut; 00107 00108 double G = G0 *Pat *ef *sqrt(p/Pat); 00109 double K = G * (2.0*(1.0+v)/(3.0*(1.0-2.0*v))); 00110 00111 // Building elasticity tensor 00112 ElasticState::ElasticStiffness = I_ijkl *(K - 2.0*G/3.0) + I4s *(2.0*G); 00113 00114 return ElasticState::ElasticStiffness; 00115 } 00116 00118 stresstensor DM04_Elastic::getStress() const 00119 { 00120 return Stress; 00121 } 00122 00123 // Get G0 00124 double DM04_Elastic::getG0(const MaterialParameter &MaterialParameter_in) const 00125 { 00126 if ( G0_index > MaterialParameter_in.getNum_Material_Parameter() || G0_index < 2) { 00127 opserr << "DM04_Elastic: Invalid Input. " << endln; 00128 exit (1); 00129 } 00130 else 00131 return MaterialParameter_in.getMaterial_Parameter(G0_index - 1); 00132 } 00133 00134 // Get v 00135 double DM04_Elastic::getv(const MaterialParameter &MaterialParameter_in) const 00136 { 00137 double v = 0.0; 00138 if ( v_index > MaterialParameter_in.getNum_Material_Parameter() || v_index < 2) { 00139 opserr << "DM04_Elastic: Invalid Input. " << endln; 00140 exit (1); 00141 } 00142 else { 00143 v = MaterialParameter_in.getMaterial_Parameter(v_index - 1); 00144 if (v >= 0.5 || v <= -1.0) { 00145 opserr << "Warning!! DM04_Elastic: Invalid possoin's ratio. " << endln; 00146 exit (1); 00147 } 00148 return v; 00149 } 00150 } 00151 00152 // Get Pat 00153 double DM04_Elastic::getPat(const MaterialParameter &MaterialParameter_in) const 00154 { 00155 if ( Pat_index > MaterialParameter_in.getNum_Material_Parameter() || Pat_index < 1) { 00156 opserr << "DM04_Elastic: Invalid Input. " << endln; 00157 exit (1); 00158 } 00159 else 00160 return MaterialParameter_in.getMaterial_Parameter(Pat_index - 1); 00161 } 00162 00163 // Get k_cut 00164 double DM04_Elastic::getk_c(const MaterialParameter &MaterialParameter_in) const 00165 { 00166 if ( k_c_index > MaterialParameter_in.getNum_Material_Parameter() || k_c_index < 1) { 00167 opserr << "DM04_Elastic: Invalid Input. " << endln; 00168 exit (1); 00169 } 00170 else 00171 return MaterialParameter_in.getMaterial_Parameter(k_c_index - 1); 00172 } 00173 00174 // Get e0 00175 double DM04_Elastic::gete0(const MaterialParameter &MaterialParameter_in) const 00176 { 00177 if ( e0_index > MaterialParameter_in.getNum_Material_Parameter() || e0_index < 2) { 00178 opserr << "DM04_Elastic: Invalid Input. " << endln; 00179 exit (1); 00180 } 00181 else 00182 return MaterialParameter_in.getMaterial_Parameter(e0_index - 1); 00183 } 00184 00185 #endif 00186
Generated on Mon Oct 23 15:05:14 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents