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