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