Modelling the Viscoelastic Response of Ceramic Materials by Commercial Finite Elements Codes

A very prominent patented technology permits to obtain bent porcelain stoneware tiles by a proper combination of machining and secondary firing in a kiln. During this firing, the ceramic materials show a viscoelastic behaviour. The viscoelasticity permits the tiles to be bent without further interventions: just using the gravity force in a so-called pyroplastic deformation. Both the viscoelastic response in general and pyroplastic deformation in particular are complex aspects of material behaviour to be modelled with accuracy. In general, the theory of viscoelasticity can be considered extremely large and precise, but its application on real cases is extremely tricky. A time-depending problem, as viscoelasticity naturally is, has to be merged with a temperature-depending situation. It means that, even if the constitutive equations could be set as general formulation, all the fundamental parameters inside these formulas change in function of both entities. Finite Elements codes could help to pass by this impasse, permitting to discretize the variability on temperature. But several cautions have to be taken into account, especially considering that few commercial codes developed proper algorithms. This paper investigates how the viscoelastic response of ceramic materials can be modelled by commercial Finite Elements codes, defining limits and proposing solutions.


INTRODUCTION
This work deals with the definition of a numerical method to study the viscous behaviour of a specific family of ceramic material, the Grès porcelain, during a specific transformation, known as pyroplastic.
A Grès porcelain stoneware tile is a ceramic building material with high bending strength.The French word "grès" means that the ceramic body of the tile is extremely "vitrified", that is to say compact, hence the exceptional great resistance, perfect for an extreme use as in the case of floor and wall tiles.In recent years, it has been increased request of specific bending elements made of grès porcelain.This kind of tiles is conventionally manufactured by iso-static pressing or slip casting and a great effort is required to accurately reproduce the aesthetic design and coloration [1].
A multi-patented technology was developed on the basis of a pyroplastic shaping of porcelain stoneware tiles [2,3].It consists of a second firing at a temperature below that of sintering (e.g.1190-1200 °C), but with a longer firing time (e.g. less than 1100 °C) and slower thermal rates in order to avoid heating and/or cooling ruptures [4,1].Pyroplastic deformation is overguided by means of refractory supports and some thin groves ground in the back of the porcelain stones, Fig. 1.In general terms, pyroplasticity is the tendency of the material to deform under its own weight, due to a low viscosity during firing.The viscosity decreases as the temperature increase in accordance with Arrhenius law.Local rearrangements and microstructural changes affect the mechanical strength and other technological properties of material [5,6,7].Pyroplasticity can be also defined as the loss of shape of the product during its firing.The measure of pyroplasticity is particularly important for bodies which must be completely sintered because during the final stage of the firing they develop an abundant vitreous phase with a sufficiently low viscosity to cause a rapid deformation of the material [5,8].Deformation is related to properties of liquid phases formed during firing [9].It occurs not only because of the flow of this viscous glass at high temperatures, but also because of the applied stress due to the weight of the product during firing [10,11].
A persuasive explanation of what happens in a ceramic material during the complete cycle of firing is provided by [12,13].From room temperature up to 800 •C the material behaviour is characterized by a density decrease due to the weight caused by the residual water evaporation and by the decomposition of organic components.At higher temperature the density increase as consequence of porosity closure.At temperatures of about 900-1000 °C the feldspar content of raw materials melts; this turns into a viscosity decrease that causes the pyroplastic deformation.On the other hand, during the thermal cycle, the sintering process involves an increase of the density and as a consequence an increase of the viscosity.Finally, the deflection stops during the cooling phase, Fig. 2. In this specific case, the pyroplastic deformation under investigation consists in the effect of bending of the grès porcelain stoneware tile caused by gravity during the secondary heat treatment.This situation is probably not so far from what wasexperimented in Fig. 3.

THEORETICAL BACKROUND
Even if in recent investigations the pyroplastic behaviour was found to depend in a complex way on microstructural and compositional features, which affect the effective viscosity [1], typically, ceramic material undergoing pyroplasticity and sintering is treated as a viscoelastic material [13].
Viscoelastic materials are quite more complex to describe than purely elastic materials [14].In fact, and in opposition to elastic materials, viscoelastic materials are characterized by having strain-rate dependence, loading-history dependence and an energy dissipative behaviour due to internal damping mechanisms.Some common phenomena in viscoelastic materials are the creep (strain increases with time for a constant stress), relaxation (stress decreases with time for a constant strain) as detailed in Fig. 4.These intrinsic material characteristics increase the constitutive modelling complexity due to the associated non-conservative effects, strain-rate dependence, loading amplitude dependence and type of loading [15].As a consequence, a very critical aspect is the definition of the material constitutive law able to replicate, in a very accurate way, the physical phenomena involved during sintering.Generally, there are many different possibilities to account for the viscous properties of materials.The simplest solution is represented by the Maxwell model consisting in a purely elastic spring and a purely viscous damper in a series connection.This Hookean spring should be visualized as representing the elastic or energetic component of the response, while the Newtonian dashpot represents the conformational or entropic component [16,17].The stress at each element is the same and equal to the imposed stress, while the total strain is the sum of the strain in each element, representing the elastic and the viscous strain contribution.
The Maxwell model, largely adopted for its ability to simplify, is far away from representing the elastic viscous behaviour of materials.For instance, it is scarcely useful when the elastic strain in negligible respect to viscous strain.In particular, Scherer pointed out the inadequacy of applying viscoelasticity model to sintering materials that often do not show linear viscoelastic behaviour [18].A model based on the purely viscous response of the sintering material was developed thanks to an analogy illustrated in [19].
For an elastic material, the behaviour is described by a constitutive equation, the proportional Hook law, along with the definition of two independent elastic constants.Typically, Poisson's ratio and Young modulus are used in equations, but Shear modulus and Bulk viscosity could also be adopted.In fact, these four constants (respectively E, v, G, K), only depending on material, are correlated by simple numerical expressions.
( ) Since these relations, Eq. ( 1) are obtained by the application of equilibrium conditions, they do not depend on properties of materials and can be used in the generalisation of the theory.
In the case of isotropic, linearly viscous, incompressible material, for instance, the constitutive equation can be obtained invoking the so called "elasticviscous analogy": the strain is replaced by the strain rate (including time), the Young's modulus is replaced by the Shear viscosity and Poisson's ratio assumes the value of ½ (strictly valid for an incompressible material).
The Scherer model represents the other extreme hypothesis, a purely viscous response, used to simplify the behaviour of materials.In most of the cases, the reality is in the middle.Additional physical models have been developed combining springs and dumpers in a crescent complexity as Kelvin-Voigt, Standard linear solid [20,21] and Burgers, Fig. 5. Generalization of Maxwell or Kelvin models to multiple elements is also considered as the most general form of the linear model for viscoelasticity taking into account that physical phenomena (as relaxation) do not occur at a single time, but at a distribution of times.

CONSTITUTIVE LAWS
The consistency of material models on practical cases have been investigated in several articles [12,25,3], also comparing theoretical and numerical predictions with experimental results.A useful generalisation is available on [12] referring to the theoretical formulations summarized in Fig. 6.In particular, the models of Olevsky [22], Riedel and Gasik-Zhang are difficult to implement in numerical commercial codes and dedicated user-subroutines are needed.Adding these models dependent on microscopic material properties, difficult to measure by experimental tests.As detailed by the same authors in [13], the viscoelastic term for Olevsky is a function of the hydrostatic and deviatoric value of the stress acting and it is correlated to microscopic properties as porosity and sintering surface energy.Due to the complexities of investigating the microscopic material properties, together with additional considerations on time calculation and conditions for convergence, these models are not commonly used.On the contrary, Scherer and Maxwell models are already implemented in numerical commercial codes as ANSYS, LS DYNA, ABAQUS [20], MARC [23], and used to replicate the viscoelastic behaviour of ceramic material during the firing process.Between the two models, Maxwell law is largely preferred as evident in [12,13,23,24,25,26,27].On the basis of an accurate experimental characterization, the parameters of viscoelastic constitutive law for materials under investigation can be easily derived.In the case of Scherer or Maxwell models, their parameters are dependent on macroscopic material properties that may be determined by conventional experimental tests.

EXPERIMENTAL TESTS
In general terms, the viscoelastic response of materials can be determined by laboratory tests at in-use conditions.Static or dynamic test methods can be used.Common tests are 1.Uniaxial tension or compression methods, as for elastic materials, but modified so as to enable observation of the time dependency of the material response in terms of Young's modulus [25].2. planar and simple shear method permitting to determine the pure shear modulus via planar shear test 3. creep test measuring the time-dependent strain resulting from the application of a steady uniaxial stress [28, 29] 4. stress relaxation test monitoring the time-dependent stress resulting from the application of a steady strain [24] 5. dynamic loading test [16] in which the stress (or strain) resulting from a sinusoidal strain (or stress) is measured, permitting to investigate the material response at shorter times (as seconds).These tests offer a set of experimental values not easy to be directly applied inside the theoretical models of viscoelasticity [30].But, they can be used to calibrate the numerical calculation.For instance, bending creep tests (e.g.two-or three-point bending configuration and optical methods for contactless measures [31]) can investigate the pyroplasticity measuring the bending rate curve.Then, a Finite Element analysis simulating the same experimental conditions can be adjusted to obtain the same results.Then, the same model is used to estimate other proprieties as the change in viscosity [13].In the same way, dilatometric tests can determine material shrinkage deriving, by simulation, the evolution of density during the sintering cycle [5,29].

DEPENDENCE ON TIME
An interesting method to include experimental data in a viscoelastic constitutive model was mathematically developed by Gaspard Riche de Prony in 1795, but had to wait for digital computers before affirming its importance.Similarly to a Fourier transformation, Prony's method is essentially a decomposition of a intricate curve by building a series of damped complex exponentials.This approximation is pertinent since several properties in viscoelasticity are decaying function of time.
Above all, instead of having constant values for share modulus (G) and bulk modulus (K) and, by extension, for elastic modulus and Poisson's ratio, these properties can be represented by the Prony series: Experimental tests can provide values for these formula, Eq. ( 2), considering that α i represents the amount of (shear or bulk) stiffness lost at time τ i and α ∞ the percentage of remaining stiffness [24].
In a relaxation test, the material is subjected to a constant strain over the duration of the test and the stress is measured over time.The initial stress is due to the elastic response of the material and it relaxes over time due to the viscous effects.The resulting stress vs. time data can be fitted with a number of equations.Only the notation changes depending of the type of strain applied: tensile, compressive, shear or bulk.For each test a Prony series can be determine fitting data with the equation by error minimization algorithms [24,25].
A creep experiment is usually easier to perform than a relaxation one.On the contrary, there is no known closed form for the creep compliance in terms of the coefficient of the Prony series.As a consequence, it is possible, but not easy to get the relaxation coefficients used in Prony series [17].

DEPENDENCE ON TEMPERATURE
In the discussion so far, the effect of temperature on the constitutive equations has not been taken into account.On the contrary, it is evident that viscoelastic materials behave as elastic solids and viscous liquids, depending on the temperature.In general, it is possible to consider this effect assuming that coefficients in equations are simultaneously a function of time and temperature.And, in principle, these unknown functions could be determined by the same tests (e.g.relaxation, creep, etc.), performed at different levels of temperature.This approach becomes very involved especially considering a large number of experiments required.Fortunately, theoretical and experimental results indicate that quite often the effect due to time and temperature can be combined into a single parameter through a "superposition principle".Practically, the formula used in viscoelasticity can be saved if a "reduce time" (τ) is introduced for the correct estimation of temperature depending modula (share, bulk... ) [16,17].While the time (t) inside the constitutive law permits to obtain a time depending solution, the "reduce time" (τ) permits to consider the parameters in function of temperature (and not time).This approach is valid under the assumption that the effect of temperature on time dependent mechanical behaviour is equivalent to a stretching or shrinking of the real time for temperature above or below the reference temperature.For instance, in the case of Maxwell's model, according to the timetemperature superposition, the constitutive law becomes as in Eq. ( 3) where K and G are expressed in function of the reduce time (τ).: Typically, there is a strict connection between the temperature and the time, like in the case under investigation.Fig 7 shows the time-depending temperature inside the kiln during the secondary firing of porcelain stoneware tiles.This time-temperature relation is sometimes considered in simulation by the use of a so called "timetemperature shift function".It allows users to define a single "master" relaxation curve, depending of the temperature, to be shifted by values also depending of temperature (according to the shift function).In other words, the behaviour of the material at high temperature and high strain rate is similar to that at low temperature and low strain rate (also called "thermorheologically simple") [16,17].
Without entering in to further details (available in [24]), the choice of using temperature-dependent constant or a shift functions mainly depends on the availability of experimental data.

RESULTS AND DISCUSSION
Following the previous analysis, several aspects deserve to be mentioned: − even if the linear viscoelasticity is the dominant model in investigation of pyroplasticity, the inelasticity is evident in our problem (the bending process is definitive in tiles) and the effect of a final viscoplastic deformation should be considered (e.g.Norton model as in [25]) − ceramic sintering, as seen in Fig. 2, involves changes of phases and other physical transformations that create a multifaceted concurrence; an unlike circumstance is the use of a single viscoelasticity model for investigating this complex phenomenon − Maxwell model, largely preferred for its general simplicity (including a formulation easy to understand by researchers) and its presence in commercial codes (limiting the efforts for its use), is often adopted without a real correspondence with the behaviour response of materials − temperature has a dramatic influence on the rate of viscoelastic response: in practical work the viscoelastic analysis has to be adjusted for varying temperature or results are useless − the use of Prony series, though representing a valid mathematical approach, is often limited to one or two elements in the series with a drastic reduction in precision − in several practical situations, with the aim of reducing the number of tests, the relaxation bulk modulus (K) is erroneously considered as constant and, sometime, even infinite.This simplification, strictly valid in limited cases (e.g.uncompressible, no hydrostatic deformation..), leads to an estimated value twice the real one [32]) − the hypothesis of "thermorheologically simple" material is often too simplistic for real cases − shift function is relatively easy to formulate, and its values determined by experiments, only when relaxation time obeys specific models (e.g.Arrhenious relation [16]) − apparently similar materials can reveal largely different behaviour in viscoelasticity − a limited number of experiments, and rarely for ceramics, provide material data suitable for the correct extrapolation of curves to be used in simulation

LESSON LEARNT
Even if the industrial problem under investigation in this paper, a ceramic tile that bends under its own weight during a secondary firing, seems, at a glance, very simple, its correct modelling is all but an easy task.A numerical simulation of that process was performed as a test-case, with the aim of validating hypotheses and alternative choices (fig.8).ANSYS Workbench Platform Ver.14.5 was used applying the following conditions [33,34]: 2D model, nonlinear analysis, large deflection [35], material isotropy, elastic stress-strain curve and experimentally determined parameters (at room temperature [8]), "pure" viscoelasticity (TB Visco), Generalized Maxwell constitutive law, shell elements (181/VISCO88), constant bulk modulus (incompressible), time-depending shear modulus, time-dependency defined by two-element Prony series, no shift function, 3 load steps, time steps and substeps carefully optimized.On purpose, settings were chosen considering a mid-level experienced user of simulation code and information available from guidelines.Material data were chosen in accordance with analogue researches or similar materials (as [1,5,8,13]).In this phase, the idea to perform a large experimental campaign was rejected since it is out of the scope.The main objective of the paper is represented by an assessment on the use of commercial FE codes for quick solving practical problems of viscoelasticity.In that sense, a complete characterisation of time and temperature-depending behaviour of materials by complex experiments was considered out of ordinary.The simulative activity highlighted the following criticalities.− unlike other fields where simulation is commonly applied, its use for investigating viscoelastic problems needs a deep familiarity with material science, constitutive models, experimental data evaluation.Without this knowledge, it is extremely difficult to double-check input parameters and results are uncertain − it is also requested by users a discrete consciousness on the applicability of particular algorithms of specific software (up to, for instance, if a particular version of code was able to solve a methodological bug).This information, when available, is difficult to be disclosed, hidden in technical guides and comparative review − viscoelasticity is not properly considered if time step is not small enough during transitions (e.g.temperature changes); this mistake is not easy to detect.

CONCLUSIONS
Mathematical models for predicting the behaviour of ceramic material during firing have been under investigation for more than half a century.Since the very complex process, involving several different mechanisms for material transport, a simple and definite approach is not possible.Sintering is a transient and highly nonlinear phenomena, both in terms of linear and material characteristics.Computer simulations have been developed in parallel with large efforts resulting in huge number of publications.Unfortunately, accuracy in prediction continues to be a real limit and not related to computational limitation.Input data are complex and algorithms not completely optimized.As a consequence, direct experiments continue to represent the best way for investigating the viscoelasticity in ceramics in spite of their complexity and costs.

Figure 4 .
Figure 4. Stress and strain vs. time with evidence of creep and recovery

Figure 7 .
Figure 7. Temperature profile during the bending process