Investigation of the Dynamic Behaviour of the ATZ Ceramic Composites Obtained by Additive Technologies

1. Introduction

Ceramic composites based on alumina and zirconia systems are a class of functional materials characterized not only by strength properties, but also by chemical and thermal stability. They are also used in biological applications, to manufacture electric, magnetic and other products [1,2]. One of the best known and actively studied classes of functional ceramic composites with transformation hardening are ATZ composites (zirconia reinforced with alumina), a matrix of zirconia. They are used not only in engineering as structural materials [3] for the manufacture of cutting tools and wear parts, but also in medicine (especially in orthopedic prostheses). When preparing such composites, the high temperature phase (tetragonal or cubic) is partially stabilized by small additional impurities of Y₂O₃, MgO and CeO oxides [4]. The tetragonal phase can undergo a martensitic phase transition to a monoclinic phase under the influence of tensile stresses in the region of microconcentrators (at the boundaries of the particles of the hardening phases, at the tops of cracks, etc.). The transition from the tetragonal to the monoclinic phase is accompanied by a 4% increase in the volume of the phase [4,5]. At a pressure of 12.5 GPa, the tetragonal phase becomes monoclinic [6]. The phase transition is accompanied by the development of shear and volume deformations. These deformations provide stress relaxation and crack surface closure. The local martensitic transformation t → m in ZrO₂ contributes to closing crack edges or collapsing micropores, thereby reducing the intensity of stress concentrators near defects. This effect, known as transformational hardening [7], stabilizes existing or newly formed microcracks while maintaining the level of external loading. Realizing the hardening effect makes it possible to achieve strength properties (crack resistance, strength) of ceramic materials comparable to those of structural alloys. It should be noted that the ATZ composites obtained by means of additive technologies are particularly interesting for study, have high crack resistance and strength under static loading conditions [8,9], and have high prospects of application in industry. Remember that ceramics are generally the most difficult materials to work with [10], and ATZ composites are no exception. Conventional bulk ceramic production techniques are very primitive and severely limited by mold geometry. In addition to the limited geometric freedom in the manufacturing process, the same problem exists in the post-heat treatment process, as it is extremely difficult to change their geometry by machining or other subtractive manufacturing methods. This geometric constraint imposed by the manufacturing process has been studied for many years.

The most promising method to mitigate this problem is additive manufacturing (AM) [11,12]. The geometric limitations inherent in traditional manufacturing methods can be eliminated by introducing additive manufacturing for the fabrication of CE-RAMIC composite structures. In recent years, extensive research efforts have been devoted to the development and demonstration of the feasibility of the use of rational geometry design to tune the properties of various structures [13,14,15,16]. On the other hand, the mechanical properties of ceramic composite materials that undergo transformation hardening under energetic influences are not well understood.

The evaluation of the dynamic behavior of such composites is a serious problem. It can be solved by computer modeling methods. This has been due to the lack of understanding of the processes involved in the evolution of the structure of composite materials under dynamic loading and the lack of adequate models of the mechanical behavior of composite materials which can take account of the characteristics of the deformation process, phase transitions, damage development and destruction of the material. The aim of this work was to develop an approach from the point of view of computational mesomechanics to study the mechanical behavior of ceramic composites with a hardened transformation matrix under dynamic effects. This approach will make it possible to create intelligent ceramic composite materials by means of additive technologies that will operate in a wide range of loading speeds.

2. Materials and Methods

In order to predict the mechanical behavior of ceramic composites with a transformation-hardenable matrix under intense dynamic effects, the computational mechanics of materials has been used in this study [17]. Numerical simulation in Lagrangian coordinates using explicit integration schemes is a commonly-utilized tool for studying features of dynamic mechanical behavior of materials, including fracture dynamics. The finite element method (FEM) [18,19] is a widely used "Lagrangian" numerical method for solving such problems. With its well-developed mathematical tools, this method enables the implementation of complex rheological models to analyze the mechanical response of materials and solve intricate problems related to the mechanics of deformable solids. It is important to note that utilization of this method for dynamic problem-solving poses significant challenges in the modeling of multiple fractures that are accompanied by intense mass transfer and mixing. In addition to the finite element method (FEM), numerical methods using discrete elements (MDE) are currently used to simulate dynamic processes involving multiple fractures and fragment transfer. These methods provide a discrete approach to describing the medium and have been actively used in [20,21,22]. Despite their well-known advantages in fracture modeling, MDEs have limitations such as reduced spatial accuracy compared to traditional FEM and insufficient mathematical formalism due to using quasi-static material models. The brittle inelastic behavior of materials is mainly associated with the accumulation of microdamage (the formation of pores and cracks of various sizes, their growth or healing, consolidation, collapse, etc.) [23,24,25]. These phenomena are discussed in [23,24,25]. The presence of defects in the crystal lattice during plastic deformation of such materials occurs only under high pressure and high temperature conditions as mentioned in [26,27]. Therefore, to explain the dynamic inelastic behavior of brittle materials, modified plasticity models are utilized. These models consider the sensitive nature of inelastic response parameters to pressure values, as well as the complicated relationship between shear and volumetric plastic deformations (unassociated flow laws). In addition, the properties of the plastic potential and the boundary surface are also important in the determination of the model's scope.

The Taylor-Chen-Kuszmaul model [28,29], Holmquist-Johnson-Cook model [30], Continuous Surface Cap (CSC) model [31], Riedel-Thom-Hermaier (RHT) model [32], and several others [33] are among the dynamic models for brittle isotropic materials. These models include a considerable number of parameters which, once accurately determined, sufficiently describe the characteristics of the dynamic deformation and fracture of brittle materials, while also taking into account the sensitivity of the properties to pressure and loading rate.

In these models, the traditional method is to take into account the sensitivity of the inelastic response of the material to changes in the local stress-strain state dynamics by incorporating the dependence of the plasticity and fracture model parameters on the strain rate. Experimental determination of the relationship between plasticity or fracture model parameters and strain rate is generally based on the average values of the entire specimen. However, the individual strain rate values may differ significantly from the overall average. Changes in the effective shear strength of the condensed phase have been determined by considering the current value of the shear strength of the damaged material and the changes in the shear strength of the condensed phase.

A hierarchical model of a structured environment can be constructed from optical, probe scanning, and electron microscopy data to provide a realistic representation of the environment. Figure 1a shows a model of a representative volume of material containing a stochastic system of Al₂O₃ particles. The cross section is 5 × 5 × 5 μm. The matrix particles and inclusions are almost tetrahedral in shape, having a characteristic size of approximately 0.5 μm. The concentration of aluminum oxide in the cell corresponds to the average concentration in the macroscopic sample. The data on the structure of the composite have been obtained by studying the properties of the composite produced by means of additive technologies [9,10].

The kinematics of the structured medium can be described by the mass velocity parameter ui and the strain rate $\dot{\epsilon }$_ij in a representative cell with dimensions exceeding 10 nm. The continuum representation can be used to describe inelastic deformations because the microparticles of the material are in the tens of nm.

Temperature and mass velocity stresses are distributed in the cell. Taking into account the distribution of the local mass density of the matrix and the reinforcing particles, equation (1) determines the mean mass density. Equations (1) to (4) determine the effective mechanical parameters. The effective tensors of strain rate and bending torsion rate are calculated using expressions from Equation (3), while the effective stress tensor is derived from Equations (4)

$$\langle \rho \rangle =\frac{1}{V}\underset{V}{\overset{ }{{\displaystyle \iiint }}}\rho \left(x_k\right)dV;$$

(1)

$$\langle u_i\rangle =\frac{1}{V}\underset{V}{\overset{ }{{\displaystyle \iiint }}}{u}_i\left(x_1,x_2,x_3\right)d{x}_{1}d{x}_{2}d{x}_3,$$

(2)

here $\langle u_i\rangle$ - the average mass velocity; u_i – the local mass velocity in a representative cell; V – the volume.

$$\langle {\dot{\epsilon }}_{ij}\rangle =\frac{1}{2}\left[\frac{\partial \langle u_i\rangle }{\partial x_j}+\frac{\partial \langle u_j\rangle }{\partial x_i}\right], {\dot{\omega }}_{ij}=\frac{1}{2}\left[\frac{\partial \langle u_i\rangle }{\partial x_j}-\frac{\partial \langle u_j\rangle }{\partial x_i}\right];$$

(3)

$$\langle {\sigma }_{ij}\rangle =\frac{1}{\langle \rho \rangle }\frac{d{E}_V}{dt\langle \dot{{\epsilon }_{ij}}\rangle },$$

(4)

here $E_V=\underset{V}{\overset{ }{{\displaystyle \iiint }}}EdV$ – the total specific internal energy in a representative cell.

Obtaining the information necessary to construct the governing equation at the macroscopic level is a major challenge when studying the deformation and fracture patterns of ceramic composites under intense dynamic loading. In this framework, information about the response of a structured medium under different forms of deformation is obtained by modifying the loading conditions of a representative cell within the composite.

For example, the boundary conditions for a typical cell deformed at the front of a planar shockwave are shown in Figure 1b.

At the initial moment of time, the medium of the computational domain is in an unstressed state at a temperature T₀ = 293 K, the strains and the velocity of the material particles of the condensed phase are zero. The computational domain is loaded on the surface S₁ by setting the velocity of material particles directed inside the sample.

$$\begin{gathered} {\sigma }_n\left(x_k,t\right)=F\left(x_k,t\right), x_kϵS_1; \\ {u}_1\left(x_k,t\right)=0, x_kϵS_2{\displaystyle \cup }{S}_3;u_2\left(x_k,t\right)=0, x_kϵS_4{\displaystyle \cup }{S}_5; \\ P=P_{ref}+\left[u_n\left(x_k,t\right)-u_{ref}\right]\left[\rho C\right], x_kϵS_6; \end{gathered}$$

(5)

here ρ - the mass density near the boundary, u_n – mass velocity in the direction normal to the boundary, the value of parameter C decreases from the longitudal sound velocity in the elastic precursor to C_B, F(x_k, t) – the function that determines the shape and duration of pulse loading, precursor to C_B, F(x_k, t) – the function that determines the shape and duration of pulse loading.

On the boundaries of the inclusion-matrix are given conditions:

$${\sigma }_{ij}{n}_j^{+}={\sigma }_{ij}{n}_j^{-}; u_i^{+}=u_i^{-},$$

(6)

here $n_J^{+}, n_j^{-}$ - components of the normal vector to the outer and inner surfaces of the interface.

The boundary conditions take the form in violation of the interphase interface of the matrix and inclusion:

$${\sigma }_n^{+}=-{\sigma }_n^{-}; {\sigma }_{\tau }=\eta P; {\sigma }_{\tau }^{+}={\sigma }_{\tau }^{-},$$

(7)

here σ_n, σ_τ - normal and tangential stresses at points at the interface, η - the coefficient of friction.

The initial conditions for the computational domain for the case under consideration are taken in the form:

$$u_i\left(x_k\right)=u_i^{\left(0\right)}, u_{ij}\left(x_k\right)=0, \rho \left(x_k\right)={\rho }^{\left(0\right)}\left(x_k\right), {\epsilon }_{ij}\left(x_k\right)=0,E\left(x_k\right)=0.$$

(8)

Effective elastic modulus was determined using (9) from the calculated longitudinal C_L and volume C_B sound velocities:

$$\langle \mu \rangle =\overline{\rho }{C}_s^2, \langle K\rangle =\overline{\rho }{C}_B^2, C_L^2=C_B^2+\frac{4}{3}{C}_s^2,$$

(9)

here $\overline{\rho }$ – 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 [34]:

$$\frac{\partial {\sigma }_{ij}}{\partial x_j}=\rho \frac{\partial u_i}{\partial t}, \frac{d{u}_i}{dt}=\rho \frac{\partial u_i}{\partial x_i}, \rho \frac{dE}{dt}={\sigma }_{ij}{\dot{\epsilon }}_{ij},$$

(10)

here σ_ij – the microstress tensor components, ρ - mass density, u_i – the components of displacement vector, ε_ij – the components of the strain rate tensor, E – 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:

$${\dot{\epsilon }}_{ij}={\dot{\epsilon }}_{ij}^e+{\dot{\epsilon }}_{ij}^p,$$

(11)

here ${\dot{\epsilon }}_{ij}$ - 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:

$${\dot{\epsilon }}_{ij}^p={[{\dot{\epsilon }}_{ij}^p]}_S+{\left[{\dot{\epsilon }}_{ij}^p\right]}_m.$$

(12)

Such a representation takes into account dilatancy 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 plasticity:

$${\dot{[\epsilon }}_{ij}^p{]}_s={[{\dot{e}}_{ij}^p]}_S, {\left[{\dot{\epsilon }}_{ij}^p\right]}_m={\left[\frac{1}{3}{\dot{\epsilon }}_{kk}^p\right]}_m{\delta }_{ij}+{\left[{\dot{e}}_{ij}^p\right]}_m,$$

(13)

here ${\delta }_{ij}$ - Kronecker symbol.

$${\dot{e}}_{ij}^p=\frac{3}{2}{\dot{\epsilon }}_u^p\left[\frac{s_{ij}}{{\sigma }_u}\right], {\sigma }_u=\sqrt{(3/2)S_{ij}{S}_{ij}} , {\dot{\epsilon }}_u^p=\sqrt{(2/3){\dot{e}}_{ij}^p{\dot{e}}_{ij}^p} , {\dot{e}}_{ij}^p={\left[{\dot{e}}_{ij}^p\right]}_S+{\left[{\dot{e}}_{ij}^p\right]}_m$$

(14)

here S_ij - stress tensor deviator components.

The components of the Kirchhoff stress tensor are represented as the sum of the pressure p and the deviator S_ij:

$${\sigma }_{ij}=-p{\delta }_{ij}+S_{ij}.$$

(15)

Pressure p in the range up to 10 GPa can be calculated using the equation of state:

$$\begin{gathered} p=K_1\theta +K_1{\theta }^2+K_1{\theta }^3+\mathsf{\Gamma }{\rho }_{0}E, \mathrm{under} \mathrm{compression} \theta >0, \\ p=K_1\theta , \mathrm{tensile} \theta \le 0, \end{gathered}$$

(16)

here K₁, K₂, K₃ - material constants, θ = (ρ/ρ₀) - 1, Γ - Gruneisen 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:

$$\frac{d{S}_{ij}}{dt}=2\mu \left({\dot{\epsilon }}_{ij}^e-\frac{1}{3}{\dot{\epsilon }}_{kk}^e{\delta }_{ij}\right), \frac{d{S}_{ij}}{dt}={\dot{S}}_{ij}-S_{ik}\dot{{\omega }_{jk}}-S_{jk}\dot{{\omega }_{ik}}, {\dot{\epsilon }}_{ij}^e={\epsilon }_{ij}-{\left[{\dot{\epsilon }}_{ij}^p\right]}_m-{\left[e_{ij}^p\right]}_S$$

(17)

The defining equation in a damaged medium can be represented as

$${\sigma }_{ij}=\left(1-D\right)\left[-P{\delta }_{ij}+S_{ij}\right], D={\displaystyle \sum }_{i=1}^m\frac{\left[\Delta {\epsilon }_{eq}^p\right]}{{\epsilon }_f},$$

(18)

here $\Delta D_i=\frac{{\left[\Delta {\epsilon }_{eq}^p\right]}_i}{{\epsilon }_f}$ - increment of the damage parameter over a period of time Δt, $\Delta {\epsilon }_{eq}^p$ - inelastic strain rates, ε_f - 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 materials, the relation proposed in [35] was used:

$${\epsilon }_f=D_1{\left(P^{*}+T^{*}\right)}^{D_2},$$

(19)

here D₁, D₂ – material coefficients; P^* = P/P_HEL, P_HEL - pressure corresponding to the Hugonios elastic limit of the main condensed phase of the ceramic; T* = (T-T_r)/(T_m-T_r), T - temperature of a material point (representative volume) on an absolute scale, T_r – room temperature, T_m - melting point or dissociation of the main condensed phase of the medium.

The current shear strength of the damaged material incorporates changes in shear resistance of the condensed phase. Relations [36] determine the shear resistance at the material point of the damaged medium:

$${\sigma }_s={\sigma }_s^c-D\left({\sigma }_s-{\sigma }_f\right),$$

(20)

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 this work, to simulate the deformation of oxide ceramic materials, we used the relations:

$${\sigma }_s^c=A_1{\left(P^{*}+T^{*}\right)}^{m_1}\left(1+C_1\ln \dot{\epsilon }\right),$$

(20)

$${\sigma }_f=B_2{\left(P^{*}\right)}^{m_2}\left(1+C_2\ln \dot{\epsilon }\right),$$

(21)

here A₁, B₂, C₁, C₂, m₁, m₂ – material coefficients, $\dot{\epsilon }$ - normalized strain rate tensor intensity.

The elastic modulus, mass velocity, and coefficients of α-Al₂O₃ and t-ZrO₂ crystalline phases' equations are dependent on the concentration of chemical impurities, presence of nanoscale inclusions, and submicron-level defects [32, 36,37,38,39,40]. The coefficients chosen for calculations are listed in Table 1 and Table 2.

The finite difference method was used to simulate a microsecond pulse load on a cell using the ANSYS-19 / AUTODYN software package.

3. Results

To investigate the effect of loading rate on the fracture toughness and resistance to cracking of t-ZrO₂-n% Al₂O₃ ceramic composites with submicron Al₂O₃ particles, rapid compression and tensile tests were performed on specimens exposed to shock waves and unloading wave interaction regions. In this case n is 5, 10, 15.

Figure 2 shows the effective plastic strain rate distribution throughout the model volume of the t-ZrO₂-15% Al₂O₃ composite during 1200 m/s shock front compression at various time intervals.

The simulation results show that the effective rate of plastic deformation is localized at the mesoscopic level in the region of reinforcing submicron alumina particles and is weakly dependent on the shape and size of these particles (within the studied range of variation). The overall picture of the evolution of the structure under the influence of high-speed loading is also not significantly affected by areas of representative volume in the pores. In this case, the reinforcing particles are completely destroyed, forming pore structures filled with the volume of the reinforcing particle material, while maintaining the overall shape of the matrix. This indicates the viscous behavior of the destruction of the representative cell.

Figure 3a shows the damage parameter distribution throughout the model volume of the t-ZrO₂-15% Al₂O₃ composite during 1200 m/s shock front compression. The zirconia dioxide matrix is hidden from view for better visualization of the material. The dependence of the specific destruction work is shown in Figure 3b.

The simulation results on Figure 3a show that when the shock front passes through a representative cell, the reinforcing particles of alumina are completely destroyed, but the shape of the domains remains unchanged. A fragile collapse of inclusions is observed simultaneously Figure 3b shows an increase in the specific work of destruction due to an increase in the contribution of the open-hearth transition, which is characteristic of zirconia.

Figure 4 shows the distribution of damage parameters in the model volume in the t-ZrO₂-15% Al₂O₃ composite upon compression in the shock front with an amplitude of 1200 m/s at successive time intervals.

The results shown in Figure 4 for successive time intervals demonstrate the formation of local ceramic matrix damage in the shock transition zone. After the increase of the time of hydrostatic compression by the shock wave front, no significant changes in the size and configuration of the local damage areas were detected. In the region of reinforcing particles, microdamages are localized at the mesoscopic level. In the region beyond the elastic precursor, microcrack regions are observed in the ceramic matrix. While the hardening particles are completely destroyed, the ceramic matrix retains its integrity until the plastic rarefaction wave.

The simulation results of the specific fracture work for ATZ ceramic composites loaded at different velocities are shown in Figure 5, Figure 6 and Figure 7.

The simulation results shown in Figure 5, Figure 6 and Figure 7a indicate an increase in the contribution to the total work of specific fracture of composite inclusions with increasing percentage of particles. At that, at the initial stages of the shock front passage through the composite material, the work of fracture is realized first of all in the material of aluminum inclusions, which corresponds to the time of 4 ns. Further passage of the impact front deep into the composite structure leads to redistribution of the specific fracture work to the zirconium dioxide composite matrix, increasing the contribution of the martensitic transition to the fracture toughness. The results of specific fracture work for zirconium dioxide presented in Figure 5, Figure 6 and Figure 7b, the sharp spikes in values correspond to the formation of cracks within the material and possible dissipation of energy associated with the martensitic phase transition and crack healing. At low loading rates, the contribution of the reinforcement particles to the total work of fracture is insignificant, but there is a slight increase in the strength of the specimens.

4. Discussion

The analysis of the results of the composite modeling showed that the volume content of the alumina particles leads to the formation of an inhomogeneous front and has an effect on the width of the transition. The strain rate in the impact transition region is affected by the value of the loading rate for the composites considered. In this case, a part of the fracture energy will be dissipated in the region of the hardening particles. In this case, the type of fracture is purely brittle. In this case, the hardening particles are completely destroyed within the composite matrix. The formation and nucleation of cracks in the matrix is carried out mainly at the boundaries of the reinforcing particles and between the reinforcing particles. A sharp increase in the specific fracture work is observed, which corresponds to the transition to the martensitic phase. It is noteworthy that the size and shape of the fragments beyond the unloading front does not change significantly as a result of the increase in shock front passage time.

Analysis of the simulated results has shown that the strain rate-dependent strength characteristics of the ceramic composite are strongly nonlinear. It generally corresponds to modelling data obtained on a similar class of brittle materials [41]. Increasing the strain rate leads to an increase in the energy dissipation rate, which is realized by accumulating microdamage at the front of the elastic precursor. Moreover, the destroyed regions continue to have a deterrent value in the development of further damage to the composite matrix for the composites considered above in terms of the front of the elastic precursor.

The overall fracture performance was not significantly affected by the presence of pores in the structure (up to 3%). This is probably due to the fact that the composite matrix acts as a restraint on the crack propagation due to the dilatancy effect inherent in the zirconia [42].

In this case, the fracture work, corresponding to the time of 4 ns, is realized first in the material of aluminum inclusions during the initial stages of the impact front passage through the composite. As the impact front continues to penetrate deeper into the composite structure, the specific fracture work is redistributed to the zirconia composite matrix. This increases the contribution of the martensitic transition to the fracture toughness. The sharp peaks in the values correspond to the formation of cracks in the material and the possible dissipation of energy associated with the martensitic phase transition and the healing of the cracks. At low loading rates, the contribution of the reinforcement particles to the total fracture work is insignificant, but there is a slight increase in specimen strength. In dynamic modeling of brittle media, the results show a qualitative dependence of the type of energy dissipation on the loading rate [43].

As a result, the dispersed Al₂O₃ particles are completely destroyed and clamped by the composite matrix while playing the role of weakly deformable particles during intense deformation in the crack nucleation and propagation zone.

The decrease in crack resistance is due to the intensive nucleation of microcracks, both in the volume of the reinforcing inclusions and in the volume of the nanostructured matrix, according to the model concepts used in this study. The nature of fracture under these conditions is not related to increasing mesoscopic crack size, but rather to merging of multiple submicron cracks formed in the bulk of the composite. Damage occurs both in the nanostructured matrix and in the submicron reinforcing particles. However, a significant portion of the energy is expended on phase transitions. As a result, a strengthening effect occurs. A quantitative comparison of fracture energy and fracture toughness of composites has also been made in some studies [44]. In this paper, we have departed from this formalism and used the concept of specific work of destruction, and the contribution of the work of destruction of each component of the composite has been shown. This approach is original and allows the evaluation of the contribution of each component of the composite to the work of destruction over a wide range of loading rates. This approach will make it possible to predict the dynamic properties of the material at the design stage of ceramic composites. It will also allow the development of "smart" materials for a wide range of applications.

Within the framework of the proposed model, an increase in the fracture toughness of composites achieved by additive technology may occur by reducing the grain size of the matrix below 500 nm. At the same time, the question of the influence of the size and shape of the brittle inclusions needs to be further investigated.