M. AMBRO@I^, A. NIKONOV: DYNAMIC BIAXIAL STRESS ANALYSIS OF FLAT LAYERED CERAMIC COMPOSITES 195–200 DYNAMIC BIAXIAL STRESS ANALYSIS OF FLAT LAYERED CERAMIC COMPOSITES ANALIZA DVOOSNIH DINAMI^NIH NAPETOSTI V RAVNIH PLASTOVITIH KERAMI^NIH KOMPOZITIH Milan Ambro`i~ 1,2* , Anatolij Nikonov 2 1 Faculty of Mathematics and Natural Sciences, University of Maribor, Koro{ka 160, 2000 Maribor, Slovenia 2 Faculty of Industrial Engineering, [egova 112, 8000 Novo mesto, Slovenia Prejem rokopisa – received: 2020-04-06; sprejem za objavo – accepted for publication: 2020-11-24 doi:10.17222/mit.2020.055 We study theoretically the biaxial bending of symmetric, flat layered ceramic composites (laminates) due to external loading. We focus on three-layered alumina/zirconia laminates. We compare the principal stresses in the samples in the case of static and harmonic dynamic loading. The dynamic equation within the Kirchhoff theory for thin homogeneous plates is first generalized to the case of multilayered plates. It is solved numerically with the relaxation method, which we have developed for this pur- pose. Keywords: layered ceramic composites, Kirchhoff theory of plates, biaxial stress, principal stresses Teoreti~no smo {tudirali dvoosni upogib simetri~nih ravnih plastovitih kerami~nih kompozitov (laminatov) zaradi zunanje obremenitve. Osredoto~ili smo se na triplastne laminate iz aluminijevega in cirkonijevega oksida. Primerjali smo lastne vrednosti napetosti za stati~no in harmoni~no dinami~no obremenitev. Dinami~no ena~bo v okviru Kirchhoffove teorije za homogene tanke plo{~e smo najprej generalizirali za primer ve~plastnih plo{~. Ena~bo smo re{evali numeri~no z relaksacijsko metodo, ki smo jo razvili v ta namen. Klju~ne besede: plastoviti kerami~ni kompoziti; Kirchhoffova teorija plo{~; dvoosna napetost; lastne vrednosti napetosti 1 INTRODUCTION Alumina (Al 2 O 3 ) based ceramics are frequently used in various applications, such as cutting tools and biomed- ical implants, because of their good mechanical proper- ties. Pure alumina has high hardness and low weight, but a relative moderate bend strength and fracture toughness. Zirconia (ZrO 2 ) can have significantly higher fracture toughness and bend strength than alumina. Therefore, composites of alumina with an appropriate volume frac- tion of tetragonal zirconia (zirconia toughened alumina, ZTA) are promising materials due to a high bend strength and fracture toughness, as well as wear resis- tance. 1–6 F. Sommer et al. 6 have reported that using 1 mol. % of yttria can raise the bend strength of ZTA ce- ramics (containing 17 vol. % of yttria stabilized zirconia) up to nearly 1200 MPa and the ISB fracture toughness up to 8.5 MPa m. 6 Multi-layered alumina/zirconia com- posites (laminates) have been studied particularly in re- gard to residual and loading stresses. 7–12 Thermal resid- ual stresses within individual layers arise upon cooling the material after the sintering process because of the mismatch of the thermal expansion coefficients of alu- mina and zirconia. Various static uniaxial or biaxial bending tests are used to measure the strength of ceramic samples. E. Carrera and A. Ciuffreda 13 compared a three-dimensional stress distribution in composite plates from different the- ories and validated their results for different static load- ings with the finite-element method (FEM). However, the ceramic and other engineering products/components are often subjected to dynamic loads in everyday use. Therefore, different aspects of the dynamics of elastic plates, shells and beams have been investigated thor- oughly. A particular interest has been devoted to the propagation of traction-free elastic waves in thin flat plates, theoretically with a finite thickness but infinite lateral dimensions. An analytical approach shows that even for homogeneous and symmetric three-layered plates the ordinary equations for dispersion relations and the corresponding strains and stresses are complex. 14–16 Different analytical or semi-analytical approaches for fi- nite plates can be used only when the problem is effec- tively one-dimensional (1D) and the boundary conditions are simple. Z. J. Ai et al. 17 studied the problem of an elastic plate on the stack of layers with the uniformly moving load on the plate, where they used the Fourier transformation and analytical layer method. An interest- ing idea for the partially analytical approach for the two-dimensional (2D) problem of multilayered plate dy- namics was given by M. Sharyjat and M. Roshanfar. 18 They expanded the deflection of the plate with sine func- tions of coordinates and with time-dependent coeffi- cients. Materiali in tehnologije / Materials and technology 55 (2021) 2, 195–200 195 UDK 521.19:666.3.015:543.58:532.6 ISSN 1580-2949 Original scientific article/Izvirni znanstveni ~lanek MTAEC9, 55(2)195(2021) *Corresponding author's e-mail: milan.ambrozic@um.si (Milan Ambro`i~) Numerical calculations for the biaxial deformations of finite plates include several variants of FEM methods. Zhang et al. have developed an efficient two-step FEM method for the calculation of the dynamics and stresses in multilayered plates: the method is composed of a fine-scale computation for parts of the system, followed by a coarse-scale computation for the whole system. 19 Among various numerical methods for the three-dimen- sional (3D) problem of multilayered plates we mention a very sophisticated sampling surface method, where the basic variables are a dense set of inner sampling sur- faces. T. Ye et al. 20 implemented this method with the spectral expansion technique where the basic orthogonal functions are Chebyshev polynomials. Using this ap- proach, they were able to reach a highly accurate numer- ical precision, even for relatively thick plates with vari- ous boundary conditions. There are several other approaches and techniques to attack the problem of mul- tilayered plates, but they are not focused directly on the simple single equation for the dynamic bending of the plate as a whole. 21–24 As regards the parabolic-type fourth-order classic equation for bending deformation of thin plates within the Kirchhoff theory, several attempts have been made to refine the theory, particularly for infi- nite plates, including Timoshenko-Reissner-Mindlin type theories. 25 In this paper we use the simple Kirchhoff theory for pure 2D bending deformation and generalize the dy- namic equation of a homogeneous plate to the case of multilayered plates. The dynamic equation is solved nu- merically on a discrete network with a relaxation method. We are particularly interested how the frequency of the external load and the composition of the sample influence the stress magnitude in the 3-layered alu- mina/zirconia sample. A benefit of our approach is that we present the vibration of the plate as a single object. The parameters of different layers enter a single dynamic equation only through the composed flexural rigidity, without the necessity of solving the coupled system of equations for different layers simultaneously. After the dynamic equation is solved, the stress tensor is calcu- lated in each layer separately. It may be interesting from the mathematical point, that the efficiency of the relax- ation method for the fourth-order equation is not signifi- cantly different from that for elliptic-type equations of the second order. 2 NUMERICAL MODEL Let us first consider the case of a homogeneous, thin square plate subject to bending deformation. Its Young’s modulus, Poisson’s ratio and mass density are E, and , respectively. The thickness of the plate is h, while its side is 2a. We choose the Cartesian coordinate system with the origin at the center of the plate, the z–axis is perpen- dicular to the plate, and the other two axes parallel to the sides of the plate. The neutral plane is at z = 0. The de- formation according to the Kirchhoff theory for a thin plate is characterized by the displacement w(x, y, t)ofthe points on the neutral plane in the direction of z–axis. The corresponding dynamic equation for the "displacement function" w is: Dwh w t p () ∇+= 22 2 2 ∂ ∂ (1) where 2 = 2 / x 2 + 2 / y 2 is the 2D Laplacian. 26,27 The quantity p on the right side of Equation (1) is the exter- nal pressure difference on both sides of the plate. In general, it can be a function of x, y and t. The flexural ri- gidity (bending stiffness) for a homogeneous plate is: D Eh = − 3 2 12 1() (2) A convenient numerical form of Equation (1) and its extension to a layered plate in dimensionless units (de- noted by * ) are given in Appendix A. We take some refer- ence Young’s modulus E 0 and density 0 ,s ot h a t E * = E/E 0 , p * = p/E 0 and * = / 0 . We introduce the ref- erence flexural rigidity as D 0 = E 0 h 3 instead of Equation (2), and the dimensionless flexural rigidity is D * = D/D 0 . The coordinates and displacement are normalized as: x * = x/a, y * = y/a, z * = z/h and w * = w/h. We define the ra- tio = h/a, which is supposed to be small. A suitable choice for the dimensionless time is t * = t/t 0 , where the characteristic time is: t E a hE a 0 0 0 2 0 0 =⋅=⋅ (3) This expression gives a correct order of magnitude of the lowest eigenfrequency =1 / t 0 . When the function w * is found, the components of the strain and stress tensors can be calculated; see the explanation in the textbook of by Landau and Lifshitz 26 and additional details in the pa- pers of V . V . Vasiliev 27 and N. I. Robinson. 28 The basic assumption in the derivation of equations for thin plates is that the stress components xz , yz and zz can be ne- glected in comparison with xx , yy and xy . The eigen- values 1 and 2 (principal stresses) of the remaining 2D stress tensor can finally be calculated. The direct rela- tionship between the stress tensor and the displacement function w * is given in Appendix A. We present in this paper the results for homogeneous pressure over the plate. For static loading the constant pressure p * 0 is applied. In the case of harmonic dynamic loading we take: p * = p * 0 sin(2 app t), with applied fre- quency app and the corresponding period T app =1 / app . We present the results for clamped edges (rigid boundary with zero values of displacement and its normal deriva- tive): (I) w = 0 and w/ x = 0 for x = a; (II) w = 0 and w/ y = 0 for y = a. M. AMBRO@I^, A. NIKONOV: DYNAMIC BIAXIAL STRESS ANALYSIS OF FLAT LAYERED CERAMIC COMPOSITES 196 Materiali in tehnologije / Materials and technology 55 (2021) 2, 195–200 3 RESULTS AND DISCUSSION 3.1 Analysis of a homogeneous plate Typical symmetrically positioned representative points of the plate at z = 0 are indicated in Figure 1. First, we consider a homogeneous "plate". We choose a homogeneous alumina sample with the following pa- rameters (Table 1): 2a = 100 mm, h =5m m( = 0.1), E = E 0 = 390 GPa, = 0.238, = 0 = 3.98 kg/dm 3 , p 0 = 1 MPa. The time dependence of the displacement of three representative points in the harmonic loading is compared in Figure 2, where we choose T app * = 10. The period of the basic eigenmode is estimated: T * = 2.32 or T A = T * t 0 = 117 s (Equation (3)); the corre- sponding eigenfrequency is A =1 / T A = 8.53 kHz (the in- dex "A" stands for alumina). To represent the distribution of stress we take a static loading with pressure p 0 =1M P a .Figure 3 shows the profile of both principal stresses along two symmetry lines: y = 0 (horizontal line through points N, P 3 ,P 0 ,P 1 and L); 2) and y = x (diagonal line through points A, P 7 , P 0 ,P 5 and C). The principal stresses 1 and 2 at the cen- tral point P 0 are equal due to symmetry (x * =0i nFig- ure 3). The lines on the top plate’s surface are taken (z = h/2) in Figure 3. We are interested in positive (tensile) stresses, which are much more detrimental for ceramic materials than compressive stresses. However, the strongly negative value 2 –120 MPa at the edges means that the corresponding principal stress is +120 MPa on the opposite surface (z =– h/2), at the point with the same coordinates x * and y * . The largest stresses thus correspond to edge points K, L, M and N in Figure 1, not to the central point P 0 . The static displacement of P 0 is w * stat = w * (0, 0) = 0.0057, or 0.03 mm in physical units. 3.2 Analysis of 3-layered plate Now, we focus on the ZTA composites and we keep the same dimensions as given above. We take a symmet- ric 3-layered composite with alumina outer layers and a ZTA inner layer. The corresponding material parameters, presented in Table 1, can be obtained from several refer- ences. 10,11 Only the value 500 MPa for the strength of pure alumina has been chosen rather arbitrarily. This is, because the data from literature vary significantly due to different fabrication procedures. 1,10,11 Table 1: Material parameters of pure alumina (A) and a particular ZTA composite: Young’s modulus, Poisson’s number, density and ap- proximate bend strength Material E (GPa) (g/cm 3 ) b (MPa) A 390 0.238 3.98 500 ZTA24 338 0.264 4.48 850 M. AMBRO@I^, A. NIKONOV: DYNAMIC BIAXIAL STRESS ANALYSIS OF FLAT LAYERED CERAMIC COMPOSITES Materiali in tehnologije / Materials and technology 55 (2021) 2, 195–200 197 Figure 3: Line profile of both principal stresses: dependence on x * = x/a. Here, 1 > 2 is chosen, regardless of the corresponding eigenvectors. The symbol (x) attached to graphs denotes the y = 0 line, and the symbol (d) the y = x diagonal line. Solid and dashed curves: 1 and 2 on y = 0 line, respectively; dash-dotted and dotted lines: 1 and 2 on y = x line, respectively Figure 1: Location of representative points on the square plate Figure 2: Time-dependent displacements of the plate at points P 0 (solid line), P 1 (dashed line) and P 5 (dash–dotted line) located at z =0 plane for harmonic pressure variation; T app * = T app /t 0 =1 0 We have varied systematically the thickness of the middle layer (keeping the total thickness as 5 mm) as well as the frequency of the pressure oscillation. Besides the pure alumina plate, we choose a ZTA24 composite (with 24 % of mass fractions of yttria stabilized zirconia) for the middle layer. 6 The lowest eigenfrequency A = 8.53 kHz for the alumina sample is the basis for the set of testing applied frequencies for all the samples. The maximum value of the larger value of the pair ( 1 , 2 ) over testing time and over the whole sample is traced. We call this "maximum stress" for briefness and denote it with max . The simulation time is several periods T app . The graphs for three compositions are presented in Fig- ure 4. It has been verified in all cases that the maximum stress appears at the edge points K, L, M and N in Fig- ure 1, both for static and dynamic pressure. The maxi- mum stress at the point P 0 is roughly twice smaller. The maximum stress increases steadily with the frequency app . It also slightly increases with the thickness of the ZTA middle layer. The exception is the largest frequency app = A , where the maximum stress of pure alumina is the largest. This is easily explained, since the calcula- tions show that when the thickness of the ZTA middle layer increases, the corresponding period T* for the first eigenmode increases from 2.32 for alumina to 2.61 for composite with 4.8 mm ZTA inner layer. Thus, for app = A the two comparable ZTA composites are slightly out of their own eigenfrequencies. Finally, we should stress another point. No damping of the plate oscillation has been considered in the pre- sented results for two reasons: (1) because of the lack of corresponding data for ceramics, (2) because the inclu- sion of the contribution of the individual layers to total damping is not straightforward. We have nevertheless tried to make a simple modification of the dynamic equa- tion in order to get some insight how the damping could influence the stresses in the material. Damping can be in- troduced by an additional term * w * / t * on the left side of the dimensionless Equation (A.1). The dimensionless parameter * represents the strength of damping, and the two terms with time derivatives on the left side of the up- graded equations compete. To test the case of moderate damping, we have chosen a particular value * = * /10, just for the pure alumina sample. As expected, the maxi- mum stress in the sample is decreased when the damping is present. For frequency app = A /8, the maximum stress i s2%l o w e rthan without damping. Even at app = A /2 the difference in maximum stress is only 4 %, while for app = A this difference becomes more significant: 20 %. We note that our criterion for the maximum stress is not the steady state, when the vibration eigenmodes cease, but we have checked the stresses from the beginning of vibration. The stresses in the steady state instead of the transient part of vibration can be simply traced, after many periods T app are elapsed. 4 CONCLUSIONS We have found that by increasing the frequency of the applied harmonic pressure from zero (static load) up to the “resonant” frequency the stresses in the ceramic composite laminate are increased by an order of magni- tude. For the frequencies about 1/8 of eigenfrequency, i.e., of the order kHz in our case, the dynamic stresses in the composite differ from the corresponding static values only by 10 %. Furthermore, while for static loading the differences in maximum stresses are relatively small when the thickness of the middle ZTA layer is increased, these differences are larger for dynamic loading. When a moderate damping is included, the results do not seem to change significantly up to half the resonant frequency, as simple calculations indicate. Some of our findings may serve researchers in the ce- ramic society, which work with composite materials. One of the goals is to optimize the composition of differ- ent layers to obtain the maximum strength of the lami- nates. 9–11 Although the loading in the uniaxial or biaxial tests is usually static and similar to point-like, the de- pendence of the stresses on frequency may be qualita- tively similar to the one presented in Figure 4. For in- stance, up to frequencies one order of magnitude smaller than the lowest eigenfrequency, the increment of stress can be safely ignored. For larger frequencies appropriate correction factors can be used: for instance, when the ap- plied frequency is half the lowest eigenfrequency, in our calculations the maximum stress has roughly twice the value corresponding to static loading. The order of mag- nitude of resonant frequency may be estimated from the characteristic time in Equation (3). We have also touched on some other aspects, e.g., boundary conditions, distribution of pressure and the elongation of the square into rectangle. If the boundary condition with the clamped edges is replaced by the con- M. AMBRO@I^, A. NIKONOV: DYNAMIC BIAXIAL STRESS ANALYSIS OF FLAT LAYERED CERAMIC COMPOSITES 198 Materiali in tehnologije / Materials and technology 55 (2021) 2, 195–200 Figure 4: Dependence of maximum stress on the applied frequency; three samples are presented, all with thickness 5 mm: (1) pure alumina (squares); (2) composite with 4 mm ZTA inner layer and two 0.5 mm alumina outer layers (triangles); (3) composite with 4.8 mm ZTA inner layer and two 0.1 mm alumina outer layers (diamonds). dition of freely (simply) supported edges, the deflection of P 0 increases as well. The increase of the deflection of P 0 is accompanied with the increase of the local stress. Furthermore, if the same force is locally distributed over the small area, the local stress is significantly increased. In general, we must keep in mind that the stress may not be the largest at the central point of the sample. When the square is elongated, so that the sides 2b in y–axis di- rection are many times larger than the sides 2a in x–axis direction, the solution for w essentially depends only on x–coordinate in the most part of the rectangle (when the y–coordinate is not near the values ±b). In this way we come to quasi-one-dimensional problem, which has ana- lytical solutions for some specific problems. For in- stance, the static solution for w* in Equation (A.1) with constant pressure in one dimension is a polynomial of the fourth order. We have used such analytical 1D solu- tions as a test for the reliability of our numerical proce- dure. Acknowledgement This research was financially supported by the ARRS (Slovenian Research Agency) in the frame of the na- tional project No. J2-9224. 5 REFERENCES 1 F. F. Lange, Transformation Toughening: Parts 1–5, Journal of Mate- rial Science 17 (1982) 5–6, 225–262 2 M. Rühle, N. Claussen, A. H. Heuer, Transformation and Microcrack Toughening as Complementary Processes in ZrO 2-Toughened Al 2O 3, Journal of the American Ceramic Society, 69 (1986) 3, 195–197 3 H. E. Lutz, N. Claussen, Duplex ceramics, Parts 1–2, Journal of the European Ceramic Society, 7 (1991) 209–226 4 S. Deville, J. Chevalie, G. Fantozzi, J. F. Bartolome, J. Requena, J. S. Moya et al., Development of advanced zirconia-toughened alumina nanocomposites for orthopaedic applications, Key Engineering Ma- terials, 264–268 (2004) 2012–2016 5 N. Mandal, B. Doloi, B. Mondal, R. Das, Optimization of flank wear using Zirconia Toughened Alumina (ZTA) cutting tool: Taguchi method and Regression analysis, Measurement, 44 (2011) 2149–2155 6 F. Sommer, R. Landfried, F. Kern, R. Gadov, Mechanical properties of zirconia toughened alumina with 10–24 vol. % 1Y–TZP reinforce- ment, Journal of the European Ceramic Society, 32 (2012) 4177–4184 7 A. V . Virkar, J. F. Jue, J. J. Hansen, and R. A. Cutler, Measurement of Residual Stresses in Oxide-ZrO 2 Three-Layer Composites, Journal of the American Ceramic Society, 71 (1988) 3, C-148–151 8 D. J. Green, P. Z. Cai, and G. L. Messing, Residual Stresses in Alu- mina–Zirconia Laminates, Journal of the European Ceramic Society, 19 (1999) 2511–2517 9 C. H. Hsueh, Modeling of Elastic Deformation of Multilayers due to Residual Stresses and External Bending, Journal of Applied Physics, 91 (2002) 12, 9652–9656 10 D. D. Barnett-Ritcey, and P. S. Nicholson, Failure Prediction Maps for a Model Al 2O 3 c-ZrO 2/Al 2O 3 Al 2O 3 Brittle Polycrystalline Trilayer Composite, Journal of the American Ceramic Society, 86 (2003) 1, 121–128 11 M. Ambro`i~, T. Kosma~, Optimization of the Bend Strength of Flat-Layered Alumina–Zirconia Composites, Journal of the Ameri- can Ceramic Society, 90 (2007) 5, 1545–1550 12 R. Bermejo, J. Pascual, T. Lube, R. Danzer, Optimal strength and toughness of Al 2O 3–ZrO 2 laminates designed with external or inter- nal compressive layers, Journal of the European Ceramic Society, 28 (2008) 1575–1583 13 E. Carrera, A. Ciuffreda, Bending of composites and sandwich plates subjected to localized lateral loadings: a comparison of various theo- ries, Composite Structures, 68 (2005) 185–202 14 K. F. Graff, Wave Motion in Elastic Solids (2012), Courier Corpora- tion 15 P. Lee, N. Chang, Harmonic in Elastic Sandwich Plates, Journal of Elasticity 9 (1979) 51–69 16 J. Kaplunov, D. A. Prikazchikov, L. A. Prikazchikova, Dispersion of Elastic Waves in a Strongly Inhomogeneous Three-Layered Plate, In- ternational Journal of Solids and Structures 113–114 (2017) 169–179, doi:10.1016/j.ijsolstr.2017.01.042 17 Z. J. Ai, C. J. Xu, G. P. Ren, Vibration of a pre-stressed plate on a transversely isotropic multilayered half-plane due to a moving load, Applied Mathematical Modeling, 59 (2018) 728–738, doi:10.1016/ j.apm.2018.02.027 18 M. Shariyat, M. Roshanfar, A new analytical solution and novel en- ergy formulations for non-linear eccentric impact analysis of com- posite multi-layer/sandwich plates resting on point supports, Thin-Walled Structures, 127 (2018) 157–168, doi:10.1016/j.tws. 2018.02.001 19 S. Zhang, J. Yin, H. W. Zhang, B. S. Chen, A two-level method for static and dynamic analysis of multilayered composite beam and plate, Finite Elements in Analysis and Design, 111 (2016) 1–18, doi:10.1016/j.finel.2015.12.001 20 T. Ye, G. Yin, Z. Su, Three-dimensional vibrational analysis of sand- wich and multilayered plates with general ply stacking sequences by a spectral sampling surface method, Composite Structures, 176 (2017) 1124–1142, doi:10.1016/j.compstruct.2017.06.008 21 W. Lacarbonara, M. Pasquali, A geometrically exact formulation for thin multi-layered laminated composite plates: Theory and experi- ment, Composite Structures, 93 (2011) 1649–1663, doi:10.1016/ j.compstruct.2010.12.005 22 L. V . Tran, S.-E. Kim, Static and free vibration analysis of multilay- ered plates by a higher-order shear and normal deformation theory and isogeometric analysis, Thin-Walled Structures, 130 (2018) 622–640, doi:10.1016/j.tws.2018.06.013 23 Y . Wang, Z. Li, W. Chen, C. Zhang, J. Zhu, Free vibration and active control of pre-streched multilayered electroactive plates, Interna- tional Journal of Solids and Structures, 180–181 (2019) 108–124, doi:10.1016/j. ijsolstr.2019.07.010 24 P. Zhang, C. Qi, H. Fang, C. Ma, Y . Huang, Semi-analytical analysis of static and dynamic responses for laminated magneto-electro-elas- tic plates, Composite Structures, 222 (2019) 110933, doi:10.1016/j.compstruct.2019.110933 25 B. Erbaº, J. Kaplunov, E. Nolde, M. Palsü, Composite wave models for elastic plates, Preceedings of Royal Society A, 474 (2018) 20180103, doi:10.1098/j.rspa.2018.0103 26 L. D. Landau, E. M. Lifshitz, Theory of Elasticity, Course of Theo- retical Physics, V ol. 7, Butterworth–Heinemann (1986) 27 V. V. Vasiliev, Modern Conceptions of Plate Theory, Composite Structures, 48 (2000), 93–84. 28 N. I. Robinson, Augmented Lévi–Michell equations for flexural plates, International Journal of Solids and Structures, 191–192 (2020) 497–508, doi:10.1016/j. ijsolstr.2019.12.021 M. AMBRO@I^, A. NIKONOV: DYNAMIC BIAXIAL STRESS ANALYSIS OF FLAT LAYERED CERAMIC COMPOSITES Materiali in tehnologije / Materials and technology 55 (2021) 2, 195–200 199 APPENDIX A: Dimensionless form of dynamic equa- tion Using the dimensionless quantities Equation (1) transforms into: Dw w t p ** ** * * * () ∇+= 22 2 24 1 ∂ ∂ (A.1) The Laplacian *2 contains derivatives with respect to renormalized coordinates x * and y * . Let’s take the symmetric composite with N layers (N is odd), denoted from "bottom" to "top" by indices i from 1t oN. Their parameters are h i * = h i /h, E i * = E i /E 0 , i * = i / 0 and i . The total thickness h is the sum of thick- nesses h i of all layers. The dimensionless flexural rigid- ity is: [] D E zz i i i N ii * * ** ()() = − − = ∑ 1 3 1 2 1 33 TB (A.2) Their symbols z iT * and z iB * mean the dimensionless z–coordinates of the top and the bottom side of the i–th layer, respectively. The neutral plane is at z * = 0 due to symmetry. We have derived Equation (A.2) in a similar manner as the derivation of flexural rigidity for homoge- neous plate in the textbook of Landau and Lifshitz. 26 The case with N = 1 (homogeneous plate) gives the expres- sion according to Equation (2) in physical units. The ef- fective density is: ** = = ∑ i i N i h 1 (A.3) The 2D stress tensor in i-th layer is: =− − ⋅ +− E z w x w y w xy i i ii * 1 1 2 2 2 2 2 2 * * * * * * * () ∂ ∂ ∂ ∂ ∂ ∂∂ () * * * * * * 1 22 2 2 2 −+ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ i * i w xy w y w x ∂ ∂∂ ∂ ∂ ∂ ∂ (A.4) We give a few notes about the numerical procedure that we have developed. A discrete quadratic network of equidistant points, typically with size 201×201, was used. Fourth-order spatial derivatives of w * , appearing in the double-Laplacian term, were derived on the basis of finite differences for neighboring points, in analogy to elliptic equations. As regards the time relaxation proce- dure, the problem of the second order in time derivative was transformed into the problem of first order by trac- ing variables w and the derivative w/ t at all discrete points in network. This enabled us to use a kind of Runge-Kutta method of second order, but simultaneously for all points in network. The numerical reliability for an appropriately chosen small time step was checked by halving the time step and comparing the results. M. AMBRO@I^, A. NIKONOV: DYNAMIC BIAXIAL STRESS ANALYSIS OF FLAT LAYERED CERAMIC COMPOSITES 200 Materiali in tehnologije / Materials and technology 55 (2021) 2, 195–200