- Email: [email protected]
Material database optimization using Monte Carlo method in engineering application
Probabilistic Engineering Mechanics 16 (2001) 349±354
www.elsevier.com/locate/probengmech
Material database optimization using Monte Carlo method in engineering application Tadeusz Niezgoda* Faculty of Mechanics, Military University of Technology, Kaliskiego Street 2, 00-908 Warsaw, Poland
Abstract Ceramic materials were simulated numerically. The practical goal of the research was to assess numerically the average residual stress value in such materials. The scienti®c goal of the work was to establish a basic analytical apparatus. One of the main aspects was to create material database, enabling simulation of polycrystalline structure with random orientations of thermal expansion orthotropy axes. Monte Carlo method was adopted for this purpose. The calculation results were analyzed with the use of statistical method, which is analogous to experimental measure analysis. However, the developed methods can be readily used to analyze both ceramics and ceramic composites. Currently, the methods can take advantage of more precise geometrical description to achieve more accurate models re¯ecting real internal structure of the material. q 2001 Published by Elsevier Science Ltd. Keywords: Finite element method; Monte Carlo method; Al2O3 ceramics; Numerical method; Residual stress
1. Introduction Ceramics and ceramic composites are materials with very interesting properties. They are often used in structures that are subject to high temperatures, friction and erosion processes caused by chemically aggressive environments or high-speed ¯ow of gases. For example, thin ceramic coatings protect elements of turbines, jet and internal combustion engines from excessive heat induced loading, and ceramic materials are often used as cutting tools. Porous ceramic materials may be used as membranes for gas ®ltration and separation. Both ceramics and ceramic composites are polycrystalline materials. Each grain in such a structure is thermally orthotropic. In cooling down polycrystalline ceramic materials, from sintering to room temperature, the random variation in crystallographic orientation of orthotropy from one grain to another induces a very complex state of constraint. Thus, after completion of the manufacturing process, a state of residual stress remains. This limits the range of acceptable loads on these ceramic materials. The complicated spatial variation of the residual stress ®eld depends on the elastic properties of grains, the grain size and the shape distribution. In some cases, the residual stress across the grain
* Tel./fax: 148-22-683-94-61. E-mail address: [email protected] (T. Niezgoda). 0266-8920/01/$ - see front matter q 2001 Published by Elsevier Science Ltd. PII: S 0266-892 0(01)00015-7
boundary can be suf®ciently large to initialize microcracking across the grain boundary. In elastic solids this residual stress s T is given by [1]
s T bEDaDT
1
where b is a coef®cient that depends on the shape of the grain and the grain orientation b [ k 2 0:5; 10:5l; E the Young's modulus, Da the difference in the singlecrystal thermal expansion coef®cients along the direction considered, DT the difference between the temperature Tc (at which the sample was stress-free) and the room temperature T0. The orthotropy of thermal expansion coef®cients is caused by building of crystallographic lattice of Al2O3 ceramic grain. For orthotropic grains aa am ± ac and hence, Da ac 2 aa for Al2O3 ceramics equals Da 0:9 £ 1026 K21 : The average residual stress value is a deterministic quantity, characteristic for a given ceramic material. Experiments to obtain these average residual stresses are expensive, so computer-aided methods are more and more frequently applied for design of new ceramic materials [3]. Mostly, computer calculations are based on the ®nite element method (FEM) [6]. Two-dimensional discrete models of polycrystalline microstructures are commonly used. This enables large FEM models to be built, and numerical simulations of microcrack propagation to be performed. In case of two-dimensional models it is possible to de®ne grain shapes on the basis of real material
350
T. Niezgoda / Probabilistic Engineering Mechanics 16 (2001) 349±354
microscopic images. Orthotropy orientations in model grains may be determined by experimental methods like piezospectroscopy [5]. In this paper, the author presents a numerical FEM-based method of average residual stress assessment in ceramic materials based on Monte Carlo simulation techniques [7]. The method is based on the following components: a threedimensional FEM model of a ceramic microstructure, a method of random material orientation assignment applying the Monte Carlo method, method of boundary condition de®nition and a statistical method of evaluation of the numerically generated results. Both ceramics and ceramic composites were taken into consideration. The ceramics were represented by corundum Al2O3 and the composite was represented by corundum with zirconia ZrO2. In the case of ceramic composites, additional factors affect the origination of residual stress state Ð the signi®cant difference in mechanical properties between components (corundum and zirconia, Young's modulus equals, respectively, 400 and 200 GPa for these materials) and the difference in thermal expansion coef®cients between the material components (for corundum and zirconia a equals, respectively, 8.9 and 10.9 £ 10 26 K 21).
Fig. 2. Model with polyhedrons as grains.
second model is based on the so-called `Kelvin polyhedrons' (Fig. 2). Here, each polyhedron models a single grain too, but it consists of 32 ®nite elements. Such models allowed us to represent ceramics only if they do not allow us to differentiate grain sizes. 2.2. Model with combined geometry
2. Microstructure modeling
The calculations were made for a prism that was assumed to be an internal part of an in®nite ceramic body. Three models of the prism have been investigated so far. The ®rst one consists of cubic grains (Fig. 1). Each cube, consisting of a single ®nite element, models a single grain. The
The model (the third one) used in the calculations presented here contains cubic as well as polyhedron grains (Fig. 2). It allows for changing of the proportion between sizes of cubic and polyhedral grains in the model (Fig. 3). Each cubic grain consists of eight ®nite elements, and the polyhedron grain consists of 32 ®nite elements. This model was used for calculations involving both ceramics and ceramic composites.
Fig. 1. Model with cubic grains (simplest one).
Fig. 3. Model of microstructure geometry consisting of cubic (in black) and polyhedron grains.
2.1. Models with homogenous geometry
T. Niezgoda / Probabilistic Engineering Mechanics 16 (2001) 349±354
351
Fig. 4. Currently generated grain and its neighbors. XYZ Ð global coordinate frame for FEM model; a, m, c Ð material coordinate frame.
3. Monte Carlo method application in material database de®nition Orthotropy axes in single monocrystals are named a, m, c (see Fig. 4). Thermal expansion coef®cients corresponding to these axes are named, respectively, a a, a m, a c and in case of orthotropy they complete a relation aa am , ac : As was mentioned earlier, combined random orientation of these axes in neighboring grains and temperature fall during manufacturing process is the cause of residual stress origination in polycrystalline materials like ceramics. The components of the residual stress are, respectively, s a, s m, s c. Because a a am , ac , tensile residual stresses occur in the volume. The average value of s c component is studied here. When modeling polycrystalline orthotropic materials the problem is to establish a method of assigning orientations of the materials property in each grain. It is enough to de®ne the c axis, the a and m axes can have any orientation in each particular grain. For the needs of these preliminary investigations and for demonstration purposes of the model, a relativity small number of grains is considered here. In such a case the neighborhood of grains with the same or close orientation of the c axis is not permitted because it may signi®cantly affect the calculation results. The orthotropy axes for realistic ceramic materials are in random order. It is possible that the adjacent grains have the same orthotropy directions. In this case the highest stress levels do not appear. Because the numerical model is limited by the number of degrees of freedom, the restrictions to the orthotropy directions were imposed, and the randomizing process applied. The author has considered pure ceramic material, and therefore, limited the number of axes to 14 in Kelvin's polyhedron. The numerical model representing more
Fig. 5. The set of permissible directions of c axes Ð the selection of the c axis direction for each particular grain de®nes the plane of orthogonal a and m axes.
complex microstructure (e.g. some kinds of hollows and insertions of other materials) will be considered in future. The model of microstructure is built automatically using automatic mesh generation software. Model grains are generated row by row, layer by layer. For each grain that is currently generated a c axis orientation is chosen at random from a set of pre-de®ned permissible directions. If the chosen direction is the same as the other in any of the neighboring grains generated earlier, the process must be repeated. It was noticed that a currently generated grain has a maximum of 13 neighbors, which were generated earlier (Fig. 4). Thus, the minimum number of permissible directions is 14. It is also a suf®cient number of possible directions of orthotropy enabling ef®cient diversity of orthotropy axes orientations in the model and should enable calculations of higher values of stresses. Permissible directions are arbitrarily pre-de®ned with the use of nodes of the ®nite element model mesh (Fig. 5). 3.1. Randomizing algorithm on the base of Monte Carlo method application It is known that the model contains m grains. An array d(m) is used as the material directions database. The array element dk contains the number associated with the direction that is attributed to the kth grain. Each grain has maximum n neighbors de®ned during the automatic generation process that is conducted in the way: row by row, plane by plane. Thus, n 13: An array g(m, x, n) is given, where its element gkl stores the number of lth neighbor of kth element. It is assumed that h is the number of arbitrary pre-de®ned
352
T. Niezgoda / Probabilistic Engineering Mechanics 16 (2001) 349±354
permissible directions. In the presented considerations h 14: An array c(h ) is given and contains numbers of these directions. Uniform number generation is used to ®nd an integer between 1 and 14, de®ning directions of orthotropy. The algorithm of direction randomizing all the m grains is as follows: Loop over m elements; do 30 k 1; m 10 continue Loop over n neighbors; do 20 l 1; n set i Ð the number of lth neighbor of the kth element; i g k; l set check parameter (1 is the truth); delta 1 choose at random direction number c(k) d(i) is the direction number of the lth neighbor of the kth grain if the current c(k) direction is the same as in the lth neighbor of kth if (c(k).eq.d(i)) delta 0 20 continue if no success Ð repeat drawing if (delta.eq. 0) go to 10 success: d(i) is the proper direction for the current grain 30 continue.
4. Boundary conditions The calculations were made for a prism that was assumed to be an internal part of an `in®nite' ceramic body. It can be any ceramic structure (in macroscopic scale) and the prism represents ®nite volume of the ceramic body in microscopic scale. The prism volume makes a small part of the in®nite body volume and shows the shapes of grains. The calculations were made with two-step method often called `global± local' method. It assumes that the ®rst calculating step is carried out for the ceramic body. In this case, it is assumed that the body is an isotropic material. It means that the same isotropic thermal expansion coef®cient a is de®ned for each axis. Because the large number of thermally orthotropic grains is situated statistically in volume body, it can be said that behavior of the ceramic body is isotropic. Then the displacements u, v, w of the prism boundaries are obtained for temperature fall DT. In the second step, calculations were performed for the prism, which consists of grains. Each grain was assumed to be orthotropic, with appropriate orthotropy of thermal expansion coef®cients a a am , ac and calculations were made based on the u, v, w displacements from the ®rst step and for the same thermal load DT (Fig. 6). As a result the residual stress state in the model is obtained.
Fig. 6. Method of the residual stress calculation.
5. Methods of results evaluation 5.1. Qualitative evaluation of results The contour plot for s c for alumina is presented (Fig. 7). They are obtained for the data in the local coordinate systems that are the same as the directions of local orthotropy in each of the grains. However, the uses of contour plots alone do not allow us to assess exactly the value of s c residual stress. In numerical calculation, the stress were performed in the integration points (called Gaussian points) and numerical model was built for solid ceramic body, the stress concentrations do not occur due to the geometry of squared corners. By using contour plots we can only determine that the most ef®cient contractions exist at grain boundaries, where orthotropy axes orientations must occur. In case of large disproportion in sizes of grains, it can also be seen that small grains surrounded by bigger ones are the stress concentrators in entire volume.
Fig. 7. An example of contour plot for s c component for the alumina ceramics.
T. Niezgoda / Probabilistic Engineering Mechanics 16 (2001) 349±354
353
Fig. 8. A solid ®nite element with nodes (in white) and integration points (in black).
As expected, for alumina in both fractions of the grains the average s c stress is a tensile one. In the case of ceramic composites based on alumina and zirconia, stress concentrations are not visible. Alumina grains are under tension and zirconia grains are under compression, for any vol% fractions of alumina. This is the effect of differences in thermal expansion coef®cients Da 0 2:0 £ 10 26 K21 and, in most, of high difference in mechanical properties (Young's modulus for alumina equals < 400 GPa and for zirconia < 200 GPa). 5.2. Quantitative evaluation of results For more precise assessment of the stress state in the model a quantitative evaluation of results must be performed. The method of quantitative evaluation of numerical calculation results is analogous to the method of evaluation of experimental data. The proposed method of evaluation of calculation results makes them comparable to the experimental ones. Stresses at all integration points (Fig. 8) are recorded, and a histogram is established as in Fig. 9. As a result it is possible to obtain the distribution of a
Fig. 10. Peak position is interpreted as searched average value of the given stress components s a, s m and s c (distributions obtained for alumina).
given stress component for any part of the model, i.e. for the fraction of cubic or polyhedron grains as well as for the entire model. This data is extremely useful, as it would be dif®cult and very expensive to generate them experimentally. The shapes of the obtained distributions of numerical data are then approximated, assuming a Gaussian distribution (Fig. 10). Taking into consideration the average values of s a, s m and s c obtained in such a manner for the entire model, it is possible to determine that the equilibrium condition given by
sa 1 sm 1 sc 0
2
is completed for this numerical model, providing justi®cation for the proposed methods (Fig. 10). 6. Conclusions Fig. 9. Distribution of the stress component s c for the entire model. It was obtained by making counts in the discrete intervals over integration points of ®nite elements.
The proposed type of modeling of ceramic microstructure enables investigations of relative grain size in¯uence on stress distribution in ceramic composites. The greater the
354
T. Niezgoda / Probabilistic Engineering Mechanics 16 (2001) 349±354
difference in size between neighboring grains, the higher the resulting stress in the small grains. The model presented here may be used not only for investigations of pure ceramics, but also for investigations of two-phase ceramic composites [6]. Using the statistic post-processing of numerical results, the quantitative assessment of stresses in the modeled microstructure is possible. The stress values calculated in integration points of ®nite elements are counted, to obtain stress distribution in the modeled volume of material. Then the peak position of the distribution is interpreted as the average stress value. This makes the numerical investigations similar to an experimental technique used for measuring average residual stress (the piezospectroscopic technique was enabled to verify numerical results). Such a method of data evaluation brings more information about stress state in the model and makes the numerical calculations results comparable with the results of the experiment. They agree with the results of our own piezospectroscopic results and those described in the works of other authors [1,2,4]. The results of calculations presented in the form of contour plots are useful for qualitative assessment of calculated stresses. They show that the signi®cant part of internal energy in the modeled volume of microstructure is concentrated on the area close to the grain boundaries. The advantage of numerical simulations is that the stress distribution may be easily accounted in any part of the model. The presented models of microstructure geometry, the
method of material database de®nition and optimization, the method of boundary conditions de®nition and the method of numerical results evaluation are the basic analytical apparatus for analyzing residual stress state in ceramic materials. The obtained values of average residual stress for ceramics complete appropriate equilibrium condition Eq. (2). This is an indicator of correctness of the models and methods proposed in the research stage. The next step is to improve on these models in order to be able to apply them to realistic materials. References [1] Blendel JE, Coble RL. Measurement of stress due to thermal expansion anisotropy in Al2O3. J Am Ceram Soc 1982;65(3):174. [2] Krell A, Teresiak A, Schlafer D. Grain size dependent residual microstress in submicron Al2O3 and ZrO2. J Eur Ceram Soc 1996;16:803. [3] Niezgoda T, Szymczyk W, Mal¤achowski J. Numerical assessment of average residual stress in orthotropic alumina±zirconia ceramic composite. Sixth Conference and Exhibition of the European Ceramic Society, 1999. [4] Ma Q, Clarke DR. Stress measurement in single-crystal and polycrystalline ceramics using their optical ¯uorescence. J Am Ceram Soc 1993;76:1433±40. [5] Ma Q, Clarke DR. Piezospectroscopic determination of residual stresses in polycrystalline alumina. J Am Ceram Soc 1994;77(2):298. [6] Dacko M, Borkowski W, DobrocinÂski S, Niezgoda T, Wieczorek M. Metoda elementoÂw skonÂczonych w mechanice konstrukcji, ARKADY (in Polish), 1994. [7] Szmelter J. Metody Komputerowe w mechanice, PWN, Warsaw, 1980. (In Polish)
www.elsevier.com/locate/probengmech
Material database optimization using Monte Carlo method in engineering application Tadeusz Niezgoda* Faculty of Mechanics, Military University of Technology, Kaliskiego Street 2, 00-908 Warsaw, Poland
Abstract Ceramic materials were simulated numerically. The practical goal of the research was to assess numerically the average residual stress value in such materials. The scienti®c goal of the work was to establish a basic analytical apparatus. One of the main aspects was to create material database, enabling simulation of polycrystalline structure with random orientations of thermal expansion orthotropy axes. Monte Carlo method was adopted for this purpose. The calculation results were analyzed with the use of statistical method, which is analogous to experimental measure analysis. However, the developed methods can be readily used to analyze both ceramics and ceramic composites. Currently, the methods can take advantage of more precise geometrical description to achieve more accurate models re¯ecting real internal structure of the material. q 2001 Published by Elsevier Science Ltd. Keywords: Finite element method; Monte Carlo method; Al2O3 ceramics; Numerical method; Residual stress
1. Introduction Ceramics and ceramic composites are materials with very interesting properties. They are often used in structures that are subject to high temperatures, friction and erosion processes caused by chemically aggressive environments or high-speed ¯ow of gases. For example, thin ceramic coatings protect elements of turbines, jet and internal combustion engines from excessive heat induced loading, and ceramic materials are often used as cutting tools. Porous ceramic materials may be used as membranes for gas ®ltration and separation. Both ceramics and ceramic composites are polycrystalline materials. Each grain in such a structure is thermally orthotropic. In cooling down polycrystalline ceramic materials, from sintering to room temperature, the random variation in crystallographic orientation of orthotropy from one grain to another induces a very complex state of constraint. Thus, after completion of the manufacturing process, a state of residual stress remains. This limits the range of acceptable loads on these ceramic materials. The complicated spatial variation of the residual stress ®eld depends on the elastic properties of grains, the grain size and the shape distribution. In some cases, the residual stress across the grain
* Tel./fax: 148-22-683-94-61. E-mail address: [email protected] (T. Niezgoda). 0266-8920/01/$ - see front matter q 2001 Published by Elsevier Science Ltd. PII: S 0266-892 0(01)00015-7
boundary can be suf®ciently large to initialize microcracking across the grain boundary. In elastic solids this residual stress s T is given by [1]
s T bEDaDT
1
where b is a coef®cient that depends on the shape of the grain and the grain orientation b [ k 2 0:5; 10:5l; E the Young's modulus, Da the difference in the singlecrystal thermal expansion coef®cients along the direction considered, DT the difference between the temperature Tc (at which the sample was stress-free) and the room temperature T0. The orthotropy of thermal expansion coef®cients is caused by building of crystallographic lattice of Al2O3 ceramic grain. For orthotropic grains aa am ± ac and hence, Da ac 2 aa for Al2O3 ceramics equals Da 0:9 £ 1026 K21 : The average residual stress value is a deterministic quantity, characteristic for a given ceramic material. Experiments to obtain these average residual stresses are expensive, so computer-aided methods are more and more frequently applied for design of new ceramic materials [3]. Mostly, computer calculations are based on the ®nite element method (FEM) [6]. Two-dimensional discrete models of polycrystalline microstructures are commonly used. This enables large FEM models to be built, and numerical simulations of microcrack propagation to be performed. In case of two-dimensional models it is possible to de®ne grain shapes on the basis of real material
350
T. Niezgoda / Probabilistic Engineering Mechanics 16 (2001) 349±354
microscopic images. Orthotropy orientations in model grains may be determined by experimental methods like piezospectroscopy [5]. In this paper, the author presents a numerical FEM-based method of average residual stress assessment in ceramic materials based on Monte Carlo simulation techniques [7]. The method is based on the following components: a threedimensional FEM model of a ceramic microstructure, a method of random material orientation assignment applying the Monte Carlo method, method of boundary condition de®nition and a statistical method of evaluation of the numerically generated results. Both ceramics and ceramic composites were taken into consideration. The ceramics were represented by corundum Al2O3 and the composite was represented by corundum with zirconia ZrO2. In the case of ceramic composites, additional factors affect the origination of residual stress state Ð the signi®cant difference in mechanical properties between components (corundum and zirconia, Young's modulus equals, respectively, 400 and 200 GPa for these materials) and the difference in thermal expansion coef®cients between the material components (for corundum and zirconia a equals, respectively, 8.9 and 10.9 £ 10 26 K 21).
Fig. 2. Model with polyhedrons as grains.
second model is based on the so-called `Kelvin polyhedrons' (Fig. 2). Here, each polyhedron models a single grain too, but it consists of 32 ®nite elements. Such models allowed us to represent ceramics only if they do not allow us to differentiate grain sizes. 2.2. Model with combined geometry
2. Microstructure modeling
The calculations were made for a prism that was assumed to be an internal part of an in®nite ceramic body. Three models of the prism have been investigated so far. The ®rst one consists of cubic grains (Fig. 1). Each cube, consisting of a single ®nite element, models a single grain. The
The model (the third one) used in the calculations presented here contains cubic as well as polyhedron grains (Fig. 2). It allows for changing of the proportion between sizes of cubic and polyhedral grains in the model (Fig. 3). Each cubic grain consists of eight ®nite elements, and the polyhedron grain consists of 32 ®nite elements. This model was used for calculations involving both ceramics and ceramic composites.
Fig. 1. Model with cubic grains (simplest one).
Fig. 3. Model of microstructure geometry consisting of cubic (in black) and polyhedron grains.
2.1. Models with homogenous geometry
T. Niezgoda / Probabilistic Engineering Mechanics 16 (2001) 349±354
351
Fig. 4. Currently generated grain and its neighbors. XYZ Ð global coordinate frame for FEM model; a, m, c Ð material coordinate frame.
3. Monte Carlo method application in material database de®nition Orthotropy axes in single monocrystals are named a, m, c (see Fig. 4). Thermal expansion coef®cients corresponding to these axes are named, respectively, a a, a m, a c and in case of orthotropy they complete a relation aa am , ac : As was mentioned earlier, combined random orientation of these axes in neighboring grains and temperature fall during manufacturing process is the cause of residual stress origination in polycrystalline materials like ceramics. The components of the residual stress are, respectively, s a, s m, s c. Because a a am , ac , tensile residual stresses occur in the volume. The average value of s c component is studied here. When modeling polycrystalline orthotropic materials the problem is to establish a method of assigning orientations of the materials property in each grain. It is enough to de®ne the c axis, the a and m axes can have any orientation in each particular grain. For the needs of these preliminary investigations and for demonstration purposes of the model, a relativity small number of grains is considered here. In such a case the neighborhood of grains with the same or close orientation of the c axis is not permitted because it may signi®cantly affect the calculation results. The orthotropy axes for realistic ceramic materials are in random order. It is possible that the adjacent grains have the same orthotropy directions. In this case the highest stress levels do not appear. Because the numerical model is limited by the number of degrees of freedom, the restrictions to the orthotropy directions were imposed, and the randomizing process applied. The author has considered pure ceramic material, and therefore, limited the number of axes to 14 in Kelvin's polyhedron. The numerical model representing more
Fig. 5. The set of permissible directions of c axes Ð the selection of the c axis direction for each particular grain de®nes the plane of orthogonal a and m axes.
complex microstructure (e.g. some kinds of hollows and insertions of other materials) will be considered in future. The model of microstructure is built automatically using automatic mesh generation software. Model grains are generated row by row, layer by layer. For each grain that is currently generated a c axis orientation is chosen at random from a set of pre-de®ned permissible directions. If the chosen direction is the same as the other in any of the neighboring grains generated earlier, the process must be repeated. It was noticed that a currently generated grain has a maximum of 13 neighbors, which were generated earlier (Fig. 4). Thus, the minimum number of permissible directions is 14. It is also a suf®cient number of possible directions of orthotropy enabling ef®cient diversity of orthotropy axes orientations in the model and should enable calculations of higher values of stresses. Permissible directions are arbitrarily pre-de®ned with the use of nodes of the ®nite element model mesh (Fig. 5). 3.1. Randomizing algorithm on the base of Monte Carlo method application It is known that the model contains m grains. An array d(m) is used as the material directions database. The array element dk contains the number associated with the direction that is attributed to the kth grain. Each grain has maximum n neighbors de®ned during the automatic generation process that is conducted in the way: row by row, plane by plane. Thus, n 13: An array g(m, x, n) is given, where its element gkl stores the number of lth neighbor of kth element. It is assumed that h is the number of arbitrary pre-de®ned
352
T. Niezgoda / Probabilistic Engineering Mechanics 16 (2001) 349±354
permissible directions. In the presented considerations h 14: An array c(h ) is given and contains numbers of these directions. Uniform number generation is used to ®nd an integer between 1 and 14, de®ning directions of orthotropy. The algorithm of direction randomizing all the m grains is as follows: Loop over m elements; do 30 k 1; m 10 continue Loop over n neighbors; do 20 l 1; n set i Ð the number of lth neighbor of the kth element; i g k; l set check parameter (1 is the truth); delta 1 choose at random direction number c(k) d(i) is the direction number of the lth neighbor of the kth grain if the current c(k) direction is the same as in the lth neighbor of kth if (c(k).eq.d(i)) delta 0 20 continue if no success Ð repeat drawing if (delta.eq. 0) go to 10 success: d(i) is the proper direction for the current grain 30 continue.
4. Boundary conditions The calculations were made for a prism that was assumed to be an internal part of an `in®nite' ceramic body. It can be any ceramic structure (in macroscopic scale) and the prism represents ®nite volume of the ceramic body in microscopic scale. The prism volume makes a small part of the in®nite body volume and shows the shapes of grains. The calculations were made with two-step method often called `global± local' method. It assumes that the ®rst calculating step is carried out for the ceramic body. In this case, it is assumed that the body is an isotropic material. It means that the same isotropic thermal expansion coef®cient a is de®ned for each axis. Because the large number of thermally orthotropic grains is situated statistically in volume body, it can be said that behavior of the ceramic body is isotropic. Then the displacements u, v, w of the prism boundaries are obtained for temperature fall DT. In the second step, calculations were performed for the prism, which consists of grains. Each grain was assumed to be orthotropic, with appropriate orthotropy of thermal expansion coef®cients a a am , ac and calculations were made based on the u, v, w displacements from the ®rst step and for the same thermal load DT (Fig. 6). As a result the residual stress state in the model is obtained.
Fig. 6. Method of the residual stress calculation.
5. Methods of results evaluation 5.1. Qualitative evaluation of results The contour plot for s c for alumina is presented (Fig. 7). They are obtained for the data in the local coordinate systems that are the same as the directions of local orthotropy in each of the grains. However, the uses of contour plots alone do not allow us to assess exactly the value of s c residual stress. In numerical calculation, the stress were performed in the integration points (called Gaussian points) and numerical model was built for solid ceramic body, the stress concentrations do not occur due to the geometry of squared corners. By using contour plots we can only determine that the most ef®cient contractions exist at grain boundaries, where orthotropy axes orientations must occur. In case of large disproportion in sizes of grains, it can also be seen that small grains surrounded by bigger ones are the stress concentrators in entire volume.
Fig. 7. An example of contour plot for s c component for the alumina ceramics.
T. Niezgoda / Probabilistic Engineering Mechanics 16 (2001) 349±354
353
Fig. 8. A solid ®nite element with nodes (in white) and integration points (in black).
As expected, for alumina in both fractions of the grains the average s c stress is a tensile one. In the case of ceramic composites based on alumina and zirconia, stress concentrations are not visible. Alumina grains are under tension and zirconia grains are under compression, for any vol% fractions of alumina. This is the effect of differences in thermal expansion coef®cients Da 0 2:0 £ 10 26 K21 and, in most, of high difference in mechanical properties (Young's modulus for alumina equals < 400 GPa and for zirconia < 200 GPa). 5.2. Quantitative evaluation of results For more precise assessment of the stress state in the model a quantitative evaluation of results must be performed. The method of quantitative evaluation of numerical calculation results is analogous to the method of evaluation of experimental data. The proposed method of evaluation of calculation results makes them comparable to the experimental ones. Stresses at all integration points (Fig. 8) are recorded, and a histogram is established as in Fig. 9. As a result it is possible to obtain the distribution of a
Fig. 10. Peak position is interpreted as searched average value of the given stress components s a, s m and s c (distributions obtained for alumina).
given stress component for any part of the model, i.e. for the fraction of cubic or polyhedron grains as well as for the entire model. This data is extremely useful, as it would be dif®cult and very expensive to generate them experimentally. The shapes of the obtained distributions of numerical data are then approximated, assuming a Gaussian distribution (Fig. 10). Taking into consideration the average values of s a, s m and s c obtained in such a manner for the entire model, it is possible to determine that the equilibrium condition given by
sa 1 sm 1 sc 0
2
is completed for this numerical model, providing justi®cation for the proposed methods (Fig. 10). 6. Conclusions Fig. 9. Distribution of the stress component s c for the entire model. It was obtained by making counts in the discrete intervals over integration points of ®nite elements.
The proposed type of modeling of ceramic microstructure enables investigations of relative grain size in¯uence on stress distribution in ceramic composites. The greater the
354
T. Niezgoda / Probabilistic Engineering Mechanics 16 (2001) 349±354
difference in size between neighboring grains, the higher the resulting stress in the small grains. The model presented here may be used not only for investigations of pure ceramics, but also for investigations of two-phase ceramic composites [6]. Using the statistic post-processing of numerical results, the quantitative assessment of stresses in the modeled microstructure is possible. The stress values calculated in integration points of ®nite elements are counted, to obtain stress distribution in the modeled volume of material. Then the peak position of the distribution is interpreted as the average stress value. This makes the numerical investigations similar to an experimental technique used for measuring average residual stress (the piezospectroscopic technique was enabled to verify numerical results). Such a method of data evaluation brings more information about stress state in the model and makes the numerical calculations results comparable with the results of the experiment. They agree with the results of our own piezospectroscopic results and those described in the works of other authors [1,2,4]. The results of calculations presented in the form of contour plots are useful for qualitative assessment of calculated stresses. They show that the signi®cant part of internal energy in the modeled volume of microstructure is concentrated on the area close to the grain boundaries. The advantage of numerical simulations is that the stress distribution may be easily accounted in any part of the model. The presented models of microstructure geometry, the
method of material database de®nition and optimization, the method of boundary conditions de®nition and the method of numerical results evaluation are the basic analytical apparatus for analyzing residual stress state in ceramic materials. The obtained values of average residual stress for ceramics complete appropriate equilibrium condition Eq. (2). This is an indicator of correctness of the models and methods proposed in the research stage. The next step is to improve on these models in order to be able to apply them to realistic materials. References [1] Blendel JE, Coble RL. Measurement of stress due to thermal expansion anisotropy in Al2O3. J Am Ceram Soc 1982;65(3):174. [2] Krell A, Teresiak A, Schlafer D. Grain size dependent residual microstress in submicron Al2O3 and ZrO2. J Eur Ceram Soc 1996;16:803. [3] Niezgoda T, Szymczyk W, Mal¤achowski J. Numerical assessment of average residual stress in orthotropic alumina±zirconia ceramic composite. Sixth Conference and Exhibition of the European Ceramic Society, 1999. [4] Ma Q, Clarke DR. Stress measurement in single-crystal and polycrystalline ceramics using their optical ¯uorescence. J Am Ceram Soc 1993;76:1433±40. [5] Ma Q, Clarke DR. Piezospectroscopic determination of residual stresses in polycrystalline alumina. J Am Ceram Soc 1994;77(2):298. [6] Dacko M, Borkowski W, DobrocinÂski S, Niezgoda T, Wieczorek M. Metoda elementoÂw skonÂczonych w mechanice konstrukcji, ARKADY (in Polish), 1994. [7] Szmelter J. Metody Komputerowe w mechanice, PWN, Warsaw, 1980. (In Polish)