BarSlipMaterial.cppGo to the documentation of this file.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.2 $ 00022 // $Date: 2004/10/06 19:21:12 $ 00023 // $Source: /usr/local/cvs/OpenSees/SRC/material/uniaxial/BarSlipMaterial.cpp,v $ 00024 00025 00026 // Written: NM (nmitra@u.washington.edu) 00027 // Created: January 2002 00028 // Updated: September 2004 00029 // 00030 // Description: This file contains the class implementation for 00031 // bar-slip material which is defined by 4 points on the positive and 00032 // negative envelopes and a bunch of damage parameters. The material accounts for 00033 // 3 types of damage rules : Strength degradation, Stiffness degradation, 00034 // unloading stiffness degradation. 00035 // The theory for the material has been provided by Prof. Lowes 00036 // Updates: new damage calculations and user interfaces 00037 00038 00039 #include <BarSlipMaterial.h> 00040 #include <math.h> 00041 #include <float.h> 00042 #include <OPS_Globals.h> 00043 00044 BarSlipMaterial::BarSlipMaterial(int tag, 00045 double f, double fs, double es, double fsu, 00046 double eh, double dbar, double ljoint, int n, double w, double d, 00047 int bsf, int typ): 00048 UniaxialMaterial(tag, MAT_TAG_BarSlip), 00049 tagMat(tag), bsflag(bsf), unit(0), type(typ), damage(1), width(w), depth(d), 00050 envlpPosStress(6), envlpPosStrain(6), envlpNegStress(6), envlpNegStrain(6), 00051 fc(f), fy(fs), Es(es), fu(fsu), Eh(eh), db(dbar), nbars(n), ld(ljoint), 00052 eP(4,2), eN(4,2), envlpPosDamgdStress(6), envlpNegDamgdStress(6), state3Stress(4), 00053 state3Strain(4), state4Stress(4), state4Strain(4) 00054 { 00055 rDispP = 0.25; rForceP = 0.25; uForceP = 0.0; 00056 rDispN = 0.25; rForceN = 0.25; uForceN = 0.0; 00057 gammaK1 = 0.3; gammaK2 = 0.0; gammaK3 = 0.1; gammaK4 = 0.0; gammaKLimit = 0.4; 00058 gammaD1 = 0.6; gammaD2 = 0.0; gammaD3 = 0.2; gammaD4 = 0.0; gammaDLimit = 0.25; 00059 gammaF1 = 0.7; gammaF2 = 0.3; gammaF3 = 0.5; gammaF4 = 0.1; gammaFLimit = 0.0; 00060 gammaE = 10.0; 00061 00062 getBondStrength(); 00063 getBarSlipEnvelope(); 00064 createMaterial(); 00065 } 00066 00067 BarSlipMaterial::BarSlipMaterial(int tag, 00068 double f, double fs, double es, double fsu, 00069 double eh, double dbar, double ljoint, int n, double w, double d, 00070 int bsf, int typ, int dam, int unt): 00071 UniaxialMaterial(tag, MAT_TAG_BarSlip), 00072 tagMat(tag), bsflag(bsf), unit(unt), type(typ), damage(dam), width(w), depth(d), 00073 envlpPosStress(6), envlpPosStrain(6), envlpNegStress(6), envlpNegStrain(6), 00074 fc(f), fy(fs), Es(es), fu(fsu), Eh(eh), db(dbar), nbars(n), ld(ljoint), 00075 eP(4,2), eN(4,2), envlpPosDamgdStress(6), envlpNegDamgdStress(6), state3Stress(4), 00076 state3Strain(4), state4Stress(4), state4Strain(4) 00077 { 00078 rDispP = 0.25; rForceP = 0.25; uForceP = 0.0; 00079 rDispN = 0.25; rForceN = 0.25; uForceN = 0.0; 00080 gammaK1 = 0.3; gammaK2 = 0.0; gammaK3 = 0.1; gammaK4 = 0.0; gammaKLimit = 0.4; 00081 gammaD1 = 0.6; gammaD2 = 0.0; gammaD3 = 0.2; gammaD4 = 0.0; gammaDLimit = 0.25; 00082 gammaF1 = 0.7; gammaF2 = 0.3; gammaF3 = 0.5; gammaF4 = 0.1; gammaFLimit = 0.0; 00083 gammaE = 10.0; 00084 00085 00086 if (damage == 0) 00087 { 00088 gammaK1 = 0.0; gammaK2 = 0.0; gammaK3 = 0.0; gammaK4 = 0.0; gammaKLimit = 0.0; 00089 gammaD1 = 0.0; gammaD2 = 0.0; gammaD3 = 0.0; gammaD4 = 0.0; gammaDLimit = 0.0; 00090 gammaF1 = 0.0; gammaF2 = 0.0; gammaF3 = 0.0; gammaF4 = 0.0; gammaFLimit = 0.0; 00091 } 00092 if (damage == 1) 00093 { 00094 // activation of linear degradation of strength 00095 gammaF1 = 0.0; gammaF2 = 0.0; gammaF3 = 0.0; gammaF4 = 0.0; gammaFLimit = 0.0; 00096 } 00097 if (damage == 2) 00098 { 00099 // activation of logarithmic degradation of strength 00100 gammaF1 = 11.8986; gammaF2 = 0.0; gammaF3 = 3.9694; gammaF4 = 0.0; gammaFLimit = 0.85; 00101 } 00102 00103 getBondStrength(); 00104 getBarSlipEnvelope(); 00105 createMaterial(); 00106 } 00107 00108 00109 BarSlipMaterial::BarSlipMaterial(): 00110 UniaxialMaterial(0, MAT_TAG_BarSlip), 00111 tagMat(0), bsflag(0), unit(0), type(0), damage(0), width(0.0), depth(0.0), 00112 envlpPosStress(6), envlpPosStrain(6), envlpNegStress(6), envlpNegStrain(6), 00113 fc(0.0), fy(0.0), Es(0.0), fu(0.0), Eh(0.0), db(0.0), nbars(0), ld(0.0), 00114 eP(4,2), eN(4,2), envlpPosDamgdStress(6), envlpNegDamgdStress(6), state3Stress(4), 00115 state3Strain(4), state4Stress(4), state4Strain(4) 00116 { 00117 // default constructor --- does nothing 00118 } 00119 00120 BarSlipMaterial::~BarSlipMaterial() 00121 { 00122 // destructor 00123 // fn->close(); 00124 // fg->close(); 00125 00126 } 00127 00128 void BarSlipMaterial::getBondStrength(void) 00129 { 00130 if (fc <= 0.0) 00131 opserr << "WARNING : BAR-SLIP -- fc should be positive entry" << endln; 00132 00133 if (unit == 0) // unit has not been specified : either psi or MPa 00134 { 00135 if (fc >= 1000.0) 00136 { 00137 unit = 2; // psi unit being used 00138 } 00139 else 00140 { 00141 unit = 1; // MPa unit being used 00142 } 00143 } 00144 00145 if (unit == 1) // "mpa" (N/sq. mm) system being used 00146 { 00147 if (bsflag == 1) // consider bond strength flag --- weak 00148 { 00149 tauYT = 0.05*pow(fc,0.5); 00150 tauET = 1.8*pow(fc,0.5); 00151 tauEC = 2.2*pow(fc,0.5); 00152 tauYC = 3.7*pow(fc,0.5); 00153 tauR = 0.15*pow(fc,0.5); 00154 } 00155 else if (bsflag == 0) // bond strength -- strong 00156 { 00157 tauYT = 0.4*pow(fc,0.5); 00158 tauET = 1.8*pow(fc,0.5); 00159 tauEC = 2.2*pow(fc,0.5); 00160 tauYC = 3.7*pow(fc,0.5); 00161 tauR = 0.15*pow(fc,0.5); 00162 } 00163 } 00164 else if (unit == 2) // "psi" (lb/sq. inch) system used 00165 { 00166 if (bsflag == 1) // bond strength -- weak 00167 { 00168 tauYT = 0.6*pow(fc,0.5); 00169 tauET = 21*pow(fc,0.5); 00170 tauEC = 26*pow(fc,0.5); 00171 tauYC = 43*pow(fc,0.5); 00172 tauR = 1.8*pow(fc,0.5); 00173 } 00174 else if (bsflag == 0) // bond strength -- strong 00175 { 00176 tauYT = 4.8*pow(fc,0.5); 00177 tauET = 21*pow(fc,0.5); 00178 tauEC = 26*pow(fc,0.5); 00179 tauYC = 43*pow(fc,0.5); 00180 tauR = 1.8*pow(fc,0.5); 00181 } 00182 } 00183 else if (unit == 3) // "Pa" (N/sq. m) system being used 00184 { 00185 if (bsflag == 1) // bond strength -- weak 00186 { 00187 tauYT = 1.58*pow(fc,0.5); 00188 tauET = 56.92*pow(fc,0.5); 00189 tauEC = 69.57*pow(fc,0.5); 00190 tauYC = 117*pow(fc,0.5); 00191 tauR = 4.74*pow(fc,0.5); 00192 } 00193 else if (bsflag == 0) // bond strength -- strong 00194 { 00195 tauYT = 12.65*pow(fc,0.5); 00196 tauET = 56.92*pow(fc,0.5); 00197 tauEC = 69.57*pow(fc,0.5); 00198 tauYC = 117*pow(fc,0.5); 00199 tauR = 4.74*pow(fc,0.5); 00200 } 00201 } 00202 else if (unit == 4) // "psf" (lb/sq.ft) system being used 00203 { 00204 if (bsflag == 1) // bond strength -- weak 00205 { 00206 tauYT = 7.2*pow(fc,0.5); 00207 tauET = 252*pow(fc,0.5); 00208 tauEC = 312*pow(fc,0.5); 00209 tauYC = 516*pow(fc,0.5); 00210 tauR = 21.6*pow(fc,0.5); 00211 } 00212 else if (bsflag == 0) // bond strength -- strong 00213 { 00214 tauYT = 57.6*pow(fc,0.5); 00215 tauET = 252*pow(fc,0.5); 00216 tauEC = 312*pow(fc,0.5); 00217 tauYC = 516*pow(fc,0.5); 00218 tauR = 21.6*pow(fc,0.5); 00219 } 00220 } 00221 else if (unit == 5) // "ksi" (kilo-lb/sq.inch) system being used 00222 { 00223 if (bsflag == 1) // bond strength -- weak 00224 { 00225 tauYT = 0.02*pow(fc,0.5); 00226 tauET = 0.66*pow(fc,0.5); 00227 tauEC = 0.82*pow(fc,0.5); 00228 tauYC = 1.36*pow(fc,0.5); 00229 tauR = 0.06*pow(fc,0.5); 00230 } 00231 else if (bsflag == 0) // bond strength -- strong 00232 { 00233 tauYT = 0.15*pow(fc,0.5); 00234 tauET = 0.66*pow(fc,0.5); 00235 tauEC = 0.82*pow(fc,0.5); 00236 tauYC = 1.36*pow(fc,0.5); 00237 tauR = 0.06*pow(fc,0.5); 00238 } 00239 } 00240 else if (unit == 6) // "ksf" (kilo-lb/sq.ft) system being used 00241 { 00242 if (bsflag == 1) // bond strength -- weak 00243 { 00244 tauYT = 0.24*pow(fc,0.5); 00245 tauET = 7.92*pow(fc,0.5); 00246 tauEC = 9.84*pow(fc,0.5); 00247 tauYC = 16.32*pow(fc,0.5); 00248 tauR = 0.72*pow(fc,0.5); 00249 } 00250 else if (bsflag == 0) // bond strength -- strong 00251 { 00252 tauYT = 1.8*pow(fc,0.5); 00253 tauET = 7.92*pow(fc,0.5); 00254 tauEC = 9.84*pow(fc,0.5); 00255 tauYC = 16.32*pow(fc,0.5); 00256 tauR = 0.72*pow(fc,0.5); 00257 } 00258 } 00259 } 00260 00261 void BarSlipMaterial::getBarSlipEnvelope(void) 00262 { 00263 //fn = new FileStream(); 00264 //fn->setFile("BarSlipEnvelope.out", APPEND); 00265 //fg = new FileStream(); 00266 //fg->setFile("BarSlipDamage.out", APPEND); 00267 00268 double delta = 0.0; 00269 double del_ult = 0.0; 00270 if (unit == 1) // MPa 00271 { 00272 delta = 3.0; 00273 del_ult = 10.0; 00274 } 00275 else if (unit == 2 || unit == 5) //psi or ksi 00276 { 00277 delta = 3.0/25.43; 00278 del_ult = 10.0/25.43; 00279 } 00280 else if (unit == 3) // Pa 00281 { 00282 delta = 0.003; 00283 del_ult = 0.01; 00284 } 00285 else if (unit == 4 || unit == 6) // psf or ksf 00286 { 00287 delta = 3.0/(25.43*12.0); 00288 del_ult = 10.0/(25.43*12.0); 00289 } 00290 00291 double Ab = PI*pow(db,2)/4; 00292 double As = Ab*nbars; 00293 00294 double frP = 0.0; 00295 double frN = 0.0; 00296 double geo = 0.0; 00297 double let = 0.0; 00298 double lec = 0.0; 00299 double lyc = 0.0; 00300 double lyt = 0.0; 00301 double k1 = 0.0; 00302 double k2 = 0.0; 00303 double k3 = 0.0; 00304 double le1 = 0.0; 00305 double le2 = 0.0; 00306 double unloadNForce = 0.0; 00307 double unloadPForce = 0.0; 00308 eP.Zero(); 00309 eN.Zero(); 00310 00311 frP = tauR*ld*PI*db*As/Ab; 00312 frN = -tauR*ld*PI*db*As/Ab; 00313 00314 geo = PI*db/Ab; 00315 let = fy/(tauET*geo); // length of elastic part in tension 00316 lyt = (fu-fy)/(tauYT*geo); // length of yielded part in tension 00317 lec = fy/(tauEC*geo); // length of elastic part in compression 00318 lyc = (fu-fy)/(tauYC*geo); // length of yielded part in compression 00319 00320 k1 = 2*Es*(tauET/fy)*geo*As; // initial slope 00321 00322 eP(0,0) = 0.5*fy*As/k1; 00323 eP(0,1) = 0.5*fy*As; 00324 eP(1,0) = fy*As/k1; 00325 eP(1,1) = fy*As; 00326 00327 if (let+lyt<ld && bsflag ==0 ) 00328 { 00329 k2 = (fu-fy)*As/(fy*lyt/Es + 0.5*geo*tauYT*pow(lyt,2)/Eh); // second slope after yield 00330 } 00331 else 00332 { 00333 le1 = fy/(tauET*geo); 00334 le2 = fy/(tauYT*geo); 00335 k2 = (fu-fy)*As/(geo*0.5*tauYT*(pow(le2,2)/Es - pow(le1,2)/Es + pow(lyt,2)/Eh) + fy*lyt/Es); 00336 } 00337 00338 00339 eP(2,0) = (delta<(fy*As/k1 + (fu-fy)*As/k2))? delta:(fy*As/k1 + (fu-fy)*As/k2); 00340 00341 if (eP(2,0) == delta) 00342 { 00343 eP(2,1) = fy*As + (delta - fy*As/k1)*k2; 00344 } 00345 else 00346 { 00347 eP(2,1) = fu*As; 00348 } 00349 00350 eP(3,0) = del_ult; 00351 eP(3,1) = eP(2,1) + (eP(2,1) - eP(1,1))*(eP(3,0) - eP(2,0))/(eP(2,0) - eP(1,0)); 00352 00353 gammaFLimit = 1.0 - frP/eP(2,1); 00354 00355 double dth = depth; 00356 double dd = 0.1*depth; 00357 double j = 0.0; 00358 double beta = 0.85; 00359 double fcc = 0.0; 00360 00361 if (unit == 2) 00362 { 00363 fcc = fc; 00364 } 00365 else if (unit == 1) 00366 { 00367 fcc = 145*fc; 00368 } 00369 else if (unit == 3) 00370 { 00371 fcc = 0.000145*fc; 00372 } 00373 else if (unit == 4) 00374 { 00375 fcc = 0.00694*fc; 00376 } 00377 else if (unit == 5) 00378 { 00379 fcc = 1000*fc; 00380 } 00381 else if (unit == 6) 00382 { 00383 fcc = 6.94*fc; 00384 } 00385 00386 double b1 = (fcc - 4000)*0.05/1000; 00387 if (b1 <= 0.0) 00388 { 00389 b1 = 0.0; 00390 } 00391 if (b1 >= 0.2) 00392 { 00393 b1 = 0.2; 00394 } 00395 00396 double beta1 = 0.85 - b1; 00397 00398 double cForce = 0.0; 00399 00400 if ((type == 0) || (type == 1)) // member type -- beam bottom and beam top (or just beam) 00401 { 00402 j = 0.85; 00403 } 00404 else if (type ==2) // member type -- column 00405 { 00406 j = 0.75; 00407 } 00408 00409 cForce = 1 + (beta*fc*dth*width*2*(1-j)/(Es*As*0.003*beta1*(1-dd*beta1/(2*dth*(1-j))))); 00410 As = As*cForce; 00411 k1 = 2*Es*(tauEC/fy)*geo*As; 00412 00413 00414 eN(0,0) = -.5*fy*As/k1; 00415 eN(0,1) = -.5*fy*As; 00416 eN(1,0) = -fy*As/k1; 00417 eN(1,1) = -fy*As; 00418 00419 if (lec+lyc < ld && bsflag ==0 ) 00420 { 00421 k2 = (fu-fy)*As/(fy*lyc/Es + 0.5*geo*tauYC*pow(lyc,2)/Eh); 00422 } 00423 else 00424 { 00425 le1 = fy/(tauEC*geo); 00426 le2 = fy/(tauYC*geo); 00427 k2 = (fu-fy)*As/(geo*0.5*tauYC*(pow(le2,2)/Es - pow(le1,2)/Es + pow(lyc,2)/Eh) + fy*lyc/Es); 00428 } 00429 00430 eN(2,0) = (delta<(fy*As/k1 + (fu-fy)*As/k2))? -delta:-(fy*As/k1 + (fu-fy)*As/k2); 00431 00432 if (eN(2,0) == -delta) 00433 { 00434 eN(2,1) = -fy*As + (-delta + fy*As/k1)*k2; 00435 } 00436 else 00437 { 00438 eN(2,1) = -fu*As; 00439 } 00440 00441 k3 = 1e-3*k1; 00442 eN(3,0) = -del_ult; 00443 eN(3,1) = eN(2,1) + k3*(eN(3,0)-eN(2,0)); 00444 00445 As = As/cForce; 00446 00447 double tP = (ld<(lec+lyc)) ? ld:(lec+lyc); 00448 double tN = (ld<(let+lyt)) ? ld:(let+lyt); 00449 00450 unloadPForce = tauR*PI*db*As*tP/Ab; 00451 unloadNForce = -tauR*PI*db*As*tN/Ab; 00452 00453 uForceP = unloadPForce/eP(2,1); 00454 uForceN = unloadNForce/eN(2,1); 00455 00456 double ratio = 2.0; //(research underway to induce more pinching in the material) 00457 rForceP = (uForceP>uForceN) ? ratio*uForceP:ratio*uForceN; 00458 rForceN = (uForceP>uForceN) ? ratio*uForceP:ratio*uForceN; 00459 // rForceP = 0.25; rForceN = 0.25; 00460 00461 } 00462 00463 void BarSlipMaterial::createMaterial(void) 00464 { 00465 /* material = new Pinching4Material (tagMat, eP(0,1), eP(0,0), eP(1,1), eP(1,0), eP(2,1), eP(2,0), eP(3,1), eP(3,0), 00466 eN(0,1), eN(0,0), eN(1,1), eN(1,0), eN(2,1), eN(2,0), eN(3,1), eN(3,0), rDispP, rForceP, uForceP, rDispN, 00467 rForceN, uForceN, gammaK1, gammaK2, gammaK3, gammaK4, gammaKLimit, gammaD1, gammaD2, gammaD3, 00468 gammaD4, gammaDLimit, gammaF1, gammaF2, gammaF3, gammaF4, gammaFLimit, gammaE, 0); 00469 */ 00470 bool error = false; 00471 00472 // check backbone parameters formed 00473 if (eP(0,0) <= 0.0 || eP(1,0) <= 0.0 || eP(2,0) <= 0.0 || eP(3,0) <= 0.0) 00474 error = true; 00475 00476 if (eN(0,0) >= 0.0 || eN(1,0) >= 0.0 || eN(2,0) >= 0.0 || eN(3,0) >= 0.0) 00477 error = true; 00478 00479 if (error) { 00480 opserr << "Error: -- input backbone not unique, BarSlipMaterial::BarSlipMaterial" << "\a"; 00481 } 00482 00483 envlpPosStress.Zero(); envlpPosStrain.Zero(); envlpNegStress.Zero(); envlpNegStrain.Zero(); 00484 energyCapacity = 0.0; kunload = 0.0; elasticStrainEnergy = 0.0; 00485 00486 // set envelope slopes 00487 SetEnvelope(); 00488 00489 envlpPosDamgdStress = envlpPosStress; envlpNegDamgdStress = envlpNegStress; 00490 state3Stress.Zero(); state3Strain.Zero(); state4Stress.Zero(); state4Strain.Zero(); 00491 // Initialize history variables 00492 revertToStart(); 00493 revertToLastCommit(); 00494 00495 00496 } 00497 //******************************************************************************* 00498 /* 00499 int BarSlipMaterial::setTrialStrain(double strain, double strainrate) 00500 { 00501 00502 return material->setTrialStrain(strain, strainrate); 00503 } 00504 00505 double BarSlipMaterial::getStrain(void) 00506 { 00507 return material->getStrain(); 00508 } 00509 00510 double BarSlipMaterial::getStress(void) 00511 { 00512 return material->getStress(); 00513 } 00514 00515 double BarSlipMaterial::getTangent(void) 00516 { 00517 return material->getTangent(); 00518 } 00519 00520 double BarSlipMaterial::getInitialTangent(void) 00521 { 00522 return material->getInitialTangent(); 00523 } 00524 int BarSlipMaterial::commitState(void) 00525 { 00526 return material->commitState(); 00527 } 00528 00529 int BarSlipMaterial::revertToLastCommit(void) 00530 { 00531 00532 return material->revertToLastCommit(); 00533 } 00534 00535 int BarSlipMaterial::revertToStart(void) 00536 { 00537 return material->revertToStart(); 00538 } 00539 */ 00540 UniaxialMaterial* BarSlipMaterial::getCopy(void) 00541 { 00542 BarSlipMaterial *theCopy = new BarSlipMaterial (this->getTag(), fc, fy, Es, 00543 fu, Eh, db, ld, nbars, width, depth, bsflag, type, damage, unit); 00544 00545 theCopy->rDispN = rDispN; 00546 theCopy->rDispP = rDispP; 00547 theCopy->rForceN = rForceN; 00548 theCopy->rForceP = rForceP; 00549 theCopy->uForceN = uForceN; 00550 theCopy->uForceP = uForceP; 00551 00552 // Trial state variables 00553 theCopy->Tstress = Tstress; 00554 theCopy->Tstrain = Tstrain; 00555 theCopy->Ttangent = Ttangent; 00556 00557 // Coverged material history parameters 00558 theCopy->Cstate = Cstate; 00559 theCopy->Cstrain = Cstrain; 00560 theCopy->Cstress = Cstress; 00561 theCopy->CstrainRate = CstrainRate; 00562 00563 theCopy->lowCstateStrain = lowCstateStrain; 00564 theCopy->lowCstateStress = lowCstateStress; 00565 theCopy->hghCstateStrain = hghCstateStrain; 00566 theCopy->hghCstateStress = hghCstateStress; 00567 theCopy->CminStrainDmnd = CminStrainDmnd; 00568 theCopy->CmaxStrainDmnd = CmaxStrainDmnd; 00569 theCopy->Cenergy = Cenergy; 00570 theCopy->CgammaK = CgammaK; 00571 theCopy->CgammaD = CgammaD; 00572 theCopy->CgammaF = CgammaF; 00573 theCopy->gammaKUsed = gammaKUsed; 00574 theCopy->gammaFUsed = gammaFUsed; 00575 00576 // trial material history parameters 00577 theCopy->Tstate = Tstate; 00578 theCopy->dstrain = dstrain; 00579 theCopy->lowTstateStrain = lowTstateStrain; 00580 theCopy->lowTstateStress = lowTstateStress; 00581 theCopy->hghTstateStrain = hghTstateStrain; 00582 theCopy->hghTstateStress = hghTstateStress; 00583 theCopy->TminStrainDmnd = TminStrainDmnd; 00584 theCopy->TmaxStrainDmnd = TmaxStrainDmnd; 00585 theCopy->Tenergy = Tenergy; 00586 theCopy->TgammaK = TgammaK; 00587 theCopy->TgammaD = TgammaD; 00588 theCopy->TgammaF = TgammaF; 00589 00590 // Strength and stiffness parameters 00591 theCopy->kElasticPos = kElasticPos; 00592 theCopy->kElasticNeg = kElasticNeg; 00593 theCopy->kElasticPosDamgd = kElasticPosDamgd; 00594 theCopy->kElasticNegDamgd = kElasticNegDamgd; 00595 theCopy->uMaxDamgd = uMaxDamgd; 00596 theCopy->uMinDamgd = uMinDamgd; 00597 00598 for (int i = 0; i<6; i++) 00599 { 00600 theCopy->envlpPosStrain(i) = envlpPosStrain(i); 00601 theCopy->envlpPosStress(i) = envlpPosStress(i); 00602 theCopy->envlpNegStrain(i) = envlpNegStrain(i); 00603 theCopy->envlpNegStress(i) = envlpNegStress(i); 00604 theCopy->envlpNegDamgdStress(i) = envlpNegDamgdStress(i); 00605 theCopy->envlpPosDamgdStress(i) = envlpPosDamgdStress(i); 00606 } 00607 00608 for (int j = 0; j<4; j++) 00609 { 00610 theCopy->state3Strain(j) = state3Strain(j); 00611 theCopy->state3Stress(j) = state3Stress(j); 00612 theCopy->state4Strain(j) = state4Strain(j); 00613 theCopy->state4Stress(j) = state4Stress(j); 00614 } 00615 00616 theCopy->energyCapacity = energyCapacity; 00617 theCopy->kunload = kunload; 00618 theCopy->elasticStrainEnergy = elasticStrainEnergy; 00619 00620 return theCopy; 00621 } 00622 00623 int BarSlipMaterial::sendSelf(int commitTag, Channel &theChannel) 00624 { 00625 return -1; 00626 } 00627 00628 int BarSlipMaterial::recvSelf(int commitTag, Channel &theChannel, 00629 FEM_ObjectBroker & theBroker) 00630 { 00631 return -1; 00632 } 00633 00634 void BarSlipMaterial::Print(OPS_Stream &s, int flag) 00635 { 00636 s << "BarSlipMaterial, tag: " << this-> getTag() << endln; 00637 s << "Positive Envelope : " << eP << endln; 00638 s << "Negative Envelope : " << eN << endln; 00639 } 00640 00641 //********************************************************************************** 00642 void BarSlipMaterial::SetEnvelope(void) 00643 { 00644 double kPos = eP(0,1)/eP(0,0); 00645 double kNeg = eN(0,1)/eN(0,0); 00646 00647 double k = (kPos>kNeg) ? kPos:kNeg; 00648 double u = (eP(0,0)>-eN(0,0)) ? 1e-4*eP(0,0):-1e-4*eN(0,0); 00649 00650 envlpPosStrain(0) = u; envlpPosStress(0) = u*k; 00651 envlpNegStrain(0) = -u; envlpNegStress(0) = -u*k; 00652 00653 for (int i1 = 1; i1 < 5; i1++) { 00654 envlpPosStrain(i1) = eP(i1-1,0); 00655 envlpPosStress(i1) = eP(i1-1,1); 00656 envlpNegStrain(i1) = eN(i1-1,0); 00657 envlpNegStress(i1) = eN(i1-1,1); 00658 } 00659 00660 double k1 = (eP(3,1) - eP(2,1))/(eP(3,0) - eP(2,0)); 00661 double k2 = (eN(3,1) - eN(2,1))/(eN(3,0) - eN(2,0)); 00662 00663 envlpPosStrain(5) = 1e+6*eP(3,0); envlpNegStrain(5) = 1e+6*eN(3,0); 00664 envlpPosStress(5) = (k1>0) ? eP(3,1) + k1*(envlpPosStrain(5)-envlpPosStrain(4)):envlpPosStress(4)*1.1; 00665 envlpNegStress(5) = (k2>0) ? eN(3,1) + k2*(envlpNegStrain(5)-envlpNegStrain(4)):envlpNegStress(4)*1.1; 00666 00667 kElasticPos = envlpPosStress(1)/envlpPosStrain(1); 00668 kElasticNeg = envlpNegStress(1)/envlpNegStrain(1); 00669 00670 double energyPos = 0.5*envlpPosStrain(0)*envlpPosStress(0); 00671 double energyNeg = 0.5*envlpNegStrain(0)*envlpNegStress(0); 00672 00673 for (int i2 = 0; i2<4; i2++) { 00674 energyPos += 0.5*(envlpPosStress(i2)+envlpPosStress(i2+1))*(envlpPosStrain(i2+1) - envlpPosStrain(i2)); 00675 } 00676 00677 for (int i3 = 0; i3<4; i3++) { 00678 energyNeg += 0.5*(envlpNegStress(i3)+envlpNegStress(i3+1))*(envlpNegStrain(i3+1) - envlpNegStrain(i3)); 00679 } 00680 00681 double max_energy = (energyPos>energyNeg) ? energyPos:energyNeg; 00682 energyCapacity = gammaE*max_energy; 00683 } 00684 00685 int BarSlipMaterial::setTrialStrain(double strain, double CstrainRate) 00686 { 00687 Tstate = Cstate; 00688 Tenergy = Cenergy; 00689 Tstrain = strain; 00690 lowTstateStrain = lowCstateStrain; 00691 hghTstateStrain = hghCstateStrain; 00692 lowTstateStress = lowCstateStress; 00693 hghTstateStress = hghCstateStress; 00694 TminStrainDmnd = CminStrainDmnd; 00695 TmaxStrainDmnd = CmaxStrainDmnd; 00696 TgammaF = CgammaF; 00697 TgammaK = CgammaK; 00698 TgammaD = CgammaD; 00699 00700 dstrain = Tstrain - Cstrain; 00701 if (dstrain<1e-12 && dstrain>-1e-12){ 00702 dstrain = 0.0; 00703 } 00704 00705 00706 // determine new state if there is a change in state 00707 getstate(Tstrain,dstrain); 00708 00709 switch (Tstate) 00710 { 00711 00712 case 0: 00713 Ttangent = envlpPosStress(0)/envlpPosStrain(0); 00714 Tstress = Ttangent*Tstrain; 00715 break; 00716 case 1: 00717 Tstress = posEnvlpStress(strain); 00718 Ttangent = posEnvlpTangent(strain); 00719 break; 00720 case 2: 00721 Ttangent = negEnvlpTangent(strain); 00722 Tstress = negEnvlpStress(strain); 00723 break; 00724 case 3: 00725 kunload = (hghTstateStrain<0.0) ? kElasticNegDamgd:kElasticPosDamgd; 00726 state3Strain(0) = lowTstateStrain; 00727 state3Strain(3) = hghTstateStrain; 00728 state3Stress(0) = lowTstateStress; 00729 state3Stress(3) = hghTstateStress; 00730 00731 getState3(state3Strain,state3Stress,kunload); 00732 Ttangent = Envlp3Tangent(state3Strain,state3Stress,strain); 00733 Tstress = Envlp3Stress(state3Strain,state3Stress,strain); 00734 00735 //Print(opserr,1); 00736 break; 00737 case 4: 00738 kunload = (lowTstateStrain<0.0) ? kElasticNegDamgd:kElasticPosDamgd; 00739 state4Strain(0) = lowTstateStrain; 00740 state4Strain(3) = hghTstateStrain; 00741 state4Stress(0) = lowTstateStress; 00742 state4Stress(3) = hghTstateStress; 00743 00744 getState4(state4Strain,state4Stress,kunload); 00745 Ttangent = Envlp4Tangent(state4Strain,state4Stress,strain); 00746 Tstress = Envlp4Stress(state4Strain,state4Stress,strain); 00747 break; 00748 } 00749 00750 double denergy = 0.5*(Tstress+Cstress)*dstrain; 00751 elasticStrainEnergy = (Tstrain>0.0) ? 0.5*Tstress/kElasticPosDamgd*Tstress:0.5*Tstress/kElasticNegDamgd*Tstress; 00752 00753 Tenergy = Cenergy + denergy; 00754 00755 updateDmg(Tstrain); 00756 return 0; 00757 } 00758 00759 double BarSlipMaterial::getStrain(void) 00760 { 00761 return Tstrain; 00762 } 00763 00764 double BarSlipMaterial::getStress(void) 00765 { 00766 return Tstress; 00767 } 00768 00769 double BarSlipMaterial::getTangent(void) 00770 { 00771 return Ttangent; 00772 } 00773 00774 double BarSlipMaterial::getInitialTangent(void) 00775 { 00776 return envlpPosStress(0)/envlpPosStrain(0); 00777 } 00778 00779 int BarSlipMaterial::commitState(void) 00780 { 00781 Cstate = Tstate; 00782 00783 if (dstrain>1e-12||dstrain<-(1e-12)) { 00784 CstrainRate = dstrain;} 00785 else { 00786 CstrainRate = TstrainRate;} 00787 00788 lowCstateStrain = lowTstateStrain; 00789 lowCstateStress = lowTstateStress; 00790 hghCstateStrain = hghTstateStrain; 00791 hghCstateStress = hghTstateStress; 00792 CminStrainDmnd = TminStrainDmnd; 00793 CmaxStrainDmnd = TmaxStrainDmnd; 00794 Cenergy = Tenergy; 00795 00796 Cstress = Tstress; 00797 Cstrain = Tstrain; 00798 00799 CgammaK = TgammaK; 00800 CgammaD = TgammaD; 00801 CgammaF = TgammaF; 00802 00803 // define adjusted strength and stiffness parameters 00804 kElasticPosDamgd = kElasticPos*(1 - gammaKUsed); 00805 kElasticNegDamgd = kElasticNeg*(1 - gammaKUsed); 00806 00807 uMaxDamgd = TmaxStrainDmnd*(1 + CgammaD); 00808 uMinDamgd = TminStrainDmnd*(1 + CgammaD); 00809 00810 envlpPosDamgdStress = envlpPosStress*(1-gammaFUsed); 00811 envlpNegDamgdStress = envlpNegStress*(1-gammaFUsed); 00812 00813 //(*fn) << tagMat << " " << envlpNegStrain(4) << " " << envlpNegStrain(3) << " " << envlpNegStrain(2) << " " << envlpNegStrain(1) << " " << envlpPosStrain(1) << " " << envlpPosStrain(2) << " " << envlpPosStrain(3) << " " << envlpPosStrain(4) << " " << envlpNegDamgdStress(4) << " " << envlpNegDamgdStress(3) << " " << envlpNegDamgdStress(2) << " " << envlpNegDamgdStress(1) << " " << envlpPosDamgdStress(1) << " " << envlpPosDamgdStress(2) << " " << envlpPosDamgdStress(3) << " " << envlpPosDamgdStress(4) << endln; 00814 //(*fg) << tagMat << " " << gammaF1 << " " << gammaF2 << " " << gammaF3 << " " << gammaF4 << " " << gammaK1 << " " << gammaK2 << " " << gammaK3 << " " << gammaK4 << " "<< gammaD1 << " " << gammaD2 << " " << gammaD3 << " " << gammaD4 << endln; 00815 //(*fg) << tagMat << " " << TgammaF << " " << TgammaD << " " << TgammaK << " " << CgammaF << " " << CgammaD << " " << CgammaK << endln; 00816 return 0; 00817 } 00818 00819 int BarSlipMaterial::revertToLastCommit(void) 00820 { 00821 Tstate = Cstate; 00822 00823 TstrainRate = CstrainRate; 00824 00825 lowTstateStrain = lowCstateStrain; 00826 lowTstateStress = lowCstateStress; 00827 hghTstateStrain = hghCstateStrain; 00828 hghTstateStress = hghCstateStress; 00829 TminStrainDmnd = CminStrainDmnd; 00830 TmaxStrainDmnd = CmaxStrainDmnd; 00831 Tenergy = Cenergy; 00832 00833 Tstrain = Cstrain; Tstress = Cstress; 00834 00835 TgammaD = CgammaD; 00836 TgammaK = CgammaK; 00837 TgammaF = CgammaF; 00838 00839 return 0; 00840 } 00841 00842 int BarSlipMaterial::revertToStart(void) 00843 { 00844 Cstate = 0; 00845 Cstrain = 0.0; 00846 Cstress = 0.0; 00847 CstrainRate = 0.0; 00848 lowCstateStrain = envlpNegStrain(0); 00849 lowCstateStress = envlpNegStress(0); 00850 hghCstateStrain = envlpPosStrain(0); 00851 hghCstateStress = envlpPosStress(0); 00852 CminStrainDmnd = envlpNegStrain(1); 00853 CmaxStrainDmnd = envlpPosStrain(1); 00854 Cenergy = 0.0; 00855 CgammaK = 0.0; 00856 CgammaD = 0.0; 00857 CgammaF = 0.0; 00858 00859 Ttangent = envlpPosStress(0)/envlpPosStrain(0); 00860 dstrain = 0.0; 00861 gammaKUsed = 0.0; 00862 gammaFUsed = 0.0; 00863 00864 kElasticPosDamgd = kElasticPos; 00865 kElasticNegDamgd = kElasticNeg; 00866 uMaxDamgd = CmaxStrainDmnd; 00867 uMinDamgd = CminStrainDmnd; 00868 00869 return 0; 00870 } 00871 00872 void BarSlipMaterial::getstate(double u,double du) 00873 { 00874 int cid = 0; 00875 int cis = 0; 00876 int newState = 0; 00877 if (du*CstrainRate<=0.0){ 00878 cid = 1; 00879 } 00880 if (u<lowTstateStrain || u>hghTstateStrain || cid) { 00881 if (Tstate == 0) { 00882 if (u>hghTstateStrain) { 00883 cis = 1; 00884 newState = 1; 00885 lowTstateStrain = envlpPosStrain(0); 00886 lowTstateStress = envlpPosStress(0); 00887 hghTstateStrain = envlpPosStrain(5); 00888 hghTstateStress = envlpPosStress(5); 00889 } 00890 else if (u<lowTstateStrain){ 00891 cis = 1; 00892 newState = 2; 00893 lowTstateStrain = envlpNegStrain(5); 00894 lowTstateStress = envlpNegStress(5); 00895 hghTstateStrain = envlpNegStrain(0); 00896 hghTstateStress = envlpNegStress(0); 00897 } 00898 } 00899 else if (Tstate==1 && du<0.0) { 00900 cis = 1; 00901 if (Cstrain>TmaxStrainDmnd) { 00902 TmaxStrainDmnd = u - du; 00903 } 00904 if (TmaxStrainDmnd<uMaxDamgd) { 00905 TmaxStrainDmnd = uMaxDamgd; 00906 } 00907 if (u<uMinDamgd) { 00908 newState = 2; 00909 gammaFUsed = CgammaF; 00910 for (int i=0; i<=5; i++) { 00911 envlpNegDamgdStress(i) = envlpNegStress(i)*(1.0-gammaFUsed); 00912 } 00913 lowTstateStrain = envlpNegStrain(5); 00914 lowTstateStress = envlpNegStress(5); 00915 hghTstateStrain = envlpNegStrain(0); 00916 hghTstateStress = envlpNegStress(0); 00917 } 00918 else { 00919 newState = 3; 00920 lowTstateStrain = uMinDamgd; 00921 gammaFUsed = CgammaF; 00922 for (int i=0; i<=5; i++) { 00923 envlpNegDamgdStress(i) = envlpNegStress(i)*(1.0-gammaFUsed); 00924 } 00925 lowTstateStress = negEnvlpStress(uMinDamgd); 00926 hghTstateStrain = Cstrain; 00927 hghTstateStress = Cstress; 00928 } 00929 gammaKUsed = CgammaK; 00930 kElasticPosDamgd = kElasticPos*(1.0-gammaKUsed); 00931 } 00932 else if (Tstate ==2 && du>0.0){ 00933 cis = 1; 00934 if (Cstrain<TminStrainDmnd) { 00935 TminStrainDmnd = Cstrain; 00936 } 00937 if (TminStrainDmnd>uMinDamgd) { 00938 TminStrainDmnd = uMinDamgd; 00939 } 00940 if (u>uMaxDamgd) { 00941 newState = 1; 00942 gammaFUsed = CgammaF; 00943 for (int i=0; i<=5; i++) { 00944 envlpPosDamgdStress(i) = envlpPosStress(i)*(1.0-gammaFUsed); 00945 } 00946 lowTstateStrain = envlpPosStrain(0); 00947 lowTstateStress = envlpPosStress(0); 00948 hghTstateStrain = envlpPosStrain(5); 00949 hghTstateStress = envlpPosStress(5); 00950 } 00951 else { 00952 newState = 4; 00953 lowTstateStrain = Cstrain; 00954 lowTstateStress = Cstress; 00955 hghTstateStrain = uMaxDamgd; 00956 gammaFUsed = CgammaF; 00957 for (int i=0; i<=5; i++) { 00958 envlpPosDamgdStress(i) = envlpPosStress(i)*(1.0-gammaFUsed); 00959 } 00960 hghTstateStress = posEnvlpStress(uMaxDamgd); 00961 } 00962 gammaKUsed = CgammaK; 00963 kElasticNegDamgd = kElasticNeg*(1.0-gammaKUsed); 00964 } 00965 else if (Tstate ==3) { 00966 if (u<lowTstateStrain){ 00967 cis = 1; 00968 newState = 2; 00969 lowTstateStrain = envlpNegStrain(5); 00970 hghTstateStrain = envlpNegStrain(0); 00971 lowTstateStress = envlpNegDamgdStress(5); 00972 hghTstateStress = envlpNegDamgdStress(0); 00973 } 00974 else if (u>uMaxDamgd && du>0.0) { 00975 cis = 1; 00976 newState = 1; 00977 lowTstateStrain = envlpPosStrain(0); 00978 lowTstateStress = envlpPosStress(0); 00979 hghTstateStrain = envlpPosStrain(5); 00980 hghTstateStress = envlpPosStress(5); 00981 } 00982 else if (du>0.0) { 00983 cis = 1; 00984 newState = 4; 00985 lowTstateStrain = Cstrain; 00986 lowTstateStress = Cstress; 00987 hghTstateStrain = uMaxDamgd; 00988 gammaFUsed = CgammaF; 00989 for (int i=0; i<=5; i++) { 00990 envlpPosDamgdStress(i) = envlpPosStress(i)*(1.0-gammaFUsed); 00991 } 00992 hghTstateStress = posEnvlpStress(uMaxDamgd); 00993 gammaKUsed = CgammaK; 00994 kElasticNegDamgd = kElasticNeg*(1.0-gammaKUsed); 00995 } 00996 } 00997 else if (Tstate == 4){ 00998 if (u>hghTstateStrain){ 00999 cis = 1; 01000 newState = 1; 01001 lowTstateStrain = envlpPosStrain(0); 01002 lowTstateStress = envlpPosDamgdStress(0); 01003 hghTstateStrain = envlpPosStrain(5); 01004 hghTstateStress = envlpPosDamgdStress(5); 01005 } 01006 else if (u<uMinDamgd && du <0.0) { 01007 cis = 1; 01008 newState = 2; 01009 lowTstateStrain = envlpNegStrain(5); 01010 lowTstateStress = envlpNegDamgdStress(5); 01011 hghTstateStrain = envlpNegStrain(0); 01012 hghTstateStress = envlpNegDamgdStress(0); 01013 } 01014 else if (du<0.0) { 01015 cis = 1; 01016 newState = 3; 01017 lowTstateStrain = uMinDamgd; 01018 gammaFUsed = CgammaF; 01019 for (int i=0; i<=5; i++) { 01020 envlpNegDamgdStress(i) = envlpNegStress(i)*(1.0-gammaFUsed); 01021 } 01022 lowTstateStress = negEnvlpStress(uMinDamgd); 01023 hghTstateStrain = Cstrain; 01024 hghTstateStress = Cstress; 01025 gammaKUsed = CgammaK; 01026 kElasticPosDamgd = kElasticPos*(1.0-gammaKUsed); 01027 } 01028 } 01029 } 01030 if (cis) { 01031 Tstate = newState; 01032 01033 } 01034 } 01035 01036 double BarSlipMaterial::posEnvlpStress(double u) 01037 { 01038 double k = 0.0; 01039 int i = 0; 01040 double f = 0.0; 01041 while (k==0.0 && i<=4){ 01042 01043 if (u<=envlpPosStrain(i+1)){ 01044 k = (envlpPosDamgdStress(i+1)-envlpPosDamgdStress(i))/(envlpPosStrain(i+1)-envlpPosStrain(i)); 01045 f = envlpPosDamgdStress(i) + (u-envlpPosStrain(i))*k; 01046 } 01047 i++; 01048 } 01049 01050 01051 if (k==0.0){ 01052 k = (envlpPosDamgdStress(5) - envlpPosDamgdStress(4))/(envlpPosStrain(5) - envlpPosStrain(4)); 01053 f = envlpPosDamgdStress(5) + k*(u-envlpPosStrain(5)); 01054 } 01055 01056 return f; 01057 01058 } 01059 01060 double BarSlipMaterial::posEnvlpTangent(double u) 01061 { 01062 double k = 0.0; 01063 int i = 0; 01064 while (k==0.0 && i<=4){ 01065 01066 if (u<=envlpPosStrain(i+1)){ 01067 k = (envlpPosDamgdStress(i+1)-envlpPosDamgdStress(i))/(envlpPosStrain(i+1)-envlpPosStrain(i)); 01068 } 01069 i++; 01070 } 01071 01072 if (k==0.0){ 01073 k = (envlpPosDamgdStress(5) - envlpPosDamgdStress(4))/(envlpPosStrain(5) - envlpPosStrain(4)); 01074 } 01075 01076 return k; 01077 01078 } 01079 01080 double BarSlipMaterial::negEnvlpStress(double u) 01081 { 01082 double k = 0.0; 01083 int i = 0; 01084 double f = 0.0; 01085 while (k==0.0 && i<=4){ 01086 if (u>=envlpNegStrain(i+1)){ 01087 k = (envlpNegDamgdStress(i)-envlpNegDamgdStress(i+1))/(envlpNegStrain(i)-envlpNegStrain(i+1)); 01088 f = envlpNegDamgdStress(i+1)+(u-envlpNegStrain(i+1))*k; 01089 } 01090 i++; 01091 } 01092 01093 if (k==0.0){ 01094 k = (envlpNegDamgdStress(4) - envlpNegDamgdStress(5))/(envlpNegStrain(4)-envlpNegStrain(5)); 01095 f = envlpNegDamgdStress(5) + k*(u-envlpNegStrain(5)); 01096 } 01097 return f; 01098 } 01099 01100 double BarSlipMaterial::negEnvlpTangent(double u) 01101 { 01102 double k = 0.0; 01103 int i = 0; 01104 while (k==0.0 && i<=4){ 01105 01106 if (u>=envlpNegStrain(i+1)){ 01107 k = (envlpNegDamgdStress(i)-envlpNegDamgdStress(i+1))/(envlpNegStrain(i)-envlpNegStrain(i+1)); 01108 } 01109 i++; 01110 } 01111 01112 if (k==0.0){ 01113 k = (envlpNegDamgdStress(4) - envlpNegDamgdStress(5))/(envlpNegStrain(4)-envlpNegStrain(5)); 01114 } 01115 return k; 01116 01117 } 01118 01119 void BarSlipMaterial::getState3(Vector& state3Strain, Vector& state3Stress, double kunload) 01120 { 01121 01122 double kmax = (kunload>kElasticNegDamgd) ? kunload:kElasticNegDamgd; 01123 01124 if (state3Strain(0)*state3Strain(3) <0.0){ 01125 // trilinear unload reload path expected, first define point for reloading 01126 state3Strain(1) = lowTstateStrain*rDispN; 01127 if (rForceN-uForceN > 1e-8) { 01128 state3Stress(1) = lowTstateStress*rForceN; 01129 } 01130 else { 01131 if (TminStrainDmnd < envlpNegStrain(3)) { 01132 double st1 = lowTstateStress*uForceN*(1.0+1e-6); 01133 double st2 = envlpNegDamgdStress(4)*(1.0+1e-6); 01134 state3Stress(1) = (st1<st2) ? st1:st2; 01135 } 01136 else { 01137 double st1 = envlpNegDamgdStress(3)*uForceN*(1.0+1e-6); 01138 double st2 = envlpNegDamgdStress(4)*(1.0+1e-6); 01139 state3Stress(1) = (st1<st2) ? st1:st2; 01140 } 01141 } 01142 // if reload stiffness exceeds unload stiffness, reduce reload stiffness to make it equal to unload stiffness 01143 if ((state3Stress(1)-state3Stress(0))/(state3Strain(1)-state3Strain(0)) > kElasticNegDamgd) { 01144 state3Strain(1) = lowTstateStrain + (state3Stress(1)-state3Stress(0))/kElasticNegDamgd; 01145 } 01146 // check that reloading point is not behind point 4 01147 if (state3Strain(1)>state3Strain(3)) { 01148 // path taken to be a straight line between points 1 and 4 01149 double du = state3Strain(3) - state3Strain(0); 01150 double df = state3Stress(3) - state3Stress(0); 01151 state3Strain(1) = state3Strain(0) + 0.33*du; 01152 state3Strain(2) = state3Strain(0) + 0.67*du; 01153 state3Stress(1) = state3Stress(0) + 0.33*df; 01154 state3Stress(2) = state3Stress(0) + 0.67*df; 01155 } 01156 else { 01157 if (TminStrainDmnd < envlpNegStrain(3)) { 01158 state3Stress(2) = uForceN*envlpNegDamgdStress(4); 01159 } 01160 else { 01161 state3Stress(2) = uForceN*envlpNegDamgdStress(3); 01162 } 01163 state3Strain(2) = hghTstateStrain - (hghTstateStress-state3Stress(2))/kunload; 01164 01165 if (state3Strain(2) > state3Strain(3)) { 01166 // point 3 should be along a line between 2 and 4 01167 double du = state3Strain(3) - state3Strain(1); 01168 double df = state3Stress(3) - state3Stress(1); 01169 state3Strain(2) = state3Strain(1) + 0.5*du; 01170 state3Stress(2) = state3Stress(1) + 0.5*df; 01171 } 01172 else if ((state3Stress(2) - state3Stress(1))/(state3Strain(2) - state3Strain(1)) > kmax) { 01173 // linear unload-reload path expected 01174 double du = state3Strain(3) - state3Strain(0); 01175 double df = state3Stress(3) - state3Stress(0); 01176 state3Strain(1) = state3Strain(0) + 0.33*du; 01177 state3Strain(2) = state3Strain(0) + 0.67*du; 01178 state3Stress(1) = state3Stress(0) + 0.33*df; 01179 state3Stress(2) = state3Stress(0) + 0.67*df; 01180 } 01181 else if ((state3Strain(2) < state3Strain(1))||((state3Stress(2)-state3Stress(1))/(state3Strain(2)-state3Strain(1))<0)) { 01182 if (state3Strain(2)<0.0) { 01183 // pt 3 should be along a line between 2 and 4 01184 double du = state3Strain(3)-state3Strain(1); 01185 double df = state3Stress(3)-state3Stress(1); 01186 state3Strain(2) = state3Strain(1) + 0.5*du; 01187 state3Stress(2) = state3Stress(1) + 0.5*df; 01188 } 01189 else if (state3Strain(1) > 0.0) { 01190 // pt 2 should be along a line between 1 and 3 01191 double du = state3Strain(2)-state3Strain(0); 01192 double df = state3Stress(2)-state3Stress(0); 01193 state3Strain(1) = state3Strain(0) + 0.5*du; 01194 state3Stress(1) = state3Stress(0) + 0.5*df; 01195 } 01196 else { 01197 double avgforce = 0.5*(state3Stress(2) + state3Stress(1)); 01198 double dfr = 0.0; 01199 if (avgforce < 0.0){ 01200 dfr = -avgforce/100; 01201 } 01202 else { 01203 dfr = avgforce/100; 01204 } 01205 double slope12 = (state3Stress(1) - state3Stress(0))/(state3Strain(1) - state3Strain(0)); 01206 double slope34 = (state3Stress(3) - state3Stress(2))/(state3Strain(3) - state3Strain(2)); 01207 state3Stress(1) = avgforce - dfr; 01208 state3Stress(2) = avgforce + dfr; 01209 state3Strain(1) = state3Strain(0) + (state3Stress(1) - state3Stress(0))/slope12; 01210 state3Strain(2) = state3Strain(3) - (state3Stress(3) - state3Stress(2))/slope34; 01211 } 01212 } 01213 } 01214 } 01215 else { 01216 // linear unload reload path is expected 01217 double du = state3Strain(3)-state3Strain(0); 01218 double df = state3Stress(3)-state3Stress(0); 01219 state3Strain(1) = state3Strain(0) + 0.33*du; 01220 state3Strain(2) = state3Strain(0) + 0.67*du; 01221 state3Stress(1) = state3Stress(0) + 0.33*df; 01222 state3Stress(2) = state3Stress(0) + 0.67*df; 01223 } 01224 01225 01226 double checkSlope = state3Stress(0)/state3Strain(0); 01227 double slope = 0.0; 01228 01229 // final check 01230 int i = 0; 01231 while (i<3) { 01232 double du = state3Strain(i+1)-state3Strain(i); 01233 double df = state3Stress(i+1)-state3Stress(i); 01234 if (du<0.0 || df<0.0) { 01235 double du = state3Strain(3)-state3Strain(0); 01236 double df = state3Stress(3)-state3Stress(0); 01237 state3Strain(1) = state3Strain(0) + 0.33*du; 01238 state3Strain(2) = state3Strain(0) + 0.67*du; 01239 state3Stress(1) = state3Stress(0) + 0.33*df; 01240 state3Stress(2) = state3Stress(0) + 0.67*df; 01241 slope = df/du; 01242 i = 3; 01243 } 01244 if (slope > 1e-8 && slope < checkSlope) { 01245 state3Strain(1) = 0.0; state3Stress(1) = 0.0; 01246 state3Strain(2) = state3Strain(3)/2; state3Stress(2) = state3Stress(3)/2; 01247 } 01248 i++; 01249 } 01250 // added jun 9th 2004 to prevent slope from getting negative 01251 if (state3Stress(2) <= state3Stress(1)) { 01252 state3Stress(1) = state3Stress(2)*1.02; 01253 } 01254 } 01255 01256 void BarSlipMaterial::getState4(Vector& state4Strain,Vector& state4Stress, double kunload) 01257 { 01258 01259 double kmax = (kunload>kElasticPosDamgd) ? kunload:kElasticPosDamgd; 01260 01261 if (state4Strain(0)*state4Strain(3) <0.0){ 01262 // trilinear unload reload path expected 01263 state4Strain(2) = hghTstateStrain*rDispP; 01264 if (uForceP==0.0){ 01265 state4Stress(2) = hghTstateStress*rForceP; 01266 } 01267 else if (rForceP-uForceP > 1e-8) { 01268 state4Stress(2) = hghTstateStress*rForceP; 01269 } 01270 else { 01271 if (TmaxStrainDmnd > envlpPosStrain(3)) { 01272 double st1 = hghTstateStress*uForceP*(1.0+1e-6); 01273 double st2 = envlpPosDamgdStress(4)*(1.0+1e-6); 01274 state4Stress(2) = (st1>st2) ? st1:st2; 01275 } 01276 else { 01277 double st1 = envlpPosDamgdStress(3)*uForceP*(1.0+1e-6); 01278 double st2 = envlpPosDamgdStress(4)*(1.0+1e-6); 01279 state4Stress(2) = (st1>st2) ? st1:st2; 01280 } 01281 } 01282 // if reload stiffness exceeds unload stiffness, reduce reload stiffness to make it equal to unload stiffness 01283 if ((state4Stress(3)-state4Stress(2))/(state4Strain(3)-state4Strain(2)) > kElasticPosDamgd) { 01284 state4Strain(2) = hghTstateStrain - (state4Stress(3)-state4Stress(2))/kElasticPosDamgd; 01285 } 01286 // check that reloading point is not behind point 1 01287 if (state4Strain(2)<state4Strain(0)) { 01288 // path taken to be a straight line between points 1 and 4 01289 double du = state4Strain(3) - state4Strain(0); 01290 double df = state4Stress(3) - state4Stress(0); 01291 state4Strain(1) = state4Strain(0) + 0.33*du; 01292 state4Strain(2) = state4Strain(0) + 0.67*du; 01293 state4Stress(1) = state4Stress(0) + 0.33*df; 01294 state4Stress(2) = state4Stress(0) + 0.67*df; 01295 } 01296 else { 01297 if (TmaxStrainDmnd > envlpPosStrain(3)) { 01298 state4Stress(1) = uForceP*envlpPosDamgdStress(4); 01299 } 01300 else { 01301 state4Stress(1) = uForceP*envlpPosDamgdStress(3); 01302 } 01303 state4Strain(1) = lowTstateStrain + (-lowTstateStress+state4Stress(1))/kunload; 01304 01305 if (state4Strain(1) < state4Strain(0)) { 01306 // point 2 should be along a line between 1 and 3 01307 double du = state4Strain(2) - state4Strain(0); 01308 double df = state4Stress(2) - state4Stress(0); 01309 state4Strain(1) = state4Strain(0) + 0.5*du; 01310 state4Stress(1) = state4Stress(0) + 0.5*df; 01311 } 01312 else if ((state4Stress(2) - state4Stress(1))/(state4Strain(2) - state4Strain(1)) > kmax) { 01313 // linear unload-reload path expected 01314 double du = state4Strain(3) - state4Strain(0); 01315 double df = state4Stress(3) - state4Stress(0); 01316 state4Strain(1) = state4Strain(0) + 0.33*du; 01317 state4Strain(2) = state4Strain(0) + 0.67*du; 01318 state4Stress(1) = state4Stress(0) + 0.33*df; 01319 state4Stress(2) = state4Stress(0) + 0.67*df; 01320 } 01321 else if ((state4Strain(2) < state4Strain(1))||((state4Stress(2)-state4Stress(1))/(state4Strain(2)-state4Strain(1))<0)) { 01322 if (state4Strain(1)>0.0) { 01323 // pt 2 should be along a line between 1 and 3 01324 double du = state4Strain(2)-state4Strain(0); 01325 double df = state4Stress(2)-state4Stress(0); 01326 state4Strain(1) = state4Strain(0) + 0.5*du; 01327 state4Stress(1) = state4Stress(0) + 0.5*df; 01328 } 01329 else if (state4Strain(2) < 0.0) { 01330 // pt 2 should be along a line between 2 and 4 01331 double du = state4Strain(3)-state4Strain(1); 01332 double df = state4Stress(3)-state4Stress(1); 01333 state4Strain(2) = state4Strain(1) + 0.5*du; 01334 state4Stress(2) = state4Stress(1) + 0.5*df; 01335 } 01336 else { 01337 double avgforce = 0.5*(state4Stress(2) + state4Stress(1)); 01338 double dfr = 0.0; 01339 if (avgforce < 0.0){ 01340 dfr = -avgforce/100; 01341 } 01342 else { 01343 dfr = avgforce/100; 01344 } 01345 double slope12 = (state4Stress(1) - state4Stress(0))/(state4Strain(1) - state4Strain(0)); 01346 double slope34 = (state4Stress(3) - state4Stress(2))/(state4Strain(3) - state4Strain(2)); 01347 state4Stress(1) = avgforce - dfr; 01348 state4Stress(2) = avgforce + dfr; 01349 state4Strain(1) = state4Strain(0) + (state4Stress(1) - state4Stress(0))/slope12; 01350 state4Strain(2) = state4Strain(3) - (state4Stress(3) - state4Stress(2))/slope34; 01351 } 01352 } 01353 } 01354 } 01355 else { 01356 // linear unload reload path is expected 01357 double du = state4Strain(3)-state4Strain(0); 01358 double df = state4Stress(3)-state4Stress(0); 01359 state4Strain(1) = state4Strain(0) + 0.33*du; 01360 state4Strain(2) = state4Strain(0) + 0.67*du; 01361 state4Stress(1) = state4Stress(0) + 0.33*df; 01362 state4Stress(2) = state4Stress(0) + 0.67*df; 01363 } 01364 01365 01366 01367 double checkSlope = state4Stress(0)/state4Strain(0); 01368 double slope = 0.0; 01369 01370 // final check 01371 int i = 0; 01372 while (i<3) { 01373 double du = state4Strain(i+1)-state4Strain(i); 01374 double df = state4Stress(i+1)-state4Stress(i); 01375 if (du<0.0 || df<0.0) { 01376 double du = state4Strain(3)-state4Strain(0); 01377 double df = state4Stress(3)-state4Stress(0); 01378 state4Strain(1) = state4Strain(0) + 0.33*du; 01379 state4Strain(2) = state4Strain(0) + 0.67*du; 01380 state4Stress(1) = state4Stress(0) + 0.33*df; 01381 state4Stress(2) = state4Stress(0) + 0.67*df; 01382 slope = df/du; 01383 i = 3; 01384 } 01385 if (slope > 1e-8 && slope < checkSlope) { 01386 state4Strain(1) = 0.0; state4Stress(1) = 0.0; 01387 state4Strain(2) = state4Strain(3)/2; state4Stress(2) = state4Stress(3)/2; 01388 } 01389 01390 i++; 01391 } 01392 // added jun 9th 2004 to prevent getting a negative slope 01393 if (state4Stress(2) <= state4Stress(1)) { 01394 state4Stress(2) = state4Stress(1)*1.02; 01395 } 01396 } 01397 01398 double BarSlipMaterial::Envlp3Tangent(Vector s3Strain, Vector s3Stress, double u) 01399 { 01400 double k = 0.0; 01401 int i = 0; 01402 while ((k==0.0||i<=2) && (i<=2)) 01403 { 01404 if (u>= s3Strain(i)) { 01405 k = (s3Stress(i+1)-s3Stress(i))/(s3Strain(i+1)-s3Strain(i)); 01406 } 01407 i++; 01408 } 01409 if (k==0.0) { 01410 if (u<s3Strain(0)) { 01411 i = 0; 01412 } 01413 else { 01414 i = 2; 01415 } 01416 k = (s3Stress(i+1)-s3Stress(i))/(s3Strain(i+1)-s3Strain(i)); 01417 01418 } 01419 return k; 01420 } 01421 01422 double BarSlipMaterial::Envlp4Tangent(Vector s4Strain, Vector s4Stress, double u) 01423 { 01424 double k = 0.0; 01425 int i = 0; 01426 while ((k==0.0||i<=2) && (i<=2)) 01427 { 01428 if (u>= s4Strain(i)) { 01429 k = (s4Stress(i+1)-s4Stress(i))/(s4Strain(i+1)-s4Strain(i)); 01430 } 01431 i++; 01432 } 01433 if (k==0.0) { 01434 if (u<s4Strain(0)) { 01435 i = 0; 01436 } 01437 else { 01438 i = 2; 01439 } 01440 k = (s4Stress(i+1)-s4Stress(i))/(s4Strain(i+1)-s4Strain(i)); 01441 01442 } 01443 return k; 01444 } 01445 01446 double BarSlipMaterial::Envlp3Stress(Vector s3Strain, Vector s3Stress, double u) 01447 { 01448 double k = 0.0; 01449 int i = 0; 01450 double f = 0.0; 01451 while ((k==0.0||i<=2) && (i<=2)) 01452 { 01453 if (u>= s3Strain(i)) { 01454 k = (s3Stress(i+1)-s3Stress(i))/(s3Strain(i+1)-s3Strain(i)); 01455 f = s3Stress(i)+(u-s3Strain(i))*k; 01456 } 01457 i++; 01458 } 01459 if (k==0.0) { 01460 if (u<s3Strain(0)) { 01461 i = 0; 01462 } 01463 else { 01464 i = 2; 01465 } 01466 k = (s3Stress(i+1)-s3Stress(i))/(s3Strain(i+1)-s3Strain(i)); 01467 f = s3Stress(i)+(u-s3Strain(i))*k; 01468 } 01469 return f; 01470 } 01471 01472 double BarSlipMaterial::Envlp4Stress(Vector s4Strain, Vector s4Stress, double u) 01473 { 01474 double k = 0.0; 01475 int i = 0; 01476 double f = 0.0; 01477 while ((k==0.0||i<=2) && (i<=2)) 01478 { 01479 if (u>= s4Strain(i)) { 01480 k = (s4Stress(i+1)-s4Stress(i))/(s4Strain(i+1)-s4Strain(i)); 01481 f = s4Stress(i)+(u-s4Strain(i))*k; 01482 } 01483 i++; 01484 } 01485 if (k==0.0) { 01486 if (u<s4Strain(0)) { 01487 i = 0; 01488 } 01489 else { 01490 i = 2; 01491 } 01492 k = (s4Stress(i+1)-s4Stress(i))/(s4Strain(i+1)-s4Strain(i)); 01493 f = s4Stress(i)+(u-s4Strain(i))*k; 01494 } 01495 return f; 01496 } 01497 01498 void BarSlipMaterial::updateDmg(double strain) 01499 { 01500 double umaxAbs = (TmaxStrainDmnd>-TminStrainDmnd) ? TmaxStrainDmnd:-TminStrainDmnd; 01501 double uultAbs = (envlpPosStrain(4)>-envlpNegStrain(4)) ? envlpPosStrain(4):-envlpNegStrain(4); 01502 if ((strain<uultAbs && strain>-uultAbs)&& Tenergy< energyCapacity) 01503 { 01504 TgammaK = gammaK1*pow((umaxAbs/uultAbs),gammaK3); 01505 TgammaD = gammaD1*pow((umaxAbs/uultAbs),gammaD3); 01506 if (damage == 2 || damage == 0) { 01507 TgammaF = gammaF1*pow((umaxAbs/uultAbs),gammaF3); 01508 } 01509 01510 if (damage == 1) { 01511 if (umaxAbs >= envlpPosStrain(3)) { 01512 double a = gammaFLimit/0.7; 01513 double b = -gammaFLimit*3.0/7.0; 01514 double x = (umaxAbs/uultAbs); 01515 TgammaF = a*x + b; 01516 } 01517 } 01518 01519 if (Tenergy>elasticStrainEnergy) { 01520 TgammaK = TgammaK + gammaK2*pow(((Tenergy-elasticStrainEnergy)/energyCapacity),gammaK4); 01521 TgammaD = TgammaD + gammaD2*pow(((Tenergy-elasticStrainEnergy)/energyCapacity),gammaD4); 01522 TgammaF = TgammaF + gammaF2*pow(((Tenergy-elasticStrainEnergy)/energyCapacity),gammaF4); 01523 } 01524 double kminP = (posEnvlpStress(TmaxStrainDmnd)/TmaxStrainDmnd); 01525 double kminN = (negEnvlpStress(TminStrainDmnd)/TminStrainDmnd); 01526 double kmin = (kminP/kElasticPos>kminN/kElasticNeg) ? kminP/kElasticPos:kminN/kElasticNeg; 01527 double gammaKLimEnv = (0.0>(1-kmin)) ? 0.0:(1-kmin); 01528 double k1 = (TgammaK<gammaKLimit) ? TgammaK:gammaKLimit; 01529 TgammaK = (k1<gammaKLimEnv) ? k1:gammaKLimEnv; 01530 TgammaD = (TgammaD<gammaDLimit) ? TgammaD:gammaDLimit; 01531 TgammaF = (TgammaF<gammaFLimit) ? TgammaF:gammaFLimit; 01532 } 01533 else if (strain<uultAbs && strain>-uultAbs) { 01534 double kminP = (posEnvlpStress(TmaxStrainDmnd)/TmaxStrainDmnd); 01535 double kminN = (negEnvlpStress(TminStrainDmnd)/TminStrainDmnd); 01536 double kmin = (kminP/kElasticPos>kminN/kElasticNeg) ? kminP/kElasticPos:kminN/kElasticNeg; 01537 double gammaKLimEnv = (0.0>(1-kmin)) ? 0.0:(1-kmin); 01538 01539 TgammaK = (gammaKLimit<gammaKLimEnv) ? gammaKLimit:gammaKLimEnv; 01540 TgammaD = gammaDLimit; 01541 TgammaF = gammaFLimit; 01542 } 01543 }
Generated on Mon Oct 23 15:05:19 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents