OpenSees » SRC » material » nD » NewTemplate3Dep

elnp_Elastic.cpp

Go 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 "e = e_ref - log(p/p_ref)"
00035 // So K = (1+e)*p/kappa (Wood, 1990, Soil Behavior and Critical State Soil Mechanics)
00036 // Parameters:
00037 // 1: kappa: unliading slope in e-log(p) relation curve
00038 // 2: v:     poissoin's ratio
00039 // 3: Kc:    cut-off bulk modulus, when K<Kc, let K = Kc
00040 // 4: e0:    initial void ratio
00041 
00042 #ifndef elnp_Elastic_CPP
00043 #define elnp_Elastic_CPP
00044 
00045 #include "elnp_Elastic.h"
00046 
00047 //BJtensor elnp_Elastic::StiffnessH(4, def_dim_4, 0.0);
00048 
00049 elnp_Elastic::elnp_Elastic(int kappa_in, 
00050                            int v_in,
00051                            int K_c_in,
00052                            int e0_in,
00053                            const stresstensor& initialStress,
00054                            const straintensor& initialStrain)
00055 : ElasticState(initialStress, initialStrain),
00056   kappa_index(kappa_in),
00057   v_index(v_in),
00058   K_c_index(K_c_in),
00059   e0_index(e0_in)
00060 {
00061 
00062 }
00063 
00064 // Create a new 
00065 ElasticState* elnp_Elastic::newObj() 
00066 {
00067     ElasticState * Els = new  elnp_Elastic(this->kappa_index, 
00068                                            this->v_index, 
00069                                            this->K_c_index,
00070                                            this->e0_index,
00071                                            this->Stress,
00072                                            this->Strain);
00073      return Els;
00074 }
00075 
00076 // Get Stiffness Tensor
00077 const BJtensor& elnp_Elastic::getElasticStiffness(const MaterialParameter &MaterialParameter_in) const
00078 {
00079     // Kronecker delta tensor
00080     BJtensor I2("I", 2, def_dim_2);
00081 
00082     BJtensor I_ijkl = I2("ij")*I2("kl");
00083     I_ijkl.null_indices();
00084     BJtensor I_ikjl = I_ijkl.transpose0110();
00085     BJtensor I_iljk = I_ijkl.transpose0111();
00086     BJtensor I4s = (I_ikjl+I_iljk)*0.5;
00087     
00088     double kappa = getkappa(MaterialParameter_in);
00089     double v = getv(MaterialParameter_in);
00090     if (v >= 0.5 || v < -1.0) {
00091       opserr << "Warning!! elnp_Elastic: Invalid Possoin's ratio. " << endln;
00092       exit (1);
00093     }
00094     double Kc = getK_c(MaterialParameter_in);
00095     double e0 = gete0(MaterialParameter_in);
00096     
00097     double p = this->getStress().p_hydrostatic();
00098     double epsilon_v = this->getStrain().Iinvariant1();
00099     double e = e0 + (1 + e0) *epsilon_v;
00100     double Kk = (1.0 + e) *p /kappa;
00101     
00102     double K = (Kk > Kc) ? Kk : Kc ;
00103     double G = K *1.5*(1.0-2.0*v)/(1.0+v);
00104     
00105     // To avoid numerical problems
00106     //if (G < (1.0e-3) *K)  G = 1.0e-3 *K;
00107        
00108     // Building elasticity tensor
00109     ElasticState::ElasticStiffness = I_ijkl *(K - 2.0*G/3.0) + I4s *(2.0*G);
00110 
00111     return ElasticState::ElasticStiffness;
00112 }
00113 
00114 // Get kappa
00115 double elnp_Elastic::getkappa(const MaterialParameter &MaterialParameter_in) const
00116 {
00117     if ( kappa_index > MaterialParameter_in.getNum_Material_Parameter() || kappa_index < 2) { 
00118         opserr << "elnp_Elastic: Invalid Input. " << endln;
00119         exit (1);
00120     }
00121     else
00122         return MaterialParameter_in.getMaterial_Parameter(kappa_index - 1); 
00123 }
00124 
00125 // Get v
00126 double elnp_Elastic::getv(const MaterialParameter &MaterialParameter_in) const
00127 {
00128     if ( v_index > MaterialParameter_in.getNum_Material_Parameter() || v_index < 2) { 
00129         opserr << "elnp_Elastic: Invalid Input. " << endln;
00130         exit (1);
00131     }
00132     else
00133       return MaterialParameter_in.getMaterial_Parameter(v_index - 1); 
00134 }
00135 
00136 // Get K_c
00137 double elnp_Elastic::getK_c(const MaterialParameter &MaterialParameter_in) const
00138 {
00139     if ( K_c_index > MaterialParameter_in.getNum_Material_Parameter() || K_c_index < 2) { 
00140         opserr << "elnp_Elastic: Invalid Input. " << endln;
00141         exit (1);
00142     }
00143     else
00144         return MaterialParameter_in.getMaterial_Parameter(K_c_index - 1); 
00145 }
00146 
00147 // Get e0
00148 double elnp_Elastic::gete0(const MaterialParameter &MaterialParameter_in) const
00149 {
00150     if ( e0_index > MaterialParameter_in.getNum_Material_Parameter() || e0_index < 2) { 
00151         opserr << "elnp_Elastic: Invalid Input. " << endln;
00152         exit (1);
00153     }
00154     else
00155         return MaterialParameter_in.getMaterial_Parameter(e0_index - 1); 
00156 }
00157 
00158 #endif

Generated on Mon Oct 23 15:05:15 2006 for OpenSees by doxygen 1.5.0