00001 #include "stdlib.h"
00002 #include "math.h"
00003 #include "homog.h"
00004 #include "matrix.h"
00005 #include "vector.h"
00006
00007
00008
00009
00010
00011 HomogData::HomogData(const double iE_m, const double inu_m, const double if_m,
00012 const double iE_i, const double inu_i, const double iHirsch)
00013
00014 {
00015 E_m = iE_m; nu_m = inu_m; f_m = if_m;
00016 E_i = iE_i; nu_i = inu_i; f_i = (1-if_m);
00017 H_chi=iHirsch;
00018
00019 ENuToKMu(E_m, nu_m, k_m, mu_m);
00020 ENuToKMu(E_i, nu_i, k_i, mu_i);
00021 }
00022
00023
00024 HomogData::HomogData(double *PhaseMatrix){
00025 f_m=*(PhaseMatrix+0);
00026 E_m=*(PhaseMatrix+1);
00027 nu_m=*(PhaseMatrix+2);
00028 f_i=*(PhaseMatrix+3);
00029 E_i=*(PhaseMatrix+4);
00030 nu_i=*(PhaseMatrix+5);
00031
00032 if((f_m+f_i)<0.99 || (f_m+f_i)>1.01){
00033 printf("volumetric fraction of phases 0,1 is not unity, terminating\n");
00034 exit(1);
00035 }
00036
00037 ENuToKMu(E_m, nu_m, k_m, mu_m);
00038 ENuToKMu(E_i, nu_i, k_i, mu_i);
00039 }
00040
00041
00042
00043 HomogData::HomogData()
00044 {
00045 }
00046
00047
00048 void HomogData::Voigt(void) {
00049
00050 k_hmg=k_m*f_m+k_i*f_i;
00051 mu_hmg=mu_m*f_m+mu_i*f_i;
00052
00053 KMuToENu(k_hmg, mu_hmg, E_hmg, nu_hmg);
00054 }
00055
00056
00057 void HomogData::Reuss(void){
00058
00059 k_hmg=1/(f_m/k_m+f_i/k_i);
00060 mu_hmg=1/(f_m/mu_m+f_i/mu_i);
00061 KMuToENu(k_hmg, mu_hmg, E_hmg, nu_hmg);
00062 }
00063
00064
00065
00066
00067 void HomogData::Hirsch(void){
00068 E_hmg=1/((1.-H_chi)*(f_m/E_m+f_i/E_i)+H_chi/(f_m*E_m+f_i*E_i));
00069
00070 }
00071
00072
00073
00074
00075 void HomogData::Hansen(void){
00076 E_hmg=E_i*((f_i*E_i+(1+f_m)*E_m)/((1.+f_m)*E_i+f_i*E_m));
00077
00078 }
00079
00080
00081
00082
00083 void HomogData::Counto(void){
00084
00085 E_hmg=1./((1.-sqrt(f_m))/E_i+1./(E_i*((1.-sqrt(f_m))/sqrt(f_m))+E_m));
00086
00087 }
00088
00089
00090
00091
00092
00093 void HomogData::MoriTanaka(void){
00094 double t, k_MT, mu_MT;
00095
00096 alpha_m = 3. * k_m / ( 3. * k_m + 4 * mu_m );
00097 beta_m = 6. * ( k_m + 2. * mu_m ) / ( 5. * ( 3. * k_m + 4. * mu_m ) );
00098
00099 t = 1. + alpha_m * ( k_i / k_m - 1. );
00100
00101 k_MT = f_m * k_m + f_i * k_i / t;
00102 k_MT /= f_m + f_i / t + ( 1. - f_m - f_i ) / ( 1. - alpha_m );
00103
00104 t = 1. + beta_m * ( mu_i / mu_m - 1. );
00105
00106 mu_MT = f_m * mu_m + f_i * mu_i / t;
00107 mu_MT /= f_m + f_i / t + ( 1. - f_m - f_i ) / ( 1. - beta_m );
00108
00109 KMuToENu( k_MT, mu_MT, E_hmg, nu_hmg );
00110
00111 k_hmg=k_MT;
00112 mu_hmg=mu_MT;
00113 }
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 void HomogData::MT_mtrx(double *PhaseMatrix, int NumRows, int MatrixRow){
00129 double f_r, E_r, nu_r, k_r, mu_r;
00130 double fr_tot = 0.;
00131 double nom_k_MT=0. , denom_k_MT=0. , nom_mu_MT=0., denom_mu_MT=0.;
00132
00133
00134
00135 f_m=*(PhaseMatrix+3*MatrixRow+0);
00136 E_m=*(PhaseMatrix+3*MatrixRow+1);
00137 nu_m=*(PhaseMatrix+3*MatrixRow+2);
00138
00139 ENuToKMu(E_m, nu_m, k_m, mu_m);
00140
00141
00142 alpha_m = 3. * k_m / ( 3. * k_m + 4 * mu_m );
00143 beta_m = 6. * ( k_m + 2. * mu_m ) / ( 5. * ( 3. * k_m + 4. * mu_m ) );
00144
00145
00146 for (int r=0; r<NumRows; r++){
00147 f_r=*(PhaseMatrix+3*r+0);
00148 E_r=*(PhaseMatrix+3*r+1);
00149 nu_r=*(PhaseMatrix+3*r+2);
00150
00151 fr_tot +=f_r;
00152 ENuToKMu(E_r, nu_r, k_r, mu_r);
00153
00154 nom_k_MT += f_r*k_r/(1.+alpha_m*(k_r/k_m-1.));
00155 denom_k_MT += f_r/(1.+alpha_m*(k_r/k_m-1.));
00156
00157 nom_mu_MT += f_r*mu_r/(1.+beta_m*(mu_r/mu_m-1.));
00158 denom_mu_MT += f_r/(1.+beta_m*(mu_r/mu_m-1.));
00159 }
00160
00161 if(fr_tot<0.99 || fr_tot>1.01){
00162 printf("volumetric fraction of phases 0-%d is not unity (= %lf), terminating (MT_mtrx)\n", NumRows, fr_tot);
00163 exit(1);
00164 }
00165
00166 k_hmg=nom_k_MT/denom_k_MT;
00167 mu_hmg=nom_mu_MT/denom_mu_MT;
00168
00169 KMuToENu(k_hmg, mu_hmg, E_hmg, nu_hmg);
00170
00171 }
00172
00173
00174
00175
00176
00177
00178
00179 void HomogData::SelfConsistent(void){
00180 double k_SCS, mu_SCS, t_m, t_i;
00181
00182
00183 k_SCS=10;
00184 mu_SCS=.3;
00185
00186
00187 for(int i=1; i<100; i++){
00188
00189 alpha_m=3.*k_SCS/(3.*k_SCS+4*mu_SCS);
00190 beta_m=6.*(k_SCS+2.*mu_SCS)/(5.*(3.*k_SCS+4.*mu_SCS));
00191
00192 t_m=1.+alpha_m*(k_m/k_SCS-1.);
00193 t_i=1.+alpha_m*(k_i/k_SCS-1.);
00194
00195 k_SCS=f_m*k_m/t_m+f_i*k_i/t_i;
00196 k_SCS/=(f_m/t_m+f_i/t_i);
00197
00198
00199 t_m=1.+beta_m*(mu_m/mu_SCS-1.);
00200 t_i=1.+beta_m*(mu_i/mu_SCS-1.);
00201
00202 mu_SCS=f_m*mu_m/t_m+f_i*mu_i/t_i;
00203 mu_SCS/=(f_m/t_m+f_i/t_i);
00204 }
00205
00206 KMuToENu( k_SCS, mu_SCS, E_hmg, nu_hmg );
00207
00208 k_hmg=k_SCS;
00209 mu_hmg=mu_SCS;
00210
00211 temp1=alpha_m;
00212 temp2=beta_m;
00213 temp3=k_SCS;
00214 temp4=mu_SCS;
00215 }
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225 void HomogData::SCS(double *PhaseMatrix, int NumRows){
00226 double f_r, E_r, nu_r, k_r, mu_r;
00227 double k_SCS, mu_SCS, nom_k_SCS, denom_k_SCS, nom_mu_SCS, denom_mu_SCS;
00228 double fr_tot;
00229
00230
00231 k_SCS=10.;
00232 mu_SCS=.3;
00233
00234
00235 for(int i=1; i<100; i++){
00236 fr_tot = 0.;
00237 nom_k_SCS=0;
00238 denom_k_SCS=0;
00239 nom_mu_SCS=0;
00240 denom_mu_SCS=0;
00241 for (int r=0; r<NumRows; r++){
00242 f_r=*(PhaseMatrix+3*r+0);
00243 E_r=*(PhaseMatrix+3*r+1);
00244 nu_r=*(PhaseMatrix+3*r+2);
00245
00246 fr_tot += f_r;
00247
00248 ENuToKMu(E_r, nu_r, k_r, mu_r);
00249
00250 alpha_m=3.*k_SCS/(3.*k_SCS+4*mu_SCS);
00251 beta_m=6.*(k_SCS+2.*mu_SCS)/(5.*(3.*k_SCS+4.*mu_SCS));
00252
00253 nom_k_SCS+=f_r*k_r/(1+alpha_m*(k_r/k_SCS-1));
00254 denom_k_SCS+=f_r/(1+alpha_m*(k_r/k_SCS-1));
00255
00256 nom_mu_SCS+=f_r*mu_r/(1+beta_m*(mu_r/mu_SCS-1));
00257 denom_mu_SCS+=f_r/(1+beta_m*(mu_r/mu_SCS-1));
00258 }
00259 k_SCS=nom_k_SCS/denom_k_SCS;
00260 mu_SCS=nom_mu_SCS/denom_mu_SCS;
00261 }
00262
00263 temp3=k_SCS;
00264 temp4=mu_SCS;
00265
00266 if(fr_tot<0.99 || fr_tot>1.01){
00267 printf("volumetric fraction of phases 0-%d is not unity (= %lf), terminating (SCS)\n", NumRows, fr_tot);
00268 exit(1);
00269 }
00270
00271 KMuToENu(k_SCS, mu_SCS, E_hmg, nu_hmg );
00272
00273 k_hmg=k_SCS;
00274 mu_hmg=mu_SCS;
00275 }
00276
00277
00278
00279 void HomogData::Voigt(double *PhaseMatrix, int NumRows){
00280
00281 E_hmg = 0.;
00282 nu_hmg = 0.;
00283 for (int r=0; r<NumRows; r++){
00284 E_hmg+=(*(PhaseMatrix+3*r+0))*(*(PhaseMatrix+3*r+1));
00285 nu_hmg+=(*(PhaseMatrix+3*r+0))*(*(PhaseMatrix+3*r+2));
00286 }
00287 ENuToKMu(E_hmg, nu_hmg, k_hmg, mu_hmg);
00288 }
00289
00290
00291 void HomogData::Reuss(double *PhaseMatrix, int NumRows){
00292
00293 E_hmg = 0.;
00294 nu_hmg = 0.;
00295 for (int r=0; r<NumRows; r++){
00296 E_hmg+=(*(PhaseMatrix+3*r+0))/(*(PhaseMatrix+3*r+1));
00297 nu_hmg+=(*(PhaseMatrix+3*r+0))/(*(PhaseMatrix+3*r+2));
00298 }
00299 E_hmg=1/E_hmg;
00300 nu_hmg=1/nu_hmg;
00301 ENuToKMu(E_hmg, nu_hmg, k_hmg, mu_hmg);
00302 }
00303
00304
00305
00306
00307
00308
00309
00310
00311 void HomogData::HashinShtrikman(void){
00312 double k_HS, mu_HS;
00313
00314 if (((k_m < k_i) && (mu_m > mu_i)) || ((k_m > k_i) && (mu_m < mu_i))){
00315 printf("Non well ordered two-phase material, use Walpole bounds\n");
00316 printf("k[0]=%f k[1]=%f mu[0]=%f mu[1]=%f\n", k_m, k_i, mu_m, mu_i);
00317 exit(0);
00318 }
00319
00320
00321 k_HS=k_m+f_i/(1./(k_i-k_m)+3.*f_m/(3.*k_m+4.*mu_m));
00322 mu_HS=mu_m+f_i/(1./(mu_i-mu_m)+6.*f_m*(k_m+2.*mu_m)/(5.*mu_m*(3.*k_m+4.*mu_m)));
00323
00324
00325
00326
00327
00328
00329 KMuToENu(k_HS, mu_HS, E_hmg, nu_hmg);
00330
00331 k_hmg=k_HS;
00332 mu_hmg=mu_HS;
00333
00334
00335 k_HS=k_i+f_m/(1./(k_m-k_i)+3.*f_i/(3.*k_i+4.*mu_i));
00336
00337 mu_HS=mu_i+f_m/(1./(mu_m-mu_i)+6.*f_i*(k_i+2.*mu_i)/(5.*mu_i*(3.*k_i+4.*mu_i)));
00338
00339 KMuToENu(k_HS, mu_HS, E_hmg_2, nu_hmg_2);
00340
00341 k_hmg_2=k_HS;
00342 mu_hmg_2=mu_HS;
00343
00344
00345
00346
00347 if (E_hmg>E_hmg_2){
00348 Swap(E_hmg,E_hmg_2);
00349 }
00350 if (nu_hmg>nu_hmg_2){
00351 Swap(nu_hmg,nu_hmg_2);
00352 }
00353
00354 if (k_hmg>k_hmg_2){
00355 Swap(k_hmg,k_hmg_2);
00356 }
00357 if (mu_hmg>mu_hmg_2){
00358 Swap(mu_hmg,mu_hmg_2);
00359 }
00360 }
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371 void HomogData::WalpoleMulti(double *PhaseMatrix, int NumRows){
00372 double k_low, mu_low, k_high, mu_high;
00373 double k_min, k_max, mu_min, mu_max;
00374 double **PhaseMatrixKMu;
00375 PhaseMatrixKMu = new double*[NumRows];
00376 long i;
00377 for(i=0; i<NumRows; i++)
00378 PhaseMatrixKMu[i]= new double[3];
00379
00380 int r;
00381 for(r=0; r<NumRows; r++){
00382
00383 PhaseMatrixKMu[r][0]=*(PhaseMatrix+3*r+0);
00384
00385 ENuToKMu(*(PhaseMatrix+3*r+1), *(PhaseMatrix+3*r+2), PhaseMatrixKMu[r][1], PhaseMatrixKMu[r][2]);
00386 }
00387
00388
00389
00390 k_min=PhaseMatrixKMu[0][1];
00391 k_max=0;
00392 mu_min=PhaseMatrixKMu[0][2];
00393 mu_max=0;
00394 for(r=0; r<NumRows; r++){
00395 if(PhaseMatrixKMu[r][1]<k_min)
00396 k_min=PhaseMatrixKMu[r][1];
00397 if(PhaseMatrixKMu[r][2]<mu_min)
00398 mu_min=PhaseMatrixKMu[r][2];
00399 if(PhaseMatrixKMu[r][1]>k_max)
00400 k_max=PhaseMatrixKMu[r][1];
00401 if(PhaseMatrixKMu[r][2]>mu_max)
00402 mu_max=PhaseMatrixKMu[r][2];
00403 }
00404
00405
00406 k_low=Lambda(*PhaseMatrixKMu, NumRows, k_min);
00407 k_high=Lambda(*PhaseMatrixKMu, NumRows, k_max);
00408 mu_low=Gamma(*PhaseMatrixKMu, NumRows, Zeta(k_min, mu_min));
00409 mu_high=Gamma(*PhaseMatrixKMu, NumRows, Zeta(k_max, mu_max));
00410
00411
00412 KMuToENu(k_low, mu_low, E_hmg, nu_hmg);
00413 KMuToENu(k_high, mu_high, E_hmg_2, nu_hmg_2);
00414
00415 k_hmg=k_low;
00416 mu_hmg=mu_low;
00417 k_hmg_2=k_high;
00418 mu_hmg_2=mu_high;
00419 }
00420
00421
00422
00423 double HomogData::Lambda(double *PhaseMatrixKMu, int NumRows, double mu){
00424 double lambda=0;
00425 for (int r=0; r<NumRows; r++){
00426 lambda+=*(PhaseMatrixKMu+3*r+0)/(*(PhaseMatrixKMu+3*r+1)+4./3.*mu);
00427 }
00428 lambda=1./lambda;
00429 lambda-=4./3.*mu;
00430 return lambda;
00431 }
00432
00433
00434 double HomogData::Gamma(double *PhaseMatrixKMu, int NumRows, double k){
00435 double gamma=0;
00436 for (int r=0; r<NumRows; r++){
00437 gamma+=*(PhaseMatrixKMu+3*r+0)/(*(PhaseMatrixKMu+3*r+2)+k);
00438 }
00439 gamma=1./gamma;
00440 gamma-=k;
00441 return gamma;
00442 }
00443
00444
00445 double HomogData::Zeta(double k, double mu){
00446 return mu/6.*(9.*k+8.*mu)/(k+2.*mu);
00447 }
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479 void HomogData::HerveZaoui(double *PhaseMatrix, int NumPhases){
00480 vector r(NumPhases-1),k(NumPhases),mu(NumPhases);
00481
00482
00483 vector Q11(NumPhases-1), Q21(NumPhases-1), F(NumPhases), G(NumPhases);
00484 matrix J(2,2), Jinv(2,2), N(2,2), Nhelp(2,2), Q(2,2);
00485
00486 double lambda_0=0.735;
00487
00488 vector A(NumPhases), B(NumPhases), C(NumPhases), D(NumPhases);
00489 vector P_v(4), Phelp(4);
00490 matrix L(4,4), Linv(4,4), M(4,4), Mhelp(4,4), Z(4,4), M_test(4,4);
00491 matrix *P_arr;
00492
00493 double gamma=2.;
00494
00495 long double a,b,c,nu_n,sqr;
00496
00497 long double ak, bk, ck, dk, ek, fk, alphak, nuk, nukp1, koef;
00498 int i, phi;
00499
00500 if(NumPhases<3){
00501 printf("Number of considered phases must be at least 3 (including hom. medium)\n");
00502 exit(1);
00503 }
00504
00505 F[NumPhases-1]=lambda_0/3.;
00506 A[NumPhases-1]=gamma;
00507
00508
00509 temp1=0.;
00510 for(i=0; i<NumPhases-1;i++){
00511 temp1+=*(PhaseMatrix+3*i+0);
00512 r[i]=1.*pow(temp1,1./3.);
00513
00514 }
00515
00516 if(r[NumPhases-2]<0.99 || r[NumPhases-2]>1.01){
00517 printf("volumetric fraction of phases 0-%d is not unity, terminating\n", NumPhases-2);
00518 exit(1);
00519 }
00520
00521
00522 for(i=0; i<NumPhases;i++){
00523 ENuToKMu(*(PhaseMatrix+3*i+1), *(PhaseMatrix+3*i+2), k[i], mu[i]);
00524
00525 }
00526
00527
00528
00529
00530
00531
00532 Nhelp[0][0]=1.;
00533 Nhelp[1][1]=1.;
00534
00535
00536 for(phi=0;phi<=NumPhases-2;phi++){
00537
00538 fillJ(J,r[phi],mu,k,phi+1);
00539 invm(J, Jinv, 1.e-10);
00540
00541 fillJ(J,r[phi],mu,k,phi);
00542
00543 mxm(Jinv,J, N);
00544
00545 mxm(N,Nhelp,Q);
00546
00547 Q11[phi]=Q[0][0];
00548 Q21[phi]=Q[1][0];
00549
00550 copym(Q, Nhelp);
00551 }
00552
00553
00554 for(phi=1;phi<=NumPhases-1;phi++){
00555 F[phi]=F[NumPhases-1]*Q11[phi-1]/Q11[NumPhases-2];
00556 G[phi]=F[NumPhases-1]*Q21[phi-1]/Q11[NumPhases-2];
00557
00558 }
00559
00560
00561 for(phi=1;phi<=NumPhases-1;phi++){
00562
00563
00564 }
00565
00566
00567
00568
00569 k_hmg=(3.*k[NumPhases-2]*pow(r[NumPhases-2],3)*Q11[NumPhases-3]-
00570 4.*mu[NumPhases-2]*Q21[NumPhases-3])/
00571 (3.*(pow(r[NumPhases-2],3)*Q11[NumPhases-3]+Q21[NumPhases-3]));
00572
00573
00574 destrv(Q11);
00575 destrv(Q21);
00576 destrv(F);
00577 destrv(G);
00578
00579 destrm(J);
00580 destrm(Jinv);
00581 destrm(N);
00582 destrm(Nhelp);
00583 destrm(Q);
00584
00585
00586
00587
00588
00589
00590
00591 Mhelp[0][0]=1.;
00592 Mhelp[1][1]=1.;
00593 Mhelp[2][2]=1.;
00594 Mhelp[3][3]=1.;
00595
00596
00597 P_arr=new matrix[NumPhases-1];
00598
00599
00600 for(phi=0;phi<=NumPhases-2;phi++){
00601
00602 fillL(L,r[phi],mu,k,phi+1);
00603 invm(L, Linv, 1.e-10);
00604
00605 fillL(L,r[phi],mu,k,phi);
00606
00607 mxm(Linv,L,M);
00608
00609 nuk=*(PhaseMatrix+3*phi+2);
00610 nukp1=*(PhaseMatrix+3*(phi+1)+2);
00611
00612
00613
00614 ak=mu[phi]/mu[phi+1]*(7.+5.*nuk)*(7.-10.*nukp1)-(7.-10*nuk)*(7.+5.*nukp1);
00615 bk=4.*(7.-10.*nuk)+mu[phi]/mu[phi+1]*(7.+5.*nuk);
00616 ck=(7.-5.*nukp1)+2.*mu[phi]/mu[phi+1]*(4.-5.*nukp1);
00617 dk=(7.+5.*nukp1)+4.*mu[phi]/mu[phi+1]*(7.-10.*nukp1);
00618 ek=2*(4.-5.*nuk)+mu[phi]/mu[phi+1]*(7.-5.*nuk);
00619 fk=(4.-5.*nuk)*(7.-5.*nukp1)-mu[phi]/mu[phi+1]*(4.-5.*nukp1)*(7.-5.*nuk);
00620 alphak=mu[phi]/mu[phi+1]-1.;
00621 koef=1./(5.*(1-nukp1));
00622
00623 M_test[0][0]=koef*ck/3.;
00624 M_test[0][1]=koef*pow(r[phi],2)*(3*bk-7*ck)/(5.*(1.-2*nuk));
00625 M_test[0][2]=koef*(-12.)*alphak/pow(r[phi],5);
00626 M_test[0][3]=koef*4.*(fk-27.*alphak)/(15.*(1.-2*nuk)*pow(r[phi],3));
00627
00628 M_test[1][0]=koef*0.;
00629 M_test[1][1]=koef*(1.-2*nukp1)*bk/(7.*(1.-2.*nuk));
00630 M_test[1][2]=koef*(-20.)*(1.-2.*nukp1)*alphak/(7*pow(r[phi],7));
00631 M_test[1][3]=koef*(-12.)*alphak*(1.-2.*nukp1)/(7.*(1.-2.*nuk)*(pow(r[phi],5)));
00632
00633 M_test[2][0]=koef*pow(r[phi],5)*alphak/2.;
00634 M_test[2][1]=koef*(-pow(r[phi],7))*(2*ak+147*alphak)/(70.*(1.-2.*nuk));
00635 M_test[2][2]=koef*dk/7.;
00636 M_test[2][3]=koef*pow(r[phi],2)*(105*(1.-nukp1)+12.*alphak*(7.-10.*nukp1)-7.*ek)/(35.*(1.-2*nuk));
00637
00638 M_test[3][0]=koef*(-5.)/6.*(1.-2*nukp1)*alphak*pow(r[phi],3);
00639 M_test[3][1]=koef*7.*(1.-2*nukp1)*alphak*pow(r[phi],5)/(2.*(1.-2.*nuk));
00640 M_test[3][2]=koef*0.;
00641 M_test[3][3]=koef*ek*(1.-2.*nukp1)/(3.*(1.-2*nuk));
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653 allocm(4,4,P_arr[phi]);
00654
00655 mxm(M,Mhelp,P_arr[phi]);
00656
00657 copym(P_arr[phi], Mhelp);
00658 }
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669 P_v[0]=P_arr[NumPhases-2].a[4+1];
00670 P_v[1]=-P_arr[NumPhases-2].a[4+0];
00671 P_v[2]=P_v[3]=0;
00672 for(phi=1;phi<=NumPhases-1;phi++){
00673 mxv(P_arr[phi-1], P_v, Phelp);
00674 cmulv(A[NumPhases-1]/(P_arr[NumPhases-2].a[4+1]*P_arr[NumPhases-2].a[0+0]-P_arr[NumPhases-2].a[0+1]*P_arr[NumPhases-2].a[4+0]),Phelp, Phelp);
00675 A[phi]=Phelp[0];
00676 B[phi]=Phelp[1];
00677 C[phi]=Phelp[2];
00678 D[phi]=Phelp[3];
00679
00680 }
00681
00682
00683 for(phi=0;phi<=NumPhases-2;phi++){
00684
00685 fillL(L,r[phi],mu,k,phi);
00686 Phelp[0]=A[phi];
00687 Phelp[1]=B[phi];
00688 Phelp[2]=C[phi];
00689 Phelp[3]=D[phi];
00690 mxv(L, Phelp, P_v);
00691
00692
00693 fillL(L,r[phi],mu,k,phi+1);
00694 Phelp[0]=A[phi+1];
00695 Phelp[1]=B[phi+1];
00696 Phelp[2]=C[phi+1];
00697 Phelp[3]=D[phi+1];
00698 mxv(L, Phelp, P_v);
00699
00700 }
00701
00702
00703
00704
00705 for(i=0;i<=3;i++){
00706 for(int j=0;j<=3;j++){
00707 Z[i][j]=P_arr[NumPhases-3].a[4*i+0]*P_arr[NumPhases-3].a[4*j+1]-
00708 P_arr[NumPhases-3].a[4*j+0]*P_arr[NumPhases-3].a[4*i+1];
00709
00710 }
00711 }
00712
00713
00714
00715
00716
00717
00718
00719 nu_n=(3.*k[NumPhases-2]-2.*mu[NumPhases-2])/(6.*k[NumPhases-2]+2.*mu[NumPhases-2]);
00720
00721 a=4.*pow(r[NumPhases-2],10)*(1.-2.*nu_n)*(7.-10.*nu_n)*Z[0][1]+
00722 20.*pow(r[NumPhases-2],7)*(7.-12.*nu_n+8.*nu_n*nu_n)*Z[3][1]+
00723 12.*pow(r[NumPhases-2],5)*(1.-2.*nu_n)*(Z[0][3]-7.*Z[1][2])+
00724 20.*pow(r[NumPhases-2],3)*(1.-2.*nu_n)*(1.-2.*nu_n)*Z[0][2]+
00725 16.*(4.-5.*nu_n)*(1.-2.*nu_n)*Z[3][2];
00726
00727 b=3.*pow(r[NumPhases-2],10)*(1.-2.*nu_n)*(15.*nu_n-7.)*Z[0][1]+
00728 60.*pow(r[NumPhases-2],7)*(nu_n-3.)*nu_n*Z[3][1]-
00729 24.*pow(r[NumPhases-2],5)*(1.-2.*nu_n)*(Z[0][3]-7.*Z[1][2])-
00730 40.*pow(r[NumPhases-2],3)*(1.-2.*nu_n)*(1.-2.*nu_n)*Z[0][2]-
00731 8.*(1.-5.*nu_n)*(1.-2.*nu_n)*Z[3][2];
00732
00733 c=-pow(r[NumPhases-2],10)*(1.-2.*nu_n)*(7.+5*nu_n)*Z[0][1]+
00734 10.*pow(r[NumPhases-2],7)*(7.-nu_n*nu_n)*Z[3][1]+
00735 12.*pow(r[NumPhases-2],5)*(1.-2.*nu_n)*(Z[0][3]-7.*Z[1][2])+
00736 20.*pow(r[NumPhases-2],3)*(1.-2.*nu_n)*(1.-2.*nu_n)*Z[0][2]-
00737 8.*(7.-5.*nu_n)*(1.-2.*nu_n)*Z[3][2];
00738
00739
00740
00741
00742 if((sqr=b*b-4.*a*c)<0){
00743 printf("shear modulus does not yield real number, terminating\n");
00744 exit(1);
00745 }
00746
00747 if((-b-pow(sqr, 0.5))>=0){
00748 printf("Two solutions for shear modulus were found, continuing with the higher value\n");
00749 }
00750
00751
00752 mu_hmg=mu[NumPhases-2]*(-b+pow(sqr, 0.5))/(2.*a);
00753
00754
00755
00756
00757
00758 KMuToENu(k_hmg, mu_hmg, E_hmg, nu_hmg);
00759
00760 destrv(r);
00761 destrv(k);
00762 destrv(mu);
00763
00764 destrv(A);
00765 destrv(B);
00766 destrv(C);
00767 destrv(D);
00768 destrv(P_v);
00769 destrv(Phelp);
00770
00771 destrm(L);
00772 destrm(Linv);
00773 destrm(M);
00774 destrm(Mhelp);
00775 destrm(Z);
00776 delete [] P_arr;
00777 }
00778
00779
00780 void HomogData::fillJ(matrix &J, double r, const vector &mu, const vector &k, int phase){
00781 J[0][0]=r;
00782 J[0][1]=1./(pow(r,2));
00783 J[1][0]=3*k[phase];
00784 J[1][1]=-4.*mu[phase]/(pow(r,3));
00785 }
00786
00787
00788 void HomogData::fillL(matrix &L, double r, const vector &mu, const vector &k, int phase){
00789 double mu_l=mu[phase];
00790 double k_l=k[phase];
00791 double nu_l=(3.*k_l-2.*mu_l)/(6.*k_l+2.*mu_l);
00792
00793 L[0][0]=r;
00794 L[0][1]=-6.*nu_l*pow(r,3)/(1.-2.*nu_l);
00795 L[0][2]=3./pow(r,4);
00796 L[0][3]=(5.-4.*nu_l)/((1.-2.*nu_l)*pow(r,2));
00797
00798 L[1][0]=r;
00799 L[1][1]=-(7.-4.*nu_l)*pow(r,3)/(1.-2.*nu_l);
00800 L[1][2]=-2./pow(r,4);
00801 L[1][3]=2./pow(r,2);
00802
00803 L[2][0]=mu_l;
00804 L[2][1]=3.*nu_l*mu_l*pow(r,2)/(1.-2.*nu_l);
00805 L[2][2]=-12.*mu_l/pow(r,5);
00806 L[2][3]=2.*(nu_l-5.)*mu_l/((1.-2.*nu_l)*pow(r,3));
00807
00808 L[3][0]=mu_l;
00809 L[3][1]=-(7.+2*nu_l)*mu_l*pow(r,2)/(1.-2.*nu_l);
00810 L[3][2]=8.*mu_l/pow(r,5);
00811 L[3][3]=2.*(1+nu_l)*mu_l/((1.-2.*nu_l)*pow(r,3));
00812 }
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824 void HomogData::KusterToksoz(void){
00825 double k_KT, mu_KT, nom, denom;
00826 nom=1+(4*mu_m*(k_i-k_m)/((3*k_i*+4*mu_m)*k_m))*f_i;
00827 denom=1-(3*(k_i-k_m)/(3*k_i+4*mu_m))*f_i;
00828 k_KT=k_m*nom/denom;
00829
00830 nom=(6*k_m+12*mu_m)*mu_i+(9*k_m+8*mu_m)*(f_m*mu_m+f_i*mu_i);
00831 denom=(9*k_m+8*mu_m)*mu_m+(6*k_m+12*mu_m)*(f_m*mu_i+f_i*mu_m);
00832 mu_KT=mu_m*nom/denom;
00833
00834 KMuToENu( k_KT, mu_KT, E_hmg, nu_hmg );
00835
00836 k_hmg=k_KT;
00837 mu_hmg=mu_KT;
00838 }
00839
00840 void HomogData::ENuToKMu(const double oE, const double onu, double &ok, double &omu )
00841 {
00842 omu = oE / ( 2. * ( 1. + onu ) );
00843 ok = oE / ( 3. * ( 1. - 2. * onu ) );
00844 }
00845
00846 void HomogData::KMuToENu(const double ok, const double omu, double &oE, double &onu )
00847 {
00848 oE = 9. * ok * omu / ( 3. * ok + omu );
00849 onu = ( 3. * ok - 2. * omu ) / ( 6. * ok + 2 * omu );
00850 }
00851
00852 void HomogData::Swap(double &A, double &B){
00853 double C;
00854 C=A;
00855 A=B;
00856 B=C;
00857 }