00001 #include "creepb.h"
00002 #include "global.h"
00003 #include "math.h"
00004 #include "globmat.h"
00005 #include "genfile.h"
00006 #include "adaptivity.h"
00007 #include "intpoints.h"
00008 #include "stdlib.h"
00009 #include "elastisomat.h"
00010 #include "therisomat.h"
00011
00012 creepb::creepb (void)
00013 {
00014 timeMax=Mp->timecon.endtime ();
00015 double m = log10(1.2*timeMax), mm; long i;
00016 if ((timeMax >= 0.0) && (timeMax < 110.0))
00017 {
00018 nRetTime=5; tncomp=40;
00019 }
00020 if ((timeMax >= 110.0) && (timeMax < 1100.0))
00021 {
00022 nRetTime=7; tncomp=40;
00023 }
00024 if ((timeMax >= 1100.0) && (timeMax < 14000.0))
00025 {
00026 nRetTime=8; tncomp=50;
00027 }
00028 allocv (nRetTime,retTime);
00029 allocv (nRetTime,ert);
00030 retTime[0]=1.0e-9;
00031
00032 retTime[1]=1.0;
00033
00034 mm=1./(nRetTime-2);
00035 for (i=2;i<nRetTime-1;i++){
00036 retTime[i]=pow(10.,(i-1)*mm*m);
00037 }
00038 retTime[nRetTime -1]=1.2*timeMax;
00039 type_h=1;
00040 type_temp=1;
00041 h_s = 0.4;
00042 desht=0.0;
00043 k_s = 1.0;
00044 k_d = 1.0;
00045 timemat=0.0;
00046 alfa=0.0;
00047 allocm(nRetTime, nRetTime, apom);
00048 ccTime = -1.0e30;
00049
00050 et=0.0;
00051 }
00052 creepb::~creepb (void)
00053 {
00054 }
00055
00056 void creepb::creepinit (long ipp,double val,nonmechquant nmq)
00057 {
00058 long nc=Mm->ip[ipp].ncompstr;
00059 if(nmq == rel_hum )Mm->ip[ipp].eqother[nc*(nRetTime+1)+1]=val;
00060 if(nmq == temperature )Mm->ip[ipp].eqother[nc*(nRetTime+1)+2]=val;
00061 }
00062
00063 void creepb::read (XFILE *in)
00064 {
00065 xfscanf (in,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %ld %ld",
00066 &tb,&t_w,&fc,&wc,&sc,&gc,&c_s,&a1,&k_d,&type_h,&type_temp);
00067
00068 if (type_h==1){xfscanf (in,"%lf ",&h_s); }
00069 if (type_temp==1){xfscanf (in,"%lf ",&temp_s); }
00070
00071
00072 }
00073
00074 long creepb::numberOfCreepb ()
00075 {
00076 return nRetTime;
00077 }
00078
00079
00080
00081
00082
00083
00084
00085
00086 double creepb::approx (vector &areacoord,vector &nodval)
00087 {
00088 double f;
00089 scprd (areacoord,nodval,f);
00090 return f;
00091 }
00092
00093
00094
00095
00096
00097
00098
00099 void creepb::inv_sym (matrix &a)
00100 {
00101 long i,j,k;
00102 double diag,an;
00103 for (k=0;k<a.n;k++){
00104 diag = a[k][k];
00105 if (diag<=0.0){
00106 fprintf (stderr,"\n Creep model: zero diagonal entry in line %ld",k);
00107 abort();
00108 }
00109 for (i=0;i<a.n;i++){
00110 an = a[i][k]/diag;
00111 if (i!=k) {
00112 for (j=i;j<a.n;j++){
00113 if(j!=k) {
00114 a[i][j] = a[i][j] - an*a[k][j];
00115 a[j][i] = a[i][j];
00116 }
00117 }
00118 }
00119 a[i][k] = an;
00120 a[k][i] = an;
00121 }
00122 a[k][k] = - 1.0/diag;
00123 }
00124
00125 for (i=0;i<a.n;i++){
00126 for (j=0;j<a.n;j++){
00127 a[i][j] = - a[i][j];
00128 }
00129 }
00130 }
00131
00132 void creepb::updateval()
00133 {
00134 if(type_h!=1&&type_temp!=1) updatevalchange();
00135 else updatevalkons();
00136
00137 }
00138 void creepb::updatevalchange()
00139 {
00140 long i,j,k;
00141 double pomt,pomt1,qq,delYi,delYj,m0,m;
00142
00143 if (ccTime != Mp->timecon.actualtime ())
00144 {
00145
00146 ddTime=Mp->timecon.actualforwtimeincr ();
00147
00148 ccTime=Mp->timecon.actualtime ();
00149 pomt=ccTime+ddTime/2.;
00150 qq=2./3.;
00151
00152 fillm(0.0, apom);
00153 if (timemat<=pomt) {
00154 pomt1=pomt+0.0001;
00155 m0 = log10(pomt1);
00156 m = log10(2.0*timeMax)*0.999;
00157 k=0;
00158 while (pomt1<=2.0*timeMax)
00159 {
00160 for (i=0;i<nRetTime;i++)
00161 {
00162 delYi=pow((pomt/retTime[i]),qq)-pow((pomt1/retTime[i]),qq);
00163 for (j=0;j<nRetTime;j++)
00164 {
00165 delYj= pow((pomt/retTime[j]),qq)-pow((pomt1/retTime[j]),qq);
00166 apom[i][j]=apom[i][j]+(1.-exp(delYi))*(1.-exp(delYj));
00167 }
00168 }
00169 k=k+1;
00170 pomt1=pow(10.,m0+(m-m0)/tncomp*k);
00171 }
00172 inv_sym(apom);
00173 }
00174
00175 }
00176 }
00177 void creepb::updatevalkons()
00178 {
00179 long i,j,k;
00180 double pomt,pomt1,qq,jt,delYi,delYj,m0,m;
00181 vector bpom(nRetTime);
00182
00183 if (ccTime != Mp->timecon.actualtime ())
00184 {
00185
00186 ddTime=Mp->timecon.actualforwtimeincr ();
00187
00188 ccTime=Mp->timecon.actualtime ();
00189 pomt=ccTime+ddTime/2.;
00190 qq=2./3.;
00191 jt=0.0;
00192
00193 if(timemat<=pomt) {
00194 et=0.;
00195 for (i=0;i<nRetTime;i++){
00196 bpom[i]=0.;
00197 for (j=0;j<nRetTime;j++){
00198 apom[i][j]=0.;
00199 }
00200 }
00201
00202 pomt1=pomt+0.0001;
00203 m0 = log10(pomt1);
00204 m = log10(2.0*timeMax)*0.999;
00205 k=0;
00206 while (pomt1<=2.0*timeMax){
00207
00208 {b3_law(jt, pomt, pomt1);}
00209
00210
00211 for (i=0;i<nRetTime;i++){
00212 delYi=pow((pomt/retTime[i]),qq)-pow((pomt1/retTime[i]),qq);
00213 bpom[i]=bpom[i]+(1.-exp(delYi))*jt;
00214 for (j=0;j<nRetTime;j++){
00215 delYj= pow((pomt/retTime[j]),qq)-pow((pomt1/retTime[j]),qq);
00216 apom[i][j]=apom[i][j]+(1.-exp(delYi))*(1.-exp(delYj));
00217 }
00218 }
00219 k=k+1;
00220 pomt1=pow(10.,m0+(m-m0)/tncomp*k);
00221
00222 }
00223
00224 inv_sym(apom);
00225 for (i=0;i<nRetTime;i++){
00226 ert[i]=0.;
00227 for (j=0;j<nRetTime;j++){
00228 ert[i]=ert[i]+apom[i][j]*bpom[j];
00229 }
00230 ert[i]=1./ert[i];
00231 if(ccTime<0.0001){
00232 delYi=pow((0.0001+ddTime)/retTime[i],qq)-pow(0.0001/retTime[i],qq);}
00233 else {
00234 delYi=pow((ccTime+ddTime)/retTime[i],qq)-pow(ccTime/retTime[i],qq);}
00235 et=et+(delYi-1.+exp(-delYi))/delYi/ert[i];
00236 }
00237 if(et>1.e-55) et=1./et;
00238 else{fprintf (Out,"\n\n Stiff is zero"); abort();}
00239
00240
00241 }
00242
00243 else {
00244
00245 et=1.;
00246 for (i=0;i<nRetTime;i++){ ert[i]=0.;}
00247 }
00248
00249 }
00250
00251 }
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262 void creepb::nlstresses (long ipp, long im, long ido)
00263 {
00264
00265 nc=Mm->ip[ipp].ncompstr;
00266
00267
00268
00269
00270
00271 ss=Mm->ip[ipp].ssst;
00272
00273
00274
00275 ddTime=Mp->timecon.actualforwtimeincr ();
00276
00277 ccTime=Mp->timecon.actualtime ();
00278
00279
00280 imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
00281
00282
00283 if(type_h!=1){
00284 h_s=Mm->givenonmechq(rel_hum, ipp);
00285 }
00286
00287
00288
00289
00290 if(type_temp!=1){
00291 temp_s=Mm->givenonmechq(temperature, ipp);
00292 }
00293
00294
00295
00296 get_h (ipp);
00297 get_temp (ipp);
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308 }
00309
00310 void creepb::phase1 (long ipp,long im, long ido)
00311 {
00312
00313
00314 long i,j,ii;
00315 double qq,dj,endTime;
00316
00317 vector epscr(nc),sig(nc);
00318 matrix d(nc,nc),screep(nc,nRetTime);
00319
00320
00321
00322
00323 endTime=ccTime+ddTime;
00324 get_desht (desht, ccTime, endTime);
00325 for (i=0; i<nc; i++){
00326 epscr[i]=0.0;
00327 }
00328 qq=2.0/3.0;
00329 for (i=0;i<nRetTime;i++){
00330 ii=nc*(i+1);
00331 dj=pow(endTime/retTime[i],qq)-pow(ccTime/retTime[i],qq);
00332 for (j=0; j<nc; j++){
00333 screep[j][i]=Mm->ip[ipp].eqother[ii+j];
00334 epscr[j]=epscr[j]+(1.-exp(-dj))*screep[j][i];
00335 }
00336 }
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348 epscr[0]=epscr[0]+desht;
00349 epscr[1]=epscr[1]+desht;
00350 if(ss==planestrain) epscr[2]=epscr[2]+desht;
00351 else if(ss==axisymm) epscr[2]=epscr[2]+desht;
00352 else if(ss==spacestress) epscr[2]=epscr[2]+desht;
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363 }
00364
00365
00366
00367 void creepb::phase2 (long ipp)
00368 {
00369
00370 long i,j,ii;
00371 double qq,dj, endTime;
00372 vector epscr(nc),deps(nc),depstot(nc),sig(nc);
00373 matrix d(nc,nc),screep(nc,nRetTime);
00374
00375 endTime=ccTime+ddTime;
00376
00377
00378
00379 for (i=0; i<nc; i++){
00380 depstot[i]=Mm->ip[ipp].strain[i]-Mm->ip[ipp].eqother[i];
00381 }
00382 qq=2.0/3.0;
00383 for (i=0;i<nRetTime;i++){
00384 ii=nc*(i+1);
00385 dj=pow(endTime/retTime[i],qq)-pow(ccTime/retTime[i],qq);
00386 for (j=0; j<nc; j++){
00387 screep[j][i]=Mm->ip[ipp].eqother[ii+j];
00388 deps[j]=depstot[j]-(1.-exp(-dj))*screep[j][i];
00389 }
00390 }
00391 get_desht (desht, ccTime, endTime);
00392
00393
00394
00395
00396
00397 deps[0]=deps[0]-desht;
00398 deps[1]=deps[1]-desht;
00399 if(ss==planestrain) deps[2]=deps[2]-desht;
00400 else if(ss==axisymm) deps[2]=deps[2]-desht;
00401 else if(ss==spacestress) deps[2]=deps[2]-desht;
00402
00403
00404
00405
00406
00407
00408
00409
00410 Mm->matstiff(d,ipp);
00411 e0=Mm->eliso[imat].e;
00412 mi=Mm->eliso[imat].nu;
00413 mxv (d,deps,sig);
00414
00415
00416 for (i=0; i<nc; i++){Mm->ip[ipp].eqother[i]=Mm->ip[ipp].strain[i];}
00417
00418 seps_time (screep,sig);
00419 for (i=0; i<nRetTime; i++){
00420 ii=nc*(i+1);
00421 for (j=0; j<nc; j++){
00422 Mm->ip[ipp].eqother[ii+j]=screep[j][i];
00423 }
00424 }
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438 Mm->ip[ipp].eqother[nc*(nRetTime+1)+0]+=desht;
00439 Mm->ip[ipp].eqother[nc*(nRetTime+1)+1]=h_s;
00440 Mm->ip[ipp].eqother[nc*(nRetTime+1)+2]=temp_s;
00441
00442
00443
00444
00445 for (i=0; i<nc; i++){
00446 Mm->ip[ipp].eqother[nc*(nRetTime+1)+3+i]+=depstot[i]-deps[i];
00447 Mm->ip[ipp].eqother[nc*(nRetTime+2)+3+i]+=sig[i];
00448 }
00449 }
00450
00451 void creepb::get_h (long ipp)
00452 {
00453 h_slast=1.0;
00454 if(type_h !=1 ){h_slast=Mm->ip[ipp].eqother[nc*(nRetTime+1)+1]; }
00455
00456
00457
00458
00459 }
00460
00461
00462 void creepb::get_temp (long ipp)
00463 {
00464 temp_slast=283.0;
00465 if(type_temp != 1 ){temp_slast=Mm->ip[ipp].eqother[nc*(nRetTime+1)+2]; }
00466
00467
00468
00469
00470 }
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481 void creepb::matstiff (matrix &d, long ipp)
00482
00483
00484 {
00485 double qq=1.;
00486
00487 nc=Mm->ip[ipp].ncompstr;
00488 ss=Mm->ip[ipp].ssst;
00489
00490 imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
00491
00492 Mm->elmatstiff (d,ipp);
00493
00494 e0=Mm->eliso[imat].e;
00495 mi=Mm->eliso[imat].nu;
00496
00497 if(type_h!=1 && type_temp!=1){
00498 matstiffchange(qq, ipp);
00499 }
00500 else{
00501 matstiffkons(qq);
00502
00503 }
00504
00505
00506
00507 cmulm( qq, d);
00508 }
00509 void creepb::matstiffkons (double &qq)
00510
00511 {
00512
00513 qq=et/e0;
00514 }
00515
00516 void creepb::matstiffchange (double &qq, long ipp)
00517
00518
00519 {
00520 long i,j,k;
00521 double jt,pomt,pomt1,q,delYi,m0,m;
00522 vector bpom(nRetTime);
00523
00524 if(Mm->ip[ipp].hmt & 1)
00525 {i=Mm->ip[ipp].idm[Mm->ip[ipp].nm-1]; alfa=Mm->tidilat[i].alpha;}
00526
00527 if(type_h!=1) {h_s=Mm->givenonmechq(rel_hum, ipp);}
00528 if(type_temp!=1) {temp_s=Mm->givenonmechq(temperature, ipp);}
00529 get_h (ipp);
00530 get_temp (ipp);
00531
00532 pomt=ccTime+ddTime/2.;
00533 q=2./3.;
00534 jt=0.0;
00535
00536 if(timemat<=pomt) {
00537 et=0.;
00538 for (i=0;i<nRetTime;i++){
00539 bpom[i]=0.;
00540
00541
00542
00543 }
00544
00545 pomt1=pomt+0.0001;
00546 m0 = log10(pomt1);
00547 m = log10(2.*timeMax)*0.999;
00548 k=0;
00549 while (pomt1<=2.*timeMax){
00550
00551 {b3_law(jt, pomt, pomt1);}
00552
00553
00554 for (i=0;i<nRetTime;i++){
00555 delYi=pow((pomt/retTime[i]),q)-pow((pomt1/retTime[i]),q);
00556 bpom[i]=bpom[i]+(1.-exp(delYi))*jt;
00557 for (j=0;j<nRetTime;j++){
00558
00559
00560 }
00561 }
00562 k=k+1;
00563 pomt1=pow(10.,m0+(m-m0)/50*k);
00564
00565 }
00566
00567
00568 for (i=0;i<nRetTime;i++){
00569 ert[i]=0.;
00570 for (j=0;j<nRetTime;j++){
00571 ert[i]=ert[i]+apom[i][j]*bpom[j];
00572 }
00573 ert[i]=1./ert[i];
00574 if(ccTime<0.0001){
00575 delYi=pow((0.0001+ddTime)/retTime[i],q)-pow(0.0001/retTime[i],q);}
00576 else {
00577 delYi=pow((ccTime+ddTime)/retTime[i],q)-pow(ccTime/retTime[i],q);}
00578 et=et+(delYi-1.+exp(-delYi))/delYi/ert[i];
00579 }
00580 if(et>1.e-50) et=1./et;
00581 else{fprintf (Out,"\n\n Stiff is zero"); abort();}
00582 }
00583
00584 else {
00585 et=1.;
00586 for (i=0;i<nRetTime;i++){ ert[i]=0.;}
00587 }
00588
00589
00590
00591
00592
00593 qq=et/e0;
00594
00595 }
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605 void creepb::seps_time (matrix &screep,vector &sig)
00606 {
00607 int i;
00608 double q,dj,pomt,pp;
00609
00610 q=2./3.;
00611 if((ccTime+ddTime)==ddTime) {
00612 }
00613 else{
00614 pomt=ccTime+ddTime;
00615 if(ss==planestress){
00616 for (i=0;i<nRetTime;i++){
00617 dj=pow(pomt/retTime[i],q)-pow(ccTime/retTime[i],q);
00618 screep[0][i]=exp(-dj)*screep[0][i]+(sig[0]-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
00619 screep[1][i]=exp(-dj)*screep[1][i]+(sig[1]-mi*sig[0])*(1.-exp(-dj))/dj/ert[i];
00620 screep[2][i]=exp(-dj)*screep[2][i]+ sig[2]*2.*(1.+mi)*(1.-exp(-dj))/dj/ert[i];
00621 screep[3][i]=0.0;
00622 }
00623 }
00624 else if(ss==planestrain){
00625 pp=1.+mi;
00626 for (i=0;i<nRetTime;i++){
00627 dj=pow(pomt/retTime[i],q)-pow(ccTime/retTime[i],q);
00628 screep[0][i]=exp(-dj)*screep[0][i]+ pp*(sig[0]*(1.-mi)-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
00629 screep[1][i]=exp(-dj)*screep[1][i]+ pp*(sig[1]*(1.-mi)-mi*sig[0])*(1.-exp(-dj))/dj/ert[i];
00630 screep[2][i]=exp(-dj)*screep[2][i]+ pp*sig[2]*2. *(1.-exp(-dj))/dj/ert[i];
00631 screep[3][i]=0.0;
00632 }
00633 }
00634 else if(ss==axisymm){
00635 for (i=0;i<nRetTime;i++){
00636 dj=pow(pomt/retTime[i],q)-pow(ccTime/retTime[i],q);
00637 screep[0][i]=exp(-dj)*screep[0][i]+(sig[0]-mi*sig[1]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
00638 screep[1][i]=exp(-dj)*screep[1][i]+(sig[1]-mi*sig[0]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
00639 screep[2][i]=exp(-dj)*screep[2][i]+(sig[2]-mi*sig[0]-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
00640 screep[3][i]=exp(-dj)*screep[3][i]+ sig[3]*(1.+mi) *(1.-exp(-dj))/dj/ert[i];
00641 }
00642 }
00643 else if(ss==spacestress){
00644 for (i=0;i<nRetTime;i++){
00645 dj=pow(pomt/retTime[i],q)-pow(ccTime/retTime[i],q);
00646 screep[0][i]=exp(-dj)*screep[0][i]+(sig[0]-mi*sig[1]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
00647 screep[1][i]=exp(-dj)*screep[1][i]+(sig[1]-mi*sig[0]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
00648 screep[2][i]=exp(-dj)*screep[2][i]+(sig[2]-mi*sig[0]-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
00649 screep[3][i]=exp(-dj)*screep[3][i]+ sig[3]*(1.+mi) *(1.-exp(-dj))/dj/ert[i];
00650 screep[4][i]=exp(-dj)*screep[4][i]+ sig[4]*(1.+mi) *(1.-exp(-dj))/dj/ert[i];
00651 screep[5][i]=exp(-dj)*screep[5][i]+ sig[5]*(1.+mi) *(1.-exp(-dj))/dj/ert[i];
00652 }
00653 }
00654 }
00655
00656 }
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667 void creepb::get_desht (double &des_hn, double t0, double t)
00668 {
00669 double esn,kh,tau,ht,ht0,dh_s,dtemp;
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684 if (tb<=t0){
00685
00686
00687 if ( type_h ==1){
00688 if ( (t0>=t_w) ){
00689
00690 tau=k_s*k_s*k_d*k_d*85000./pow(t,0.08)/pow(fc,0.25);
00691 esn=-a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
00692 ht0 =tanh(sqrt((t0-t_w)/tau));
00693 ht =tanh(sqrt((t -t_w)/tau));
00694 if(h_s<=0.98){kh=(1.0-pow(h_s,3.0));}
00695 else {kh=-10.0*(1.-h_s);}
00696 des_hn = esn*kh/1000000.0*(ht-ht0);
00697 }
00698 else {
00699 des_hn=0.0;
00700 }
00701 }
00702 else {
00703
00704 dh_s = h_s-h_slast;
00705 esn = a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
00706 kh = 3.0*h_s*h_s;
00707 des_hn = esn*kh*dh_s*1.e-6;
00708 }
00709 if ( type_temp ==1 ){
00710
00711 }
00712 else {
00713
00714 dtemp = temp_s-temp_slast;
00715 esn=alfa*1.e6;
00716 des_hn=des_hn+esn*dtemp*1.e-6;
00717 }
00718 }
00719
00720 else {
00721 des_hn=0.0;
00722 }
00723 }
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733 void creepb::b3_law (double &jt, double t0, double t)
00734 {
00735 double n,c0,cd, ac,ag,m,x,y,q,q1,q2,q3,q4,q5,esn,kh,tau,ht,ht0,dh_s,dtemp;
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756 if (tb<t0){
00757 if (t_w<tb) t_w = tb;
00758 if (t_w<=0.0) t_w = 1.0;
00759 ac=sc+gc; ag=ac/gc; m=0.5;
00760
00761 x=(2.1*ac/pow(sc,1.4)+0.1*pow((fc/6895),1.5)*pow(wc,(1.0/3.0))*pow(ag,2.2))*a1-4;
00762 if (x<=4){n=0.12;}
00763 else {n=0.12+(0.07*pow(x,6.0)/(5130.0+pow(x,6.0)));}
00764 q1=600000.0/4734.0/sqrt(fc);
00765 x=1.7*pow(t0,0.12)+8.0;
00766 y=(0.086*pow(t0,(2.0/9.0))+1.21*pow(t0,(4.0/9.0)));
00767 q=1.0/y/pow( ( 1.0+ pow( (pow(t0,m)/log(1.+pow((t-t0),n))/y),x ) ),(1./x));
00768 q2=185.4*sqrt(c_s)/pow(fc,0.9);
00769 q3=0.29*pow(wc,4.0)*q2;
00770 q4=20.3/pow(ac,0.7);
00771
00772 c0=q2*q+q3*log(1.0+pow((t-t0),n))+q4*log(t/t0);
00773
00774
00775 if ( type_h ==1){
00776 if ( (t0>=t_w) ){
00777
00778 tau=k_s*k_s*k_d*k_d*85000./pow(t,0.08)/pow(fc,0.25);
00779 esn=a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
00780 ht0 =tanh(sqrt((t0-t_w)/tau));
00781 ht =tanh(sqrt((t -t_w)/tau));
00782 if(h_s<=0.98){kh=(1.0-pow(h_s,3.0));}
00783 else {kh=-10.0*(1.-h_s);}
00784
00785 q5 = 0.757/fc/pow(esn,0.6);
00786 cd = q5*sqrt( exp(-8+8.*(1.0-h_s)*ht)-exp(-8.+8.*(1.0-h_s)*ht0) );
00787 }
00788 else {
00789 cd = 0.0;
00790 }
00791 }
00792 else {
00793
00794 dh_s = h_s-h_slast;
00795 esn = a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
00796 kh = 3.0*h_s*h_s;
00797 q5 = 0.35/fc;
00798 cd = fabs(q5*esn*kh*dh_s);
00799 }
00800 if ( type_temp ==1 ){
00801
00802 }
00803 else {
00804
00805 dtemp = temp_s-temp_slast;
00806 esn=alfa*1.e6;
00807 q5= 1.5/fc;
00808 cd = cd+fabs(q5*esn*dtemp);
00809 }
00810 }
00811
00812 else {
00813 q1 = 0.0;
00814 c0 = 0.0;
00815 cd = 0.0;
00816 }
00817
00818
00819
00820 jt=(q1+c0+cd)*1.e-6;
00821
00822
00823 }
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836 void creepb::giveirrstrains (long ipp, long ido, vector &epscr)
00837 {
00838 long i;
00839 long ncompstr = Mm->ip[ipp].ncompstr;
00840 long id=ncompstr+ncompstr*nRetTime+3;
00841 for (i=0;i<epscr.n;i++)
00842 epscr[i] = Mm->ip[ipp].eqother[ido+id+i];
00843 }