00001 #include <math.h>
00002 #include <stdlib.h>
00003 #include "microM4.h"
00004 #include "global.h"
00005 #include "intpoints.h"
00006
00007
00008
00009
00010 microM4::microM4 (void)
00011 {
00012 c1 = 6.20e-1;
00013 c2 = 2.76;
00014 c4 = 70.;
00015 c5 = 2.50;
00016 c6 = 1.30;
00017 c7 = 50.;
00018 c8 = 8.0;
00019 c9 = 1.30;
00020 c10 = 0.73;
00021 c11 = 0.2;
00022 c12 = 7000.;
00023 c13 = 0.23;
00024 c14 = 0.8;
00025 c15 = 1.0;
00026 c16 = 0.02;
00027 c17 = 0.01;
00028 c18 = 1.0;
00029 c19 = 0.4;
00030 mu = 1.0;
00031
00032 kronecker.a = new double [6];
00033 kronecker.n = 6;
00034
00035
00036 kronecker [0]= 1.;
00037 kronecker [1]= 1.;
00038 kronecker [2]= 1.;
00039 kronecker [3]= 0.;
00040 kronecker [4]= 0.;
00041 kronecker [5]= 0.;
00042 }
00043 microM4::~microM4 (void)
00044 {
00045 }
00046
00047 void microM4::read(XFILE *in)
00048 {
00049 xfscanf(in,"%ld %lf %lf %lf %lf %lf %lf %lf %lf",
00050 &numberOfMicroplanes,&k1,&k2,&k3,&k4,&c3,&c20,&e,&nu);
00051
00052 ev = e/(1-2*nu);
00053 ed = 5*e/(2+3*mu)/(1+nu);
00054 et = mu*ed;
00055
00056 initializeData (numberOfMicroplanes);
00057 }
00058
00059 void microM4::initializeData (long numberOfMicroplanes )
00060 {
00061
00062 long i,i4,mPlane;
00063 long ij[6][2] = {{1,1},{2,2},{3,3},{2,3},{3,1},{1,2}};
00064 matrix microplaneNormals(numberOfMicroplanes,3);
00065
00066 projL.a = new double [numberOfMicroplanes*6];
00067 projL.m = numberOfMicroplanes;
00068 projL.n = 6;
00069
00070 projM.a = new double [numberOfMicroplanes*6];
00071 projM.m = numberOfMicroplanes;
00072 projM.n = 6;
00073
00074 projN.a = new double [numberOfMicroplanes*6];
00075 projN.m = numberOfMicroplanes;
00076 projN.n = 6;
00077
00078 microplaneWeights.a = new double [numberOfMicroplanes];
00079 microplaneWeights.n = numberOfMicroplanes;
00080
00081
00082
00083
00084 if (numberOfMicroplanes == 28) {
00085 for (i=0; i<4; i++) {
00086 microplaneWeights[i] = 0.0160714276E0;
00087 microplaneWeights[i+4] = 0.0204744730E0;
00088 microplaneWeights[i+8] = 0.0204744730E0;
00089 microplaneWeights[i+12] = 0.0204744730E0;
00090 microplaneWeights[i+16] = 0.0158350505E0;
00091 microplaneWeights[i+20] = 0.0158350505E0;
00092 microplaneWeights[i+24] = 0.0158350505E0;
00093 }
00094 double help[7][3]={{0.577350259E0, 0.577350259E0, 0.577350259E0},
00095 {0.935113132E0, 0.250562787E0, 0.250562787E0},
00096 {0.250562787E0, 0.935113132E0, 0.250562787E0},
00097 {0.250562787E0, 0.250562787E0, 0.935113132E0},
00098 {0.186156720E0, 0.694746614E0, 0.694746614E0},
00099 {0.694746614E0, 0.186156720E0, 0.694746614E0},
00100 {0.694746614E0, 0.694746614E0, 0.186156720E0}};
00101 for (i=0; i< 7; i++) {
00102 i4 = i * 4;
00103 microplaneNormals[i4 ][0] = help[i][0];
00104 microplaneNormals[i4 ][1] = help[i][1];
00105 microplaneNormals[i4 ][2] = help[i][2];
00106 microplaneNormals[i4+1][0] = help[i][0];
00107 microplaneNormals[i4+1][1] = help[i][1];
00108 microplaneNormals[i4+1][2] =-help[i][2];
00109 microplaneNormals[i4+2][0] = help[i][0];
00110 microplaneNormals[i4+2][1] =-help[i][1];
00111 microplaneNormals[i4+2][2] = help[i][2];
00112 microplaneNormals[i4+3][0] = help[i][0];
00113 microplaneNormals[i4+3][1] =-help[i][1];
00114 microplaneNormals[i4+3][2] =-help[i][2];
00115 }
00116 }
00117 else if (numberOfMicroplanes == 21){
00118 for (i=0; i<3; i++) microplaneWeights[i] = 0.02652141274;
00119 for (i=3; i<9; i++) microplaneWeights[i] = 0.01993014153;
00120 for (i=9; i<21; i++) microplaneWeights[i] = 0.02507124272;
00121
00122 microplaneNormals[0][0] = microplaneNormals[1][1] = microplaneNormals[2][2] = 1. ;
00123 microplaneNormals[0][1] = microplaneNormals[0][2] = microplaneNormals[1][0] =
00124 microplaneNormals[1][2] = microplaneNormals[2][0] = microplaneNormals[2][1] = 0.;
00125 microplaneNormals[3][0] = microplaneNormals[3][1] = microplaneNormals[4][0] =
00126 microplaneNormals[5][0] = microplaneNormals[5][2] = microplaneNormals[6][0] =
00127 microplaneNormals[7][1] = microplaneNormals[7][2] = microplaneNormals[8][1] =
00128 0.7071067812;
00129 microplaneNormals[4][1] = microplaneNormals[6][2] = microplaneNormals[8][2] =
00130 -0.7071067812;
00131 microplaneNormals[3][2] = microplaneNormals[4][2] = microplaneNormals[5][1] =
00132 microplaneNormals[6][1] = microplaneNormals[7][0] = microplaneNormals[8][0] = 0.;
00133
00134
00135 double help [3][3] = {{0.3879072746, 0.3879072746, 0.8360956240},
00136 {0.3879072746, 0.8360956240, 0.3879072746},
00137 {0.8360956240, 0.3879072746, 0.3879072746}};
00138
00139 for(i=0; i<3;i++){
00140 i4 = i * 4;
00141 microplaneNormals[9+i4][0] = help [i][0];
00142 microplaneNormals[9+i4][1] = help [i][1];
00143 microplaneNormals[9+i4][2] = help [i][2];
00144
00145 microplaneNormals[10+i4][0] = help [i][0];
00146 microplaneNormals[10+i4][1] = help [i][1];
00147 microplaneNormals[10+i4][2] = -help [i][2];
00148
00149 microplaneNormals[11+i4][0] = help [i][0];
00150 microplaneNormals[11+i4][1] = -help [i][1];
00151 microplaneNormals[11+i4][2] = help [i][2];
00152
00153 microplaneNormals[12+i4][0] = help [i][0];
00154 microplaneNormals[12+i4][1] = -help [i][1];
00155 microplaneNormals[12+i4][2] = -help [i][2];
00156 }
00157 }
00158 else if(numberOfMicroplanes == 61){
00159
00160 double help [61][4]={
00161 {1.000000000000, 0.000000000000, 0.000000000000, 0.00795844204678},
00162 {0.745355992500, 0.0 , 0.666666666667, 0.00795844204678},
00163 {0.745355992500,-0.577350269190,-0.333333333333, 0.00795844204678},
00164 {0.745355992500, 0.577350269190,-0.333333333333, 0.00795844204678},
00165 {0.333333333333, 0.577350269190, 0.745355992500, 0.00795844204678},
00166 {0.333333333333,-0.577350269190, 0.745355992500, 0.00795844204678},
00167 {0.333333333333,-0.934172358963, 0.127322003750, 0.00795844204678},
00168 {0.333333333333,-0.356822089773,-0.872677996250, 0.00795844204678},
00169 {0.333333333333, 0.356822089773,-0.872677996250, 0.00795844204678},
00170 {0.333333333333, 0.934172358963, 0.127322003750, 0.00795844204678},
00171 {0.794654472292,-0.525731112119, 0.303530999103, 0.0105155242892},
00172 {0.794654472292, 0.0 ,-0.607061998207, 0.0105155242892},
00173 {0.794654472292, 0.525731112119, 0.303530999103, 0.0105155242892},
00174 {0.187592474085, 0.0 , 0.982246946377, 0.0105155242892},
00175 {0.187592474085,-0.850650808352,-0.491123473188, 0.0105155242892},
00176 {0.187592474085, 0.850650808352,-0.491123473188, 0.0105155242892},
00177 {0.934172358963, 0.0 , 0.356822089773, 0.0100119364272},
00178 {0.934172358963,-0.309016994375,-0.178411044887, 0.0100119364272},
00179 {0.934172358963, 0.309016994375,-0.178411044887, 0.0100119364272},
00180 {0.577350269190, 0.309016994375, 0.755761314076, 0.0100119364272},
00181 {0.577350269190,-0.309016994375, 0.755761314076, 0.0100119364272},
00182 {0.577350269190,-0.809016994375,-0.110264089708, 0.0100119364272},
00183 {0.577350269190,-0.5 ,-0.645497224368, 0.0100119364272},
00184 {0.577350269190, 0.5 ,-0.645497224368, 0.0100119364272},
00185 {0.577350269190, 0.809016994375,-0.110264089708, 0.0100119364272},
00186 {0.356822089773,-0.809016994375, 0.467086179481, 0.0100119364272},
00187 {0.356822089773, 0.0 ,-0.934172358963, 0.0100119364272},
00188 {0.356822089773, 0.809016994375, 0.467086179481, 0.0100119364272},
00189 {0.0 , 0.5 , 0.866025403784, 0.0100119364272},
00190 {0.0 ,-1.0 , 0.0 , 0.0100119364272},
00191 {0.0 , 0.5 ,-0.866025403784, 0.0100119364272},
00192 {0.947273580412,-0.277496978165, 0.160212955043, 0.00690477957966},
00193 {0.812864676392,-0.277496978165, 0.512100034157, 0.00690477957966},
00194 {0.595386501297,-0.582240127941, 0.553634669695, 0.00690477957966},
00195 {0.595386501297,-0.770581752342, 0.227417407053, 0.00690477957966},
00196 {0.812864676392,-0.582240127941,-0.015730584514, 0.00690477957966},
00197 {0.492438766306,-0.753742692223,-0.435173546254, 0.00690477957966},
00198 {0.274960591212,-0.942084316623,-0.192025554687, 0.00690477957966},
00199 {-0.076926487903,-0.942084316623,-0.326434458707, 0.00690477957966},
00200 {-0.076926487903,-0.753742692223,-0.652651721349, 0.00690477957966},
00201 {0.274960591212,-0.637341166847,-0.719856173359, 0.00690477957966},
00202 {0.947273580412, 0.0 ,-0.320425910085, 0.00690477957966},
00203 {0.812864676392,-0.304743149777,-0.496369449643, 0.00690477957966},
00204 {0.595386501297,-0.188341624401,-0.781052076747, 0.00690477957966},
00205 {0.595386501297, 0.188341624401,-0.781052076747, 0.00690477957966},
00206 {0.812864676392, 0.304743149777,-0.496369449643, 0.00690477957966},
00207 {0.492438766306, 0.753742692223,-0.435173546254, 0.00690477957966},
00208 {0.274960591212, 0.637341166847,-0.719856173359, 0.00690477957966},
00209 {-0.076926487903, 0.753742692223,-0.652651721349, 0.00690477957966},
00210 {-0.076926487903, 0.942084316623,-0.326434458707, 0.00690477957966},
00211 {0.274960591212, 0.942084316623,-0.192025554687, 0.00690477957966},
00212 {0.947273580412, 0.277496978165, 0.160212955043, 0.00690477957966},
00213 {0.812864676392, 0.582240127941,-0.015730584514, 0.00690477957966},
00214 {0.595386501297, 0.770581752342, 0.227417407053, 0.00690477957966},
00215 {0.595386501297, 0.582240127941, 0.553634669695, 0.00690477957966},
00216 {0.812864676392, 0.277496978165, 0.512100034157, 0.00690477957966},
00217 {0.492438766306, 0.0 , 0.870347092509, 0.00690477957966},
00218 {0.274960591212, 0.304743149777, 0.911881728046, 0.00690477957966},
00219 {-0.076926487903, 0.188341624401, 0.979086180056, 0.00690477957966},
00220 {-0.076926487903,-0.188341624401, 0.979086180056, 0.00690477957966},
00221 {0.274960591212,-0.304743149777, 0.911881728046, 0.00690477957966}};
00222
00223 for(i=0; i<numberOfMicroplanes; i++){
00224 microplaneNormals[i][0]= help[i][0];
00225 microplaneNormals[i][1]= help[i][1];
00226 microplaneNormals[i][2]= help[i][2];
00227 microplaneWeights[i] = help[i][3];
00228 }
00229 }else{
00230 fprintf (stderr,"microM4::initializeData Unsupported number of microplanes");
00231 }
00232
00233
00234 long ii,jj;
00235 vector n(3), m(3), l(3);
00236 double aux;
00237
00238 for (mPlane=0; mPlane < numberOfMicroplanes; mPlane++) {
00239 n[0] = microplaneNormals[mPlane][0];
00240 n[1] = microplaneNormals[mPlane][1];
00241 n[2] = microplaneNormals[mPlane][2];
00242
00243 if ((mPlane+1)%3 == 0) {
00244 aux = sqrt (n[0]*n[0] + n[1]*n[1]);
00245 if(fabs(aux) > 1.0e-12){
00246 m[0] = n[1]/aux;
00247 m[1] = -n[0]/aux;
00248 m[2] = 0.0;
00249 }
00250 else {
00251 m[0] = 1.0;
00252 m[1] = 0.0;
00253 m[2] = 0.0;
00254 }
00255 } else if ((mPlane+1)%3 == 1) {
00256 aux = sqrt (n[1]*n[1] + n[2]*n[2]);
00257 if(fabs(aux) > 1.0e-12){
00258 m[0] = 0.0;
00259 m[1] = n[2]/aux;
00260 m[2] = -n[1]/aux;
00261 }
00262 else {
00263 m[0] = 0.0;
00264 m[1] = 1.0;
00265 m[2] = 0.0;
00266 }
00267 } else {
00268 aux = sqrt (n[0]*n[0] + n[2]*n[2]);
00269 if(fabs(aux) > 1.0e-12){
00270 m[0] = n[2]/aux;
00271 m[1] = 0.0;
00272 m[2] = -n[0]/aux;
00273 }
00274 else {
00275 m[0] = 0.0;
00276 m[1] = 0.0;
00277 m[2] = 1.0;
00278 }
00279 }
00280
00281
00282 crprd(m,n,l);
00283
00284 for (i=0; i<6; i++) {
00285 ii = ij[i][0]-1;
00286 jj = ij[i][1]-1;
00287
00288 projN[mPlane][i] = n(ii)*n(jj);
00289 projM[mPlane][i] = 0.5*(m(ii)*n(jj) + m(jj)*n(ii)) ;
00290 projL[mPlane][i] = 0.5*(l(ii)*n(jj) + l(jj)*n(ii)) ;
00291 }
00292 }
00293 return;
00294 }
00295
00296 void microM4::matstiff (matrix &d,long ipp,long ido)
00297
00298
00299
00300
00301
00302 {
00303 Mm->elmatstiff (d,ipp);
00304 }
00305
00306 void microM4::nlstresses (long ipp, long ido)
00307 {
00308 long i, i_other,mPlaneIndex;
00309 double previousSigmaV,previousSigmaL,previousSigmaM,previousSigmaN,previousSigmaD;
00310 double sel,sem,sev,sed,f;
00311 double cv,cd,meanN=0.;
00312 double sigmaN,sigmaL,sigmaM,sigmaV,sigmaD;
00313 double epsV,depsV, epsD, epsN, depsD, depsM, depsL,depsN;
00314 vector macroStrain(6), macroStrainIncrement(6), macroSigma(6);
00315
00316
00317 for (i=0;i<6;i++){
00318 macroStrain[i] = Mm->ip[ipp].strain[i];
00319 macroStrainIncrement[i] = Mm->ip[ipp].strain[i] - Mm->ip[ipp].eqother[ido+i+numberOfMicroplanes*3+1];
00320
00321
00322 Mm->ip[ipp].other[ido+i+numberOfMicroplanes*3+1]=Mm->ip[ipp].strain[i];
00323 }
00324
00325
00326 epsV = (macroStrain[0]+macroStrain[1]+macroStrain[2])/3.0;
00327 depsV= (macroStrainIncrement[0]+macroStrainIncrement[1]+macroStrainIncrement[2])/3.0;
00328
00329
00330 previousSigmaV=Mm->ip[ipp].eqother[ido+0];
00331
00332
00333 for (mPlaneIndex=0; mPlaneIndex < numberOfMicroplanes; mPlaneIndex++) {
00334
00335
00336 epsN=0., depsL=0., depsM=0.; depsN=0.;
00337 for (i=0; i<6; i++) {
00338 epsN += projN [mPlaneIndex][i]*macroStrain[i];
00339 depsL += projL [mPlaneIndex][i]*macroStrainIncrement[i];
00340 depsM += projM [mPlaneIndex][i]*macroStrainIncrement[i];
00341 depsN += projN [mPlaneIndex][i]*macroStrainIncrement[i];
00342 }
00343 epsD = epsN - epsV;
00344 depsD = depsN - depsV;
00345
00346
00347 i_other=ido+(mPlaneIndex+1)*3-2;
00348 previousSigmaL=Mm->ip[ipp].eqother[i_other];
00349 previousSigmaM=Mm->ip[ipp].eqother[i_other+1];
00350 previousSigmaN=Mm->ip[ipp].eqother[i_other+2];
00351 previousSigmaD=previousSigmaN-previousSigmaV;
00352
00353
00354
00355
00356 cv = ev;
00357 cd = ed;
00358
00359
00360 sev=previousSigmaV + cv*depsV;
00361 sigmaV=minim( maxim(sev,FVminus(epsV)), FVplus(epsV));
00362
00363 sed=previousSigmaD + cd*depsD;
00364 sigmaD=minim(maxim(sed,FDminus(epsD)),FDplus(epsD));
00365
00366 sigmaN=minim( (sigmaV+sigmaD) ,FN(epsN,previousSigmaV));
00367
00368 sel=previousSigmaL + et*depsL;
00369 sem=previousSigmaM + et*depsM;
00370 f=FT(sigmaN,epsV);
00371 if (sel>f) sigmaL=f; else if (sel<-f) sigmaL=-f; else sigmaL=sel;
00372 if (sem>f) sigmaM=f; else if (sem<-f) sigmaM=-f; else sigmaM=sem;
00373
00374
00375 meanN += sigmaN* microplaneWeights[mPlaneIndex]*6.0;
00376
00377
00378 Mm->ip[ipp].other[i_other]=sigmaL;
00379 Mm->ip[ipp].other[i_other+1]=sigmaM;
00380 Mm->ip[ipp].other[i_other+2]=sigmaN;
00381
00382 }
00383
00384
00385
00386
00387 if (sigmaV > (meanN/3.)) {sigmaV = meanN/3.;}
00388
00389
00390 Mm->ip[ipp].other[ido+0]=sigmaV;
00391
00392 fillv(0,macroSigma);
00393
00394
00395 for (mPlaneIndex=0; mPlaneIndex < numberOfMicroplanes; mPlaneIndex++) {
00396 i_other=ido+(mPlaneIndex+1)*3-2;
00397 sigmaL=Mm->ip[ipp].other[i_other];
00398 sigmaM=Mm->ip[ipp].other[i_other+1];
00399 sigmaN=Mm->ip[ipp].other[i_other+2];
00400
00401 sigmaD = sigmaN - sigmaV;
00402
00403
00404 for (i=0; i<6; i++) {
00405 macroSigma[i] +=((projN[mPlaneIndex][i]-kronecker[i]/3.)*sigmaD +
00406 projL[mPlaneIndex][i]*sigmaL +
00407 projM[mPlaneIndex][i]*sigmaM)* microplaneWeights[mPlaneIndex]*6.0;
00408
00409
00410 }
00411
00412 }
00413
00414
00415 macroSigma[0]+=sigmaV;
00416 macroSigma[1]+=sigmaV;
00417 macroSigma[2]+=sigmaV;
00418
00419
00420 Mm->storestress(0,ipp,macroSigma);
00421
00422 }
00423
00424
00425 void microM4::nonloc_nlstresses (long ipp, long ido)
00426 {
00427 long i, j;
00428 matrix e(6,6);
00429 double sigma_elast;
00430 vector sigma_nonloc(6);
00431
00432
00433 Mm->elmatstiff(e,ipp);
00434
00435 for (i=0;i<6;i++){
00436 sigma_elast=0;
00437 for (j=0;j<6;j++){
00438 sigma_elast += e[i][j] * (Mm->ip[ipp].strain[j]);
00439 }
00440
00441 sigma_nonloc[i]=Mm->ip[ipp].nonloc[i] + sigma_elast;
00442 }
00443
00444
00445 Mm->storestress(0,ipp,sigma_nonloc);
00446 }
00447
00448 void microM4::updateval(long ipp, long im, long ido)
00449 {
00450 long i,n = Mm->givencompother(ipp, im);
00451
00452 for (i=0;i<n;i++){
00453 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00454 }
00455 }
00456
00457 void microM4::changeparam (atsel &atm,vector &val)
00458 {
00459 long i;
00460
00461 for (i=0;i<atm.num;i++){
00462 switch (atm.atrib[i]){
00463 case 0:{
00464 k1=val[i];
00465 break;
00466 }
00467 case 1:{
00468 k2=val[i];
00469 break;
00470 }
00471 case 2:{
00472 k3=val[i];
00473 break;
00474 }
00475 case 3:{
00476 k4=val[i];
00477 break;
00478 }
00479 case 4:{
00480 c3=val[i];
00481 break;
00482 }
00483 case 5:{
00484 c20=val[i];
00485 break;
00486 }
00487 case 6:{
00488 e=val[i];
00489 break;
00490 }
00491 case 7:{
00492 nu=val[i];
00493 break;
00494 }
00495 default:{
00496 fprintf (stderr,"\n\n wrong number of atribute in function changeparam (%s, line %d).\n",__FILE__,__LINE__);
00497 }
00498 }
00499 }
00500
00501 ev = e/(1-2*nu);
00502 ed = 5*e/(2+3*mu)/(1+nu);
00503 et = mu*ed;
00504 }
00505
00506 inline double
00507 microM4::macbra(double x)
00508 {
00509 return(maxim(x,0.0));
00510 }
00511
00512 inline double
00513 microM4::FVplus (double epsV)
00514
00515 {
00516 return(ev*k1*c13/(1+(c14/k1)*macbra(epsV-c13*c15*k1)));
00517 }
00518
00519 inline double
00520 microM4::FVminus (double epsV)
00521
00522 {
00523 return(-e*k1*k3*exp(-epsV/(k1*k4)));
00524 }
00525
00526 inline double
00527 microM4::FDminus(double epsD)
00528
00529 {
00530 double a;
00531 a=macbra(-epsD-c8*c9*k1)/(k1*c7);
00532 return(-e*k1*c8/(1+a*a));
00533 }
00534
00535 inline double
00536 microM4::FDplus(double epsD)
00537
00538 {
00539 double a;
00540 a=(macbra(epsD-c5*c6*k1)/(k1*c20*c7));
00541 return (e*k1*c5/(1+a*a));
00542 }
00543
00544 inline double
00545 microM4::FN(double epsN,double sigmaV)
00546
00547 {
00548 return(e*k1*c1*exp(-macbra(epsN-c1*c2*k1)/(k1*c3+macbra(-c4*(sigmaV/ev)))));
00549 }
00550
00551 inline double
00552 microM4::FT (double sigmaN,double epsV)
00553
00554 {
00555 double a, sn0;
00556 sn0= et*k1*c11/(1+c12*fabs(epsV));
00557 a=macbra(-sigmaN+sn0);
00558 return(et*k1*k2*c10*a/(et*k1*k2+c10*a));
00559 }
00560
00561 inline double
00562 microM4::maxim (double a,double b)
00563 {
00564 return(a>b ? a:b);
00565 }
00566
00567 inline double
00568 microM4::minim (double a,double b)
00569 {
00570 return(a<b ? a:b);
00571 }