OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
homogenize.C
Go to the documentation of this file.
1 /*
2  *
3  * ##### ##### ###### ###### ### ###
4  * ## ## ## ## ## ## ## ### ##
5  * ## ## ## ## #### #### ## # ##
6  * ## ## ## ## ## ## ## ##
7  * ## ## ## ## ## ## ## ##
8  * ##### ##### ## ###### ## ##
9  *
10  *
11  * OOFEM : Object Oriented Finite Element Code
12  *
13  * Copyright (C) 1993 - 2013 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 
35 #include "homogenize.h"
36 #include "mathfem.h"
37 
38 namespace oofem {
40 { }
41 
42 //parallel configuration of Voigt
43 void Homogenize :: voigt(FloatMatrix &PhaseMatrix)
44 {
45  double k, mu;
46  double volTot = 0.;
47  int r;
48  int NumPhases = PhaseMatrix.giveNumberOfRows();
49 
50  checkVolFraction(PhaseMatrix);
51 
52  k_hmg = 0.;
53  mu_hmg = 0.;
54  for ( r = 0; r < NumPhases; r++ ) {
55  ENuToKMu(PhaseMatrix(r, 1), PhaseMatrix(r, 2), k, mu);
56  k_hmg += k * PhaseMatrix(r, 0);
57  mu_hmg += mu * PhaseMatrix(r, 0);
58  volTot += PhaseMatrix(r, 0);
59  }
60 
62 }
63 
64 
65 //serial configuration of Reuss
66 void Homogenize :: reuss(FloatMatrix &PhaseMatrix)
67 {
68  double k, mu;
69  int r;
70  int NumPhases = PhaseMatrix.giveNumberOfRows();
71 
72  checkVolFraction(PhaseMatrix);
73 
74  k_hmg = 0.;
75  mu_hmg = 0.;
76  for ( r = 0; r < NumPhases; r++ ) {
77  ENuToKMu(PhaseMatrix(r, 1), PhaseMatrix(r, 2), k, mu);
78  k_hmg += PhaseMatrix(r, 0) / k;
79  mu_hmg += PhaseMatrix(r, 0) / mu;
80  }
81 
82  k_hmg = 1. / k_hmg;
83  mu_hmg = 1. / mu_hmg;
85 }
86 
87 /*
88  * Hashin-Shtrikman-Walpole lower and upper bounds, arbitrary number of phases
89  * If phases satisfy ascending order of bulk or shear moduli simultaneously, i.e. k_0<k_1<...k_N together
90  * with mu_0<mu_1<...mu_N, Hashin-Shtrikman bounds are obtained
91  * If the order of bulk and shear moduli is not satisfied, wider Hashin-Shtrikman-Walpole bounds are computed
92  * Implementation according to J.G. Berryman: Mixture Theories for Rock Properties, Physics and Phase Relations,
93  * A Handbook of Physical Constants, Am. Geophys, 1995 or download directly from the web
94  */
96 {
97  int phase = 0, i, r;
98  int counter = 0;
99  int numRows = PhaseMatrix.giveNumberOfRows();
100  double k, k_min, mu, mu_phase = 0., muMin, muMax;
101  double dummy;
102  FloatMatrix PhaseMatrix1;
103  PhaseMatrix1 = PhaseMatrix;
104  FloatMatrix SortedPhaseMatrix(numRows, 3);
105 
106 
107  checkVolFraction(PhaseMatrix1);
108 
109  //Sort phases in ascending order of bulk moduli
110  for ( i = 0; i < numRows; i++ ) {
111  k_min = 1.e+40;
112  for ( r = 0; r < numRows; r++ ) {
113  if ( PhaseMatrix1(r, 0) > 0. + 1.e-40 ) {
114  ENuToKMu(PhaseMatrix1(r, 1), PhaseMatrix1(r, 2), k, mu);
115  if ( k < k_min ) {
116  phase = r;
117  k_min = k;
118  mu_phase = mu;
119  }
120  }
121  }
122 
123  SortedPhaseMatrix(counter, 0) = PhaseMatrix1(phase, 0);
124  SortedPhaseMatrix(counter, 1) = k_min;
125  SortedPhaseMatrix(counter, 2) = mu_phase;
126 
127  PhaseMatrix1(phase, 0) = 0.;
128  counter++;
129  }
130 
131  //SortedPhaseMatrix contains (vol. fraction, bulk modulus, shear modulus), sorted in ascending order of bulk modulus
132  //Find phase numbers for maximum and minimum shear moduli
133  muMin = 1.e+40;
134  muMax = 0.;
135  for ( i = 0; i < numRows; i++ ) {
136  if ( SortedPhaseMatrix(i, 2) > muMax ) {
137  muMax = SortedPhaseMatrix(i, 2);
138  }
139 
140  if ( SortedPhaseMatrix(i, 2) < muMin ) {
141  muMin = SortedPhaseMatrix(i, 2);
142  }
143  }
144 
145  // find bounds for bulk moduli
146  k_hmg = 0; //lower bound
147  k_hmg_2 = 0; //upper bound
148  for ( i = 0; i < numRows; i++ ) {
149  k_hmg += SortedPhaseMatrix(i, 0) / ( SortedPhaseMatrix(i, 1) + 4. / 3. * muMin );
150  k_hmg_2 += SortedPhaseMatrix(i, 0) / ( SortedPhaseMatrix(i, 1) + 4. / 3. * muMax );
151  }
152 
153  k_hmg = 1. / k_hmg;
154  k_hmg -= 4. / 3. * muMin;
155  k_hmg_2 = 1. / k_hmg_2;
156  k_hmg_2 -= 4. / 3. * muMax;
157 
158  //find bounds for shear moduli
159  mu_hmg = gamma( SortedPhaseMatrix, zeta(SortedPhaseMatrix(0, 1), muMin) ); //lower bound
160  mu_hmg_2 = gamma( SortedPhaseMatrix, zeta(SortedPhaseMatrix(numRows - 1, 1), muMax) ); //upper bound
161 
162  //Young's moduli are determined from k_min, mu_min and k_max an mu_max
163  kMuToENu(k_hmg, mu_hmg, E_hmg, dummy);
164  kMuToENu(k_hmg_2, mu_hmg_2, E_hmg_2, dummy);
165  //Poisson's ratii are determined from k_min, mu_max and k_max, mu_min
166  kMuToENu(k_hmg, mu_hmg_2, dummy, nu_hmg);
167  kMuToENu(k_hmg_2, mu_hmg, dummy, nu_hmg_2);
168 }
169 
170 /*
171  * Mori-Tanaka homogenization method
172  * isotropic strain concentration tensor = 3*k_S[r]*J + 2*mu_S[r]*I
173  * refRow is the reference phase of this scheme (0 is the first phase)
174  * for derivation see Bernard et al.: A multiscale micromechanics-hydration model... CCR,pp. 1293-1309, 2003
175  */
176 void Homogenize :: moriTanaka(FloatMatrix &PhaseMatrix, int refRow)
177 {
178  double f_r, E_r, nu_r, k_r, mu_r;
179  double E_m, nu_m, k_m, mu_m;
180  double fr_tot = 0.;
181  double nom_k_MT = 0., denom_k_MT = 0., nom_mu_MT = 0., denom_mu_MT = 0.;
182  double k_S_denom = 0., mu_S_denom = 0.;
183  double alpha_m, beta_m;
184  int numRows = PhaseMatrix.giveNumberOfRows();
185  int r;
186  FloatArray k_S(numRows), mu_S(numRows);
187 
188  checkVolFraction(PhaseMatrix);
189 
190  //determine matrix parameters
191  E_m = PhaseMatrix(refRow, 1);
192  nu_m = PhaseMatrix(refRow, 2);
193 
194  ENuToKMu(E_m, nu_m, k_m, mu_m);
195 
196  //auxiliary parameters
197  alpha_m = 3. * k_m / ( 3. * k_m + 4 * mu_m );
198  beta_m = 6. * ( k_m + 2. * mu_m ) / ( 5. * ( 3. * k_m + 4. * mu_m ) );
199 
200  //cycle through all phases, including the matrix phase
201  for ( r = 0; r < numRows; r++ ) {
202  f_r = PhaseMatrix(r, 0);
203  E_r = PhaseMatrix(r, 1);
204  nu_r = PhaseMatrix(r, 2);
205 
206  fr_tot += f_r;
207  ENuToKMu(E_r, nu_r, k_r, mu_r);
208 
209  nom_k_MT += f_r * k_r / ( 1. + alpha_m * ( k_r / k_m - 1. ) );
210  denom_k_MT += f_r / ( 1. + alpha_m * ( k_r / k_m - 1. ) );
211 
212  nom_mu_MT += f_r * mu_r / ( 1. + beta_m * ( mu_r / mu_m - 1. ) );
213  denom_mu_MT += f_r / ( 1. + beta_m * ( mu_r / mu_m - 1. ) );
214 
215  k_S [ r ] = 1. / ( 1. + alpha_m * ( k_r / k_m - 1. ) );
216  mu_S [ r ] = 1. / ( 1. + beta_m * ( mu_r / mu_m - 1. ) );
217  }
218 
219 
220  k_hmg = nom_k_MT / denom_k_MT;
221  mu_hmg = nom_mu_MT / denom_mu_MT;
222 
224 
225  //calculate denominator of strain concentration tensor (in terms of k_S and mu_S) and final values
226  for ( r = 0; r < numRows; r++ ) {
227  f_r = PhaseMatrix(r, 0);
228  E_r = PhaseMatrix(r, 1);
229  nu_r = PhaseMatrix(r, 2);
230  ENuToKMu(E_r, nu_r, k_r, mu_r);
231  k_S_denom += f_r * 1. / ( 1. + alpha_m * ( k_r / k_m - 1. ) );
232  mu_S_denom += f_r * 1. / ( 1. + beta_m * ( mu_r / mu_m - 1. ) );
233  }
234 
235  for ( r = 0; r < numRows; r++ ) {
236  E_r = PhaseMatrix(r, 1);
237  nu_r = PhaseMatrix(r, 2);
238  ENuToKMu(E_r, nu_r, k_r, mu_r);
239  // These are pointless operations, and will not do anything as far as I can see / Mikael
240  k_S [ r ] /= k_S_denom;
241  mu_S [ r ] /= mu_S_denom;
242  }
243 }
244 
245 // for derivation see Bernard et al.:A multiscale micromechanics-hydration model... CCR,pp. 1293-1309, 2003
247 {
248  double f_r, E_r, nu_r, k_r, mu_r;
249  double k_SCS, mu_SCS, nom_k_SCS, denom_k_SCS, nom_mu_SCS, denom_mu_SCS;
250  double fr_tot;
251  double alpha_m = 0., beta_m = 0.;
252  double k_S_denom = 0., mu_S_denom = 0.;
253  int numRows = PhaseMatrix.giveNumberOfRows();
254  FloatArray k_S(numRows), mu_S(numRows);
255 
256  checkVolFraction(PhaseMatrix);
257 
258  /*first estimation*/
259  k_SCS = 10.;
260  mu_SCS = .3;
261 
262  /*iteration for nonlinear equations*/
263  for ( int i = 1; i < 100; i++ ) {
264  fr_tot = 0.;
265  nom_k_SCS = 0;
266  denom_k_SCS = 0;
267  nom_mu_SCS = 0;
268  denom_mu_SCS = 0;
269  for ( int r = 0; r < numRows; r++ ) {
270  f_r = PhaseMatrix(r, 0);
271  E_r = PhaseMatrix(r, 1);
272  nu_r = PhaseMatrix(r, 2);
273 
274  fr_tot += f_r;
275 
276  ENuToKMu(E_r, nu_r, k_r, mu_r);
277 
278  alpha_m = 3. * k_SCS / ( 3. * k_SCS + 4 * mu_SCS );
279  beta_m = 6. * ( k_SCS + 2. * mu_SCS ) / ( 5. * ( 3. * k_SCS + 4. * mu_SCS ) );
280 
281  nom_k_SCS += f_r * k_r / ( 1 + alpha_m * ( k_r / k_SCS - 1 ) );
282  denom_k_SCS += f_r / ( 1 + alpha_m * ( k_r / k_SCS - 1 ) );
283 
284  nom_mu_SCS += f_r * mu_r / ( 1 + beta_m * ( mu_r / mu_SCS - 1 ) );
285  denom_mu_SCS += f_r / ( 1 + beta_m * ( mu_r / mu_SCS - 1 ) );
286  }
287 
288  k_SCS = nom_k_SCS / denom_k_SCS;
289  mu_SCS = nom_mu_SCS / denom_mu_SCS;
290  }
291 
292  k_hmg = k_SCS;
293  mu_hmg = mu_SCS;
294 
296 
297  //calculate nominator and denominator of strain concentration tensor (in terms of k_S and mu_S) and final values
298  for ( int r = 0; r < numRows; r++ ) {
299  f_r = PhaseMatrix(r, 0);
300  E_r = PhaseMatrix(r, 1);
301  nu_r = PhaseMatrix(r, 2);
302  ENuToKMu(E_r, nu_r, k_r, mu_r);
303  k_S [ r ] = 1. / ( 1. + alpha_m * ( k_r / k_hmg - 1. ) );
304  mu_S [ r ] = 1. / ( 1. + beta_m * ( mu_r / mu_hmg - 1. ) );
305  k_S_denom += f_r * 1. / ( 1. + alpha_m * ( k_r / k_hmg - 1. ) );
306  mu_S_denom += f_r * 1. / ( 1. + beta_m * ( mu_r / mu_hmg - 1. ) );
307  }
308 
309  for ( int r = 0; r < numRows; r++ ) {
310  E_r = PhaseMatrix(r, 1);
311  nu_r = PhaseMatrix(r, 2);
312  ENuToKMu(E_r, nu_r, k_r, mu_r);
313  k_S [ r ] /= k_S_denom;
314  mu_S [ r ] /= mu_S_denom;
315  }
316 }
317 
318 /*
319  * Homogenization scheme of Herve and Zaoui for n-spherical isotropic domains
320  * see Herve E., Zaoui A. n-layered Inclusion-based Micromechanical Modelling,
321  * Int. J. Engng Sci. 31, pp 1-10, 1993.
322  *
323  * The sum of n phases must yield a unity, (n+1) may have arbitrary radius
324  *
325  * The scheme is derived from stress and displacement equivalence at boundaries among phases
326  * homogenized bulk and shear response is calculated
327  * NumPhases is the number of phases including homogeneous medium (n+1)
328  * In homogenized values, (n+1) phase may have arbitrary values - if it has the right
329  * homogenized values of the assembly, the equations (49) and (51) are satisfied exactly without solution
330  *
331  * In case of n=2, the Hashin-Shtrikman bounds are calculated if the well-ordered moduli are used
332  * Homogenized bulk modulus is 100 % OK (corresponds to HS bounds for 2 phases)
333  *
334  * Shear modulus is with small difference when compared to HS bound for 2 phases
335  * typical usage for mortars and concrete is
336  * phase 0 - aggregate
337  * phase 1 - ITZ
338  * phase 2 (n) - cement paste
339  * phase 3 (n+1) - equivalent concrete medium
340  *
341  * The bulk modulus of this scheme coincides with HS bound. This is not true for shear modulus, see Christensen: Two Theoretical Elasticity Micromechanics Models, Journal of Elasticity 50: 15–25, 1998.
342  */
343 
345 {
346  int i, j, phi;
347  int NumPhases = PhaseMatrix.giveNumberOfRows() + 1; //need to extend for an equivalent homogeneous medium
348  FloatMatrix PhaseMatrix1( NumPhases, PhaseMatrix.giveNumberOfColumns() );
349  checkVolFraction(PhaseMatrix);
350  PhaseMatrix1.setSubMatrix(PhaseMatrix, 1, 1);
351  //add arbitrary values for a reference medium
352  PhaseMatrix1(NumPhases - 1, 0) = 0.1;
353  PhaseMatrix1(NumPhases - 1, 1) = 10.;
354  PhaseMatrix1(NumPhases - 1, 2) = 0.3;
355 
356  // a)HYDROSTATIC pressure - conditions of displacement and radial stress compatibility between adjacent phases
357  FloatArray r(NumPhases - 1), k(NumPhases), mu(NumPhases);
358  FloatArray Q11(NumPhases - 1), Q21(NumPhases - 1), F(NumPhases), G(NumPhases);
359  FloatMatrix J(2, 2), Jinv(2, 2), N(2, 2), Nhelp(2, 2), Q(2, 2);
360  //set lambda_0 for strain control - arbitrary value for homogenization
361  double lambda_0 = 0.735;
362  double temp1;
363 
364  if ( NumPhases < 3 ) {
365  OOFEM_ERROR("Number of considered phases must be at least 3 (including hom. medium)");
366  }
367 
368  F(NumPhases - 1) = lambda_0 / 3.;
369 
370  //radii of spheres, assume 1 as the maximal radius
371  temp1 = 0.;
372  for ( i = 0; i < NumPhases - 1; i++ ) {
373  temp1 += PhaseMatrix1(i, 0);
374  r(i) = 1. * cbrt(temp1);
375  //printf( "r(%d)=%f\n", i, r(i) );
376  }
377 
378  for ( i = 0; i < NumPhases; i++ ) {
379  ENuToKMu( PhaseMatrix1(i, 1), PhaseMatrix1(i, 2), k(i), mu(i) );
380  //printf("k(%d)=%f mu(%d)=%f\n", i, k(i), i, mu(i));
381  }
382 
383  //create unit matrix for the first multiplication
384  Nhelp(0, 0) = 1.;
385  Nhelp(1, 1) = 1.;
386 
387  //calculate Q
388  for ( phi = 0; phi <= NumPhases - 2; phi++ ) {
389  //calculate inverse of J_{phi+1}
390  fillJ(J, r(phi), mu, k, phi + 1);
391  Jinv.beInverseOf(J);
392  //calculate J_{phi}
393  fillJ(J, r(phi), mu, k, phi);
394  //calculate N^{phi}=J_{phi+1}^{-1}*J_{phi}
395  N.beProductOf(Jinv, J);
396  //calculate a part of Q^{phi}
397  Q.beProductOf(N, Nhelp);
398  //store part of matrix Q in vector Q11 and Q21
399  Q11(phi) = Q(0, 0);
400  Q21(phi) = Q(1, 0);
401  //copy Q to Nhelp and use in the next cycle
402  Nhelp = Q;
403  }
404 
405  //calculate V (contains vector components F,G), F(0) and G(0) make no sense
406  for ( phi = 1; phi <= NumPhases - 1; phi++ ) {
407  F(phi) = F(NumPhases - 1) * Q11(phi - 1) / Q11(NumPhases - 2);
408  G(phi) = F(NumPhases - 1) * Q21(phi - 1) / Q11(NumPhases - 2);
409  //printf("F(%d)=%f G(%d)=%f\n", phi, F(phi), phi, G(phi));
410  }
411 
412  //check displacement and stress continuity, phi is the more distant phase
413  /*
414  * for ( int phi = 1; phi <= NumPhases - 1; phi++ ) {
415  * printf("r=%f displ(in)=%f displ(out)=%f\n", r(phi-1), F(phi-1)*r(phi-1)+G(phi-1)/(r(phi-1)*r(phi-1)), F(phi)*r(phi-1)+G(phi)/(r(phi-1)*r(phi-1)));
416  * printf("r=%f stress(in)=%f stress(out)=%f\n", r(phi-1), 3.*k(phi-1)*F(phi-1)-4.*mu(phi-1)*G(phi-1)/(r(phi-1)*r(phi-1)*r(phi-1)), 3.*k(phi)*F(phi)-4.*mu(phi)*G(phi)/(r(phi-1)*r(phi-1)*r(phi-1)));
417  * }
418  */
419 
420  /*replace n+1 layer with equivalent homogeneous media calculate homogenized k
421  * checked with Mori-Tanaka method for two phases, bulk modulus is exact
422  */
423  k_hmg = ( 3. * k(NumPhases - 2) * pow(r(NumPhases - 2), 3.) * Q11(NumPhases - 3) - 4. * mu(NumPhases - 2) * Q21(NumPhases - 3) ) / ( 3. * ( pow(r(NumPhases - 2), 3.) * Q11(NumPhases - 3) + Q21(NumPhases - 3) ) );
424 
425  //printf("Homogenized k=%f\n", k_hmg);
426 
427  //b)SIMPLE SHEAR - continuity of radial and phi displacements and two shear stresses
428  FloatArray B(NumPhases), C(NumPhases), D(NumPhases), A(NumPhases);
429  FloatArray P_v(4), Phelp(4);
430  FloatMatrix L(4, 4), Linv(4, 4), M(4, 4), Mhelp(4, 4), Z(4, 4), M_test(4, 4);
431  std :: vector< FloatMatrix > P_arr(NumPhases); //array to store individual P matrices, allocate memory for P matrices, store there P11, P12, P21, P22
432  //set gamma for displacement approach - arbitrary value for homogenization
433  //variables for homogenized values
434  double a, b, c, nu_n, sqr = 0.;
435  double gamma = 10.; //no effect for homogenziation
436 
437  A.zero();
438  A(NumPhases - 1) = gamma;
439 
440  //create unit matrix for the first multiplication
441  Mhelp.zero();
442  Mhelp.beUnitMatrix();
443 
444  //calculate P
445  for ( phi = 0; phi < NumPhases - 1; phi++ ) {
446  //calculate inverse of L_{phi+1}
447  fillL(L, r(phi), mu, k, phi + 1);
448  //printf("%f", L(0,1));
449  Linv.beInverseOf(L);
450  //calculate L_{phi}
451  fillL(L, r(phi), mu, k, phi);
452  //calulate M^{phi}=L_{phi+1}^{-1}*L_{phi}
453  M.beProductOf(Linv, L);
454 
455  /*test routine directly for M
456  * this part calculate inverse of M directly
457  * M=M_test so the L, Linv and M assembly are OK
458  *
459  * double ak, bk, ck, dk, ek, fk, alphak, nuk, nukp1, koef;
460  *
461  * nuk=PhaseMatrix1(phi,2);
462  * nukp1=PhaseMatrix1(phi+1,2);
463  *
464  * ak=mu(phi)/mu(phi+1)*(7.+5.*nuk)*(7.-10.*nukp1)-(7.-10*nuk)*(7.+5.*nukp1);
465  * bk=4.*(7.-10.*nuk)+mu(phi)/mu(phi+1)*(7.+5.*nuk);
466  * ck=(7.-5.*nukp1)+2.*mu(phi)/mu(phi+1)*(4.-5.*nukp1);
467  * dk=(7.+5.*nukp1)+4.*mu(phi)/mu(phi+1)*(7.-10.*nukp1);
468  * ek=2*(4.-5.*nuk)+mu(phi)/mu(phi+1)*(7.-5.*nuk);
469  * fk=(4.-5.*nuk)*(7.-5.*nukp1)-mu(phi)/mu(phi+1)*(4.-5.*nukp1)*(7.-5.*nuk);
470  * alphak=mu(phi)/mu(phi+1)-1.;
471  * koef=1./(5.*(1-nukp1));
472  *
473  * M_test(0,0)=koef*ck/3.;
474  * M_test(0,1)=koef*pow(r(phi),2)*(3*bk-7*ck)/(5.*(1.-2*nuk));
475  * M_test(0,2)=koef*(-12.)*alphak/pow(r(phi),5);
476  * M_test(0,3)=koef*4.*(fk-27.*alphak)/(15.*(1.-2*nuk)*pow(r(phi),3));
477  *
478  * M_test(1,0)=koef*0.;
479  * M_test(1,1)=koef*(1.-2*nukp1)*bk/(7.*(1.-2.*nuk));
480  * M_test(1,2)=koef*(-20.)*(1.-2.*nukp1)*alphak/(7*pow(r(phi),7));
481  * M_test(1,3)=koef*(-12.)*alphak*(1.-2.*nukp1)/(7.*(1.-2.*nuk)*(pow(r(phi),5)));
482  *
483  * M_test(2,0)=koef*pow(r(phi),5)*alphak/2.;
484  * M_test(2,1)=koef*(-pow(r(phi),7))*(2*ak+147*alphak)/(70.*(1.-2.*nuk));
485  * M_test(2,2)=koef*dk/7.;
486  * M_test(2,3)=koef*pow(r(phi),2)*(105*(1.-nukp1)+12.*alphak*(7.-10.*nukp1)-7.*ek)/(35.*(1.-2*nuk));
487  *
488  * M_test(3,0)=koef*(-5.)/6.*(1.-2*nukp1)*alphak*pow(r(phi),3);
489  * M_test(3,1)=koef*7.*(1.-2*nukp1)*alphak*pow(r(phi),5)/(2.*(1.-2.*nuk));
490  * M_test(3,2)=koef*0.;
491  * M_test(3,3)=koef*ek*(1.-2.*nukp1)/(3.*(1.-2*nuk));
492  *
493  * for(int ii=0;ii<4;ii++)
494  * for(int jj=0;jj<4;jj++)
495  * M(ii,jj)=M_test(ii,jj);
496  */
497  P_arr [ phi ].resize(4, 4);
498  P_arr [ phi ].beProductOf(M, Mhelp);
499  Mhelp = P_arr [ phi ];
500  }
501 
502  //printf("PPP %f\n", P_arr[1](3,0)*P_arr[1](1,1) - P_arr[1](3,1)*P_arr[1](1,0));
503 
504  //calculate vector W (contains scalar components A,B,C,D for each step)
505 
506  P_v(0) = P_arr [ NumPhases - 2 ](1, 1); //matrix stored column wise (P11, P21, P31....), P_22
507  P_v(1) = -P_arr [ NumPhases - 2 ](1, 0); // P_21
508  P_v(2) = P_v(3) = 0.;
509 
510  for ( phi = 1; phi <= NumPhases - 1; phi++ ) {
511  Phelp.beProductOf(P_arr [ phi - 1 ], P_v);
512  Phelp.times( A(NumPhases - 1) / sqrt(2.) / ( P_arr [ NumPhases - 2 ](0, 0) * P_arr [ NumPhases - 2 ](1, 1) - P_arr [ NumPhases - 2 ](0, 1) * P_arr [ NumPhases - 2 ](1, 0) ) );
513  A(phi) = Phelp(0);
514  B(phi) = Phelp(1);
515  C(phi) = Phelp(2);
516  D(phi) = Phelp(3);
517  //printf("A(%d]=%f B(%d)=%f C(%d)=%f D(%d)=%f\n", phi, A(phi), phi, B(phi), phi, C(phi), phi, D(phi));
518  }
519 
520  //check displacement and stress continuity, phi is the less distant phase
521  //this cycle for for checking purposes only
522  for ( phi = 0; phi <= NumPhases - 2; phi++ ) {
523  //phase closer to the center
524  fillL(L, r(phi), mu, k, phi);
525  Phelp(0) = A(phi);
526  Phelp(1) = B(phi);
527  Phelp(2) = C(phi);
528  Phelp(3) = D(phi);
529  P_v.beProductOf(L, Phelp);
530  //printf("(in) n=%d r=%f ur=%f uphi=%f srp=%f srf=%f\n",phi, r(phi), P_v(0),P_v(1), P_v(2), P_v(3));
531 
532  fillL(L, r(phi), mu, k, phi + 1);
533  Phelp(0) = A(phi + 1);
534  Phelp(1) = B(phi + 1);
535  Phelp(2) = C(phi + 1);
536  Phelp(3) = D(phi + 1);
537  P_v.beProductOf(L, Phelp);
538  //printf("(out) n=%d r=%f ur=%f uphi=%f srp=%f srf=%f\n", phi,r(phi), P_v(0),P_v(1), P_v(2), P_v(3));
539  }
540 
541  // replace n+1 layer with equivalent homogeneous media, calculate homogenized mu
542  for ( i = 0; i <= 3; i++ ) { //alpha, rows in Z
543  for ( j = 0; j <= 3; j++ ) { //beta, columns in Z
544  Z(i, j) = P_arr [ NumPhases - 3 ](i, 0) * P_arr [ NumPhases - 3 ](j, 1) - P_arr [ NumPhases - 3 ](j, 0) * P_arr [ NumPhases - 3 ](i, 1);
545  // printf("Z=%f ", Z(i,j));
546  }
547  }
548 
549  nu_n = ( 3. * k(NumPhases - 2) - 2. * mu(NumPhases - 2) ) / ( 6. * k(NumPhases - 2) + 2. * mu(NumPhases - 2) );
550 
551  a = 4. * pow(r(NumPhases - 2), 10.) * ( 1. - 2. * nu_n ) * ( 7. - 10. * nu_n ) * Z(0, 1) +
552  20. * pow(r(NumPhases - 2), 7.) * ( 7. - 12. * nu_n + 8. * nu_n * nu_n ) * Z(3, 1) +
553  12. * pow(r(NumPhases - 2), 5.) * ( 1. - 2. * nu_n ) * ( Z(0, 3) - 7. * Z(1, 2) ) +
554  20. * pow(r(NumPhases - 2), 3.) * ( 1. - 2. * nu_n ) * ( 1. - 2. * nu_n ) * Z(0, 2) +
555  16. * ( 4. - 5. * nu_n ) * ( 1. - 2. * nu_n ) * Z(3, 2);
556 
557  b = 3. * pow(r(NumPhases - 2), 10.) * ( 1. - 2. * nu_n ) * ( 15. * nu_n - 7. ) * Z(0, 1) +
558  60. * pow(r(NumPhases - 2), 7.) * ( nu_n - 3. ) * nu_n * Z(3, 1) -
559  24. * pow(r(NumPhases - 2), 5.) * ( 1. - 2. * nu_n ) * ( Z(0, 3) - 7. * Z(1, 2) ) -
560  40. * pow(r(NumPhases - 2), 3.) * ( 1. - 2. * nu_n ) * ( 1. - 2. * nu_n ) * Z(0, 2) -
561  8. * ( 1. - 5. * nu_n ) * ( 1. - 2. * nu_n ) * Z(3, 2);
562 
563  c = -pow(r(NumPhases - 2), 10.) * ( 1. - 2. * nu_n ) * ( 7. + 5 * nu_n ) * Z(0, 1) +
564  10. * pow(r(NumPhases - 2), 7.) * ( 7. - nu_n * nu_n ) * Z(3, 1) +
565  12. * pow(r(NumPhases - 2), 5.) * ( 1. - 2. * nu_n ) * ( Z(0, 3) - 7. * Z(1, 2) ) +
566  20. * pow(r(NumPhases - 2), 3.) * ( 1. - 2. * nu_n ) * ( 1. - 2. * nu_n ) * Z(0, 2) -
567  8. * ( 7. - 5. * nu_n ) * ( 1. - 2. * nu_n ) * Z(3, 2);
568 
569  //solve quadratic equation, report higher number
570  sqr = b * b - 4. * a * c;
571  if ( sqr < 0 ) {
572  OOFEM_ERROR("Shear modulus does not yield a real number");
573  }
574 
575  if ( ( -b - pow(sqr, 0.5) ) >= 0 ) {
576  OOFEM_WARNING("Two solutions for the shear modulus were found, continuing with the higher value");
577  }
578 
579  //usually this higher value is reported
580  mu_hmg = mu(NumPhases - 2) * ( -b + pow(sqr, 0.5) ) / ( 2. * a );
581  //uncomment if lower value of mu should be used, only when positive
582  //mu_hmg=mu(NumPhases-2)*(-b+pow(sqr, 0.5))/(2.*a);
583 
584  //printf("Homogenized mu=%f\n", mu_hmg);
585 
587 }
588 
589 /*Hirsh model as a combination of Voigt and Reuss (parameter chi)*/
590 void Homogenize :: hirsch(FloatMatrix &PhaseMatrix, double chi)
591 {
592  double k_Voigt, mu_Voigt, k_Reuss, mu_Reuss;
593  voigt(PhaseMatrix);
594  k_Voigt = k_hmg;
595  mu_Voigt = mu_hmg;
596  reuss(PhaseMatrix);
597  k_Reuss = k_hmg;
598  mu_Reuss = mu_hmg;
599  k_hmg = ( 1. - chi ) * k_Reuss + chi * k_Voigt;
600  mu_hmg = ( 1. - chi ) * mu_Reuss + chi * mu_Voigt;
602 }
603 
604 
606 {
607  double f_i, E_i, f_m, E_m;
608 
609  if ( PhaseMatrix.giveNumberOfRows() != 2 ) {
610  OOFEM_ERROR("Only two phases are allowed");
611  }
612 
613  checkVolFraction(PhaseMatrix);
614 
615  f_i = PhaseMatrix(0, 0);
616  E_i = PhaseMatrix(0, 1);
617  f_m = PhaseMatrix(1, 0);
618  E_m = PhaseMatrix(1, 1);
619 
620  E_hmg = E_i * ( ( f_i * E_i + ( 1 + f_m ) * E_m ) / ( ( 1. + f_m ) * E_i + f_i * E_m ) );
621  nu_hmg = 0.2;
622  ENuToKMu(E_hmg, nu_hmg, k_hmg, mu_hmg);
623 }
624 
625 
627 {
628  double E_i, f_m, E_m, nu_m;
629  if ( PhaseMatrix.giveNumberOfRows() != 2 ) {
630  OOFEM_ERROR("Only two phases are allowed");
631  }
632 
633  checkVolFraction(PhaseMatrix);
634 
635  E_i = PhaseMatrix(0, 1);
636  f_m = PhaseMatrix(1, 0);
637  E_m = PhaseMatrix(1, 1);
638  nu_m = PhaseMatrix(1, 2);
639 
640  E_hmg = 1. / ( ( 1. - sqrt(f_m) ) / E_i + 1. / ( E_i * ( ( 1. - sqrt(f_m) ) / sqrt(f_m) ) + E_m ) );
641  nu_hmg = nu_m;
643 }
644 
645 /* Kuster-Toksoz model from Monteiro: Elastic Moduli Lightweight Aggregate Model,
646  * Cement and Concrete Research, 1995, No.2, pp. 276-280.
647  */
649 {
650  double k_KT, mu_KT, nom, denom;
651  double f_i, E_i, nu_i, k_i, mu_i, f_m, E_m, nu_m, mu_m, k_m;
652  if ( PhaseMatrix.giveNumberOfRows() != 2 ) {
653  OOFEM_ERROR("Only two phases are allowed");
654  }
655 
656  checkVolFraction(PhaseMatrix);
657 
658  f_i = PhaseMatrix(0, 0);
659  E_i = PhaseMatrix(0, 1);
660  nu_i = PhaseMatrix(0, 2);
661  ENuToKMu(E_i, nu_i, k_i, mu_i);
662  f_m = PhaseMatrix(1, 0);
663  E_m = PhaseMatrix(1, 1);
664  nu_m = PhaseMatrix(1, 2);
665  ENuToKMu(E_m, nu_m, k_m, mu_m);
666 
667  nom = 1 + ( 4 * mu_m * ( k_i - k_m ) / ( ( 3 * k_i * +4 * mu_m ) * k_m ) ) * f_i;
668  denom = 1 - ( 3 * ( k_i - k_m ) / ( 3 * k_i + 4 * mu_m ) ) * f_i;
669  k_KT = k_m * nom / denom;
670 
671  nom = ( 6 * k_m + 12 * mu_m ) * mu_i + ( 9 * k_m + 8 * mu_m ) * ( f_m * mu_m + f_i * mu_i );
672  denom = ( 9 * k_m + 8 * mu_m ) * mu_m + ( 6 * k_m + 12 * mu_m ) * ( f_m * mu_i + f_i * mu_m );
673  mu_KT = mu_m * nom / denom;
674 
675  k_hmg = k_KT;
676  mu_hmg = mu_KT;
677 
679 }
680 
682 {
683  if ( out == 0 ) {
684  printf("E %.4e\t nu %.4e\t k %.4e\t mu %.4e\n", this->E_hmg, this->nu_hmg, this->k_hmg, this->mu_hmg);
685  } else {
686  printf("E %.4e\t nu %.4e\t k %.4e\t mu %.4e\n", this->E_hmg_2, this->nu_hmg_2, this->k_hmg_2, this->mu_hmg_2);
687  }
688 }
689 
690 
691 void Homogenize :: ENuToKMu(const double E, const double nu, double &k, double &mu)
692 {
693  mu = E / ( 2. * ( 1. + nu ) );
694  k = E / ( 3. * ( 1. - 2. * nu ) );
695 }
696 
697 void Homogenize :: kMuToENu(const double k, const double mu, double &E, double &nu)
698 {
699  E = 9. * k * mu / ( 3. * k + mu );
700  nu = ( 3. * k - 2. * mu ) / ( 6. * k + 2 * mu );
701 }
702 
703 //Auxiliary function for Hashin-Shtrikman-Walpole shear bounds
704 double Homogenize :: zeta(double k, double mu) {
705  return mu / 6. * ( 9. * k + 8. * mu ) / ( k + 2. * mu );
706 }
707 
708 //Auxiliary function for Hashin-Shtrikman-Walpole shear bounds
709 double Homogenize :: gamma(FloatMatrix &SortedPhaseMatrix, double zeta)
710 {
711  double gamma = 0;
712  int r;
713  int NumPhases = SortedPhaseMatrix.giveNumberOfRows();
714 
715  for ( r = 0; r < NumPhases; r++ ) {
716  gamma += SortedPhaseMatrix(r, 0) / ( SortedPhaseMatrix(r, 2) + zeta );
717  }
718 
719  gamma = 1. / gamma;
720  gamma -= zeta;
721  return gamma;
722 }
723 
725 {
726  double volTot = 0.;
727  int NumPhases = PhaseMatrix.giveNumberOfRows();
728  int r;
729  for ( r = 0; r < NumPhases; r++ ) {
730  volTot += PhaseMatrix(r, 0);
731  }
732 
733  if ( volTot < 0.99 || volTot > 1.01 ) {
734  OOFEM_ERROR("Volumetric fraction of phases 0-%d is %f, which not the unity\n", NumPhases, volTot);
735  }
736 }
737 
738 //help function for filling the J matrix
739 void Homogenize :: fillJ(FloatMatrix &J, double r, const FloatArray &mu, const FloatArray &k, int phase)
740 {
741  J(0, 0) = r;
742  J(0, 1) = 1. / ( pow(r, 2.) );
743  J(1, 0) = 3 * k(phase);
744  J(1, 1) = -4. * mu(phase) / ( pow(r, 3.) );
745 }
746 
747 //help function for filling the L matrix
748 void Homogenize :: fillL(FloatMatrix &L, double r, const FloatArray &mu, const FloatArray &k, int phase)
749 {
750  double mu_l = mu(phase);
751  double k_l = k(phase);
752  double nu_l = ( 3. * k_l - 2. * mu_l ) / ( 6. * k_l + 2. * mu_l ); //OK
753 
754  L(0, 0) = r;
755  L(0, 1) = -6. *nu_l *pow(r, 3.) / ( 1. - 2. * nu_l );
756  L(0, 2) = 3. / pow(r, 4.);
757  L(0, 3) = ( 5. - 4. * nu_l ) / ( ( 1. - 2. * nu_l ) * pow(r, 2.) );
758 
759  L(1, 0) = r;
760  L(1, 1) = -( 7. - 4. * nu_l ) * pow(r, 3.) / ( 1. - 2. * nu_l );
761  L(1, 2) = -2. / pow(r, 4.);
762  L(1, 3) = 2. / pow(r, 2.);
763 
764  L(2, 0) = mu_l;
765  L(2, 1) = 3. *nu_l *mu_l *pow(r, 2.) / ( 1. - 2. * nu_l );
766  L(2, 2) = -12. * mu_l / pow(r, 5.);
767  L(2, 3) = 2. * ( nu_l - 5. ) * mu_l / ( ( 1. - 2. * nu_l ) * pow(r, 3.) );
768 
769  L(3, 0) = mu_l;
770  L(3, 1) = -( 7. + 2 * nu_l ) * mu_l * pow(r, 2.) / ( 1. - 2. * nu_l );
771  L(3, 2) = 8. * mu_l / pow(r, 5.);
772  L(3, 3) = 2. * ( 1 + nu_l ) * mu_l / ( ( 1. - 2. * nu_l ) * pow(r, 3.) );
773 }
774 } // end namespace oofem
int giveNumberOfColumns() const
Returns number of columns of receiver.
Definition: floatmatrix.h:158
void hirsch(FloatMatrix &PhaseMatrix, double chi)
Hirsch&#39;s scheme combining Voigt and Reuss bounds.
Definition: homogenize.C:590
void kusterToksoz(FloatMatrix &PhaseMatrix)
Kuster-Toksoz model for two phases.
Definition: homogenize.C:648
void checkVolFraction(FloatMatrix &PhaseMatrix)
Check that the total volumetric fraction is the unity.
Definition: homogenize.C:724
double sqr(double x)
Definition: mathfem.h:112
Homogenize()
Constructor.
Definition: homogenize.C:39
#define E(p)
Definition: mdm.C:368
#define OOFEM_ERROR(...)
Definition: error.h:61
void ENuToKMu(const double E, const double nu, double &k, double &mu)
Convert Young&#39;s modulus and Poisson&#39;s ratio to bulk and shear moduli, isotropic material.
Definition: homogenize.C:691
double gamma(FloatMatrix &SortedPhaseMatrix, double zeta)
Auxiliary function.
Definition: homogenize.C:709
void fillL(FloatMatrix &L, double r, const FloatArray &mu, const FloatArray &k, int phase)
Auxiliary function for Herve-Zaoui scheme.
Definition: homogenize.C:748
void counto(FloatMatrix &PhaseMatrix)
Counto&#39;s model for a prismatic inclusion in a prismatic matrix, two phases.
Definition: homogenize.C:626
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
#define N(p, q)
Definition: mdm.C:367
double mu_hmg_2
Upper shear modulus if applicable.
Definition: homogenize.h:146
double k_hmg_2
Upper bound of bulk modulus if applicable.
Definition: homogenize.h:140
Class representing vector of real numbers.
Definition: floatarray.h:82
void selfConsistent(FloatMatrix &PhaseMatrix)
Self-consistent homogenization method of Hill and Budiansky for spherical isotropic inclusions...
Definition: homogenize.C:246
void hansen(FloatMatrix &PhaseMatrix)
Hansen&#39;s model for a spherical inclusion in a spherical matrix, for Poisson&#39;s ratio 0...
Definition: homogenize.C:605
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
double nu_hmg
Effective Poisson&#39;s ratio.
Definition: homogenize.h:131
double E_hmg
Effective Young&#39;s modulus or the lower bound.
Definition: homogenize.h:125
double cbrt(double x)
Returns the cubic root of x.
Definition: mathfem.h:109
void reuss(FloatMatrix &PhaseMatrix)
Serial scheme of Reuss for any number of isotropic phases.
Definition: homogenize.C:66
void moriTanaka(FloatMatrix &PhaseMatrix, int refRow)
Mori-Tanaka homogenization method for spherical isotropic inclusions.
Definition: homogenize.C:176
double nu_hmg_2
Effective Poisson&#39;s ratio or the lower bound.
Definition: homogenize.h:134
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
void hashinShtrikmanWalpole(FloatMatrix &PhaseMatrix)
Hashin-Shtrikman-Walpole lower and upper bounds for arbitrary number of isotropic phases...
Definition: homogenize.C:95
void printYourself(int out=0)
Print homogenized output.
Definition: homogenize.C:681
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
void fillJ(FloatMatrix &J, double r, const FloatArray &mu, const FloatArray &k, int phase)
Auxiliary function for Herve-Zaoui scheme.
Definition: homogenize.C:739
the oofem namespace is to define a context or scope in which all oofem names are defined.
void kMuToENu(const double k, const double mu, double &E, double &nu)
Convert bulk and shear moduli to Young&#39;s modulus and Poisson&#39;s ratio, isotropic material.
Definition: homogenize.C:697
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
void voigt(FloatMatrix &PhaseMatrix)
Parallel scheme of Voigt for any number of isotropic phases.
Definition: homogenize.C:43
double k_hmg
Effective bulk modulus or the lower bound.
Definition: homogenize.h:137
#define OOFEM_WARNING(...)
Definition: error.h:62
void herveZaoui(FloatMatrix &PhaseMatrix)
Herve and Zaoui&#39;s homogenization scheme for n-spherical isotropic domains and arbitrary number of iso...
Definition: homogenize.C:344
double zeta(double k, double mu)
Auxiliary function for Hashin-Shtrikman bounds.
Definition: homogenize.C:704
double E_hmg_2
Upper bound of Young&#39;s modulus if applicable.
Definition: homogenize.h:128
double mu_hmg
Effective shear modulus or the lower bound.
Definition: homogenize.h:143

This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:29 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011