OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
eurocode2creep.C
Go to the documentation of this file.
1 /*
2  *
3  * ##### ##### ###### ###### ### ###
4  * ## ## ## ## ## ## ## ### ##
5  * ## ## ## ## #### #### ## # ##
6  * ## ## ## ## ## ## ## ##
7  * ## ## ## ## ## ## ## ##
8  * ##### ##### ## ###### ## ##
9  *
10  *
11  * OOFEM : Object Oriented Finite Element Code
12  *
13  * Copyright (C) 1993 - 2013 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 
35 #include "eurocode2creep.h"
36 #include "mathfem.h"
37 #include "gausspoint.h"
38 #include "timestep.h"
39 #include "contextioerr.h"
40 #include "datastream.h"
41 #include "classfactory.h"
42 #include "engngm.h"
43 
44 namespace oofem {
45 REGISTER_Material(Eurocode2CreepMaterial);
46 
49 {
50  IRResultType result; // Required by IR_GIVE_FIELD macro
51 
53 
54  /* scaling factor transforming PREDICTED stiffness Ecm
55  * e.g. if the stiffness should be in MPa, then stiffnessFactor = 1.e6
56  * e.g. if the stiffness should be in Pa, then stiffnessFactor = 1. */
57  this->stiffnessFactor = 1.e6;
59 
60  IR_GIVE_FIELD(ir, t0, _IFT_Eurocode2CreepMaterial_t0); // duration of curing (in days)
61  IR_GIVE_FIELD(ir, h0, _IFT_Eurocode2CreepMaterial_h0); // equivalent thickness, in mm !!!
62 
66  }
67 
68  this->temperatureDependent = false;
70  this->temperatureDependent = true;
71  }
72 
73 
74  int cemType;
76 
77  double henv = 1.;
79 
80  // read shrinkage type
81  int sht = 0;
83  if ( sht > 3 ) {
84  OOFEM_ERROR("unsupported shrinkage type");
85  }
86  this->shType = ( ec2ShrinkageType ) sht;
87 
88  this->computeElasticityStrengthParams(cemType);
89  this->computeShrinkageParams(cemType, henv);
90  this->computeCreepParams(cemType, henv);
91 
92  // the initialization of the upper class has to be done here
93  // because RheoChainMat calls at initialization computeCharTimes
94  // but it does not know which method should be chosen (retardation
95  // spectrum switch is initialized here)
96  // initialize Kelvin chain aging material
98 
99  return IRRT_OK;
100 }
101 
102 
103 
104 void
106 {
107  // Ecm28 in GPa (Eurocode formula given in Table 3.1
108  this->Ecm28 = 22. * pow( ( this->fcm28 * ( this->stiffnessFactor / 1.e6 ) / 10. ), 0.3 );
109  this->Ecm28 *= 1.e9 / this->stiffnessFactor; // converted to stiffness units of the analysis
110 
111  // "s" is given in 3.2 in section 3.2.1
112  if ( cemType == 1 ) { // class R
113  this->s = 0.2;
114  } else if ( cemType == 2 ) { // class N
115  this->s = 0.25;
116  } else if ( cemType == 3 ) { // class S
117  this->s = 0.38;
118  } else {
119  OOFEM_ERROR("unsupported value of cement type");
120  }
121 }
122 
123 
124 void
126 {
127  // drying shrinkage (used in B.11)
128  double alpha_ds1, alpha_ds2;
129 
130  if ( cemType == 1 ) { // class R
131  alpha_ds1 = 6;
132  alpha_ds2 = 0.11;
133  } else if ( cemType == 2 ) { // class N
134  alpha_ds1 = 4;
135  alpha_ds2 = 0.12;
136  } else if ( cemType == 3 ) { // class S
137  alpha_ds1 = 3;
138  alpha_ds2 = 0.13;
139  } else {
140  OOFEM_ERROR("unsupported value of cement type");
141  }
142 
143 
144  if ( ( this->shType == EC2_TotalShrinkage ) || ( this->shType == EC2_DryingShrinkage ) ) {
145  // kh from table 3.3
146  if ( this->h0 >= 500 ) {
147  this->kh = 0.7;
148  } else if ( this->h0 >= 300 ) {
149  this->kh = 0.75 - 0.05 * ( h0 - 300. ) / 200.;
150  } else if ( this->h0 >= 200 ) {
151  this->kh = 0.85 - 0.10 * ( h0 - 200. ) / 100.;
152  } else if ( this->h0 >= 100 ) {
153  this->kh = 1.00 - 0.15 * ( h0 - 100. ) / 100.;
154  } else {
155  this->kh = 1.;
156  }
157 
158 
159  double beta_RH = 1.55 * ( 1. - henv * henv * henv );
160  this->eps_cd_0 = -0.85e-6 * ( ( 220. + 110. * alpha_ds1 ) * exp(-alpha_ds2 * this->fcm28 * ( this->stiffnessFactor / 1.e6 ) / 10.) * beta_RH );
161  }
162 
163  if ( ( this->shType == EC2_TotalShrinkage ) || ( this->shType == EC2_AutogenousShrinkage ) ) {
164  // characteristic strength in MPa
165  double fck = this->fcm28 * ( this->stiffnessFactor / 1.e6 ) - 8.;
166  // below 10 MPa no autogenous shrinkage
167  fck = max(fck, 10.);
168 
169  this->eps_ca_infty = -2.5e-6 * ( fck - 10. ); // 3.12
170  }
171 
172  return;
173 }
174 
175 void
177 {
178  // "alpha_T_cement" influences the equivalent age (B.9)
179 
180  if ( cemType == 1 ) { // class R
181  this->alpha_T_cement = 1.;
182  } else if ( cemType == 2 ) { // class N
183  this->alpha_T_cement = 0.;
184  } else if ( cemType == 3 ) { // class S
185  this->alpha_T_cement = -1.;
186  } else {
187  OOFEM_ERROR("unsupported value of cement type");
188  }
189 
190  if ( this->fcm28 * ( this->stiffnessFactor / 1.e6 ) > 35. ) {
191  double alpha_1, alpha_2, alpha_3;
192 
193  // B.8c
194  alpha_1 = pow( ( 35. / ( this->fcm28 * ( this->stiffnessFactor / 1.e6 ) ) ), 0.7 );
195  alpha_2 = pow( ( 35. / ( this->fcm28 * ( this->stiffnessFactor / 1.e6 ) ) ), 0.2 );
196  alpha_3 = pow( ( 35. / ( this->fcm28 * ( this->stiffnessFactor / 1.e6 ) ) ), 0.5 );
197 
198  // influence of relative humidity on drying creep
199  // B.3b
200  this->phi_RH = ( 1. + alpha_1 * ( 1. - henv ) / ( 0.1 * pow(this->h0, 1. / 3.) ) ) * alpha_2;
201 
202  // B.8b
203  // in EC there is 0.012 * RH, it is equivalent to 1.2 * henv
204  this->beta_H = 1.5 * this->h0 * ( 1. + pow( ( 1.2 * henv ), 18. ) ) + 250. * alpha_3;
205  this->beta_H = min(this->beta_H, 1500. * alpha_3);
206 
207  // convert to proper time units
208  //this->beta_H *= this->timeFactor;
209  // this can't be done here because timeFactor has not been initialized yet! (done in RheoChM)
210  // beta_H remains in days
211  } else {
212  // influence of relative humidity on drying creep
213  // B.3a
214  this->phi_RH = 1. + ( 1. - henv ) / ( 0.1 * pow(this->h0, 1. / 3.) );
215 
216  // B.8b
217  // in EC there is 0.012 * RH, it is equivalent to 1.2 * henv
218  this->beta_H = 1.5 * this->h0 * ( 1. + pow( ( 1.2 * henv ), 18. ) ) + 250.;
219  this->beta_H = min(this->beta_H, 1500.);
220 
221  // convert to proper time units
222  // this->beta_H *= this->timeFactor;
223  // this can't be done here because timeFactor has not been initialized yet! (done in RheoChM)
224  }
225 
226  // B.4
227  this->beta_fcm = 16.8 / sqrt( this->fcm28 * ( this->stiffnessFactor / 1.e6 ) );
228 
229  return;
230 }
231 
232 
233 double
235 {
236  double beta = exp( this->s * ( 1. - sqrt( 28. / ( age / this->timeFactor ) ) ) );
237 
238  return beta * this->fcm28;
239 }
240 
241 double
243 {
244  double fcm_at_age;
245  double Ecm_at_age;
246 
247  fcm_at_age = this->computeConcreteStrengthAtAge(age);
248  Ecm_at_age = pow( ( fcm_at_age / this->fcm28 ), 0.3 ) * this->Ecm28;
249 
250  return Ecm_at_age;
251 }
252 
253 double
255 // compute the value of the creep function of concrete
256 // corresponding to the given load duration (in days)
257 // t-t_prime = duration of loading
258 
259 {
260  double J; //return value
261 
262 
263  // elastic modulus isn not temperature adjusted too
264  J = 1. / this->computeMeanElasticModulusAtAge(t_prime);
265 
266  J += this->computeCreepCoefficient(t, t_prime, gp, tStep) / ( 1.05 * this->Ecm28 );
267 
268  return J;
269 }
270 
271 
272 double
274 // compute the value of the creep coefficient
275 // corresponding to the given load duration (in days)
276 // t-t_prime = duration of loading
277 
278 {
279  double phi; //return value
280 
281  double beta_t0, beta_c;
282  double t_prime_equiv;
283 
284  if ( temperatureDependent ) {
285  t_prime_equiv = this->computeEquivalentAge(gp, tStep);
286  } else {
287  t_prime_equiv = t_prime;
288  }
289 
290  // B.5
291  beta_t0 = 1. / ( 0.1 + pow(t_prime_equiv / this->timeFactor, 0.2) );
292 
293  // B.7
294  beta_c = pow( ( ( t - t_prime ) / ( this->beta_H * this->timeFactor + t - t_prime ) ), 0.3 );
295 
296  // t-t'_prime should not be temperature-adjusted (so says the EC)
297 
298  // B.1 + B.2
299  phi = this->phi_RH * this->beta_fcm * beta_t0 * beta_c;
300 
301 
302  return phi;
303 }
304 
305 
306 
307 // implements B.10
308 double
310  Eurocode2CreepMaterialStatus *status = static_cast< Eurocode2CreepMaterialStatus * >( this->giveStatus(gp) );
311 
312  FloatArray et;
313 
314  double initMaturity, incrMaturity;
315  double averageTemperature;
316 
317  // Element *elem = gp->giveElement();
318  StructuralElement *selem = dynamic_cast< StructuralElement * >( gp->giveElement() );
319 
320  selem->computeResultingIPTemperatureAt(et, tStep, gp, VM_Total);
321 
322  if ( tStep->isTheFirstStep() || ( tStep->giveIntrinsicTime() - tStep->giveTimeIncrement() - this->castingTime ) < 0. ) {
323  initMaturity = this->relMatAge;
324  averageTemperature = et.at(1);
325  } else {
326  initMaturity = status->giveConcreteMaturity();
327  averageTemperature = ( et.at(1) + status->giveTemperature() ) / 2.;
328  }
329 
330  averageTemperature = min(averageTemperature, 80.);
331  averageTemperature = max(averageTemperature, 0.);
332 
333  incrMaturity = exp( -1. * ( 4000. / ( 273. + averageTemperature ) - 13.65 ) ) * tStep->giveTimeIncrement();
334 
335  status->setTempConcreteMaturity(initMaturity + incrMaturity); // full increment - for updating
336  status->setTempTemperature( et.at(1) );
337 
338  // middle of the step
339  return initMaturity + incrMaturity / 2.;
340 }
341 
342 
343 // implements B.9
344 double
346 {
347  double maturity;
348  double age;
349 
350  maturity = this->computeEquivalentMaturity(gp, tStep);
351 
352  age = maturity * pow( ( 9. / ( 2. + pow(maturity, 1.2) ) + 1. ), this->alpha_T_cement );
353  age = max(age, 0.5);
354 
355  return age;
356 }
357 
358 
359 
360 double
362 {
363  /*
364  * This function returns the incremental modulus for the given time increment.
365  * The modulus may also depend on the specimen geometry (gp - dependence).
366  *
367  * It is stored as "Einc" for further expected requests from other gaussPoints that correspond to the same material.
368  *
369  * Note: time -1 refers to the previous time.
370  */
371 
372  // !!! chartime exponents are assumed to be equal to 1 !!!
373 
374  double chainStiffness = 0.;
375 
376  if ( !this->isActivated(tStep) ) {
377  return 1.; // stresses are cancelled in giveRealStressVector;
378  }
379 
380  // updateEparModuli is called in KelvinChainMaterial
381  chainStiffness = KelvinChainMaterial :: giveEModulus(gp, tStep);
382 
383  if ( retardationSpectrumApproximation ) { //retardation spectrum used
384  double sum;
385  double t_halfstep;
386 
387  sum = 1. / chainStiffness; // convert stiffness into compliance
388 
389  sum += 1 / this->EspringVal; // add zeroth unit
390 
391  t_halfstep = this->relMatAge + ( tStep->giveTargetTime() - 0.5 * tStep->giveTimeIncrement() );
392 
393  if ( t_halfstep <= 0. ) {
394  OOFEM_ERROR("attempt to evaluate material stiffness at negative age");
395  }
396 
397  sum += 1. / this->computeMeanElasticModulusAtAge(t_halfstep); // add initial compliance
398 
399  // convert to stiffness
400  chainStiffness = 1. / sum;
401  }
402 
403  return chainStiffness;
404 }
405 
406 
407 void
409 {
410  /*
411  * This function generates discrete characteristic times
412  * according to rules which guarantee a good approximation
413  * of the relaxation or creep function by the Dirichlet series
414  */
415 
416  int j;
417 
418  if ( this->begOfTimeOfInterest == -1. ) {
419  this->begOfTimeOfInterest = 0.001 * this->timeFactor;
420  }
421 
422  if ( this->endOfTimeOfInterest == -1. ) {
423  this->endOfTimeOfInterest = 10000. * this->timeFactor;
424  }
425 
426 
427  // perhaps different values should be used for moduli computed by the least-squares method
428  // this gives really big mistakes for times near the interest boundary
429 
430  if ( retardationSpectrumApproximation ) { //retardation spectrum used
431  if ( this->begOfTimeOfInterest > 1.0 * this->timeFactor ) {
432  OOFEM_WARNING("begOfTimeOfInterest was chosen bigger than 1 days, reseting its value to 1 days (could have lead to big errors in the numerical integration of the stiffness of the zeroth Kelvin unit (the retardation spectrum is very steep)");
433  this->begOfTimeOfInterest = 1.0 * this->timeFactor;
434  }
435 
436  this->tau1 = this->begOfTimeOfInterest;
437  } else {
438  this->tau1 = this->begOfTimeOfInterest / 10.;
439  }
440 
441  //first retardation can be treated as equal to begOfTimeOfInterest
442 
443 
444 
445  if ( retardationSpectrumApproximation ) { //retardation spectrum used
446  //the last retardation time has to be bigger than approx sqrt(10) * endOfTimeOfInterest
447  // e.g. if endoftimeofinterest is 10^2, then the last retardation time has to be at least 10^2.5
448  if ( this->endOfTimeOfInterest > pow(10., 5.5) * timeFactor ) {
449  OOFEM_WARNING("endOfTimeOfInterest was chosen bigger than 10.000 days, reseting to 10.000 days (the retardation spectrum is almost zero afterwards)");
450  this->endOfTimeOfInterest = pow(10., 5.5) * timeFactor;
451  }
452  }
453 
454  j = 1;
455  while ( this->endOfTimeOfInterest >= this->tau1 * pow(10.0, ( double ) j - 1. - 0.5) ) {
456  j++;
457  }
458 
459  this->nUnits = j;
460 
461  this->charTimes.resize(this->nUnits);
462  this->charTimes.zero();
463 
464  for ( int mu = 1; mu <= this->nUnits; mu++ ) {
465  this->charTimes.at(mu) = this->tau1 * pow(10., mu - 1);
466 
467  if ( retardationSpectrumApproximation ) { // apply correction factors
468  this->charTimes.at(mu) *= this->computeRetardationTimeCorrection(mu);
469  }
470  }
471 }
472 
473 
474 
475 double
477 {
478  return 1. + 0.555 * exp( -4. * pow( ( this->tau1 * pow( 10., double( mu - 1 ) ) / ( this->beta_H * this->timeFactor ) ), 2. ) );
479 
480 
481  //return 1.;
482 }
483 
484 
485 double
487 {
488  // second-order approximation of the retardation spectrum
489  double L;
490 
491  double t = 2 * tau;
492  double n = 0.3; // exponent in the model
493 
494  double f, df, ddf; // function "t/(beta+t)" and its derivatives
495  double ddPhi; // second derivative of (f(t))^n
496 
497  f = t / ( this->beta_H * this->timeFactor + t );
498 
499  df = this->beta_H * this->timeFactor / ( ( this->beta_H * this->timeFactor + t ) * ( this->beta_H * this->timeFactor + t ) );
500 
501  ddf = -2. * this->beta_H * this->timeFactor / ( ( this->beta_H * this->timeFactor + t ) * ( this->beta_H * this->timeFactor + t ) * ( this->beta_H * this->timeFactor + t ) );
502 
503  ddPhi = n * ( n - 1. ) * pow( f, ( n - 2. ) ) * df * df + n *pow( f, ( n - 1. ) ) * ddf;
504 
505  L = -4. * tau * tau * ddPhi;
506 
507  return L;
508 }
509 
510 
511 void
513 {
514  /*
515  * If retardationSpectrumApproximation == true then analysis of continuous retardation spectrum is used for
516  * computing characteristic coefficients (moduli) of Kelvin chain
517  * Else least-squares method is used
518  */
520  // all moduli must be multiplied by g(t') / c - see equation (46) in Jirasek's retardation spectrum paper
521 
522  double coefficient;
523  double equivalentAge;
524 
525  if ( temperatureDependent ) {
526  equivalentAge = this->computeEquivalentAge(gp, tStep);
527  } else {
528  equivalentAge = atTime;
529  }
530 
531  coefficient = 1.05 * this->Ecm28 * ( 0.1 + pow(equivalentAge / this->timeFactor, 0.2) ) / ( this->phi_RH * this->beta_fcm );
532 
533 
534  // evaluate stiffness of the zero-th unit of the Kelvin chain
535  // (aging elastic spring with retardation time = 0)
536  // this is done employing Simpson's rule. the begOfTimeOfInterest cannot exceed 0.1 day
537  // E0 = EspringVal = int( L, 0, tau1/sqrt(10) )
538 
539  double tau0 = this->tau1 / sqrt(10.0); // upper bound of the integral
540 
541  this->EspringVal = 1. / ( ( log(10.) / 3. ) * (
542  this->evaluateSpectrumAt(tau0 * 1.e-8) + 4. * this->evaluateSpectrumAt(tau0 * 1.e-7) +
543  2. * this->evaluateSpectrumAt(tau0 * 1.e-6) + 4. * this->evaluateSpectrumAt(tau0 * 1.e-5) +
544  2. * this->evaluateSpectrumAt(tau0 * 1.e-4) + 4. * this->evaluateSpectrumAt(tau0 * 1.e-3) +
545  2. * this->evaluateSpectrumAt(tau0 * 1.e-2) + 4. * this->evaluateSpectrumAt(tau0 * 1.e-1) +
546  this->evaluateSpectrumAt(tau0) ) );
547 
548  this->EspringVal *= coefficient;
549 
550  // process remaining units
551  answer.resize(nUnits);
552  answer.zero();
553 
554  double tauMu;
555 
556  for ( int mu = 1; mu <= this->nUnits; mu++ ) {
557  tauMu = this->tau1 * pow( 10., double( mu - 1 ) );
558 
559  answer.at(mu) = 2. / ( log(10.) * ( this->evaluateSpectrumAt( tauMu * pow(10., -sqrt(3.) / 6.) ) + this->evaluateSpectrumAt( tauMu * pow(10., sqrt(3.) / 6.) ) ) );
560  }
561 
562  answer.times(coefficient);
563  } else { // moduli computed using the least-squares method
564  KelvinChainMaterial :: computeCharCoefficients(answer, atTime, gp, tStep);
565  }
566 }
567 
568 
569 void
571  GaussPoint *gp,
572  TimeStep *tStep,
573  ValueModeType mode)
574 {
576  answer.zero();
577 
578  if ( ( this->shType == EC2_NoShrinkage ) || ( !this->isActivated(tStep) ) ) {
579  return;
580  }
581 
582  if ( ( mode != VM_Total ) && ( mode != VM_Incremental ) ) {
583  OOFEM_ERROR("unsupported mode");
584  }
585 
586  FloatArray eps_sh;
587 
588  // the relative matereial age can be negative to properly reflect the casting time
589  double dryingTimeNow = ( tStep->giveTargetTime() + this->relMatAge - this->t0 ) / timeFactor;
590  double autoShrTimeNow = ( this->relMatAge + tStep->giveTargetTime() ) / this->timeFactor;
591 
592  if ( mode == VM_Incremental ) {
593  double dt = tStep->giveTimeIncrement() / this->timeFactor;
594 
595  if ( dryingTimeNow > 0. ) {
596  if ( ( this->shType == EC2_TotalShrinkage ) || ( this->shType == EC2_DryingShrinkage ) ) {
597  this->computeIncrementOfDryingShrinkageVector(eps_sh, gp, dryingTimeNow, dryingTimeNow - dt);
598  answer.add(eps_sh);
599  }
600  }
601 
602  if ( autoShrTimeNow > ( this->castingTime + this->relMatAge ) ) {
603  if ( ( this->shType == EC2_TotalShrinkage ) || ( this->shType == EC2_AutogenousShrinkage ) ) {
604  this->computeIncrementOfAutogenousShrinkageVector(eps_sh, gp, autoShrTimeNow, autoShrTimeNow - dt);
605  answer.add(eps_sh);
606  }
607  }
608  } else { // total
609  if ( dryingTimeNow > 0. ) {
610  if ( ( this->shType == EC2_TotalShrinkage ) || ( this->shType == EC2_DryingShrinkage ) ) {
611  this->computeIncrementOfDryingShrinkageVector( eps_sh, gp, dryingTimeNow, max(0., ( this->relMatAge - this->t0 ) / this->timeFactor) );
612  answer.add(eps_sh);
613  }
614  }
615 
616  if ( autoShrTimeNow > ( this->castingTime + this->relMatAge ) ) {
617  if ( ( this->shType == EC2_TotalShrinkage ) || ( this->shType == EC2_AutogenousShrinkage ) ) {
618  if (this->castingTime > 0.) {
619  this->computeIncrementOfAutogenousShrinkageVector(eps_sh, gp, autoShrTimeNow, (this->relMatAge+this->castingTime) / this->timeFactor);
620  } else {
621  this->computeIncrementOfAutogenousShrinkageVector(eps_sh, gp, autoShrTimeNow, (this->relMatAge) / this->timeFactor);
622  }
623  answer.add(eps_sh);
624  }
625  }
626  }
627 
628  if ( answer.at(1) != answer.at(1) ) {
629  OOFEM_ERROR("shrinkage is NaN: %f", answer.at(1));
630  }
631 }
632 
633 void
635 {
636  int size;
637  FloatArray fullAnswer;
638  MaterialMode mode = gp->giveMaterialMode();
639 
640  if ( ( mode == _3dShell ) || ( mode == _3dBeam ) || ( mode == _2dPlate ) || ( mode == _2dBeam ) ) {
641  size = 12;
642  } else {
643  size = 6;
644  }
645 
646  fullAnswer.resize(size);
647  fullAnswer.zero();
648 
649  if ( tNow > tThen ) {
650  double dEpsSh;
651 
652  // B.10
653  double beta_ds_now, beta_ds_then;
654 
655  beta_ds_now = tNow / ( tNow + 0.04 * pow(this->h0, 3. / 2.) );
656  beta_ds_then = tThen / ( tThen + 0.04 * pow(this->h0, 3. / 2.) );
657 
658  dEpsSh = ( beta_ds_now - beta_ds_then ) * this->kh * this->eps_cd_0;
659 
660  fullAnswer.at(1) = fullAnswer.at(2) = fullAnswer.at(3) = dEpsSh;
661  }
662 
664 }
665 
666 
667 void
669 {
670  int size;
671  FloatArray fullAnswer;
672  MaterialMode mode = gp->giveMaterialMode();
673 
674  if ( ( mode == _3dShell ) || ( mode == _3dBeam ) || ( mode == _2dPlate ) || ( mode == _2dBeam ) ) {
675  size = 12;
676  } else {
677  size = 6;
678  }
679 
680  fullAnswer.resize(size);
681  fullAnswer.zero();
682 
683 
684  if ( tNow > tThen ) {
685  double dEpsAu;
686 
687  // B.13
688  double beta_as_now, beta_as_then;
689 
690  beta_as_now = 1. - exp( -0.2 * sqrt(tNow) );
691  beta_as_then = 1. - exp( -0.2 * sqrt(tThen) );
692 
693  dEpsAu = ( beta_as_now - beta_as_then ) * this->eps_ca_infty;
694 
695  fullAnswer.at(1) = fullAnswer.at(2) = fullAnswer.at(3) = dEpsAu;
696  }
697 
699 }
700 
701 
704 {
705  return new Eurocode2CreepMaterialStatus(1, this->giveDomain(), gp, nUnits);
706 }
707 
708 
709 void
711 {
712  KelvinChainMaterial :: giveRealStressVector(answer, gp, reducedStrain, tStep);
713 
714  if ( !this->isActivated(tStep) ) {
715  return;
716  }
717 }
718 
719 
720 /****************************************************************************************/
721 /********** Eurocode2CreepMaterialStatus - HUMIDITY ****************************************/
722 
723 
725  KelvinChainMaterialStatus(n, d, g, nunits) {
726  tempMaturity = 0.;
728 
729  tempTemperature = 0.;
731 }
732 
733 void
735 {
737  tempMaturity = 0.;
738 
740  tempTemperature = 0.;
741 
743 }
744 
745 
748 //
749 // saves full information stored in this Status
750 //
751 {
752  contextIOResultType iores;
753 
754  if ( ( iores = KelvinChainMaterialStatus :: saveContext(stream, mode, obj) ) != CIO_OK ) {
755  THROW_CIOERR(iores);
756  }
757 
758  // write maturity value
759  if ( !stream.write(maturity) ) {
761  }
762 
763  // write temperature value
764  if ( !stream.write(temperature) ) {
766  }
767 
768  return CIO_OK;
769 }
770 
773 //
774 // restore the state variables from a stream
775 //
776 {
777  contextIOResultType iores;
778  if ( ( iores = KelvinChainMaterialStatus :: restoreContext(stream, mode, obj) ) != CIO_OK ) {
779  THROW_CIOERR(iores);
780  }
781 
782  // read maturity value
783  if ( !stream.read(maturity) ) {
784  return CIO_IOERR;
785  }
786 
787  // read temperature value
788  if ( !stream.read(temperature) ) {
789  return CIO_IOERR;
790  }
791 
792 
793  return CIO_OK;
794 }
795 } // end namespace oofem
static int giveSizeOfVoigtSymVector(MaterialMode mmode)
Returns the size of symmetric part of a reduced stress/strain vector according to given mode...
double h0
effective thickness [mm]
virtual double giveEModulus(GaussPoint *gp, TimeStep *tStep)
Evaluation of the incremental modulus.
MaterialMode giveMaterialMode()
Returns corresponding material mode of receiver.
Definition: gausspoint.h:191
double relMatAge
Physical age of the material at simulation time = 0.
Definition: rheoChM.h:138
ec2ShrinkageType
shrinkage option
virtual void computeResultingIPTemperatureAt(FloatArray &answer, TimeStep *tStep, GaussPoint *gp, ValueModeType mode)
Computes at given time (tStep) the the resulting temperature component array.
#define _IFT_Eurocode2CreepMaterial_shType
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
Definition: material.C:244
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
virtual bool isActivated(TimeStep *tStep)
Definition: rheoChM.h:292
Class and object Domain.
Definition: domain.h:115
virtual void giveShrinkageStrainVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep, ValueModeType mode)
Computes, for the given integration point, the strain vector induced by stress-independent shrinkage...
double t0
duration of curing [day, sec, ...]
void computeElasticityStrengthParams(int)
sets parameters for elasticity and strength according to formulas from EC
double computeRetardationTimeCorrection(int mu)
computes correction factor which multiplies the retardation times
virtual void computeCharCoefficients(FloatArray &answer, double tPrime, GaussPoint *gp, TimeStep *tStep)
Evaluation of characteristic moduli of the Kelvin chain.
void computeIncrementOfDryingShrinkageVector(FloatArray &answer, GaussPoint *gp, double tNow, double tThen)
computes increment of drying shrinkage - the shrinkage strain is isotropic
double s
parameter determined by cement type
The purpose of DataStream abstract class is to allow to store/restore context to different streams...
Definition: datastream.h:54
double eps_cd_0
asymptotic value of drying shrinkage at zero relative humidity, B.11 in EC2
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
int max(int i, int j)
Returns bigger value form two given decimals.
Definition: mathfem.h:71
double evaluateSpectrumAt(double tau)
evaluates retardation spectrum at given time (t-t&#39;)
ValueModeType
Type representing the mode of UnknownType or CharType, or similar types.
Definition: valuemodetype.h:78
double begOfTimeOfInterest
Time from which the model should give a good approximation. Optional field. Default value is 0...
Definition: rheoChM.h:149
double alpha_T_cement
influence of cement type on concrete equivalent age, B.9 in EC2
double tempMaturity
temperature-dependent equivalent age, maturity (temporary value)
virtual bool hasField(InputFieldType id)=0
Returns true if record contains field identified by idString keyword.
#define _IFT_Eurocode2CreepMaterial_h0
double giveTargetTime()
Returns target time.
Definition: timestep.h:146
#define _IFT_Eurocode2CreepMaterial_temperatureDependent
Element * giveElement()
Returns corresponding element to receiver.
Definition: gausspoint.h:188
General IO error.
double maturity
temperature-dependent equivalent age, maturity (equilibrated value)
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
MaterialMode
Type representing material mode of integration point.
Definition: materialmode.h:89
virtual double giveEModulus(GaussPoint *gp, TimeStep *tStep)
Evaluation of the incremental modulus.
Definition: kelvinChM.C:113
void computeShrinkageParams(int, double)
sets parameters for shrinkage according to formulas from EC
double Ecm28
Young&#39;s modulus at 28 days default [MPa].
#define _IFT_Eurocode2CreepMaterial_spectrum
virtual int read(int *data, int count)=0
Reads count integer values into array pointed by data.
virtual void computeCharTimes()
computes retardation times of the aging Kelvin chain
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
FloatArray charTimes
Characteristic times of individual units (relaxation or retardation times).
Definition: rheoChM.h:158
double giveTimeIncrement()
Returns solution step associated time increment.
Definition: timestep.h:150
bool temperatureDependent
switch for temperature dependence of concrete maturity (default option is off)
bool isTheFirstStep()
Check if receiver is first step.
Definition: timestep.C:134
#define _IFT_Eurocode2CreepMaterial_fcm28
virtual void computeCharCoefficients(FloatArray &answer, double tPrime, GaussPoint *gp, TimeStep *tStep)
Evaluation of the moduli of individual units.
Definition: kelvinChM.C:48
virtual double computeConcreteStrengthAtAge(double age)
evaluates concrete strength at given age
virtual int write(const int *data, int count)=0
Writes count integer values from array pointed by data.
double timeFactor
Scaling factor transforming the simulation time units into days (gives the number of simulation time ...
Definition: rheoChM.h:167
Abstract base class for all "structural" finite elements.
virtual void giveRealStressVector(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrain, TimeStep *tStep)
Computes the real stress vector for given total strain and integration point.
double endOfTimeOfInterest
Time (age???) up to which the model should give a good approximation.
Definition: rheoChM.h:151
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: kelvinChM.C:314
static void giveReducedSymVectorForm(FloatArray &answer, const FloatArray &vec, MaterialMode matMode)
Converts the full unsymmetric Voigt vector (2nd order tensor) to reduced form.
#define OOFEM_ERROR(...)
Definition: error.h:61
#define _IFT_Eurocode2CreepMaterial_t0
double eps_ca_infty
asymptotic value of autogenous shrinakge, 3.12 in EC2
virtual MaterialStatus * CreateStatus(GaussPoint *gp) const
Creates new copy of associated status and inserts it into given integration point.
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
double tau1
fixed retardation time of the first unit
virtual double computeMeanElasticModulusAtAge(double age)
evaluates concrete mean elastic modulus at given age
double fcm28
mean compressive strength at 28 days default - to be specified in units of the analysis (e...
This class implements associated Material Status to KelvinChainMaterial.
Definition: kelvinChM.h:44
double beta_fcm
drying creep coefficient
double giveIntrinsicTime()
Returns intrinsic time, e.g. time in which constitutive model is evaluated.
Definition: timestep.h:148
virtual double computeCreepCoefficient(double t, double t_prime, GaussPoint *gp, TimeStep *tStep)
Evaluation of the compliance function (according to appendix B from the EC)
double stiffnessFactor
factor unifying stiffnesses (Ecm is predicted from fcm...)
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
Definition: kelvinChM.C:302
Abstract base class representing a material status information.
Definition: matstatus.h:84
double phi_RH
drying creep coefficient
void computeIncrementOfAutogenousShrinkageVector(FloatArray &answer, GaussPoint *gp, double tNow, double tThen)
computes increment of autogenous shrinkage - the shrinkage strain is isotropic
Class representing vector of real numbers.
Definition: floatarray.h:82
double beta_H
drying creep coefficient
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
double temperature
temperature (equilibrated value)
bool retardationSpectrumApproximation
If true, analysis of retardation spectrum is used for evaluation of Kelvin units moduli If false...
double EspringVal
stiffness of the zeroth Kelvin unit
Class representing the general Input Record.
Definition: inputrecord.h:101
#define _IFT_Eurocode2CreepMaterial_cemType
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: kelvinChM.C:327
virtual void giveRealStressVector(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrain, TimeStep *tStep)
Computes the real stress vector for given total strain and integration point.
Definition: kelvinChM.C:198
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual double computeEquivalentAge(GaussPoint *gp, TimeStep *tStep)
implements B.9
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
long ContextMode
Context mode (mask), defining the type of information written/read to/from context.
Definition: contextmode.h:43
This class implements associated Material Status to Eurocode2CreepMaterial.
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: kelvinChM.C:279
Domain * giveDomain() const
Definition: femcmpnn.h:100
enum oofem::Eurocode2CreepMaterial::ec2ShrinkageType shType
int min(int i, int j)
Returns smaller value from two given decimals.
Definition: mathfem.h:59
#define _IFT_Eurocode2CreepMaterial_henv
REGISTER_Material(DummyMaterial)
double castingTime
Casting time.
Definition: material.h:113
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
Eurocode2CreepMaterialStatus(int n, Domain *d, GaussPoint *g, int nunits)
virtual double computeEquivalentMaturity(GaussPoint *gp, TimeStep *tStep)
implements B.10
the oofem namespace is to define a context or scope in which all oofem names are defined.
double tempTemperature
temperature (temporary value)
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
virtual double computeCreepFunction(double t, double t_prime, GaussPoint *gp, TimeStep *tStep)
Evaluation of the compliance function.
int nUnits
Number of (Maxwell or Kelvin) units in the rheologic chain.
Definition: rheoChM.h:136
Class representing integration point in finite element program.
Definition: gausspoint.h:93
#define OOFEM_WARNING(...)
Definition: error.h:62
double kh
drying shrinkage coefficient
Class representing solution step.
Definition: timestep.h:80
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
void computeCreepParams(int, double)
sets parameters for creep according to formulas from EC
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
#define _IFT_Eurocode2CreepMaterial_stiffnessFactor

This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:28 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011