00001 #include "microSIM.h"
00002 #include "global.h"
00003 #include <math.h>
00004 #include <stdlib.h>
00005
00006
00007
00008 microSIM::microSIM (void)
00009 {
00010 c1 = 6.20e-1;
00011 c2 = 2.76;
00012 c4 = 70.;
00013 c5 = 2.50;
00014 c6 = 1.30;
00015 c7 = 50.;
00016 c8 = 8.0;
00017 c9 = 1.30;
00018 c10 = 0.73;
00019 c11 = 0.2;
00020 c12 = 7000.;
00021 c13 = 0.23;
00022 c14 = 0.8;
00023 c15 = 1.0;
00024 c16 = 0.02;
00025 c17 = 0.01;
00026 c18 = 1.0;
00027 c19 = 0.4;
00028 mu = 1.0;
00029
00030 kronecker.a = new double [6];
00031 kronecker.n = 6;
00032
00033
00034 kronecker [0]= 1.;
00035 kronecker [1]= 1.;
00036 kronecker [2]= 1.;
00037 kronecker [3]= 0.;
00038 kronecker [4]= 0.;
00039 kronecker [5]= 0.;
00040 }
00041 microSIM::~microSIM (void)
00042 {
00043 }
00044
00045 void microSIM::read(XFILE *in)
00046 {
00047 xfscanf(in,"%ld %lf %lf %lf %lf %lf %lf %lf %lf",
00048 &numberOfMicroplanes,&k1,&k2,&k3,&k4,&c3,&c20,&e,&nu);
00049
00050 ev = e/(1-2*nu);
00051 ed = 5*e/(2+3*mu)/(1+nu);
00052 et = mu*ed;
00053
00054 initializeData (numberOfMicroplanes);
00055 }
00056
00057 void microSIM::initializeData (long numberOfMicroplanes )
00058 {
00059
00060 long i,i4,mPlane;
00061 long ij[6][2] = {{1,1},{2,2},{3,3},{2,3},{3,1},{1,2}};
00062 matrix microplaneNormals(numberOfMicroplanes,3);
00063
00064 projN.a = new double [numberOfMicroplanes*6];
00065 projN.m = numberOfMicroplanes;
00066 projN.n = 6;
00067
00068 microplaneWeights.a = new double [numberOfMicroplanes];
00069 microplaneWeights.n = numberOfMicroplanes;
00070
00071
00072
00073
00074 if (numberOfMicroplanes == 28) {
00075 for (i=0; i<4; i++) {
00076 microplaneWeights[i] = 0.0160714276E0;
00077 microplaneWeights[i+4] = 0.0204744730E0;
00078 microplaneWeights[i+8] = 0.0204744730E0;
00079 microplaneWeights[i+12] = 0.0204744730E0;
00080 microplaneWeights[i+16] = 0.0158350505E0;
00081 microplaneWeights[i+20] = 0.0158350505E0;
00082 microplaneWeights[i+24] = 0.0158350505E0;
00083 }
00084 double help[7][3]={{0.577350259E0, 0.577350259E0, 0.577350259E0},
00085 {0.935113132E0, 0.250562787E0, 0.250562787E0},
00086 {0.250562787E0, 0.935113132E0, 0.250562787E0},
00087 {0.250562787E0, 0.250562787E0, 0.935113132E0},
00088 {0.186156720E0, 0.694746614E0, 0.694746614E0},
00089 {0.694746614E0, 0.186156720E0, 0.694746614E0},
00090 {0.694746614E0, 0.694746614E0, 0.186156720E0}};
00091 for (i=0; i< 7; i++) {
00092 i4 = i * 4;
00093 microplaneNormals[i4 ][0] = help[i][0];
00094 microplaneNormals[i4 ][1] = help[i][1];
00095 microplaneNormals[i4 ][2] = help[i][2];
00096 microplaneNormals[i4+1][0] = help[i][0];
00097 microplaneNormals[i4+1][1] = help[i][1];
00098 microplaneNormals[i4+1][2] =-help[i][2];
00099 microplaneNormals[i4+2][0] = help[i][0];
00100 microplaneNormals[i4+2][1] =-help[i][1];
00101 microplaneNormals[i4+2][2] = help[i][2];
00102 microplaneNormals[i4+3][0] = help[i][0];
00103 microplaneNormals[i4+3][1] =-help[i][1];
00104 microplaneNormals[i4+3][2] =-help[i][2];
00105 }
00106 }
00107 else if (numberOfMicroplanes == 21){
00108 for (i=0; i<3; i++) microplaneWeights[i] = 0.02652141274;
00109 for (i=3; i<9; i++) microplaneWeights[i] = 0.01993014153;
00110 for (i=9; i<21; i++) microplaneWeights[i] = 0.02507124272;
00111
00112 microplaneNormals[0][0] = microplaneNormals[1][1] = microplaneNormals[2][2] = 1. ;
00113 microplaneNormals[0][1] = microplaneNormals[0][2] = microplaneNormals[1][0] =
00114 microplaneNormals[1][2] = microplaneNormals[2][0] = microplaneNormals[2][1] = 0.;
00115 microplaneNormals[3][0] = microplaneNormals[3][1] = microplaneNormals[4][0] =
00116 microplaneNormals[5][0] = microplaneNormals[5][2] = microplaneNormals[6][0] =
00117 microplaneNormals[7][1] = microplaneNormals[7][2] = microplaneNormals[8][1] =
00118 0.7071067812;
00119 microplaneNormals[4][1] = microplaneNormals[6][2] = microplaneNormals[8][2] =
00120 -0.7071067812;
00121 microplaneNormals[3][2] = microplaneNormals[4][2] = microplaneNormals[5][1] =
00122 microplaneNormals[6][1] = microplaneNormals[7][0] = microplaneNormals[8][0] = 0.;
00123
00124
00125 double help [3][3] = {{0.3879072746, 0.3879072746, 0.8360956240},
00126 {0.3879072746, 0.8360956240, 0.3879072746},
00127 {0.8360956240, 0.3879072746, 0.3879072746}};
00128
00129 for(i=0; i<3;i++){
00130 i4 = i * 4;
00131 microplaneNormals[9+i4][0] = help [i][0];
00132 microplaneNormals[9+i4][1] = help [i][1];
00133 microplaneNormals[9+i4][2] = help [i][2];
00134
00135 microplaneNormals[10+i4][0] = help [i][0];
00136 microplaneNormals[10+i4][1] = help [i][1];
00137 microplaneNormals[10+i4][2] = -help [i][2];
00138
00139 microplaneNormals[11+i4][0] = help [i][0];
00140 microplaneNormals[11+i4][1] = -help [i][1];
00141 microplaneNormals[11+i4][2] = help [i][2];
00142
00143 microplaneNormals[12+i4][0] = help [i][0];
00144 microplaneNormals[12+i4][1] = -help [i][1];
00145 microplaneNormals[12+i4][2] = -help [i][2];
00146 }
00147 }
00148 else if(numberOfMicroplanes == 61){
00149
00150 double help [61][4]={
00151 {1.000000000000, 0.000000000000, 0.000000000000, 0.00795844204678},
00152 {0.745355992500, 0.0 , 0.666666666667, 0.00795844204678},
00153 {0.745355992500,-0.577350269190,-0.333333333333, 0.00795844204678},
00154 {0.745355992500, 0.577350269190,-0.333333333333, 0.00795844204678},
00155 {0.333333333333, 0.577350269190, 0.745355992500, 0.00795844204678},
00156 {0.333333333333,-0.577350269190, 0.745355992500, 0.00795844204678},
00157 {0.333333333333,-0.934172358963, 0.127322003750, 0.00795844204678},
00158 {0.333333333333,-0.356822089773,-0.872677996250, 0.00795844204678},
00159 {0.333333333333, 0.356822089773,-0.872677996250, 0.00795844204678},
00160 {0.333333333333, 0.934172358963, 0.127322003750, 0.00795844204678},
00161 {0.794654472292,-0.525731112119, 0.303530999103, 0.0105155242892},
00162 {0.794654472292, 0.0 ,-0.607061998207, 0.0105155242892},
00163 {0.794654472292, 0.525731112119, 0.303530999103, 0.0105155242892},
00164 {0.187592474085, 0.0 , 0.982246946377, 0.0105155242892},
00165 {0.187592474085,-0.850650808352,-0.491123473188, 0.0105155242892},
00166 {0.187592474085, 0.850650808352,-0.491123473188, 0.0105155242892},
00167 {0.934172358963, 0.0 , 0.356822089773, 0.0100119364272},
00168 {0.934172358963,-0.309016994375,-0.178411044887, 0.0100119364272},
00169 {0.934172358963, 0.309016994375,-0.178411044887, 0.0100119364272},
00170 {0.577350269190, 0.309016994375, 0.755761314076, 0.0100119364272},
00171 {0.577350269190,-0.309016994375, 0.755761314076, 0.0100119364272},
00172 {0.577350269190,-0.809016994375,-0.110264089708, 0.0100119364272},
00173 {0.577350269190,-0.5 ,-0.645497224368, 0.0100119364272},
00174 {0.577350269190, 0.5 ,-0.645497224368, 0.0100119364272},
00175 {0.577350269190, 0.809016994375,-0.110264089708, 0.0100119364272},
00176 {0.356822089773,-0.809016994375, 0.467086179481, 0.0100119364272},
00177 {0.356822089773, 0.0 ,-0.934172358963, 0.0100119364272},
00178 {0.356822089773, 0.809016994375, 0.467086179481, 0.0100119364272},
00179 {0.0 , 0.5 , 0.866025403784, 0.0100119364272},
00180 {0.0 ,-1.0 , 0.0 , 0.0100119364272},
00181 {0.0 , 0.5 ,-0.866025403784, 0.0100119364272},
00182 {0.947273580412,-0.277496978165, 0.160212955043, 0.00690477957966},
00183 {0.812864676392,-0.277496978165, 0.512100034157, 0.00690477957966},
00184 {0.595386501297,-0.582240127941, 0.553634669695, 0.00690477957966},
00185 {0.595386501297,-0.770581752342, 0.227417407053, 0.00690477957966},
00186 {0.812864676392,-0.582240127941,-0.015730584514, 0.00690477957966},
00187 {0.492438766306,-0.753742692223,-0.435173546254, 0.00690477957966},
00188 {0.274960591212,-0.942084316623,-0.192025554687, 0.00690477957966},
00189 {-0.076926487903,-0.942084316623,-0.326434458707, 0.00690477957966},
00190 {-0.076926487903,-0.753742692223,-0.652651721349, 0.00690477957966},
00191 {0.274960591212,-0.637341166847,-0.719856173359, 0.00690477957966},
00192 {0.947273580412, 0.0 ,-0.320425910085, 0.00690477957966},
00193 {0.812864676392,-0.304743149777,-0.496369449643, 0.00690477957966},
00194 {0.595386501297,-0.188341624401,-0.781052076747, 0.00690477957966},
00195 {0.595386501297, 0.188341624401,-0.781052076747, 0.00690477957966},
00196 {0.812864676392, 0.304743149777,-0.496369449643, 0.00690477957966},
00197 {0.492438766306, 0.753742692223,-0.435173546254, 0.00690477957966},
00198 {0.274960591212, 0.637341166847,-0.719856173359, 0.00690477957966},
00199 {-0.076926487903, 0.753742692223,-0.652651721349, 0.00690477957966},
00200 {-0.076926487903, 0.942084316623,-0.326434458707, 0.00690477957966},
00201 {0.274960591212, 0.942084316623,-0.192025554687, 0.00690477957966},
00202 {0.947273580412, 0.277496978165, 0.160212955043, 0.00690477957966},
00203 {0.812864676392, 0.582240127941,-0.015730584514, 0.00690477957966},
00204 {0.595386501297, 0.770581752342, 0.227417407053, 0.00690477957966},
00205 {0.595386501297, 0.582240127941, 0.553634669695, 0.00690477957966},
00206 {0.812864676392, 0.277496978165, 0.512100034157, 0.00690477957966},
00207 {0.492438766306, 0.0 , 0.870347092509, 0.00690477957966},
00208 {0.274960591212, 0.304743149777, 0.911881728046, 0.00690477957966},
00209 {-0.076926487903, 0.188341624401, 0.979086180056, 0.00690477957966},
00210 {-0.076926487903,-0.188341624401, 0.979086180056, 0.00690477957966},
00211 {0.274960591212,-0.304743149777, 0.911881728046, 0.00690477957966}};
00212
00213 for(i=0; i<numberOfMicroplanes; i++){
00214 microplaneNormals[i][0]= help[i][0];
00215 microplaneNormals[i][1]= help[i][1];
00216 microplaneNormals[i][2]= help[i][2];
00217 microplaneWeights[i] = help[i][3];
00218 }
00219 }else{
00220 fprintf (stderr,"microSIM::initializeData Unsupported number of microplanes");
00221 }
00222
00223
00224 long ii,jj;
00225 vector n(3), m(3), l(3);
00226
00227 for (mPlane=0; mPlane < numberOfMicroplanes; mPlane++) {
00228 n[0] = microplaneNormals[mPlane][0];
00229 n[1] = microplaneNormals[mPlane][1];
00230 n[2] = microplaneNormals[mPlane][2];
00231
00232 for (i=0; i<6; i++) {
00233 ii = ij[i][0]-1;
00234 jj = ij[i][1]-1;
00235
00236 projN[mPlane][i] = n(ii)*n(jj);
00237 }
00238 }
00239 return;
00240 }
00241
00242 void microSIM::matstiff (matrix &d,long ipp,long ido)
00243
00244
00245 {
00246 Mm->elmatstiff (d,ipp);
00247 }
00248
00249 void microSIM::nlstresses (long ipp,long ido)
00250 {
00251 long i, i_other,mPlaneIndex;
00252 double previousSigmaV,previousSigmaN,previousSigmaD;
00253 double sev,sed;
00254 double cv,cd,meanN=0.;
00255 double sigmaN,sigmaV,sigmaD;
00256 double epsV,depsV, epsD, epsN, depsD,depsN;
00257 vector macroStrain(6), macroStrainIncrement(6), macroSigma(6);
00258
00259
00260 for (i=0;i<6;i++){
00261 macroStrain[i] = Mm->ip[ipp].strain[i];
00262 macroStrainIncrement[i] = Mm->ip[ipp].strain[i] - Mm->ip[ipp].eqother[ido+i+numberOfMicroplanes+1];
00263
00264
00265 Mm->ip[ipp].other[ido+i+numberOfMicroplanes+1]=Mm->ip[ipp].strain[i];
00266 }
00267
00268
00269 epsV = (macroStrain[0]+macroStrain[1]+macroStrain[2])/3.0;
00270 depsV= (macroStrainIncrement[0]+macroStrainIncrement[1]+macroStrainIncrement[2])/3.0;
00271
00272
00273 previousSigmaV=Mm->ip[ipp].eqother[ido+0];
00274
00275
00276 for (mPlaneIndex=0; mPlaneIndex < numberOfMicroplanes; mPlaneIndex++) {
00277
00278
00279 epsN=0., depsN=0.;
00280 for (i=0; i<6; i++) {
00281 epsN += projN [mPlaneIndex][i]*macroStrain[i];
00282 depsN += projN [mPlaneIndex][i]*macroStrainIncrement[i];
00283 }
00284 epsD = epsN - epsV;
00285 depsD = depsN - depsV;
00286
00287
00288 i_other=ido+mPlaneIndex+1;
00289 previousSigmaN=Mm->ip[ipp].eqother[i_other];
00290 previousSigmaD=previousSigmaN-previousSigmaV;
00291
00292
00293
00294
00295 cv = ev;
00296 cd = ed;
00297
00298
00299 sev=previousSigmaV + cv*depsV;
00300 sigmaV=minim( maxim(sev,FVminus(epsV)), FVplus(epsV));
00301
00302 sed=previousSigmaD + cd*depsD;
00303 sigmaD=minim(maxim(sed,FDminus(epsD)),FDplus(epsD));
00304
00305 sigmaN=minim( (sigmaV+sigmaD) ,FN(epsN,previousSigmaV));
00306
00307
00308 meanN += sigmaN* microplaneWeights[mPlaneIndex]*6.0;
00309
00310
00311 Mm->ip[ipp].other[i_other]=sigmaN;
00312
00313 }
00314
00315
00316
00317
00318 if (sigmaV > (meanN/3.)) {sigmaV = meanN/3.;}
00319
00320
00321 Mm->ip[ipp].other[ido+0]=sigmaV;
00322
00323 fillv(0,macroSigma);
00324
00325
00326 for (mPlaneIndex=0; mPlaneIndex < numberOfMicroplanes; mPlaneIndex++) {
00327 i_other=ido+mPlaneIndex+1;
00328 sigmaN=Mm->ip[ipp].other[i_other];
00329
00330 sigmaD = sigmaN - sigmaV;
00331
00332
00333 for (i=0; i<6; i++) {
00334 macroSigma[i]+=(projN[mPlaneIndex][i]-kronecker[i]/3.)*sigmaD*microplaneWeights[mPlaneIndex]*6.0;
00335 }
00336
00337 }
00338
00339
00340 macroSigma[0]+=sigmaV;
00341 macroSigma[1]+=sigmaV;
00342 macroSigma[2]+=sigmaV;
00343
00344
00345 Mm->storestress(0,ipp,macroSigma);
00346
00347 }
00348
00349 void microSIM::updateval(long ipp,long im,long ido)
00350 {
00351 long i,n = Mm->givencompother(ipp, im);
00352
00353 for (i=0;i<n;i++){
00354 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00355 }
00356 }
00357
00358 inline double
00359 microSIM::macbra(double x)
00360 {
00361 return(maxim(x,0.0));
00362 }
00363
00364 inline double
00365 microSIM::FVplus (double epsV)
00366
00367 {
00368 return(ev*k1*c13/(1+(c14/k1)*macbra(epsV-c13*c15*k1)));
00369 }
00370
00371 inline double
00372 microSIM::FVminus (double epsV)
00373
00374 {
00375 return(-e*k1*k3*exp(-epsV/(k1*k4)));
00376 }
00377
00378 inline double
00379 microSIM::FDminus(double epsD)
00380
00381 {
00382 double a;
00383 a=macbra(-epsD-c8*c9*k1)/(k1*c7);
00384 return(-e*k1*c8/(1+a*a));
00385 }
00386
00387 inline double
00388 microSIM::FDplus(double epsD)
00389
00390 {
00391 double a;
00392 a=(macbra(epsD-c5*c6*k1)/(k1*c20*c7));
00393 return (e*k1*c5/(1+a*a));
00394 }
00395
00396 inline double
00397 microSIM::FN(double epsN,double sigmaV)
00398
00399 {
00400 return(e*k1*c1*exp(-macbra(epsN-c1*c2*k1)/(k1*c3+macbra(-c4*(sigmaV/ev)))));
00401 }
00402
00403 inline double
00404 microSIM::maxim (double a,double b)
00405 {
00406 return(a>b ? a:b);
00407 }
00408
00409 inline double
00410 microSIM::minim (double a,double b)
00411 {
00412 return(a<b ? a:b);
00413 }