My title

00001 /* ****************************************************************** **
00002 **    OpenSees - Open System for Earthquake Engineering Simulation    **
00003 **          Pacific Earthquake Engineering Research Center            **
00004 **                                                                    **
00005 **                                                                    **
00006 ** (C) Copyright 1999, The Regents of the University of California    **
00007 ** All Rights Reserved.                                               **
00008 **                                                                    **
00009 ** Commercial use of this program without express permission of the   **
00010 ** University of California, Berkeley, is strictly prohibited.  See   **
00011 ** file 'COPYRIGHT'  in main directory for information on usage and   **
00012 ** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
00013 **                                                                    **
00014 ** Developed by:                                                      **
00015 **   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
00016 **   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
00017 **   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
00018 **                                                                    **
00019 ** ****************************************************************** */
00020                                                                         
00021 // $Revision: 1.4 $
00022 // $Date: 2001/06/20 04:37:09 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/material/uniaxial/HystereticMaterial.cpp,v $
00024 
00025 // Written: MHS
00026 // Created: July 2000
00027 //
00028 // Description: This file contains the implementation of 
00029 // HystereticMaterial.  HystereticMaterial is
00030 // a one-dimensional hysteretic model with pinching of both
00031 // force and deformation, damage due to deformation and energy, and
00032 // degraded unloading stiffness based on maximum ductility.  This
00033 // is a modified implementation of Hyster2.f90 by Filippou.
00034 
00035 #include <HystereticMaterial.h>
00036 #include <G3Globals.h>
00037 #include <math.h>
00038 #include <float.h>
00039 
00040 HystereticMaterial::HystereticMaterial(int tag,
00041    double m1p, double r1p, double m2p, double r2p, double m3p, double r3p,
00042    double m1n, double r1n, double m2n, double r2n, double m3n, double r3n,
00043    double px, double py, double d1, double d2, double b):
00044 UniaxialMaterial(tag, MAT_TAG_Hysteretic),
00045 mom1p(m1p), rot1p(r1p), mom2p(m2p), rot2p(r2p), mom3p(m3p), rot3p(r3p),
00046 mom1n(m1n), rot1n(r1n), mom2n(m2n), rot2n(r2n), mom3n(m3n), rot3n(r3n),
00047 pinchX(px), pinchY(py), damfc1(d1), damfc2(d2), beta(0.0)
00048 {
00049     if (b != 0.0)
00050         g3ErrorHandler->warning("%s -- beta parameter has been temporarily disabled",
00051             "HystereticMaterial::HystereticMaterial");
00052 
00053  bool error = false;
00054  
00055  // Positive backbone parameters
00056  if (rot1p <= 0.0)
00057   error = true;
00058 
00059  if (rot2p <= rot1p)
00060   error = true;
00061 
00062  if (rot3p <= rot2p)
00063   error = true;
00064 
00065  // Negative backbone parameters
00066  if (rot1n >= 0.0)
00067   error = true;
00068 
00069  if (rot2n >= rot1n)
00070   error = true;
00071 
00072  if (rot3n >= rot2n)
00073   error = true;
00074  
00075  if (error)
00076   g3ErrorHandler->fatal("%s -- input backbone is not unique (one-to-one)",
00077   "HystereticMaterial::HystereticMaterial");
00078 
00079  energyA = 0.5 * (rot1p*mom1p + (rot2p-rot1p)*(mom2p+mom1p) + (rot3p-rot2p)*(mom3p+mom2p) +
00080   rot1n*mom1n + (rot2n-rot1n)*(mom2n+mom1n) * (rot3n-rot2n)*(mom3n+mom2n));
00081 
00082  // Set envelope slopes
00083  this->setEnvelope();
00084 
00085  // Initialize history variables
00086  this->revertToStart();
00087  this->revertToLastCommit();
00088 }
00089 
00090 HystereticMaterial::HystereticMaterial(int tag,
00091    double m1p, double r1p, double m2p, double r2p,
00092    double m1n, double r1n, double m2n, double r2n,
00093    double px, double py, double d1, double d2, double b):
00094 UniaxialMaterial(tag, MAT_TAG_Hysteretic),
00095 mom1p(m1p), rot1p(r1p), mom3p(m2p), rot3p(r2p),
00096 mom1n(m1n), rot1n(r1n), mom3n(m2n), rot3n(r2n),
00097 pinchX(px), pinchY(py), damfc1(d1), damfc2(d2), beta(0.0)
00098 {
00099     if (b != 0.0)
00100         g3ErrorHandler->warning("%s -- beta parameter has been temporarily disabled",
00101             "HystereticMaterial::HystereticMaterial");
00102 
00103  bool error = false;
00104  
00105  // Positive backbone parameters
00106  if (rot1p <= 0.0)
00107   error = true;
00108 
00109  if (rot3p <= rot1p)
00110   error = true;
00111 
00112  // Negative backbone parameters
00113  if (rot1n >= 0.0)
00114   error = true;
00115 
00116  if (rot3n >= rot1n)
00117   error = true;
00118 
00119  if (error)
00120   g3ErrorHandler->fatal("%s -- input backbone is not unique (one-to-one)",
00121   "HystereticMaterial::HystereticMaterial");
00122 
00123  energyA = 0.5 * (rot1p*mom1p + (rot3p-rot1p)*(mom3p+mom1p) +
00124   rot1n*mom1n + (rot3n-rot1n)*(mom3n+mom1n));
00125 
00126  mom2p = 0.5*(mom1p+mom3p);
00127  mom2n = 0.5*(mom1n+mom3n);
00128 
00129  rot2p = 0.5*(rot1p+rot3p);
00130  rot2n = 0.5*(rot1n+rot3n);
00131 
00132  // Set envelope slopes
00133  this->setEnvelope();
00134 
00135  // Initialize history variables
00136  this->revertToStart();
00137  this->revertToLastCommit();
00138 }
00139 
00140 HystereticMaterial::HystereticMaterial():
00141 UniaxialMaterial(0, MAT_TAG_Hysteretic),
00142 mom1p(0.0), rot1p(0.0), mom2p(0.0), rot2p(0.0), mom3p(0.0), rot3p(0.0),
00143 mom1n(0.0), rot1n(0.0), mom2n(0.0), rot2n(0.0), mom3n(0.0), rot3n(0.0),
00144 pinchX(0.0), pinchY(0.0), damfc1(0.0), damfc2(0.0), beta(0.0)
00145 {
00146 
00147 }
00148 
00149 HystereticMaterial::~HystereticMaterial()
00150 {
00151 
00152 }
00153 
00154 int
00155 HystereticMaterial::setTrialStrain(double strain, double strainRate)
00156 {
00157  TrotMax = CrotMax;
00158  TrotMin = CrotMin;
00159  TenergyD = CenergyD;
00160  TrotPu = CrotPu;
00161  TrotNu = CrotNu;
00162 
00163  Tstrain = strain;
00164  double dStrain = Tstrain - Cstrain;
00165 
00166  TloadIndicator = CloadIndicator;
00167  
00168  if (TloadIndicator == 0)
00169   TloadIndicator = (dStrain < 0.0) ? 2 : 1;
00170 
00171  if (Tstrain >= CrotMax) {
00172   TrotMax = Tstrain;
00173   Ttangent = posEnvlpTangent(Tstrain);
00174   Tstress = posEnvlpStress(Tstrain);
00175  }
00176  else if (Tstrain <= CrotMin) {
00177   TrotMin = Tstrain;
00178   Ttangent = negEnvlpTangent(Tstrain);
00179   Tstress = negEnvlpStress(Tstrain);
00180  }
00181  else {
00182    if (dStrain < 0.0)
00183      negativeIncrement(dStrain);
00184    else if (dStrain > 0.0)
00185      positiveIncrement(dStrain);
00186  }
00187 
00188  TenergyD = CenergyD + 0.5*(Cstress+Tstress)*dStrain;
00189 
00190  return 0;
00191 }
00192 
00193 
00194 double
00195 HystereticMaterial::getStrain(void)
00196 {
00197  return Tstrain;
00198 }
00199 
00200 double
00201 HystereticMaterial::getStress(void)
00202 {
00203  return Tstress;
00204 }
00205 
00206 double
00207 HystereticMaterial::getTangent(void)
00208 {
00209  return Ttangent;
00210 }
00211 
00212 void
00213 HystereticMaterial::positiveIncrement(double dStrain)
00214 {
00215  double k = pow((CrotMin+1.0e-6)/rot1n,beta);
00216  k = (k < 1.0) ? 1.0 : 1.0/k;
00217 
00218  if (TloadIndicator == 2) {
00219   TloadIndicator = 1;
00220   if (Cstress <= 0.0) {
00221    TrotNu = Cstrain - Cstress/(E1n*k);
00222    double energy = CenergyD - 0.5*Cstress/(E1n*k)*Cstress;
00223    double damfc = 1.0;
00224    if (CrotMin < rot1n) {
00225     damfc += damfc2*energy/energyA;
00226 
00227     if (Cstrain == CrotMin) {
00228      damfc += damfc1*(CrotMax/rot1p-1.0);
00229     }
00230    }
00231 
00232    TrotMax = CrotMax * damfc;
00233   }
00234  }
00235 
00236   TloadIndicator = 1;
00237 
00238  TrotMax = (TrotMax > rot1p) ? TrotMax : rot1p;
00239 
00240  double maxmom = posEnvlpStress(TrotMax);
00241  double rotlim = negEnvlpRotlim(CrotMin);
00242  double rotrel = (rotlim > TrotNu) ? rotlim : TrotNu;
00243  rotrel = TrotNu;
00244  if (negEnvlpStress(CrotMin) >= 0.0)
00245    rotrel = rotlim;
00246 
00247  double rotmp1 = rotrel + pinchY*(TrotMax-rotrel);
00248  double rotmp2 = TrotMax - (1.0-pinchY)*maxmom/(E1n*k);
00249  double rotch = rotmp1 + (rotmp2-rotmp1)*pinchX;
00250 
00251  double tmpmo1;
00252  double tmpmo2;
00253 
00254  if (Tstrain < TrotNu) {
00255   Ttangent = E1n*k;
00256   Tstress = Cstress + Ttangent*dStrain;
00257   if (Tstress >= 0.0) {
00258    Tstress = 0.0;
00259    Ttangent = E1n*1.0e-9;
00260   }
00261  }
00262 
00263  else if (Tstrain >= TrotNu && Tstrain < rotch) {
00264   if (Tstrain <= rotrel) {
00265    Tstress = 0.0;
00266    Ttangent = E1n*1.0e-9;
00267   }
00268   else {
00269    Ttangent = maxmom*pinchY/(rotch-rotrel);
00270    tmpmo1 = Cstress + E1n*dStrain;
00271    tmpmo2 = (Tstrain-rotrel)*Ttangent;
00272    if (tmpmo1 < tmpmo2) {
00273     Tstress = tmpmo1;
00274     Ttangent = E1n;
00275    }
00276    else
00277     Tstress = tmpmo2;
00278   }
00279  }
00280 
00281  else {
00282   Ttangent = (1.0-pinchY)*maxmom/(TrotMax-rotch);
00283   tmpmo1 = Cstress + E1n*dStrain;
00284   tmpmo2 = pinchY*maxmom + (Tstrain-rotch)*Ttangent;
00285   if (tmpmo1 < tmpmo2) {
00286    Tstress = tmpmo1;
00287    Ttangent = E1n;
00288   }
00289   else
00290    Tstress = tmpmo2;
00291  }
00292 }
00293 
00294 void
00295 HystereticMaterial::negativeIncrement(double dStrain)
00296 {
00297  double k = pow((CrotMax-1.0e-6)/rot1p,beta);
00298  k = (k < 1.0) ? 1.0 : 1.0/k;
00299 
00300  if (TloadIndicator == 1) {
00301   TloadIndicator = 2;
00302   if (Cstress >= 0.0) {
00303    TrotPu = Cstrain - Cstress/(E1p*k);
00304    double energy = CenergyD - 0.5*Cstress/(E1p*k)*Cstress;
00305    double damfc = 1.0;
00306    if (CrotMax > rot1p) {
00307     damfc += damfc2*energy/energyA;
00308 
00309     if (Cstrain == CrotMax) {
00310      damfc += damfc1*(CrotMin/rot1n-1.0);
00311     }
00312    }
00313 
00314    TrotMin = CrotMin * damfc;
00315   }
00316  }
00317 
00318   TloadIndicator = 2;
00319 
00320  TrotMin = (TrotMin < rot1n) ? TrotMin : rot1n;
00321 
00322  double minmom = negEnvlpStress(TrotMin);
00323  double rotlim = posEnvlpRotlim(CrotMax);
00324  double rotrel = (rotlim < TrotPu) ? rotlim : TrotPu;
00325  rotrel = TrotPu;
00326  if (posEnvlpStress(CrotMax) <= 0.0)
00327    rotrel = rotlim;
00328 
00329  double rotmp1 = rotrel + pinchY*(TrotMin-rotrel);
00330  double rotmp2 = TrotMin - (1.0-pinchY)*minmom/(E1p*k);
00331  double rotch = rotmp1 + (rotmp2-rotmp1)*pinchX;
00332 
00333  double tmpmo1;
00334  double tmpmo2;
00335 
00336  if (Tstrain > TrotPu) {
00337   Ttangent = E1p*k;
00338   Tstress = Cstress + Ttangent*dStrain;
00339   if (Tstress <= 0.0) {
00340    Tstress = 0.0;
00341    Ttangent = E1p*1.0e-9;
00342   }
00343  }
00344 
00345  else if (Tstrain <= TrotPu && Tstrain > rotch) {
00346   if (Tstrain >= rotrel) {
00347    Tstress = 0.0;
00348    Ttangent = E1p*1.0e-9;
00349   }
00350   else {
00351    Ttangent = minmom*pinchY/(rotch-rotrel);
00352    tmpmo1 = Cstress + E1p*dStrain;
00353    tmpmo2 = (Tstrain-rotrel)*Ttangent;
00354    if (tmpmo1 > tmpmo2) {
00355     Tstress = tmpmo1;
00356     Ttangent = E1p;
00357    }
00358    else
00359     Tstress = tmpmo2;
00360   }
00361  }
00362 
00363  else {
00364   Ttangent = (1.0-pinchY)*minmom/(TrotMin-rotch);
00365   tmpmo1 = Cstress + E1p*dStrain;
00366   tmpmo2 = pinchY*minmom + (Tstrain-rotch)*Ttangent;
00367   if (tmpmo1 > tmpmo2) {
00368    Tstress = tmpmo1;
00369    Ttangent = E1p;
00370   }
00371   else
00372    Tstress = tmpmo2;
00373  }
00374 }
00375 
00376 int
00377 HystereticMaterial::commitState(void)
00378 {
00379  CrotMax = TrotMax;
00380  CrotMin = TrotMin;
00381  CrotPu = TrotPu;
00382  CrotNu = TrotNu;
00383  CenergyD = TenergyD;
00384  CloadIndicator = TloadIndicator;
00385 
00386  Cstress = Tstress;
00387  Cstrain = Tstrain;
00388 
00389  return 0;
00390 }
00391 
00392 int
00393 HystereticMaterial::revertToLastCommit(void)
00394 {
00395  TrotMax = CrotMax;
00396  TrotMin = CrotMin;
00397  TrotPu = CrotPu;
00398  TrotNu = CrotNu;
00399  TenergyD = CenergyD;
00400  TloadIndicator = CloadIndicator;
00401 
00402  Tstress = Cstress;
00403  Tstrain = Cstrain;
00404 
00405  return 0;
00406 }
00407 
00408 int
00409 HystereticMaterial::revertToStart(void)
00410 {
00411  CrotMax = 0.0;
00412  CrotMin = 0.0;
00413  CrotPu = 0.0;
00414  CrotNu = 0.0;
00415  CenergyD = 0.0;
00416  CloadIndicator = 0;
00417 
00418  Cstress = 0.0;
00419  Cstrain = 0.0;
00420 
00421  Ttangent = E1p;
00422 
00423  return 0;
00424 }
00425 
00426 UniaxialMaterial*
00427 HystereticMaterial::getCopy(void)
00428 {
00429  HystereticMaterial *theCopy = new HystereticMaterial (this->getTag(),
00430   mom1p, rot1p, mom2p, rot2p, mom3p, rot3p,
00431   mom1n, rot1n, mom2n, rot2n, mom3n, rot3n,
00432   pinchX, pinchY, damfc1, damfc2, beta);
00433 
00434  theCopy->CrotMax = CrotMax;
00435  theCopy->CrotMin = CrotMin;
00436  theCopy->CrotPu = CrotPu;
00437  theCopy->CrotNu = CrotNu;
00438  theCopy->CenergyD = CenergyD;
00439  theCopy->CloadIndicator = CloadIndicator;
00440  theCopy->Cstress = Cstress;
00441  theCopy->Cstrain = Cstrain;
00442 
00443  theCopy->Ttangent = Ttangent;
00444 
00445  return theCopy;
00446 }
00447 
00448 int
00449 HystereticMaterial::sendSelf(int commitTag, Channel &theChannel)
00450 {
00451  return -1;
00452 }
00453 
00454 int
00455 HystereticMaterial::recvSelf(int commitTag, Channel &theChannel, 
00456    FEM_ObjectBroker &theBroker)
00457 {
00458  return -1;
00459 }
00460     
00461 void
00462 HystereticMaterial::Print(ostream &s, int flag)
00463 {
00464  s << "Hysteretic Material, tag: " << this->getTag() << endl;
00465  s << "mom1p: " << mom1p << endl;
00466  s << "rot1p: " << rot1p << endl;
00467  s << "E1p: " << E1p << endl;
00468  s << "mom2p: " << mom2p << endl;
00469  s << "rot2p: " << rot2p << endl;
00470  s << "E2p: " << E2p << endl;
00471  s << "mom3p: " << mom3p << endl;
00472  s << "rot3p: " << rot3p << endl;
00473  s << "E3p: " << E3p << endl;
00474 
00475  s << "mom1n: " << mom1n << endl;
00476  s << "rot1n: " << rot1n << endl;
00477  s << "E1n: " << E1n << endl;
00478  s << "mom2n: " << mom2n << endl;
00479  s << "rot2n: " << rot2n << endl;
00480  s << "E2n: " << E2n << endl;
00481  s << "mom3n: " << mom3n << endl;
00482  s << "rot3n: " << rot3n << endl;
00483  s << "E3n: " << E3n << endl;
00484 
00485  s << "pinchX: " << pinchX << endl;
00486  s << "pinchY: " << pinchY << endl;
00487  s << "damfc1: " << damfc1 << endl;
00488  s << "damfc2: " << damfc2 << endl;
00489  s << "energyA: " << energyA << endl;
00490  s << "beta: " << beta << endl;
00491 }
00492 
00493 void
00494 HystereticMaterial::setEnvelope(void)
00495 {
00496  E1p = mom1p/rot1p;
00497  E2p = (mom2p-mom1p)/(rot2p-rot1p);
00498  E3p = (mom3p-mom2p)/(rot3p-rot2p);
00499 
00500  E1n = mom1n/rot1n;
00501  E2n = (mom2n-mom1n)/(rot2n-rot1n);
00502  E3n = (mom3n-mom2n)/(rot3n-rot2n);
00503 }
00504 
00505 double
00506 HystereticMaterial::posEnvlpStress(double strain)
00507 {
00508  if (strain <= 0.0)
00509   return 0.0;
00510  else if (strain <= rot1p)
00511   return E1p*strain;
00512  else if (strain <= rot2p)
00513   return mom1p + E2p*(strain-rot1p);
00514  else if (strain <= rot3p || E3p > 0.0)
00515   return mom2p + E3p*(strain-rot2p);
00516  else
00517   return mom3p;
00518 }
00519 
00520 double
00521 HystereticMaterial::negEnvlpStress(double strain)
00522 {
00523  if (strain >= 0.0)
00524   return 0.0;
00525  else if (strain >= rot1n)
00526   return E1n*strain;
00527  else if (strain >= rot2n)
00528   return mom1n + E2n*(strain-rot1n);
00529  else if (strain >= rot3n || E3n > 0.0)
00530   return mom2n + E3n*(strain-rot2n);
00531  else
00532   return mom3n;
00533 }
00534 
00535 double
00536 HystereticMaterial::posEnvlpTangent(double strain)
00537 {
00538  if (strain < 0.0)
00539   return E1p*1.0e-9;
00540  else if (strain <= rot1p)
00541   return E1p;
00542  else if (strain <= rot2p)
00543   return E2p;
00544  else if (strain <= rot3p || E3p > 0.0)
00545   return E3p;
00546  else
00547   return E1p*1.0e-9;
00548 }
00549 
00550 double
00551 HystereticMaterial::negEnvlpTangent(double strain)
00552 {
00553  if (strain > 0.0)
00554   return E1n*1.0e-9;
00555  else if (strain >= rot1n)
00556   return E1n;
00557  else if (strain >= rot2n)
00558   return E2n;
00559  else if (strain >= rot3n || E3n > 0.0)
00560   return E3n;
00561  else
00562   return E1n*1.0e-9;
00563 }
00564 
00565 double
00566 HystereticMaterial::posEnvlpRotlim(double strain)
00567 {
00568   double strainLimit = POS_INF_STRAIN;
00569 
00570   if (strain <= rot1p)
00571     return POS_INF_STRAIN;
00572   if (strain > rot1p && strain <= rot2p && E2p < 0.0)
00573     strainLimit = rot1p - mom1p/E2p;
00574   if (strain > rot2p && E3p < 0.0)
00575     strainLimit = rot2p - mom2p/E3p;
00576 
00577   if (strainLimit == POS_INF_STRAIN)
00578     return POS_INF_STRAIN;
00579   else if (posEnvlpStress(strainLimit) > 0)
00580     return POS_INF_STRAIN;
00581   else
00582     return strainLimit;
00583 }
00584 
00585 double
00586 HystereticMaterial::negEnvlpRotlim(double strain)
00587 {
00588   double strainLimit = NEG_INF_STRAIN;
00589 
00590   if (strain >= rot1n)
00591     return NEG_INF_STRAIN;
00592   if (strain < rot1n && strain >= rot2n && E2n < 0.0)
00593     strainLimit = rot1n - mom1n/E2n;
00594   if (strain < rot2n && E3n < 0.0)
00595     strainLimit = rot2n - mom2n/E3n;
00596 
00597   if (strainLimit == NEG_INF_STRAIN)
00598     return NEG_INF_STRAIN;
00599   else if (negEnvlpStress(strainLimit) < 0)
00600     return NEG_INF_STRAIN;
00601   else
00602     return strainLimit;
00603 }