Computational Materials Science 27 (2003) 446–460 www.elsevier.com/locate/commatsci
Wavelet-based homogenization of unidirectional multiscale composites Marcin Kami nski
*
Division of Mechanics of Materials, Technical University of Ł odz, Al. Politechniki 6, 93-590 Ł odz, Poland Received 15 April 2002; received in revised form 6 October 2002; accepted 6 January 2003
Abstract The main aim is to present a homogenization algorithm for the multiscale heterogeneous (composite) materials, which is based on the wavelet representation of material properties and the relevant multiscale reduction. It is shown that classical homogenization method used before for two-scale composites (with micro and macro scales) is a special case of general multiresolutional strategy, where a single scale parameter tends to 0. The approach presented is applied to unidirectional wavelet-based homogenization of linear elasticity heterogeneous problem and to wave propagation, which may be applied in conjunction with various discrete numerical methods for efficient modeling of heterogeneous solids, fluids and multiphase media. Ó 2003 Elsevier Science B.V. All rights reserved. Keywords: Multiresolutional analysis; Composite materials; Homogenization method; Symbolic computations; Elastostatics; Wave propagation
1. Introduction The homogenization method is still one of the most efficient approaches to computational modeling of composite systems. Usually it is assumed that there exists some scale relation between composite components and the entire system–– two-scales are introduced, which are related by a scale parameter being a small real value (tending most frequently to 0). Essential disadvantage of such techniques is impossibility of introducing *
Tel./fax: +48-42-63-13551. E-mail address:
[email protected] (M. Kami nski). URL: http://kmm.p.lodz.pl/pracownicy/Marcin_Kaminski/ index.html.
more than two various scales in composite and, further, the sensitivity of composite homogenized characteristics with respect to the geometrical scalesÕ interrelations. Multiscale methodology based on wavelet analysis [1,2,5–7,9], being a very modern and extensively developed numerical technique in signals theory, makes it possible to analyze the composite systems with multiple geometrical scales, which is more realistic for most of engineering composites (the scales of microdefects, interface, reinforcement and the entire structure). Wavelet analysis is especially promising tool in the area of composite materials, where it enables: (1) construction of the multiscale heterogeneous structures using particular wavelets what reflects perfectly the manufacturing process, for instance,
0927-0256/03/$ - see front matter Ó 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0927-0256(03)00046-6
M. Kaminski / Computational Materials Science 27 (2003) 446–460
and (2) multidimensional decomposition of the spatial distribution of composites material and physical properties by the use of the wavelets of various types introduced in different scales (heat conductivity or Young modulus. The first opportunity corresponds to the analysis of experimental results (image analysis of composite morphology), while the second one corresponds to the theoretical and computational analysis. Let us notice that the wavelet analysis introduces a new meaning for the term ÔcompositeÕ. In view of the analysis presented below we can dis-
447
tinguish homogeneous materials from composites using the following definition: Ôthe composite material and/or structure is such a continuum, whose material or physical properties are related in macro and microscales by at least a single Haar wavelet’. This definition extends traditional, engineering approach to composites, where laminated or fiber-reinforced structures were considered (with partially constant material characteristics) to the media with the 1D sinusoidal variability of these properties at least (see Figs. 1–6 below). Fig. 1 shows the spatial variability of the
Fig. 1. Distribution of the Young modulus in the real composite.
Fig. 2. Zeroth-order wavelet approximation of Young modulus in zeroth scale.
448
M. Kaminski / Computational Materials Science 27 (2003) 446–460
Fig. 3. First-order wavelet approximation of Young modulus in zeroth scale.
Fig. 4. Second-order wavelet approximation of Young modulus in zeroth scale.
Fig. 5. Second-order wavelet approximation of Young modulus in––first scale.
M. Kaminski / Computational Materials Science 27 (2003) 446–460
449
Fig. 6. Third-order wavelet approximation of Young modulus in––first scale.
Young modulus for the following wavelet function: 2 px 10 px eðxÞ ¼ e0 þ sin þ 0; 1 sin l l 4 10 px þ 0; 1 sin ; where l ¼ 10: l The next figures present the contributions of various scales to the macroscale elastic modulus of the entire composite structure. As it is demonstrated by Figs. 7 and 8, this modulus can be represented using some special
algebraic combinations of the basic wavelets (Haar, Mexican hat and/or sinusoidal waves). Spatial variability of Young modulus for the two component composite with and without some interphase can be computationally simulated using a theoretical approximation of spatial variability of this modulus and symbolic computations mathematical package MAPLE, for instance. We consider for illustration the Representative Volume Element (RVE) of two-component composite with the following elastic characteristics: e1 ¼ 209e9 and e2 ¼ 209e8 with the RVE length l ¼ 1:0 and
Fig. 7. Wavelet approximation of the elastic properties of two-component composite.
450
M. Kaminski / Computational Materials Science 27 (2003) 446–460
Fig. 8. Wavelet approximation of the elastic properties of two-component composite with interface defects.
volume fractions being equal for both components. The following wavelet functions are proposed to model the multiresolutional character of Young modulus spatial variability in the RVE (without the interface defects, cf. Fig. 7): eðxÞ ¼ hðxÞ þ 0:2 1010 sinð5 101 xÞ þ 0:2 1010 sinð5 104 xÞ for hðxÞ being the Haar wavelet basis. It can be noticed that thanks to the multiscale character of this function the illustration of composite Young modulus shows the randomness on its microscale as well, however the character of spatial variability of this modulus is still deterministic. Furthermore, much more complicated and sometimes more realistic effects can be modeled by this way in composites––the RVE can be almost damaged at the interface (according to ageing and fatigue processes) and the spatial distribution of elastic properties can vary along the specimen heterogeneity main axis. It is approximated by a combination of Haar, some sinusoidal and the socalled Mexican hat wavelets (see Fig. 8) in the following manner: eðxÞ ¼ hðxÞ þ 0:2 1010 sinð5 101 xÞ þ 0:2 1010 sinð5 104 xÞ þ 0:6 1010 0:76597 1011 expð8:00 x2 Þ ð16:0 x2 1Þ:
Algebraic structure of this wavelet is relatively complicated, however (1) it illustrates the opportunities for wavelet-based approximation of mechanical and physical properties of the real composites, (2) it can be used together with structural image analysis tools for the relevant analyses of composites and (3) it enables direct symbolic computational homogenization of such media. As long as this composite is unidirectional, some classical homogenization closed-form equations can be used to construct an equivalent medium using the corresponding differential equilibrium equations directly. Structural sensitivity analysis can be carried out with respect to the variabilities of material properties in quite different scales of composite; it can be performed similarly to the considerations presented in [6]. The situation complicates significantly in case of planar distribution of material tensors, where the composite cell problems are to be solved by wavelet decomposition and construction to determine the effective behavior of the entire composite. As it is mathematically proved in the paper, when the structure is heterogeneous in many scales, the effective elastic modulus differs from that obtained for the corresponding classical two-scale and twocomponent composite beam. Another disadvantage of the wavelet-based analysis of composite materials is the assumption
M. Kaminski / Computational Materials Science 27 (2003) 446–460
of very arbitrary character, which means that the physical model and the accompanying equations of thermodynamical equilibrium have exactly the same form in each scale of considered medium, what follows a purely mathematical nature of the wavelet transform. It eliminates the opportunity of the physical transition from the particle scale (atomistic forces and Van der Waals interactions) through chemical interface reactions in various composites to the global scale of the entire engineering structure, which is consistent with the concept that the transition between the corresponding medium scale must strongly depend on the physical scale we are operating with. The approach discussed is compared theoretically and numerically with spatial averaging and classical now asymptotic homogenization. The classical result obtained through the asymptotic homogenization theory is given for deterministic composites exhibiting two separate geometrical scales linked by the scale parameter e, which is the weakest point of this method. Sometimes e is treated as a positive real number tending to 0 (practically infinite number of the RVEs in the composite) and, alternatively, as some small positive parameter. Now, this parameter is treated as some real function introduced as the wavelet function relating two or more separate geometrical scale of the composite. As it is demonstrated, the essential differences are observed in these two models. In contrast to the classical approach to the homogenization problem, the multiresolutional technique uses the algebraic transformation between scales provided by the wavelet analysis. It is done to determine the fine-scale behavior and introduce it explicitly to the macroscopic equation and, therefore, the coefficients may vary on arbitrary many scales. The chain of subspaces . . . V2 V1 V0 V1 V2 . . .
ð1Þ
defines the hierarchy of scales used in the multiresolutional scheme. This chain of subspaces is defined in such a way than the space Vj is ‘‘finer’’ than the space Vjþ1 in the sense that (1) all of Vjþ1 is contained in Vj , and (2) the component of Vj , which is not in Vjþ1 . It consists of the functions, that resolve features on a scale finer than any function in
451
Vjþ1 may resolve. The difference between successive spaces in this chain is captured by the so-called wavelet function Wjþ1 , defined as the orthogonal complement of Vjþ1 in Vj . An orthogonal basis for the wavelet space Wjþ1 is constructed, which has vanishing moments, i.e. the basis elements are L2 orthogonal to low-degree polynomials. The existence of orthogonal wavelet-bases with vanishing moments distinguishes the multiresolutional approach from typical multi-scale discretizations provided by finite-element or hierarchical bases. If we are considering a multiresolutional analysis defined on a bounded domain, then the hierarchy of scales defined above has the coarsest scale (which is called V0 ) and it can be written as V0 V1 V2 . . .
ð2Þ
Such an assumption is consistent with traditional understanding of composites, where each time the smallest scale is introduced.
2. Elements of the multiscale reduction Let us review the multiresolution strategy for the reduction and homogenization of linear and unidirectional problems. Let us consider to this purpose a bounded linear operator Sj : Vj ! Vj . Since Vj is spanned by translations of the function /ð2j x kÞ, then it can be written as a matrix and as long as the multiresolutional analysis is defined on a bounded domain, this matrix is finite. Let us consider the equation Sj x ¼ f :
ð3Þ
The decomposition Vj ¼ Vjþ1 Wjþ1 allows us to split the operator Sj into four pieces (Wjþ1 is called the wavelet space and is the ‘‘detail’’ or finescale component of Vj ) and to write ASj BSj df dx ¼ ; ð4Þ CSj TSj sf sx where ASj : Wjþ1 ! Wjþ1 ; BSj : Vjþ1 ! Wjþ1 ; CSj : Wjþ1 ! Vjþ1 ; TSj : Vjþ1 ! Vjþ1 ;
ð5Þ
452
M. Kaminski / Computational Materials Science 27 (2003) 446–460
and dx , dy 2 Wjþ1 , sx , sy 2 Vjþ1 are the L2 -orthogonal projections of x and f on the Wjþ1 and Vjþ1 spaces. The projection sx is thus the coarse-scale component of the solution x and dx is its fine-scale component. Formally, elimination of dx from Eq. (4) by substituting dx ¼ A1 Sj ðdf BSj sx Þ yields 1 ðTSj CSj A1 Sj BSj Þsx ¼ sf CSj ASj df :
ð6Þ
This equation is called the reduced equation, while the operator RSj ¼ TSj CSj A1 Sj BSj
ð7Þ
is one step reduction of the operator Sj known also as the Schur complement of the block matrix ASj BSj . Let us note that the solution sx of the CSj TSj reduced equation is exactly Pjþ1 x, where Pjþ1 is the projection on Vjþ1 and x is the solution of Eq. (3). Let us note that the reduced equation is not the same as the averaged equation, which is given by TSj ~sx ¼ sf :
ð8Þ
Once we have obtained the reduced equation, it may formally be reduced again to produce an equation on Vjþ2 and the solution of this equation is just the Vjþ2 -component of the solution for Eq. (4) and etc. These equations can be reduced recursively n times (assuming that if the multiresolutional analysis is on a bounded domain, then j þ n 6 0) to produce an equation on Vjþn , the solution of which is the projection of the solution of Eq. (3) on Vjþn . It should be mentioned that in the finitedimensional case, if multiresolutional analysis defined on a domain in R is considered, the reduced Eq. (4) has half as many unknowns as the original equation (3) (for R2 the reduced equations have one-fourth as many unknowns as the original equation). A reduction, therefore, preserves the coarse-scale behavior of solutions significantly reducing the number of unknowns. In order to iterate the reduction step over many scales, the form of the equation is to be preserved as a way of deriving a recurrence relation. In Eqs. (3) and (4), both Sj and RSj are matrices, and thus the procedure may be repeated. However, identi-
fication of the matrix structure is usually not sufficient. In particular, even though the typical matrix A for ordinary differential equations (ODEs) and PDEs is sparse, the component A1 Sj term may become dense, changing the equation from a local one to a global one. It is important to know under what, if any, circumstances the local nature of the differential operator may be (approximately) preserved. Furthermore, if the equation is of the form of d d eðxÞ uðxÞ ¼ f ðxÞ ð9Þ dx dx or some other variable–coefficient differential equation, it is necessary to verify if the reduction procedure preserves this form, so that effective coefficients of the equation on the coarse-scale can be found. This is the basic goal of homogenization techniques and it extracts information from the reduced equation based on the original equation. Thus, a reduction and homogenization are closely related in the multiresolution approach, however have different goals: homogenization––to find effective equations and their coefficients on the coarse-scale, whereas reduction––to find a coarsescale version of a given system of equilibrium equations.
3. Multiresolutional homogenization equations The multiresolutional homogenization procedure is applied to all the ODEs, which may be written in the form Bx þ q þ k ¼ KðAx þ pÞ:
ð10Þ
Further, the equations of the following form are analyzed now: ðI þ BðtÞÞxðtÞ þ qðtÞ þ k Z t ¼ ðAðsÞxðsÞ þ pðsÞÞ ds;
t 2 ð0; 1Þ
ð11Þ
0
on the domain L2 ð0; 1Þ, where BðtÞ and AðtÞ are the n-rank matrix-valued functions, pðtÞ and qðtÞ are vector-valued forcing terms, while xðtÞ is the so-
M. Kaminski / Computational Materials Science 27 (2003) 446–460
lution vector. As a differential equation it can be formulated as d ððI þ BðtÞÞxðtÞ þ qðtÞÞ ¼ AðtÞxðtÞ þ pðtÞ dt
ð12Þ
with the initial conditions ðI þ Bð0ÞÞxð0Þ ¼ qð0Þ k. The projection of Eq. (11) on Vj , j < 0, is written as
453
for Bj and Aj their ith diagonal blocks by ðBj Þi and ðAj Þi , respectively. The matrices are given by the Haar coefficients of the n-rank matrix-valued functions BðxÞ and AðxÞ on the scale Vj i¼2j 1
Bj ¼ diagfI þ ðBj Þi gi¼0
ð17Þ
and j
Bj xj þ qj þ k ¼ Kj ðAj xj þ pj Þ
ð13Þ
i¼2 1 ; Aj ¼ diagfðAj Þi gi¼0
ð18Þ
where
or
0
Sj xj ¼ fj ;
ð14Þ
where Sj ¼ Bj Kj Aj ;
fj ¼ Kj pj qj k:
ð15Þ
After a single reduction, the equation on Vjþ1 is derived in the form ðjÞ ðjÞ ðjÞ ðjÞ ðjÞ ðjÞ Bjþ1 xjþ1 þ qjþ1 þ k ¼ Kjþ1 Ajþ1 xjþ1 þ pjþ1 ; ð16Þ ðjÞ
ðjÞ
where xjþ1 ¼ Pjþ1 xj , Bjþ1 ¼ Pjþ1 Bj etc. The notaðjÞ tion B1 is used to indicate that the operator is first projected to scale Vj , and then the reduction procedure is applied 1 j times to obtain an equation on V1 . Therefore, this notation indicates that the formula (16) is obtained by a single reduction of the same form equation on Vj to produce an equation on coarser scale Vjþ1 . Generally, it enables one to establish a recurrence relation for k ¼ j; j þ 1; . . . ; 0 between the ðjÞ ðjÞ ðjÞ ðjÞ operators and forcing terms Bk , xk , pk , qk on Vkþ1 . The recurrence relations are simplified significantly if multiresolutional analysis is used with the basis functions having non-overlapping support; the Haar basis is proposed here, but any multiwavelet basis may be used instead if only higher order elements are needed. In the Haar basis, the operators Bj , Aj and Kj derived from Eq. (13) have a simple form and are of Nj n rank, where Nj ¼ 2j is the number of unknowns on the scale Vj and n denotes the number of equations in the original system. Furthermore, Bj and Aj are both block-diagonal and n-rank matrices and, therefore, there are Nj diagonal blocks, each of which is of n-rank. Let us denote
1 I 2
B I B Kj ¼ dj B B I @... I
0 1 I 2 I
0 0 1 I 2
... ... ... I
...
1 0 0 C C ...C C A 1 I 2
ð19Þ
for dj ¼ 2j ; I is the n n identity matrix, ðBj Þi and ðAj Þi are the ith Haar coefficients on scale Vj of the n-rank matrix-values functions BðxÞ as well as AðxÞ. The following recursive relations can be written for Eq. (16) dk ðjÞ 1 Akþ1 ¼ ðSA Þi ðDA Þi F ðDB Þi þ ðSA Þi ; i 2 ð20Þ dk ðjÞ Bkþ1 ¼ ðSB Þi ðDA Þi ðDB Þi i 2 dk dk 1 ðSA Þi F ðDB Þi þ ðSA Þi ; ð21Þ 2 2 dk ðjÞ pkþ1 ¼ ðSp Þi ðDA Þi F 1 ððDq Þi þ ðSp Þi Þ; ð22Þ i 2 dk dk ðjÞ qkþ1 ¼ ðSq Þi ðDp Þi ððDB Þi ðSA Þi ÞF 1 i 2 2 ððDq Þi ðSp Þi Þ; ð23Þ where ðjÞ
ðjÞ
ðSA Þi ¼ 12ððAk Þ2i þ ðAk Þ2iþ1 Þ; ðjÞ
ðjÞ
ðDA Þi ¼ 12ððAk Þ2i ðAk Þ2iþ1 Þ; ðjÞ
ð24Þ
ðjÞ
ðSB Þi ¼ 12ððBk Þ2i þ ðBk Þ2iþ1 Þ; ðjÞ
ðjÞ
ðDB Þi ¼ 12ððBk Þ2i ðBk Þ2iþ1 Þ;
ð25Þ
454
M. Kaminski / Computational Materials Science 27 (2003) 446–460 ðjÞ
ð1Þ
ðjÞ
ðSp Þi ¼ p1ffiffi2ððpk Þ2i þ ðpk Þ2iþ1 Þ; ðjÞ
ðjÞ
ðDp Þi ¼ p1ffiffi2ððpk Þ2i ðpk Þ2iþ1 Þ; ðjÞ
ð26Þ
ð1Þ
Ah ¼ A0
ðjÞ
ðSq Þi ¼ p1ffiffi2ððqk Þ2i þ ðqk Þ2iþ1 Þ; ðjÞ
ðjÞ
ðDq Þi ¼ p1ffiffi2ððqk Þ2i ðqk Þ2iþ1 Þ;
ð27Þ
F ¼ I þ ðSB Þi þ
dk ðDA Þi : 2
ð28Þ
The recurrence relations are local and can be carried out at as many scales as necessary (assuming the existence of F 1 at each scale). Starting from Eq. (16) on Vj and, reducing it j times, it yields on V0 ðjÞ ðjÞ ðjÞ ðjÞ ðjÞ ðjÞ B0 x0 þ q0 þ k ¼ K0 A0 x0 þ p0 ; ð29Þ where the recurrence relations is applied j times to ðjÞ ðjÞ ðjÞ ðjÞ compute B0 , A0 , p0 , q0 . Multiresolutional homogenization theorem is obtained as a limit of relation (29) for j ! 1 ð1Þ ð1Þ ð1Þ ð1Þ ð1Þ ð1Þ B0 x0 þ q0 k ¼ K0 A0 x0 þ p0 ; ð30Þ which enables us to eliminate infinite number of ð1Þ ð1Þ the geometrical scales; the matrices B0 , A0 are the reduced coefficients of Eq. (16). Then, the operators and forcing terms Bh ðtÞ, Ah ðtÞ, ph ðtÞ, qh ðtÞ are to be found, such that the equation ðI þ Bh ðtÞÞxðtÞ þ qh ðtÞ þ k Z t ¼ ðAh ðsÞxðsÞ þ ph ðsÞÞ ds;
t 2 ð0; 1Þ
ð31Þ
0
subjected to the same reduction and limit procedure as (11), results for V0 with the same equation as Eq. (30). For the relation (11) we usually require that Bh , Ah , ph , qh be constant. The results of homogenization in this case are summarized as follows. Homogenization theorem. If the limits defining the ð1Þ ð1Þ matrices B0 and A0 exist, then there exist h constant matrices B , Ah and forcing terms ph , qh , such that the reduced coefficients and forcing terms
ð1Þ
ð1Þ
ð32Þ
;
e 1 I; Bh ¼ Ah A
ð33Þ
ð1Þ
ð34Þ
p h ¼ p0
and, finally
ð1Þ
of Eq. (31) are given by B0 , A0 , q0 , p0 . The homogenized coefficients Bh , Ah and forcing terms ph , qh are defined by
ð1Þ
qh ¼ q0
;
eA e ðexpð A e IÞ Ah Þ ph Þ; þ ðI 12 A 1
1
ð35Þ where e ¼ logðI þ ðI þ Bð1Þ 1Ah ÞI Ah Þ: A 0 2
ð36Þ
Proof. The recursive relations (20) and (21) simplify for the constant coefficients to Ahkþ1 ¼ Ahk ;
ð37Þ
Bhkþ1 ¼ Bhk þ
d2k h 1 A ðI þ Bhk Þ Ahk : 4 k
ð38Þ
The forcing terms can be rewritten as h ¼ pkh ; pkþ1
ð39Þ
dk h 1 A ðI þ Bhk Þ pkh : ð40Þ 2 k Since Ah is unchanged by the reduction, Ah ¼ ð1Þ A0 . The reduction does not change ph , hence ð1Þ ph ¼ p0 . The matrices Bh and qh are determined analytically from the solution of Eq. (31). In case ð1Þ ð1Þ p0 ¼ 0, there holds qh ¼ q0 . Therefore, the solution of Eq. (31) is given by
qhkþ1 ¼ qhk
e tÞ~ xðtÞ ¼ expð A q;
ð41Þ
e ¼ ðI þ Bh Þ1 Ah , q~ ¼ ðI þ Bh Þ1 ðqh þ kÞ. where A The average of this solution must also solve Eq. (30) (equation for the average value of the solution from definition). The average value of xðtÞ in Eq. (31) on the interval [0,1] is given by Z 1 e e ÞÞ A e 1 q~: hxi ¼ expð A tÞ dt q~ ¼ ðI expð A 0
ð42Þ
M. Kaminski / Computational Materials Science 27 (2003) 446–460
The solution to Eq. (30) is given by ð1Þ
x0
1 1 ð1Þ ð1Þ ð1Þ ¼ I þ B0 A0 q0 þk : 2 ð43Þ
The right hand sides of Eqs. (42) and (43) are demonstrated to be equal for all k. The solution given in the statement of the theorem can be obtained by setting k ¼ 0 and solving for Bh ; a soð1Þ lution process in case of p0 6¼ 0 proceeds quite similarly. Solutions of Eq. (31) have the same ‘‘average’’ or coarse-scale behavior as the solutions of Eq. (11). The main point is that this homogenization procedure allows coefficients to vary on arbitrarily many intermediate scales. This is a contrast to the classical homogenization examples, which did not allow for any intermediate scales. As it has been formulated above, the multiresolution approach to homogenization requires the ð1Þ ð1Þ computation of A0 and B0 , i.e. the limits over infinitely many scales. The typical practice is to ðJ Þ ðJ Þ compute successive A0 and B0 terms until finer approximations vary by less than some specified tolerance, and to use these matrices as approxið1Þ ð1Þ mations to A0 and B0 . Apart from establishing the general framework for multiresolution reduction and homogenization, it is observed that the usage the Haar basis (or a multiwavelet basis) for systems of linear ODEs, provides a technical advantage. Since the Haar functions on a fixed scale do not have overlapping supports, the recurrence relations for the operators and forcing terms in the equation may be written as local relations and solved explicitly. Thus, an explicit local reduction and homogenization procedure is possible for the ordinary differential equations; Eq. (9), for instance, may be rewritten as the coupled first-order system 8 d > > < dx vðxÞ ¼ f ðxÞ; > > : d uðxÞ ¼ eðxÞ1 vðxÞ: dx
ð44Þ
455
By writing in an integral form one can obtain ! Z x 1 uðtÞ uðxÞ uð0Þ 0 eðtÞ ¼ vðtÞ vðxÞ vð0Þ 0 0 0 ! 0 dt ð45Þ þ f ðtÞ
Thus, in Eq. (11), there BðxÞ ¼ 0, holds 1 uð0Þ 0 eðxÞ AðxÞ ¼ , k¼ and qðtÞ ¼ 0 0 0 vð0Þ 0 as well as pðtÞ ¼ . Using the reduction f ðtÞ procedure based on the Haar wavelets for a system of linear differential equations, the goal is to find constants Bh , Ah , ph , qh such that uðxÞ ðI þ Bh Þ þ qh þ k vðxÞ Z x uðtÞ ¼ þ ph dt Ah ð46Þ vðtÞ 0 which, after reduction to the scale V0 , will be the same as Eq. (45) reduced to that scale. This is accomplished by solving the recursion relations between the operators in the reduced equations explicitly, element-by-element in each matrix, which is possible because of the non-overlapping supports of the Haar basis functions on a fixed scale. The result for the first two coefficients is 0 0 0 M1 2M2 Bh ¼ ; Ah ¼ ; ð47Þ 0 0 0 0 where M1 ¼
Z
1 0
1 dt; eðtÞ
M2 ¼
Z
1 0
t 12 dt: eðtÞ
ð48Þ
Similar expressions for ph and qh can be found; let us note that ph ¼ qh ¼ 0 if f ðxÞ ¼ 0 identically. Furthermore, in general Bh , Ah do not depend on p and q. As the first-order system of ODEs, the homogenized equation yields 8 d > < vðxÞ ¼ f h ðxÞ; dx ð49Þ > : d uðxÞ ¼ ðM1 2M2 ÞvðxÞ ¼ 1 vðxÞ; dx eh what differs from the classical result. This difference follows the fact that the multiresolution homogenization procedure allows the coefficients
456
M. Kaminski / Computational Materials Science 27 (2003) 446–460
eðxÞ to vary on arbitrary many scales, whereas the classical approach allows only for coefficients of the form eðx=eÞ. In the multiresolution context this amounts to restricting the coefficients to an asymptotically fine-scale. Let us apply the same limit in the preceding section to the coefficients appearing in the multiresolution approach. We start with the coefficients of the form eðx=eÞ. Applying this homogenization scheme to the elliptic equation with these coefficients yields two terms, M1 ðeÞ, M2 ðeÞ. It is found for e ! 0 that lim M1 ðeÞ ¼ M1 e!0
ð50Þ Fig. 9. Wavelet-based definition of elastic modulus in the RVE.
and lim M2 ðeÞ ¼ 0: e!0
ð51Þ
Thus, the factor M2 is present in the multiresolution context, but it does not appear in the classical approach, and it is zero when the limit found in the classical method is applied to the result of the multiresolution methodology. Let us note that the formula for the homogenized parameter eðhÞ introduces new, closer bounds on the wavelet function defining material parameters than formulation: inR L it is done byRclassical L tegrals 0 dx=eðxÞ and 0 x dx=eðxÞ must be real and eðxÞ must be positively defined to assure homogenizability of the problem. The counterexample is the family of sinusoidal wavelets of the form eðxÞ ¼ e0 þ a sinðpx=LÞ, where e0 , a, L 2 R. Taking for example L ¼ 10, e0 ¼ 20R and a ¼ 0:1, a L symbolic integration returns ðx dx=eðxÞÞ ¼ 0 2:99242 19:07199i. Classical small parameter homogenization method should be applied in that case; otherwise another wavelet decomposition of the real composite material properties is to be proposed. The formulas presented above are implemented in the program MAPLE together with the spatial averaging method in order to compare the homogenized modulus computed by various ways (spatial averaging, classical and multiresolutional) for the same composite. Fig. 9 illustrates the variability of this modulus along the RVE, where the function eðxÞ is extracted from the following Haar and Mexican hat wavelets:
hðxÞ ¼
20:0E9; 0 6 x 6 0:5; 2:0E9; 0:5 < x 6 1;
1 x2 exp mðxÞ ¼ 2 þ pffiffiffiffiffiffi 2 3 2pr r 1
ð52Þ x2 ; 2r2
r ¼ 0:4; ð53Þ
as eðxÞ ¼ 10:0 hðxÞ þ 2:0E9 mðxÞ:
ð54Þ
The final form of these functions is established on the basis of the mathematical conditions for homogenizability analyzed before as well as to obtain final variability of composite properties similar to the traditional multicomponent structures (classical definition of periodic composite material properties contained the piecewise constant characteristics only). The following homogenized material properties are obtained from this input: eðavÞ ¼ 114:548E9 > eðeff;wavÞ ¼ 60:217E9 > eðeffÞ ¼ 35:437E9, which means that for this particular example, the highest value is obtained for spatial averaging method, then––for multiresolutional approach and at last––for classical asymptotic homogenization method.
4. Multiscale homogenization of the wave equation Let us consider the following ODE corresponding to unidirectional acoustic waves propagation in a multiscale medium with uniaxial and
M. Kaminski / Computational Materials Science 27 (2003) 446–460
multiresolutional distribution of non-homogeneities [2]:
ð1Þ ð1J Þ ¼ lim Aj ; A1 j
J !1
457
ð1Þ
B1
j
ð1J Þ
¼ lim Bj J !1
:
ð64Þ
d uðxÞ ¼ ixMðxÞuðxÞ; dx
x 2 ½0; 1;
ð55Þ
where physical coefficients MðxÞ for both composite layers used in the Haar wavelet approximation are defined by M0 ; 0 6 x < 12 ; MðxÞ ¼ ð56Þ M1 ; 12 6 x 6 1:
After some algebra it is obtained that [2] ð1Þ A1 ¼ Aj ; j ! 1 Aj Aj 1 ð1Þ exp I B1 ¼ þ I I; j 2 2 2 j ¼ 1; 2;
ð65Þ
It should be mentioned that the equations collected are derived for equal volume ratios of both layers. Otherwise, they should be complemented with the volume parameters of composite components c1 and c2 . Then, the corresponding homogenized equation can be rewritten for the deterministic system as
where the following expansion is used: Aj ix exp Mj ¼ exp 2 2
d uðxÞ ¼ K ðeffÞ uðxÞ: dx
Taking into account that the coefficients B and A represent physical properties of the composite components with the total number of various scales tending to infinity, it is possible to determine an analogous definition of the homogenized coefficient for a composite with some finite scales number; analogous considerations dealing with the probabilistic effective characteristics of various engineering structures are demonstrated in [5]. Computational experiments are performed using the symbolic computations system MAPLE 7 and additional implementation of the multiresolution homogenization analysis. Basic computations are carried out with respect to interrelation between physical constants of both layers as well as expansion order. Finally, let us observe that a homogenized system, both in terms of deterministic or stochastic effective coefficients, can be analyzed numerically using classical finite element method (FEM) [10] for instance, or by application of various stochastic numerical methods (random simulation, perturbation-based or spectral analysis [3]). Homogenization-based numerical approach will considerably speed up the process of computational modeling of composites and, in case of very complex multiscale heterogeneous media, it can be the only available method. Further, real and imaginary parts of K ðeffÞ are computed according to Eqs. (58)–(60). The
ð57Þ
It can be demonstrated that the homogenized coefficient K ðeffÞ is equal to ð1Þ
K ðeffÞ ¼ logðI þ ðI þ B0
ð1Þ 1
12A0
ð1Þ
Þ A0
Þ ð58Þ
for ð1Þ
B0
¼ SB0 14D0A ðD0B 14SA0 ÞF 1 ðD0B þ 14SA0 Þ ð59Þ
and ð1Þ
A0
¼ SA0 D0A F 1 ðD0B þ 14SA0 Þ:
ð60Þ
The right hand side coefficients denote ð1Þ ð1Þ SA0 ¼ 12 ðA1 Þ0 þ ðA1 Þ1 ; ð1Þ ð1Þ D0A ¼ 12 ðA1 Þ0 ðA1 Þ1 ;
ð61Þ
ð1Þ ð1Þ SB0 ¼ 12 ðB1 Þ0 þ ðB1 Þ1 ; ð1Þ ð1Þ D0B ¼ 12 ðB1 Þ0 ðB1 Þ1 ;
ð62Þ
F ¼ I þ SB0 þ 14D0A ;
ð63Þ
where
ix x2 ix3 3 Mj Mj2 þ M 2 8 48 j þ þ Oðxn Þ:
¼Iþ
ð66Þ
458
M. Kaminski / Computational Materials Science 27 (2003) 446–460
following data are adopted: M0 ¼ 10:0, M1 ¼ 1:0 with c1 ¼ c2 ¼ 0:5, where the parameter x ! 0 (cf. Figs. 10 and 11) and x ! 1 (see Figs. 12 and 13). As it can be observed, the real and imaginary parts tend to 0 in both cases, which finally gives K ðeffÞ ! 0, too. Further, such a combination of input parameters results in a minimum of K ðeffÞ real part for x 1:15. On the other hand, the singularity of ImðK ðeffÞ Þ is obtained with x 0:75 (Fig. 11). Next, the effective parameter in its real and imaginary part is determined as a function of x Fig. 10. Real part of K ðeffÞ near 0.
Fig. 13. Imaginary part of K ðeffÞ in x domain.
Fig. 11. Imaginary part of K ðeffÞ near 0.
Fig. 12. Real part of K ðeffÞ in x domain.
Fig. 14. Real part of K ðeffÞ .
M. Kaminski / Computational Materials Science 27 (2003) 446–460
Fig. 15. Imaginary part of K ðeffÞ .
value and the ratio relating material parameters of the composite components M0 ¼ 2; . . . ; 20. The results are presented in Figs. 14 and 15. As it can be seen from the comparison with Figs. 10 and 11, the material parameters interrelation influences significantly the effective parameters in the same range as x values. Analogous limiting values in real and imaginary parts of the homogenized parameter as well as imaginary part singularities are noticed as nonlinear functions of both design parameters of the study. Probabilistic moments of real and imaginary surfaces are expected in the probabilistic case, however a more important problem (from the physical point of view) is to determine the relations for homogenized coefficients in terms of volume fractions of the layers as well as to extend the homogenization method to the heterogeneous multiscale media with a more general periodic geometry of the RVE. The entire methodology can be adopted without any changes in the stochastic reliability studies, where probabilistic coefficients of the effective properties or the state functions can be used in the reliability index computations.
5. Concluding remarks 1. As it was demonstrated above, the waveletbased multiresolutional computational techniques can be very efficient, considering the
459
capability of making heterogeneity analysis for extremely different geometrical scales in the same time. Such phenomena appear frequently in engineering composites––at the interface between the components, in microscale connected with the periodicity cell, for window at mesoscale for a couple of reinforcing fibers or particles as well as for the macroscale connected with the global composite structure. As it can be observed, the wavelet-based numerical methods [4] (FEM, especially) can be successfully used even for the heterogeneous media with random or stochastic microstructure thanks to implementation of a randomization method (simple algebra, PDF integration, Monte-Carlo simulation, stochastic perturbation or even spectral analyses). 2. The homogenization method discussed in the paper enables us to apply an alternative approach, where the effective material parameters (or their probabilistic moments) are determined first and then the entire composite is analyzed using traditional computational techniques. Wavelet-based multiresolutional approach to the homogenization problem should however be formulated to introduce the components characteristics on many scales into the final effective structural parameters. As it was demonstrated in the mathematical considerations, homogenized properties in multiscale analysis and classical macro-micro passage are essentially different, even in deterministic formulation, which was observed previously in three scale MonteCarlo simulation based homogenization studies for the fiber-reinforced composites [8]. 3. Finally, let us note that due to the character of the homogenized modulus for the 1D elastostatic problem [5], the results obtained can be applied without any further modifications in heat transfer [5] and/or free vibrations [7] for a composite with exactly the same multiscale internal structure as well as for any linear field problem. The real and imaginary parts of the effective coefficient for wave equation can be used in acoustic waves propagation for random media. It is observed that for wave propagation, homogenized coefficients strongly depend in the same range on angular velocity and the
460
M. Kaminski / Computational Materials Science 27 (2003) 446–460
interrelation of material properties of the layered medium components. [6]
References [7] [1] M.A. Christon, D.W. Roach, The numerical performance of wavelets for PDEs: the multi-scale finite element, Comput. Mech. 25 (2000) 230–244. [2] N.A. Coult, A Multiresolutional Strategy for Homogenization of Partial Differential Equations, Ph.D. Thesis, Colorado University, 1997. [3] R.G. Ghanem, P.D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer-Verlag, 1991. [4] L. Jameson et al., Wavelet-based numerical methods, Int. J. Comput. Fluid Dyn. 10 (1998) 267–280. [5] M. Kami nski, Multiresolutional homogenization technique in transient heat transfer for unidirectional composites, in:
[8]
[9] [10]
B.H.V. Topping, Z. Bittnar (Eds.), Proc. Sixth Int. Conf. Comput. Struct. Technology, Prague, Civil-Comp Press, 2002. M. Kami nski, Multiresolutional wavelet-based homogenization of random composites, in: Proc. Second ECCOMAS Conference, Krak ow, 2001. M. Kami nski, Wavelet-based finite element elastodynamic analysis of composite beams, in: H.A. Mang, F.G. Rammerstorfer, J. Eberhardsteiner (Eds.), Proc. Fifth World Congress on Computational Mechanics, Vienna, ISBN 39501554-0-6, 2002. Available from
. M. Kami nski, M. Kleiber, Numerical homogenization of n-phase composites including stochastic interface defects, Int. J. Numer. Meth. Engrg. 47 (2000) 1001–1025. R.K. Young, Wavelet Theory and Its Applications, Kluwer, 1993. O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, fifth ed., McGraw-Hill, London, 1997.