Simulation of the Dynamic Behaviour of the ZTA Composites Obtained by Additive Technologies

This paper presents a physical and mathematical model that has been developed in the framework of the approach used in the computational mechanics of materials. The model is designed to enable the study of the patterns of deformation and fracture of ceramic composites with a transformation-hardened matrix that are obtained by additive technologies at the mesoscopic and macroscopic levels under intense dynamic loading. The influence of the loading rate on the formation of the fracture and energy dissipation fronts for composite materials, based on the Al2O3-20%ZrO2 system, is shown. Nonlinear effects under intense dynamic loading in the considered composites are associated with the processes of self-organization of structural fragments at the mesoscopic level, as well as the occurrence of martensitic phase transformations in matrix volumes adjacent to the strengthening particles.


Introduction
Functional ceramics composites represent a class of materials that exhibit characteristic properties, in addition to those inherent to ceramics, such as chemical and thermal stability. Functional ceramics composites usually exhibit one or more unique properties in the biologic, elastic, electric, magnetic, thermal or chemical regimes [1,2]. Among the most prominent type of functional ceramics are ZTA ceramic composites. Ceramic materials with transformation hardening are widely used in engineering practice as structural materials [3]. The vast majority of such ceramic materials were developed on the basis of zirconium dioxide, which was partially stabilized in the tetragonal phase [4]. Under the action of tensile stresses in the region of microconcentrators (at the boundaries of the particles of the strengthening phases, the tops of the cracks, etc.), the tetragonal phase is capable of experiencing a martensitic phase transition to the monoclinic phase. The phase transition is accompanied by the development of shear and volumetric deformations, which provide stress relaxation and the closure of crack surfaces [5]. The realisation of the hardening effect makes it possible to achieve strength characteristics (crack resistance, strength) in ceramic materials that are comparable with structural alloys. In order to increase the resistance of transformation-hardened ceramic materials to local contact loads and the formation of surface damage, the materials are further strengthened by the introduction of high-strength particles of carbides, nitrides, borides, and metal oxides, as well as metal (Ni, Nb) microparticles and fibres [6]. Ceramic composites in general are the most challenging materials to fabricate due to their brittle nature [7]. ZTA composites exhibit the same challenges regarding their fabrication. The conventional route to fabricate bulk ceramics is rather rudimentary and this fabrication process greatly limits the geometries. In addition to the limited geometric freedom during fabrication, post-fabrication following heat treatment presents the same problem considering that it is extremely difficult to modify their geometry through machining or other subtractive manufacturing techniques. This geometry limitation due to the fabrication process was a topic of interest for years. The most promising technique to mitigate and alleviate this challenge is additive manufacturing (AM) [8,9]. The implementation of additive manufacturing for the fabrication of ceramic composite structures will eliminate the geometry limitations from conventional fabrication methods. Extensive research efforts were made in recent years to design and demonstrate the viability of using a rational design of geometries to tune the properties of different structures [10][11][12][13].
Composite materials, obtained by additive technologies, that underwent the strengthening process have high crack resistance and strength under static loading conditions [6,9]. At the same time, the mechanical properties of transformation-hardened ceramic composite materials under highenergy influences are poorly studied [8]. Describing the mechanical behaviour of transformationhardened ceramic composites under dynamic loads is a significant challenge and can be solved by computer simulation methods. This is due to insufficient knowledge of the processes of evolution of the structure of composites under dynamic loading and the lack of adequate models of the mechanical behaviour of composites, which would allow taking into account the features of the course of deformation, the phase transformations, and the development of damage and destruction of materials.
The aim of this work was to develop an approach from the perspective of computational mesomechanics for studying the mechanical behaviour of composite ceramic materials with a transformation-hardened matrix under dynamic effects. This approach will provide the opportunity to create smart ceramic composites using additive technologies operating in a wide range of loading speeds.

Materials and Methods
Due to the small amount of additional impurities of the oxides Y2O3, MgO, and CeO, it is possible to partially stabilise the high-temperature phases (tetragonal or cubic) with different physicomechanical properties at low temperatures [3]. Under the action of an external load, martensitic phase transformations can occur in the ceramic matrix, which are accompanied not only by the development of shears, but also by a significant change in the volume of the material. The transition from the tetragonal to the monoclinic phase is accompanied by a 4% increase in the phase volume [4,5,14].
The local occurrence of martensitic transformations t → m in ZrO2 facilitates the closure of the edges of cracks or the collapse of micropores, thereby reducing the intensity of stress concentrators near defects. As a result, existing or newly formed microcracks become stable while maintaining the level of external load; an effect known as transformation hardening [15].
Currently, the present authors are undertaking intensive research to create high-strength ceramic materials with a zirconia matrix that is obtained by using additive technologies and is hardened by dispersed particles of metal borides and carbides.
To predict the mechanical behaviour of ceramic composites with a transformation-hardenable matrix under intense dynamic effects, in the present study the approach of computational mechanics of materials was used [16]. A common tool used to study in detail the features of the dynamic mechanical behavior of materials, including fracture dynamics, is numerical simulation in Lagrangian coordinates using explicit integration schemes. The most famous "Lagrangian" numerical method used to solve such problems is the finite element method (FEM) [17,18]. The welldeveloped mathematical apparatus of this method allows you to implement complex rheological models of the mechanical response of materials and solve complex problems of mechanics of a deformable solid. Note that when using this method for solving dynamic problems, well-known difficulties arise associated with modeling of multiple fracture, accompanied by intense mass transfer and mixing. In addition to the FEM, at present, numerical methods of the discrete approach to the description of the medium, in particular, the methods of discrete elements (MDE) [19][20][21] are actively used to model dynamic processes of multiple fracture and fragment transfer. Despite the well-known advantages in fracture modeling, MDEs have a number of limitations, which usually include a lower order of accuracy in space in comparison with the traditionally used finite element realizations and insufficient development of the mathematical formalism limited by the use of quasistatic mechanical models of materials. The inelastic behavior of brittle materials is mainly associated with the accumulation of microdamages (the formation of pores and cracks of various sizes, their growth or healing, consolidation, collapse, etc.) [22][23][24]. The involvement of crystal lattice defects in the process of deformation of such materials takes place only at high pressures and temperatures [25,26]. Therefore, to describe the dynamic inelastic behavior of brittle materials, modified plasticity models are used, which take into account the high sensitivity of the inelastic response parameters to the pressure value, as well as the complex nature of the relationship between shear and volumetric plastic deformations (unassociated flow laws). Moreover, the features of the used plastic potential and the limiting surface determine the scope of the model.
The dynamic models of brittle isotropic materials, such as the Taylor -Chen -Kuszmaul model [27,28], Holmquist -Johnson -Cook model [29], "Continuous Surface Cap" (CSC) model [30], the Riedel -Thom -Hermaier model (RHT) [31] and several others [32]. These models include a significant number of model parameters and, when correctly determined, provide an adequate description of the features of dynamic deformation and fracture of brittle materials, including the sensitivity of the properties to pressure and loading speed.
In the described models, a traditional approach is used that takes into account the sensitivity of the inelastic response of the material to the dynamics of changes in the local stress-strain state by introducing the dependences of the parameters of the plasticity and fracture models on the strain rate. As a rule, in the experimental determination of the dependences of the parameters of the applied model of plasticity or fracture on the strain rate, integral values are used for the entire sample. Moreover, the local values of strain rates can differ significantly from the integral value. Changes in the effective shear strength of the condensed phase was found from the current value of the shear strength of the damaged material, taking into account changes in the shear strength of the condensed phase.
Within the framework of this approach, a hierarchical model of a structured environment can be built an optical, probe scanning, and electron microscopy data can be used to create a realistic model. A model representative volume of a material with a cross section of 5 × 5 × 5 μm containing a stochastic system of ZrO2 particles is shown in Figure 1a. The matrix particles and inclusions have a shape that is close to tetrahedral with a characteristic size of ~ 0.5 μm. The concentration of ZrO2 in the cell corresponds to the average values of the corresponding concentration in the macroscopic sample of the material. Data on the structure of the composite were obtained by studying the properties of the composite that was obtained by using additive technologies [6,8,9].
The kinematics of a structured medium in a representative cell with dimensions exceeding 10 nm can be described by parameters of mass velocity ui and strain rate ij. The dimensions of the microparticles of the material in this case are tens of nm, which means that continuum representations can be used to describe inelastic deformations.
The stresses arising in the cell due to temperature and mass velocities have distributions. The average mass density is determined by relation (1), taking into account the distribution of the local mass density of the matrix and the strengthening particles. Effective mechanical parameters are determined by Equations (1) - (4). The components of the effective tensors of the strain rate and the rate of bending-torsion are determined by formulas in Equation (3), and the components of the effective stress tensor are determined by Equation (4): where < > is the average mass velocity; ui is the local mass velocity in a representative cell; V is the volume.
where = ∭ is the total specific internal energy in a representative cell. When studying the patterns of deformation and fracture of ceramic composites under intense dynamic effects, one of the key problems is obtaining information that allows us to construct the governing equation for the macroscopic level. Within the framework of the approach used, varying the conditions of exposure to a representative cell of the composite allows one to obtain data on the response of a structured medium under various conditions of deformation.
For example, for a representative cell that is deformed in the front of a plane shock wave, the boundary conditions can be specified in accordance with the diagram in Figure 1b At the initial moment of time, the medium of the computational domain is in an unstressed state at the temperature T0 = 293 K, and the strains and the velocity of the material particles of the condensed phase are zero. The computational domain is loaded on the surface S1 by setting the velocity of material particles that are directed into the sample.
where  is the mass density near the boundary, un is mass velocity in the direction normal to the boundary, the value of parameter C decreases from the longitudinal sound velocity in the elastic precursor to CB, and F(xk, t) is the function that determines the shape and duration of pulse loading.
On the boundaries of the inclusion-matrix the following conditions are given: where + , − are components of the normal vector to the outer and inner surfaces of the interface. In violation of the interphase interface of the matrix and inclusion, the boundary conditions take the form: where n,  are normal and tangential stresses at points at the interface, and η is the coefficient of friction. The initial conditions for the computational domain for the case under consideration are taken in the form: The effective elastic modulus was determined using Equation (9) from the calculated longitudinal CL and volume CB sound velocities: where ̅ is the average mass density in the loaded region of the cell behind the front of the corresponding wave. The dynamics of a structured medium in a cell is described in the Lagrangian reference frame by the equations of conservation of mass, momentum and energy: where ij is the microstress tensor components,  is mass density, ui represents the components of the displacement vector, ij represents the components of the strain rate tensor, and E is the local specific internal energy per unit mass. The increments of the components of the strain tensor are represented as the sum of the elastic and inelastic components: where ̇ are the components of the strain rate tensor. The components of the inelastic strain rate tensor are presented in the form of components related to dislocation plasticity mechanisms and mechanically activated mechanisms of martensitic phase transformations: Such a representation takes into account the dilatancy that is caused only by the increment of the bulk inelastic deformation of the martensitic transformation, and the deviator of the rate of inelastic deformation is determined by the contributions of shears from both plasticities: where is the Kronecker symbol.
where Sij are stress tensor deviator components. The components of the Kirchhoff stress tensor are represented as the sum of the pressure p and the deviator Sij: Pressure p in the range up to 10 GPa can be calculated using the equation of state: = 1 + 1 2 + 1 3 + Γ 0 , under compression > 0, where K1, K2, K3 are the material constants,  = (/0) -1, and  is the Grüneisen coefficient. The values of the constants for the mixture of phases of zirconium dioxide can be, to a first approximation, obtained within the framework of the mixed model. The stress tensor deviator is calculated from the solution of the relaxation Equation (17): The defining equation in a damaged medium can be represented as where ∆ = [∆ ] is the increment of the damage parameter over a period of time t, ∆ are the inelastic strain rates, and f is the ultimate strain at the time of macroscopic fracture. The ultimate strain at the time of failure for a particular medium is considered as a function of pressure and temperature. To describe the ultimate strain at the moment of destruction of the material point of a condensed medium imitating oxide ceramic material, the relation proposed in [34] was used: where D1, D2 are the material coefficients; P * = P/PHEL and PHEL is the pressure corresponding to the Hugoniot elastic limit of the main condensed phase of the ceramic; T* = (T-Tr)/(Tm-Tr), where T is the temperature of a material point (representative volume) on an absolute scale, Tr is the room temperature, and Tm is the melting point or dissociation of the main condensed phase of the medium. The current value of the shear strength of the damaged material is found taking into account changes in the shear resistance of the condensed phase. The shear resistance at the material point of the damaged medium is determined by relations [33]: where -variable shear strength of condensed phases with damage at a lower level. When solving dynamic problems and problems associated with thermomechanical effects on ceramic structured materials, it is necessary to take into account the influence of pressure, temperature and strain rate on the effective shear strength of the material. In the present work, to simulate the deformation of oxide ceramic materials, we used the relations [34,35]: = 2 ( * ) 2 (1 + 2 ln), where A1, B2, C1, C2, m1, m2 are the material coefficients; and ̇ is normalised strain rate tensor intensity.
The numerical values of the elastic modulus, mass velocity and the coefficients of the equations for the crystalline phases of α-Al2O3, t-ZrO2, depend on the concentration of chemical impurities, the presence of nanoscale inclusions, and defects at the submicron level [36][37][38][39][40][41][42][43][44][45][46]. The values of the coefficients selected for calculations are presented in the Tables 1 and 2.

Results
To study the effect of the loading rate on the fracture toughness and crack resistance of Al2O3 -20%t-ZrO2 ceramic composites with a mass content of submicron particles of t-ZrO2, on the highspeed compression of model samples in shock waves, and the high-speed tension in the region of interaction of unloading waves. Figure 2 shows the distribution in the model volume of the effective plastic strain rate in the Al2O3-20%t-ZrO2 composite upon compression in the shock front with an amplitude of 33 GPa at successive time intervals. The simulation results show that an increase in the effective rate of plastic deformation is localised at the mesoscopic level in the region of reinforcing submicron particles and weakly depends on the shape and size of those particles (in the studied range of variation). In this case, regions with a plastic flow velocity are not observed behind the wave front of the elastic precursor, indicating the fragile behaviour of a representative cell. Figure 3 shows the distribution of the specific work of fracture in the model volume of the Al2O3-20%t-ZrO2 composite upon compression in the shock front with an amplitude of 33 GPa at successive time intervals. The simulation results show that the specific work of fracture is realized behind the front of the elastically precursor and is localized at the mesoscopic level in the region of reinforcing particles, at the same time, not exceeding the maximum value of 8.5 kJ/Kg. Dynamic fracture toughness changes weakly for the entire range of loading rates.    The results shown in Figure 4 for successive time intervals testify the formation of local damages of the ceramic matrix in the zone of shock transition. There was no considerable alteration detected regarding the size and configuration of local damage areas after the increase in time of the hydrostatic compression by the shockwave front. Microdamages are localised at the mesoscopic level in the areas of reinforcing particles. In the region beyond the elastic precursor, the regions of microcrack development in the ceramic matrix are observed. However, the ceramic matrix retains its integrity, up to the plastic rarefaction wave, while the hardening particles are completely destroyed. Figure 5 shows the results of modelling the loading by the various speeds of a representative cell. Free surface velocity profiles and dissipated energy velocity are shown. The energy dissipation rate increases with increasing loading speed. However, at speeds above 1000 m/s, a sharp change in the nature of the dissipation curves is observed, instead of a monotonic increase, the curve is divided into 2 parts. This effect is a consequence of the inertia of the samples, which begins to play a significant role in high-speed material deformation. To solve this problem, it is necessary to introduce additional parameters into the design scheme that take into account high strain rates as a function of relaxation time. In the ideal case, this information should be obtained on the basis of a synthesis of experiments on the dynamic deformation of samples of the material under consideration of the required scale using various types of stress state. The determination of these parameters can be implemented in dynamic loading schemes using split Hopkinson rods. Comparison of simulation results with experimental data showed good agreement in the speed of the free surface [8].
(a) (b) Figure 5. Results of modelling loading using various speeds: (a) free surface velocity profiles; (b) dissipated energy rate. Figure 6 shows the simulation results of the specific work of the fracture, which went into energy dissipation, for different loading rates. The simulation results show that as the loading rate increases, the contribution of the martensitic phase transition to the dynamic viscosity gradually increases. When the loading speed exceeds 1000 m/s, the fracture toughness increases with increasing concentration of hardening particles of zirconium dioxide. Analysis of the fracture work curves showed that for the aluminum matrix, a purely brittle fracture pattern is observed with a classical elastic precursor, and a fracture wave following it. At low loading rates, the contribution of the reinforcing particles to the overall work of fracture is insignificant, but there is a slight increase in the strength of the samples. With an increase in the loading rate, the contribution of reinforcing particles to the overall work of fracture increases. Moreover, the nature of the accumulation of loading energy has a similar character, with a specific time of about 20 ns.

Discussion
An analysis of the results of the modelling of the composites showed that the volume content of transformation-strengthening particles leads to the formation of an inhomogeneous front and affects the width of the transition. The magnitude of the loading rate for the considered composites affects the strain rate in the region of the shock transition. In this case, the main deformation energy is dissipated in the region of strengthening particles. The formation and nucleation of mesocracks is accompanied by the development of plastic deformations in the region of the particles that undergo martensitic phase transformations with the dilatancy effect.
Analysis of the simulation results has shown that the dependence of the strength characteristics of ceramic composite on strain rate is strongly non-linear. It is generally corresponding to the modelling data, obtained on a similar class of brittle materials [47]. An increase in the loading rate leads to an increase in the rate of energy dissipation, which is realised in the accumulation of microdamage at the front of the elastic precursor. Moreover, for the composites that are considered in the above regarding the front of the elastic precursor, the destroyed regions continue to have a deterrent value in the development of further damages of the composite matrix. If pores are present in the structure of a composite, then damage primarily arises precisely at the boundaries between the pore and the grain [48]. In this simulation, we used a non-porous medium, therefore, the effect of the presence of pores was not taken into account.
The aluminium oxide matrix does not undergo plastic flows and a purely brittle fracture pattern is observed. At low loading rates, a purely brittle fracture effect is realised, mainly with an aluminium oxide matrix. In this case, the contribution of strengthening particles to the strength characteristics is minimal. However, at high loading rates, there are clear energy jumps that are associated with the effect of a change in the volume of ZrO2 particles during the phase transition t → Ortho I and t → Ortho II, or from the tetragonal to monoclinic phase. The results show a qualitative dependence of the nature of energy dissipation on the loading rate during dynamic modeling of brittle media [49].
It is noteworthy that the magnitude of the strengthening effect from filling the nanostructured oxide-aluminium matrix with submicron inclusions of t-ZrO2 depends on the particle size. For particle sizes less than 400 nm, the occurrence of mechanically activated martensitic phase transitions t → m is difficult. As a result, during intense deformation in the crack initiation and propagation zone, dispersed particles of t-ZrO2 do not undergo phase transformations, but they do play the role of weakly deformable particles.
In the framework of the model concepts that were used in the present study, the decrease in crack resistance is due to the intensive nucleation of microcracks, both in the volume of reinforcing inclusions and in the volume of the nanostructured matrix. The nature of the fracture under these conditions is not associated with an increase in the size of mesoscopic cracks, but rather with the merging of multiple submicron cracks formed in the bulk of the composite. Damage originates both in the nanostructured matrix and in the submicron reinforcing particles. However, a significant part of the energy is spent on phase transitions, due to which a strengthening effect occurs. In addition, some studies provide a quantitative comparison of the fracture energy and fracture toughness of composites [50]. In this paper, we departed from this formalism and used the concept of the specific work of destruction, and the contribution of the work of destruction of each component of the composite was shown. This approach is original and allows us to evaluate the contribution of each component of the composite to the work of failure in a wide range of loading rates. This approach will allow, at the design stage of ceramic composites, to predict the dynamic characteristics of the material and create «smart» materials for a wide range of applications.
Within the framework of the presented model, it can be said that an increase in the fracture toughness of composite materials that is obtained by additive technologies can be achieved by reducing the grain size of the matrix below 500 nm, and by filling the oxide matrix with reinforcing particles with sizes less than 500 nm.

Conclusions
The simulation results show that the nonlinear effects of the mechanical behaviour of ZrO2 -Al2O3 ceramic composites with a transformation-hardened matrix, that is obtained by additive technologies, are manifested at shock loading amplitudes that are close to or exceeding the Hugoniot elastic limit σHEL.
Nonlinear effects under intense dynamic effects on the considered composites are associated with the processes of self-organisation of deformation modes at the mesoscopic level (the formation of block substructures), as well as the occurrence of martensitic phase transformations in matrix volumes adjacent to the strengthening particles.
The presence of submicron reinforcing particles leads to a change in the shape of the local damage area, while the increase in the relative number of damaged particles can be defined as the formation of shear mesocracks in the volume of the composite.
The modeling approach presented in this paper can be used to determine the dynamic characteristics of ceramic composites up to impact loads of 1000 m/s. With an increase in the loading speed, in order to achieve an adequate operation of the design scheme, it is necessary to introduce additional relaxation equations that take into account the formation and accumulation of damage over time. This non-trivial task requires additional theoretical and experimental studies of composites obtained by additive technologies.