00001 #include "consol.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 "winpast.h"
00009 #include "stdlib.h"
00010 #include "elastisomat.h"
00011
00012 consol::consol (void)
00013 {
00014 nRetTime=7;
00015 }
00016
00017 consol::~consol (void)
00018 {
00019
00020 }
00021
00022 void consol::consolinit (long mie)
00023 {
00024 }
00025
00026 void consol::read (XFILE *in)
00027 {
00028 xfscanf (in,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",&cn,&gama,&vv,&ni,&e0,&m,&hpodl,&hpodlmax,&rr,&nf);
00029
00030 }
00031
00032 long consol::numberOfConsol ()
00033 {
00034 return nRetTime;
00035 }
00036
00037
00038
00039
00040
00041
00042
00043 void consol::nlstressesincr (long ipp)
00044 {
00045 nc=Mm->ip[ipp].ncompstr;
00046 if(nc==6) ncc=2;
00047 if(nc==3) ncc=3;
00048
00049 ddTime=Mp->timecon.actualforwtimeincr ();
00050 ccTime=Mp->timecon.actualtime ();
00051 imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
00052 phase1(ipp);
00053
00054 }
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066 void consol::nlstresses (long ipp)
00067 {
00068 nc=Mm->ip[ipp].ncompstr;
00069 if(nc==6) ncc=2;
00070 if(nc==3) ncc=3;
00071
00072 ddTime=Mp->timecon.actualforwtimeincr ();
00073 ccTime=Mp->timecon.actualtime ();
00074 imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
00075 phase2(ipp);
00076 }
00077
00078
00079
00080 void consol::phase1 (long ipp)
00081 {
00082 long i,ii,liche;
00083 double delY,dj;
00084 vector epscr(nc),deps(nc),sig(nc);
00085 matrix d(nc,nc),screep(ncc,nRetTime);
00086
00087
00088
00089 for (i=0; i<nc; i++){epscr[i]=0.0;}
00090 get_hc ( nc, ipp);
00091 delY=2.4674*ncn*cn*ddTime/hpodl/nh/hpodl/nh;
00092
00093 for (i=0;i<nRetTime;i++){
00094 ii=nc+ncc*i;
00095 liche=2*i+1;
00096 dj=delY*liche*liche;
00097 if(nc==6){
00098 screep[0][i]=Mm->ip[ipp].other[ii];
00099 epscr[2]=epscr[2]+(1.-exp(-dj))*screep[0][i];
00100 screep[1][i]=Mm->ip[ipp].other[ii+1];
00101 epscr[4]=epscr[4]+(1.-exp(-dj))*screep[1][i];
00102 }
00103
00104 else{
00105 screep[0][i]=Mm->ip[ipp].other[ii];
00106 epscr[0]=epscr[0]+(1.-exp(-dj))*screep[0][i];
00107 screep[1][i]=Mm->ip[ipp].other[ii+1];
00108 epscr[1]=epscr[1]+(1.-exp(-dj))*screep[1][i];
00109 screep[2][i]=Mm->ip[ipp].other[ii+2];
00110 epscr[2]=epscr[2]+(1.-exp(-dj))*screep[2][i];
00111 }
00112
00113 }
00114
00115
00116
00117
00118
00119
00120
00121
00122 Mm->storestress (0,ipp,epscr);
00123
00124
00125 }
00126
00127 void consol::phase2 (long ipp)
00128 {
00129 long i,j,ii,liche;
00130 double delY,dj;
00131 vector epscr(nc),deps(nc),sig(nc);
00132 matrix d(nc,nc),screep(ncc,nRetTime);
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142 for (i=0; i<nc; i++){deps[i]=Mm->ip[ipp].strain[i] - Mm->ip[ipp].other[i];}
00143 get_hc ( nc,ipp);
00144 delY=2.4674*ncn*cn*ddTime/hpodl/nh/hpodl/nh;
00145
00146 for (i=0;i<nRetTime;i++){
00147 ii=nc+ncc*i;
00148 liche=2*i+1;
00149 dj=delY*liche*liche;
00150 if(nc==6){
00151 screep[0][i]=Mm->ip[ipp].other[ii];
00152 deps[2]=deps[2]-(1.-exp(-dj))*screep[0][i];
00153 screep[1][i]=Mm->ip[ipp].other[ii+1];
00154 deps[4]=deps[4]-(1.-exp(-dj))*screep[1][i];
00155 }
00156 else{
00157 screep[0][i]=Mm->ip[ipp].other[ii];
00158 deps[0]=deps[0]-(1.-exp(-dj))*screep[0][i];
00159 screep[1][i]=Mm->ip[ipp].other[ii+1];
00160 deps[1]=deps[1]-(1.-exp(-dj))*screep[1][i];
00161 screep[2][i]=Mm->ip[ipp].other[ii+2];
00162 deps[2]=deps[2]-(1.-exp(-dj))*screep[2][i];
00163 }
00164 }
00165
00166 for (i=0; i<nc; i++){
00167 Mm->ip[ipp].other[i]=Mm->ip[ipp].strain[i];
00168 sig[i]=deps[i];
00169 }
00170
00171 seps_time (screep,sig,nc,ncc,ipp);
00172 for (i=0; i<nRetTime; i++){
00173 ii=nc+ncc*i;
00174 for (j=0; j<ncc; j++){
00175 Mm->ip[ipp].other[ii+j]=screep[j][i];
00176 }
00177 }
00178
00179
00180
00181
00182
00183
00184
00185 ii=nc+ncc*nRetTime;
00186 if(nc==6){
00187 Mm->ip[ipp].other[ii]+=deps[2];
00188 }
00189 else{
00190 Mm->ip[ipp].other[ii]+=deps[0];
00191 }
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203 }
00204
00205
00206
00207
00208
00209
00210
00211
00212 void consol::matstiff (matrix &d,long ipp)
00213
00214 {
00215 long iwink,j,liche,nc=Mm->ip[ipp].ncompstr;
00216 double delY,dj;
00217 if(Mm->ip[ipp].gemid()==2)
00218 {
00219
00220 iwink=Mm->ip[ipp].idm[1];
00221 Mm->wpast[iwink].matstiff (d,ipp);
00222 }
00223
00224
00225
00226
00227 ddTime=Mp->timecon.actualforwtimeincr ();
00228 vlivTCSum=0.0;
00229 get_hc ( d.n,ipp);
00230 delY=2.4674*ncn*cn*ddTime/hpodl/nh/hpodl/nh;
00231 for (j=0;j<nRetTime;j++){
00232 liche=2*j+1;
00233 dj=delY*liche*liche;
00234 vlivTCSum=vlivTCSum+0.81057*(1.-(1.-exp(-dj))/dj)/liche/liche;
00235 }
00236
00237
00238
00239
00240
00241 for (j=0;j<d.n;j++){
00242 d[j][j]=d[j][j]/vlivTCSum;
00243 }
00244
00245 if(nc==6){
00246 for (j=0;j<d.n/2;j++){
00247 d[j][j]=d[j][j]*nc1;
00248 }
00249 for (j=d.n/2;j<d.n;j++){
00250 d[j][j]=d[j][j]*nc2;
00251 }
00252 }
00253 else{
00254 d[0][0]=d[0][0]*nc1;
00255 for (j=1;j<d.n;j++){
00256 d[j][j]=d[j][j]*nc2;
00257 }
00258 }
00259
00260
00261
00262
00263
00264
00265 }
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275 void consol::seps_time (matrix &screep,vector &sig,long nc,long ncc,long ipp)
00276 {
00277 int j,i,liche;
00278 double dj,delY;
00279 vector rp(ncc);
00280
00281
00282
00283 ddTime=Mp->timecon.actualforwtimeincr ();
00284 vlivTCSum=0.0;
00285 get_hc ( nc,ipp);
00286 delY=2.4674*ncn*cn*ddTime/hpodl/nh/hpodl/nh;
00287 for (j=0;j<nRetTime;j++){
00288 liche=2*j+1;
00289 dj=delY*liche*liche;
00290 vlivTCSum=vlivTCSum+0.81057*(1.-(1.-exp(-dj))/dj)/liche/liche;
00291 }
00292
00293 if(nc==6){
00294 rp[0]=sig[2]/vlivTCSum;
00295 rp[1]=sig[4]/vlivTCSum;
00296
00297
00298 }
00299 else {
00300 rp[0]=sig[0]/vlivTCSum;
00301 rp[1]=sig[1]/vlivTCSum;
00302 rp[2]=sig[2]/vlivTCSum;
00303 }
00304
00305
00306
00307 delY=2.4674*ncn*cn*ddTime/hpodl/nh/hpodl/nh;
00308 for (j=0;j<nRetTime;j++){
00309 liche=2*j+1;
00310 dj=delY*liche*liche;
00311 for (i=0;i<ncc;i++){
00312 screep[i][j]=exp(-dj)*screep[i][j]+0.81057*rp[i]*(1.-exp(-dj))/dj/liche/liche;
00313 }
00314 }
00315
00316
00317
00318 }
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329 void consol::get_hc (long nc,long ipp)
00330 {
00331 int ii, ncc;
00332 double f,smax,ee,p,pp,hactual,*sumf,*sizexyz;
00333 vector deps(nc);
00334
00335
00336 if(nf<1.e-6)
00337 {
00338 nh=1.; nc1=1.; nc2=1.; ncn=1.;
00339 return;
00340 }
00341
00342 sizexyz = new double [3];
00343 sumf = new double [nc];
00344 Gtm->give_domain_sizes (sizexyz);
00345 Mb->give_comp_sum (sumf);
00346
00347 if(nc==6){
00348 if(rr<1.e-6)
00349 {if(sizexyz[0]>sizexyz[1])
00350 {smax=sizexyz[0]/2.;}
00351 else
00352 {smax=sizexyz[1]/2.;}
00353 }
00354 else
00355 {smax=rr;}
00356
00357 f=nf*sumf[2]/smax/2.;
00358 ncc=2;
00359 }
00360 else {
00361 if(rr<1.e-6)
00362 {smax=(sizexyz[0]+sizexyz[1])/4.;
00363 }
00364 else
00365 {smax=rr;}
00366
00367 f=nf*sumf[0]/smax/smax/3.14;
00368 ncc=3;
00369 }
00370
00371 p=0.5*(1.-gama*vv/(-f));
00372 if(p>=0.5)
00373 pp=0.0001;
00374 if(p>=0.4)
00375 pp=(0.5*0.4/(0.5-0.4)+0.5)-0.5/(0.5-0.4)*p;
00376 if(p>=0.24)
00377 pp=(0.5*0.24/(0.4-0.24)+1.)-0.5/(0.4-0.24)*p;
00378 else if(p>=0.06)
00379 pp=(0.06/(0.24-0.06)+2.)-1./(0.24-0.06)*p;
00380 else if(p>=0.02)
00381 pp=(0.02/(0.06-0.02)+3.)-1./(0.06-0.02)*p;
00382 else if(p>=0.005)
00383 pp=(0.005/(0.02-0.005)+4.)-1./(0.02-0.005)*p;
00384 else
00385 pp=5.;
00386
00387 hactual=smax/pp*sqrt((2.-2.*ni)/(1.-2.*ni));
00388 if(hpodlmax>=0.0)
00389 if(hactual>hpodlmax)
00390 hactual=hpodlmax;
00391
00392 nh=hactual/hpodl;
00393
00394 nc1=(0.051/nh/nh+0.888/nh+0.0606);
00395
00396 nc2=(1.2639/nh/nh-4.7236/nh+4.46);
00397
00398
00399 if(e0<1.e-6)
00400 {
00401 ncn=1.;
00402 }
00403 else
00404 {
00405 ii=nc+ncc*nRetTime;
00406 ee = Mm->ip[ipp].other[ii]/nh/hpodl;
00407 ee = -ee*(1.+e0)+e0;
00408 ncn=pow((ee/e0),m);
00409 }
00410 delete [] sizexyz;
00411 delete [] sumf;
00412
00413 }
00414
00415 void consol::updateval (long ipp, long im, long ido)
00416 {
00417 long i,n=Mm->givencompeqother(ipp,im);
00418 for (i=0;i<n;i++){
00419 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00420 }
00421 }
00422
00423
00424 void consol::give_dstresses_eqother(vector &dsigma,long ipp,long ido)
00425 {
00426 long i;
00427 long nc;
00428
00429 nc=Mm->ip[ipp].ncompstr;
00430
00431 for (i=0;i<nc;i++){
00432 dsigma[i] = Mm->ip[ipp].eqother[ido+(nc+nc+nc)+i];
00433 }
00434 }