A Lamination Model for Pressure-Assisted Sintering of Multilayered Porous Structures

Abstract

This work describes a lamination model for pressure-assisted sintering of thin, multilayered, and porous structures based on the linear viscous constitutive theory of sintering and the classical laminated plate theory of continuum mechanics. A constant out-of-plane normal stress is assumed in the constitutive relation. The lamination relations between the force/moment resultants and the strain/curvature rates are presented. Numerical simulations were performed for a symmetric tri-layer laminate consisting of a 10% gadolinia doped ceria (Ce_0.9Gd_0.1O_1.95-δ) composite structure, where porous layers were adhered to the top and bottom of a denser layer under uniaxially-applied pressures and the sinter forging conditions. The numerical results show that, compared with free sintering, the applied pressure can significantly reduce the sintering time required to achieve given layer thicknesses and porosities. Unlike free sintering, which results in a monotonic decrease of the laminate in-plane dimension, pressure-assisted sintering may produce an in-plane dimension increase or decrease, depending on the applied pressure and sintering time. Finally, the individual layers in the laminate exhibit different stress characteristics under pressure-assisted sintering.

1. Introduction

Ceramic-based multilayers are of importance for a wide range of applications, including solid oxide fuel cells, batteries, capacitors, sensors, thermo-electrics, and windows [1,2,3,4,5]. Integrated printing and sintering is a new processing technology for manufacturing multilayered ceramic materials. This integrated method is a fast and cost-effective way to enhance the functional and mechanical properties of materials compared with traditional material processing approaches, for example, vapor deposition. However, differential densifications in the individual layers during sintering will influence sintering kinetics, and can cause stresses resulting in defects, cracks, and macro-structural distortion. Modelling and understanding sintering kinetics, as well as the shrinkage and distortion, is thus critically important for the efficient manufacturing of multilayers and their applications to provide quality parts with limited defects [6,7,8].

In processing multilayers with powder-based materials, the sintering of these materials is typically done once the material is printed, casted, or pressed [9,10,11]. When sintering thin ceramic material on top of another material of greater or lesser density, or even on a substrate, the sintering behavior of the multilayers is much different from free sintering of a part, because the strain rates, or rate of densification, of the two materials are different, causing one layer to peel off the other or the warping of the multilayer as a whole [12]. The phenomena that causes this distortion is known as constrained sintering [8], and continuum mechanics models have been developed to understand the difference in stress states of the material as a macroscopic viscoelastic body [13,14]. Pressure-assisted sintering has been used to lower the sintering time and temperature, and to minimize grain growth [15]. The reason for these kinetic enhancements is from enhanced diffusion from the creep mechanisms, such as Coble creep and Nabarro-Herring creep, resulting in mechanisms like grain boundary sliding [16,17]. Cramer et al. [18] assessed the constraints and strain rates with pressure-assisted sintering on porous oxide bodies with controlled porosity. Continuum-based models for free and pressure-assisted sintering of multilayered ceramic structures are essential to the design of geometries or configurations in which no warping or distortion occurs, except for the inherent uniform shrinkage.

There has also been some work with 3D printed and tape-casted parts utilizing multilayer and multi-material approaches for making bulk structures [19,20,21,22,23], but a better understanding of the sintering and processing with and without an applied mechanical load is needed to achieve the best densities and properties. Fundamentally, 3D printing of materials needing sintering may have the same constraints as multilayers with different densities because of gaps, defects, and layering effects [22,23], but the focus of this research does not discuss this further.

In modeling shrinkage and distortion behavior of sintered multilayered porous structures, reasonable assumptions on the strain distribution through the thickness can often be used without compromising the accuracy of the models. An assumption of a uniform strain rate through the thickness of each layer was used by Chang et al. [6], Chiang et al. [7], and Ni et al. [24]. The uniform strain rate assumption may be used for a laminate consisting of extremely thin layers or a laminate that is symmetric about its geometrical mid-plane. For general multilayered structures, better accuracy may be achieved using a linear strain distribution assumption, as suggested by Kanters et al. [25], who simulated densification and warpage during sintering of a bi-layered nanocrystalline zirconia structure. Most recently, Molla et al. [26] modeled shrinkage and distortion of a bi-layered structure using a continuum, linear viscous theory of sintering applied to Ce_0.9Gd_0.1O_(1.95−δ) (CGO). Olevsky et al. [27] and Ni et al. [28] developed a constitutive model for free sintering of multilayered porous structures, and applied the model to sintering of bi-layered CGO structures, showing how all of the parameters and state variables effect the stress, porosity, and distortion. They also developed some conditions for tri-layered structures, but the pressure effects were not considered. Besides structural mechanics-based models as reviewed above, Molla et al. [29] studied shape distortion of bilayers with different densities using a multiscale model. Shabana et al. [30] modeled stress evolutions in a metal–ceramic functionally graded material during free sintering using a viscoplastic constitutive model. We point out that pressure-assisted sintering of multilayered structures was not considered in the studies of Olevsky and co-workers [26,27,28,29]. Continuum-based models for pressure-assisted sintering of multilayers have not been available in the literature to the best knowledge of the authors.

The purpose of this work is to present a lamination model for pressure-assisted sintering of thin multilayered porous structures. The model is based on the linear viscous constitutive theory of sintering [31] and a laminated plate theory of continuum mechanics [32]. The remainder of the paper is organized as follows. Section 2 reviews the constitutive relations of a linear viscous theory of sintering [31]. In Section 3, the basic equations for sintering of multilayered structures are derived based on the Kirchhoff hypotheses of classical laminated plate theory, considering the applied out-of-plane pressure. In Section 4, the lamination model is applied to pressure-assisted sintering of a symmetric, tri-layer structure. A sinter forging condition is assumed, i.e., the tri-layer structure is not confined in the in-plane directions. Numerical examples are given in Section 5 for a CGO_P/CGO_D/CGO_P tri-layer laminate, where CGO_P is the porous layer and CGO_D is a denser layer. Finally, concluding remarks and discussions of the limitations of the lamination model are provided in Section 6.

2. Constitutive Equations of Sintering

In this work, we employ the following linear viscous theory of sintering [31]:

$${\sigma }_{ij}=2{\eta }_0\left[\varphi {\dot{\epsilon }}_{ij}+\left(\psi -\frac{1}{3}\varphi \right){\dot{\epsilon }}_{kk}{\delta }_{ij}\right]+P_L{\delta }_{ij}$$

(1)

where σ_ij is the stresses, ${\dot{\epsilon }}_{ij}$ is the strain rates, P_L is the effective sintering stress, δ_ij is the Kronecker delta, η₀ is the shear viscosity of the fully dense material, φ is the normalized shear viscosity, ψ is the normalized bulk viscosity, the indices i, j, and k have the range of 1 to 3 (or x, y, and z), and a repeated index implies summation over the range of the index. According to Olevsky [31], φ, ψ, and P_L are related to the porosity by:

$$\begin{array}{l}\varphi ={(1-\theta )}^2,\\ \psi =\frac{2{(1-\theta )}^3}{3\theta },\\ P_L=\frac{3\alpha }{2G}{(1-\theta )}^2\end{array}$$

(2)

where θ is the porosity of the body being sintered, G is the grain size (average particle radius), and α is the surface energy per unit area (surface energy density). The evolution equation of porosity is based on the continuity equation, and is given by:

$$\frac{\dot{\theta }}{1-\theta }={\dot{\epsilon }}_{kk}$$

(3)

where the superimposed dot denotes the derivative with respect to time.

The total strain rate is the sum of the stress-induced creep strain rate and the free sintering strain rate, i.e.,

$${\dot{\epsilon }}_{ij}={\dot{\epsilon }}_{ij}^c+{\dot{\epsilon }}_f{\delta }_{ij}$$

(4)

where ${\dot{\epsilon }}_{ij}^c$ is the creep strain rates, and ${\dot{\epsilon }}_f$ is the free linear sintering strain rate given by [31].

$${\dot{\epsilon }}_f=-\frac{P_L}{6{\eta }_0\psi }$$

(5)

Substituting Equations (4) and (5) into Equation (1) yields:

$${\sigma }_{ij}=2{\eta }_0\left[\varphi {\dot{\epsilon }}_{ij}^c+\left(\psi -\frac{1}{3}\varphi \right){\dot{\epsilon }}_{kk}^c{\delta }_{ij}\right]$$

(6)

3. A Lamination Model for Pressure-Assisted Sintering of Multilayered Porous Structures

Consider an N-layered porous structure, as shown in Figure 1, where h_i (i = 1, 2, …, N) is the thickness of the ith layer, and H is the total thickness of the laminate, i.e., H = h₁ + h₂ + … + h_N. The rectangular coordinates (x, y, z) are selected so that the x–y plane is the mid-plane of the laminate and z is normal to the laminate plane.

3.1. Stress-Strain Rate Relation

For thin multilayered structures, the normal stress in the z direction may be assumed as the applied pressure (a constant) in pressure-assisted sintering, and the shear stresses associated with the z direction may be assumed as zero, i.e., ${\sigma }_z={\sigma }_z^0$ (${\sigma }_z^0=0$ for free sintering, as in Olevsky et al. [27]), and σ_xz = σ_yz = 0. Now Equation (6) reduces to:

$$\begin{array}{l}{\sigma }_{xx}=\frac{2{\eta }_0\varphi }{\psi +\frac{2}{3}\varphi }\left[\left(2\psi +\frac{1}{3}\varphi \right){\dot{\epsilon }}_{xx}^c+\left(\psi -\frac{1}{3}\varphi \right){\dot{\epsilon }}_{yy}^c\right]+\frac{\psi -\frac{1}{3}\varphi }{\psi +\frac{2}{3}\varphi }{\sigma }_z^0,\\ {\sigma }_{yy}=\frac{2{\eta }_0\varphi }{\psi +\frac{2}{3}\varphi }\left[\left(2\psi +\frac{1}{3}\varphi \right){\dot{\epsilon }}_{yy}^c+\left(\psi -\frac{1}{3}\varphi \right){\dot{\epsilon }}_{xx}^c\right]+\frac{\psi -\frac{1}{3}\varphi }{\psi +\frac{2}{3}\varphi }{\sigma }_z^0,\\ {\sigma }_{xy}=2{\eta }_0\varphi {\dot{\epsilon }}_{xy}^c=2{\eta }_0\varphi {\dot{\epsilon }}_{xy}\end{array}$$

(7)

The inverse form of the above relation is:

$$\begin{array}{l}{\dot{\epsilon }}_{xx}^c=\frac{1}{6{\eta }_0\varphi \psi }\left[\left(2\psi +\frac{1}{3}\varphi \right){\sigma }_{xx}-\left(\psi -\frac{1}{3}\varphi \right){\sigma }_{yy}\right]+\left[\left(\psi -\frac{1}{3}\varphi \right)\frac{1}{3\psi }-\frac{1}{2}\right]\frac{\psi -\frac{1}{3}\varphi }{\psi +\frac{2}{3}\varphi }\frac{{\sigma }_z^0}{\varphi {\eta }_0},\\ {\dot{\epsilon }}_{yy}^c=\frac{1}{6{\eta }_0\varphi \psi }\left[\left(2\psi +\frac{1}{3}\varphi \right){\sigma }_{yy}-\left(\psi -\frac{1}{3}\varphi \right){\sigma }_{xx}\right]+\left[\left(\psi -\frac{1}{3}\varphi \right)\frac{1}{3\psi }-\frac{1}{2}\right]\frac{\psi -\frac{1}{3}\varphi }{\psi +\frac{2}{3}\varphi }\frac{{\sigma }_z^0}{\varphi {\eta }_0},\\ {\dot{\epsilon }}_{xy}^c={\dot{\epsilon }}_{xy}=\frac{{\sigma }_{xy}}{2{\eta }_0\varphi }\end{array}$$

(8)

The normal creep strain rate in the thickness direction of the multilayered structures can be derived as follows:

$$\begin{array}{ll} {\dot{\epsilon }}_{zz}^c& =-\frac{1}{6{\eta }_0\varphi \psi }\left(\psi -\frac{1}{3}\varphi \right)({\sigma }_{xx}+{\sigma }_{yy})+\left[\frac{1}{2}+{\left(\psi -\frac{1}{3}\varphi \right)}^2\frac{1}{3\psi \varphi }\right]\frac{{\sigma }_z^0/{\eta }_0}{\psi +\frac{2}{3}\varphi }\\ & =-\frac{1}{\psi +\frac{2}{3}\varphi }\left(\psi -\frac{1}{3}\varphi \right)({\dot{\epsilon }}_{xx}^c+{\dot{\epsilon }}_{yy}^c)+\frac{{\sigma }_z^0/(2{\eta }_0)}{\psi +\frac{2}{3}\varphi }\end{array}$$

(9)

The relations in Equations (7)–(9) are applicable to individual layers in the laminate, and the layers may have different properties.

3.2. Strain Rates

The Kirchhoff hypotheses of the classical plate theory may be used for sintering of thin multilayered structures [26,27]. Based on the Kirchhoff hypotheses [32], the in-plane strain rates in the laminate can be expressed as:

$$\begin{array}{l}{\dot{\epsilon }}_{xx}={\dot{\epsilon }}_x^0(x,y)+z{\dot{\kappa }}_x^0(x,y),\\ {\dot{\epsilon }}_{yy}={\dot{\epsilon }}_y^0(x,y)+z{\dot{\kappa }}_y^0(x,y),\\ {\dot{\gamma }}_{xy}=2{\dot{\epsilon }}_{xy}(x,y)={\dot{\gamma }}_{xy}^0(x,y)+z{\dot{\kappa }}_{xy}^0(x,y)\end{array}$$

(10)

where ${\dot{\epsilon }}_x^0,{\dot{\epsilon }}_y^0$^, and ${\dot{\gamma }}_{xy}^0$ are the strain rates of the geometrical mid-plane (z = 0), and ${\dot{\kappa }}_x^0,{\dot{\kappa }}_y^0$^, and ${\dot{\kappa }}_{xy}^0$ are the curvature rates of the mid-plane. The creep strain rates are given by:

$$\begin{array}{l}{\dot{\epsilon }}_{xx}^c={\dot{\epsilon }}_x^0+z{\dot{\kappa }}_x^0-{\dot{\epsilon }}_f,\\ {\dot{\epsilon }}_{yy}^c={\dot{\epsilon }}_y^0+z{\dot{\kappa }}_y^0-{\dot{\epsilon }}_f,\\ {\dot{\gamma }}_{xy}^c={\dot{\gamma }}_{xy}={\dot{\gamma }}_{xy}^0+z{\dot{\kappa }}_{xy}^0\end{array}$$

(11)

Substituting the creep strain rates from Equation (11) into Equation (7), we have the following stress expressions:

$$\begin{array}{l}{\sigma }_{xx}=\frac{2{\eta }_0\varphi }{\psi +\frac{2}{3}\varphi }\left[\left(2\psi +\frac{1}{3}\varphi \right)({\dot{\epsilon }}_x^0+z{\dot{\kappa }}_x^0)+\left(\psi -\frac{1}{3}\varphi \right)({\dot{\epsilon }}_y^0+z{\dot{\kappa }}_y^0)-3\psi {\dot{\epsilon }}_f\right]+\frac{\psi -\frac{1}{3}\varphi }{\psi +\frac{2}{3}\varphi }{\sigma }_z^0\\ {\sigma }_{yy}=\frac{2{\eta }_0\varphi }{\psi +\frac{2}{3}\varphi }\left[\left(2\psi +\frac{1}{3}\varphi \right)({\dot{\epsilon }}_y^0+z{\dot{\kappa }}_y^0)+\left(\psi -\frac{1}{3}\varphi \right)({\dot{\epsilon }}_x^0+z{\dot{\kappa }}_x^0)-3\psi {\dot{\epsilon }}_f\right]+\frac{\psi -\frac{1}{3}\varphi }{\psi +\frac{2}{3}\varphi }{\sigma }_z^0\\ {\sigma }_{xy}={\eta }_0\varphi \left[{\dot{\gamma }}_{xy}^0+z{\dot{\kappa }}_{xy}^0\right]\end{array}$$

(12)

The above expressions are valid for stresses in every layer of the laminate.

3.3. Force and Moment Resultants

Define the following force and moment resultants [32] as:

$$\begin{array}{l}{N}_x={\displaystyle \underset{-H/2}{\overset{H/2}{\int }}{\sigma }_{xx}}dz,N_y={\displaystyle \underset{-H/2}{\overset{H/2}{\int }}{\sigma }_{yy}}dz,N_{xy}={\displaystyle \underset{-H/2}{\overset{H/2}{\int }}{\sigma }_{xy}}dz,\\ M_x={\displaystyle \underset{-H/2}{\overset{H/2}{\int }}{\sigma }_{xx}}zdz,M_y={\displaystyle \underset{-H/2}{\overset{H/2}{\int }}{\sigma }_{yy}}zdz,M_{xy}={\displaystyle \underset{-H/2}{\overset{H/2}{\int }}{\sigma }_{xy}}zdz\end{array}$$

(13)

Substituting the stresses in Equation (12) into the above equation, we obtain the following relations between the force/moment resultants and the strain/curvature rates:

$$\begin{array}{l}{N}_x=A_{11}{\dot{\epsilon }}_x^0+A_{12}{\dot{\epsilon }}_y^0+B_{11}{\dot{\kappa }}_x^0+B_{12}{\dot{\kappa }}_y^0-N^f+N_z^0\\ N_y=A_{12}{\dot{\epsilon }}_x^0+A_{11}{\dot{\epsilon }}_y^0+B_{12}{\dot{\kappa }}_x^0+B_{11}{\dot{\kappa }}_y^0-N^f+N_z^0\\ N_{xy}=A_{66}{\dot{\gamma }}_{xy}^0+B_{66}{\dot{\kappa }}_{xy}^0\end{array}$$

(14)

and

$$\begin{array}{l}{M}_x=B_{11}{\dot{\epsilon }}_x^0+B_{12}{\dot{\epsilon }}_y^0+D_{11}{\dot{\kappa }}_x^0+D_{12}{\dot{\kappa }}_y^0-M^f+M_z^0\\ M_y=B_{12}{\dot{\epsilon }}_x^0+B_{11}{\dot{\epsilon }}_y^0+D_{12}{\dot{\kappa }}_x^0+D_{11}{\dot{\kappa }}_y^0-M^f+M_z^0\\ M_{xy}=B_{66}{\dot{\gamma }}_{xy}^0+D_{66}{\dot{\kappa }}_{xy}^0\end{array}$$

(15)

where

$$\begin{array}{l}{A}_{11}={\displaystyle \sum _{i=1}^N}\frac{2{\eta }_{0i}{\varphi }_i}{{\psi }_i+\frac{2}{3}{\varphi }_i}\left(2{\psi }_i+\frac{1}{3}{\varphi }_i\right)(z_i-z_{i-1}),\\ A_{12}={\displaystyle \sum _{i=1}^N}\frac{2{\eta }_{0i}{\varphi }_i}{{\psi }_i+\frac{2}{3}{\varphi }_i}\left({\psi }_i-\frac{1}{3}{\varphi }_i\right)(z_i-z_{i-1}),\\ A_{66}={\displaystyle \sum _{i=1}^N}{\eta }_{0i}{\varphi }_i(z_i-z_{i-1})\end{array}$$

(16a)

$$\begin{array}{l}{B}_{11}={\displaystyle \sum _{i=1}^N}\frac{2{\eta }_{0i}{\varphi }_i}{{\psi }_i+\frac{2}{3}{\varphi }_i}\left(2{\psi }_i+\frac{1}{3}{\varphi }_i\right)\frac{1}{2}(z_i^2-z_{i-1}^2),\\ B_{12}={\displaystyle \sum _{i=1}^N}\frac{2{\eta }_{0i}{\varphi }_i}{{\psi }_i+\frac{2}{3}{\varphi }_i}\left({\psi }_i-\frac{1}{3}{\varphi }_i\right)\frac{1}{2}(z_i^2-z_{i-1}^2),\\ B_{66}={\displaystyle \sum _{i=1}^N}{\eta }_{0i}{\varphi }_i\frac{1}{2}(z_i^2-z_{i-1}^2) \end{array}$$

(16b)

$$\begin{array}{l}{D}_{11}={\displaystyle \sum _{i=1}^N}\frac{2{\eta }_{0i}{\varphi }_i}{{\psi }_i+\frac{2}{3}{\varphi }_i}\left(2{\psi }_i+\frac{1}{3}{\varphi }_i\right)\frac{1}{3}(z_i^3-z_{i-1}^3),\\ D_{12}={\displaystyle \sum _{i=1}^N}\frac{2{\eta }_{0i}{\varphi }_i}{{\psi }_i+\frac{2}{3}{\varphi }_i}\left({\psi }_i-\frac{1}{3}{\varphi }_i\right)\frac{1}{3}(z_i^3-z_{i-1}^3),\\ D_{66}={\displaystyle \sum _{i=1}^N}{\eta }_{0i}{\varphi }_i\frac{1}{3}(z_i^3-z_{i-1}^3) \end{array}$$

(16c)

and

$$\begin{array}{l}{N}^f={\displaystyle \sum _{i=1}^N}\frac{6{\eta }_{0i}{\varphi }_i{\psi }_i{\dot{\epsilon }}_{fi}}{{\psi }_i+\frac{2}{3}{\varphi }_i}(z_i-z_{i-1}),M^f={\displaystyle \sum _{i=1}^N}\frac{6{\eta }_{0i}{\varphi }_i{\psi }_i{\dot{\epsilon }}_{fi}}{{\psi }_i+\frac{2}{3}{\varphi }_i}\frac{1}{2}(z_i^2-z_{i-1}^2),\\ N_z^0={\sigma }_z^0{\displaystyle \sum _{i=1}^N}\frac{{\psi }_i-\frac{1}{3}{\varphi }_i}{{\psi }_i+\frac{2}{3}{\varphi }_i}(z_i-z_{i-1}),M_z^0={\sigma }_z^0{\displaystyle \sum _{i=1}^N}\frac{{\psi }_i-\frac{1}{3}{\varphi }_i}{{\psi }_i+\frac{2}{3}{\varphi }_i}\frac{1}{2}(z_i^2-z_{i-1}^2)\end{array}$$

(17)

In Equations (16a)–(17):

$$z_i=-\frac{H}{2}+{\displaystyle \sum _{j=1}^i{h}_j} (i=1,2,\dots ,N), z_0=-\frac{H}{2}=-\frac{1}{2}{\displaystyle \sum _{j=1}^N{h}_j}$$

(18)

The above relations between the force/moment resultants and the strain/curvature rates are an extension of those introduced in Olevsky et al. [27] for free sintering with iso-stretch in the laminate plane, i.e., ${\dot{\epsilon }}_x={\dot{\epsilon }}_y$. Moreover, the $N_z^0$ and $M_z^0$ terms are zero for free sintering.

4. Pressure-Assisted Sintering of a Symmetric Tri-Layer Structure

In this section, we apply the lamination model presented in Section 3 to pressure-assisted sintering of a symmetric tri-layer structure. A sinter forging condition is assumed, i.e., the tri-layer structure is not confined in the in-plane directions. Free sintering of this laminate was studied by Olevsky et al. [27] and Ni et al. [28], but their stress and stress resultants need to be developed further. The outside layers, i.e., layer #1 and #3, have the same material properties and thickness, i.e., η₀₁ = η₀₃, α₁ = α₃, G₁ = G₃, φ₁ = φ₃, ψ₁ = ψ₃, and h₁ = h₃. Bending and twisting thus do not occur due to the symmetry, i.e., ${\dot{\kappa }}_x^0={\dot{\kappa }}_y^0={\dot{\kappa }}_{xy}^0=0$. Moreover, B₁₁ = B₁₂ = B₆₆ = 0 and $M^f=M_z^0=0$ in Equations (14)–(17). For sinter forging with the applied pressure in the z direction, all stress resultants in Equation (14) are zero. The stresses in individual layers, however, do exist because of the material mismatch between the outside and middle layers. Because no shear deformation exists, the lamination Equation (14) now reduces to:

$$(A_{11}+A_{12}){\dot{\epsilon }}^0=N^f-N_z^0$$

(19)

where ${\dot{\epsilon }}^0$ is the normal strain rate in the x and y directions, i.e.,

$${\dot{\epsilon }}^0={\dot{\epsilon }}_x^0={\dot{\epsilon }}_y^0$$

(20)

and

$$\begin{array}{l}{A}_{11}+A_{12}=2\frac{6{\eta }_{01}{\varphi }_1{\psi }_1}{{\psi }_1+\frac{2}{3}{\varphi }_1}{h}_1+\frac{6{\eta }_{02}{\varphi }_2{\psi }_2}{{\psi }_2+\frac{2}{3}{\varphi }_2}{h}_2,\\ N^f=2\frac{6{\eta }_{01}{\varphi }_1{\psi }_1{\dot{\epsilon }}_{f1}}{{\psi }_1+\frac{2}{3}{\varphi }_1}{h}_1+\frac{6{\eta }_{02}{\varphi }_2{\psi }_2{\dot{\epsilon }}_{f2}}{{\psi }_2+\frac{2}{3}{\varphi }_2}{h}_2,\\ N_z^0={\sigma }_z^0\left(2\frac{{\psi }_1-\frac{1}{3}{\varphi }_1}{{\psi }_1+\frac{2}{3}{\varphi }_1}{h}_1+\frac{{\psi }_2-\frac{1}{3}{\varphi }_2}{{\psi }_2+\frac{2}{3}{\varphi }_2}{h}_2\right)\end{array}$$

(21)

Submitting Equation (21) into Equation (19) yields the in-plane strain rate as follows:

$$\begin{array}{ll}{\dot{\epsilon }}^0& =\frac{2{\varphi }_1{\psi }_1\left({\psi }_2+\frac{2}{3}{\varphi }_2\right){\dot{\epsilon }}_{f1}+({\eta }_{02}/{\eta }_{01}){\varphi }_2{\psi }_2\left({\psi }_1+\frac{2}{3}{\varphi }_1\right){\dot{\epsilon }}_{f2}(h_2/h_1)}{2{\varphi }_1{\psi }_1\left({\psi }_2+\frac{2}{3}{\varphi }_2\right)+({\eta }_{02}/{\eta }_{01}){\varphi }_2{\psi }_2\left({\psi }_1+\frac{2}{3}{\varphi }_1\right)(h_2/h_1)}\\ & -\frac{{\sigma }_z^0}{6{\eta }_{01}}\left(2\frac{{\psi }_1-\frac{1}{3}{\varphi }_1}{{\psi }_1+\frac{2}{3}{\varphi }_1}+\frac{{\psi }_2-\frac{1}{3}{\varphi }_2}{{\psi }_2+\frac{2}{3}{\varphi }_2}\frac{h_2}{h_1}\right)/\left(\frac{2{\varphi }_1{\psi }_1}{{\psi }_1+\frac{2}{3}{\varphi }_1}+\frac{({\eta }_{02}/{\eta }_{01}){\varphi }_2{\psi }_2}{{\psi }_2+\frac{2}{3}{\varphi }_2}\frac{h_2}{h_1}\right)\end{array}$$

(22)

which equals $\dot{L}/L$, i.e., the in-plane dimension kinetics, in which L is an in-plane dimension.

The governing equations for the evolution of porosities in the individual layers can be obtained from Equations (2), (3), (5), (8), (9), (11), and (20) as follows:

$$\begin{array}{l}\frac{{\dot{\theta }}_1}{1-{\theta }_1}=2{\dot{\epsilon }}^0\frac{{\varphi }_1}{{\psi }_1+\frac{2}{3}{\varphi }_1}-\frac{9{\theta }_1}{8(1-{\theta }_1)}\frac{{\alpha }_1}{G_1{\eta }_{01}}+\frac{{\sigma }_z^0/(2{\eta }_{01})}{{\psi }_1+\frac{2}{3}{\varphi }_1},\\ \frac{{\dot{\theta }}_2}{1-{\theta }_2}=2{\dot{\epsilon }}^0\frac{{\varphi }_2}{{\psi }_2+\frac{2}{3}{\varphi }_2}-\frac{9{\theta }_2}{8(1-{\theta }_2)}\frac{{\alpha }_2}{G_2{\eta }_{02}}+\frac{{\sigma }_z^0/(2{\eta }_{02})}{{\psi }_2+\frac{2}{3}{\varphi }_2}\end{array}$$

(23)

The governing equations for the thickness shrinkage rates of the individual layers can be obtained from Equations (4), (9), (11), and (20), with the results:

$$\begin{array}{l}\frac{{\dot{h}}_1}{h_1}={\dot{\epsilon }}_{z1}={\dot{\epsilon }}_{f1}-\frac{1}{{\psi }_1+\frac{2}{3}{\varphi }_1}\left({\psi }_1-\frac{1}{3}{\varphi }_1\right)(2{\dot{\epsilon }}^0-2{\dot{\epsilon }}_{f1})+\frac{{\sigma }_z^0/(2{\eta }_{01})}{{\psi }_1+\frac{2}{3}{\varphi }_1},\\ \frac{{\dot{h}}_2}{h_2}={\dot{\epsilon }}_{z2}={\dot{\epsilon }}_{f2}-\frac{1}{{\psi }_2+\frac{2}{3}{\varphi }_2}\left({\psi }_2-\frac{1}{3}{\varphi }_2\right)(2{\dot{\epsilon }}^0-2{\dot{\epsilon }}_{f2})+\frac{{\sigma }_z^0/(2{\eta }_{02})}{{\psi }_2+\frac{2}{3}{\varphi }_2}\end{array}$$

(24)

Equations (23) and (24) are four simultaneous first-order differential equations to determine the evolutions of thicknesses and porosities h₁, h₂, θ₁, and θ₂ with ${\dot{\epsilon }}^0$ given in Equation (22). Moreover, the free sintering strains ${\dot{\epsilon }}_{fi}$ are given by:

$${\dot{\epsilon }}_{fi} =-\frac{3{\theta }_i}{8(1-{\theta }_i)}\frac{{\alpha }_i}{G_i{\eta }_{0i}},i=1,2$$

(25)

Following Olevsky [31], we use the following dimensionless thicknesses and specific dimensionless sintering time:

$$\begin{array}{l}{h}_1^{*}=\frac{h_1}{h_{01}},h_2^{*}=\frac{h_2}{h_{02}}, \\ {\tau }_s=3{\displaystyle \underset{0}{\overset{t}{\int }}}\frac{{\alpha }_1}{G_1{\eta }_{01}}dt\end{array}$$

(26)

where h₀₁ and h₀₂ are the initial thicknesses. With the above dimensionless quantities, Equations (23) and (24) take the following normalized forms:

$$\begin{array}{l}\frac{d{\theta }_1}{d{\tau }_s}=2\left({\dot{\epsilon }}^0\frac{G_1{\eta }_{01}}{3{\alpha }_1}\right)\frac{{\varphi }_1(1-{\theta }_1)}{{\psi }_1+\frac{2}{3}{\varphi }_1}-\frac{3{\theta }_1}{8}+\frac{(1-{\theta }_1)}{{\psi }_1+\frac{2}{3}{\varphi }_1}\frac{{\sigma }_z^0}{2{\eta }_{01}}\frac{G_1{\eta }_{01}}{3{\alpha }_1},\\ \frac{d{\theta }_2}{d{\tau }_s}=2\left({\dot{\epsilon }}^0\frac{G_1{\eta }_{01}}{3{\alpha }_1}\right)\frac{{\varphi }_2(1-{\theta }_2)}{{\psi }_2+\frac{2}{3}{\varphi }_2}-\frac{3{\theta }_2}{8}\frac{G_1{\eta }_{01}}{G_2{\eta }_{02}}\frac{{\alpha }_2}{{\alpha }_1}+\frac{(1-{\theta }_2)}{{\psi }_2+\frac{2}{3}{\varphi }_2}\frac{{\sigma }_z^0}{2{\eta }_{02}}\frac{G_1{\eta }_{01}}{3{\alpha }_1}\end{array}$$

(27)

$$\begin{array}{ll}\frac{d{h}_1^{*}}{d{\tau }_s}& =-\frac{{\theta }_1}{8(1-{\theta }_1)}{h}_1^{*}-\frac{2h_1^{*}}{{\psi }_1+\frac{2}{3}{\varphi }_1}\left({\psi }_1-\frac{1}{3}{\varphi }_1\right)\left[{\dot{\epsilon }}^0\frac{G_1{\eta }_{01}}{3{\alpha }_1}+\frac{{\theta }_1}{8(1-{\theta }_1)}\right]+\frac{h_1^{*}}{{\psi }_1+\frac{2}{3}{\varphi }_1}\frac{{\sigma }_z^0}{2{\eta }_{01}}\frac{G_1{\eta }_{01}}{3{\alpha }_1},\\ \frac{d{h}_2^{*}}{d{\tau }_s}& =-\frac{{\theta }_2}{8(1-{\theta }_2)}\frac{{\alpha }_2}{{\alpha }_1}\frac{G_1{\eta }_{01}}{G_2{\eta }_{02}}{h}_2^{*}-\frac{2h_2^{*}}{{\psi }_2+\frac{2}{3}{\varphi }_2}\left({\psi }_2-\frac{1}{3}{\varphi }_2\right)\left[{\dot{\epsilon }}^0\frac{G_1{\eta }_{01}}{3{\alpha }_1}+\frac{{\theta }_2}{8(1-{\theta }_2)}\frac{{\alpha }_2}{{\alpha }_1}\frac{G_1{\eta }_{01}}{G_2{\eta }_{02}}\right]\\ & +\frac{h_2^{*}}{{\psi }_2+\frac{2}{3}{\varphi }_2}\frac{{\sigma }_z^0}{2{\eta }_{02}}\frac{G_1{\eta }_{01}}{3{\alpha }_1}\end{array}$$

(28)

The initial conditions now are:

$$\begin{array}{l}{\theta }_1={\theta }_{01},{\tau }_{\mathrm{s}}=0\\ {\theta }_2={\theta }_{02},{\tau }_{\mathrm{s}}=0\\ h_1^{*}=h_2^{*}=1,{\tau }_{\mathrm{s}}=0\end{array}$$

(29)

The normalized stresses in each layer can be obtained from Equations (2), (12), and (20) as follows:

$$\begin{array}{l}\frac{{\sigma }_{x1}}{P_{L1}}=\frac{{\sigma }_{y1}}{P_{L1}}=\frac{6{\varphi }_1{\psi }_1}{{\psi }_1+\frac{2}{3}{\varphi }_1}\left[\frac{G_1{\eta }_{01}}{{\alpha }_1}{\dot{\epsilon }}^0-\frac{G_1{\eta }_{01}}{{\alpha }_1}{\dot{\epsilon }}_{f1}\right]\frac{2}{3{(1-{\theta }_1)}^2}+\frac{{\psi }_1-\frac{1}{3}{\varphi }_1}{{\psi }_1+\frac{2}{3}{\varphi }_1}\frac{2G_1{\sigma }_z^0/{\alpha }_1}{3{(1-{\theta }_1)}^2},\\ \frac{{\sigma }_{x2}}{P_{L1}}=\frac{{\sigma }_{y2}}{P_{L1}}=\frac{2}{3{(1-{\theta }_1)}^2}\frac{6({\eta }_{02}/{\eta }_{01}){\varphi }_2{\psi }_2}{{\psi }_2+\frac{2}{3}{\varphi }_2}\left[\frac{G_1{\eta }_{01}}{{\alpha }_1}{\dot{\epsilon }}^0-\frac{G_1{\eta }_{01}}{{\alpha }_1}{\dot{\epsilon }}_{f2}\right]+\frac{{\psi }_2-\frac{1}{3}{\varphi }_2}{{\psi }_2+\frac{2}{3}{\varphi }_2}\frac{2G_1{\sigma }_z^0/{\alpha }_1}{3{(1-{\theta }_1)}^2}\end{array}$$

(30)

5. Numerical Results

This section presents numerical results of porosity evolution, thickness shrinkage, in-plane dimension change, and stresses for sinter forging of a symmetric CGO_P/CGO_D/CGO_P laminate, where CGO stands for gadolinium 10% doped ceria Ce_0.9Gd_0.1O_1.95−δ. CGO_P and CGO_D stand for a porous and a dense (less porous) CGO, respectively. Ni et al. [28] manufactured a CGO_P/CGO_D/CGO_P laminate under free sintering without applying pressures. The outside layers (i.e., layers 1 and 3, CGO_P) and mid-layer (i.e., layer 2, CGO_D) are the same material under the fully dense condition. Hence, we assume that they have the same material properties, i.e., the same shear viscosity of fully dense material η₀, and the same surface energy per unit area α. Moreover, we assume they also have the same grain size, G. Therefore, the property ratios are η₀₂/η₀₁ = 1, α₂/α₁ = 1, and G₂/G₁ = 1. Following Ni et al. [28], the initial porosities are assumed as θ₀₁ = 0.65 and θ₀₂ = 0.45, respectively, and the initial thicknesses are h₀₁ = h₀₃ = 0.4 mm and h₀₂ = 0.03 mm, respectively. These initial thicknesses correspond to the thicknesses after a debinder cycle before sintering [28]. Grain growth is not considered. The grain sizes are G₁ = 1.0 μm and G₂ = 1.0 μm for the CGO_P and CGO_D, respectively. The surface energies per unit area are α₁ = 1.0 J/m² and α₂ = 1.0 J/m² for the CGO_P and CGO_D, respectively.

Figure 2a,b shows the porosity evolutions in the CGO_P (layers 1 and 3) and CGO_D (layer 2) layers, respectively. Two applied axial pressures are considered, i.e., P = $-{\sigma }_z^0$ = 1 MPa and 5 MPa. The porosities monotonically decrease as time evolves. The porosity is significantly reduced by the applied pressure at a given sintering time. The middle CGO_D layer becomes fully dense (i.e., zero porosity) at specific sintering times of 2.75, 2.1, and 1.05 under free sintering, with 1 MPa of pressure and 5 MPa of pressure, respectively. The porosity of the outside CGO_P layers reduces to 0.15, 0.094, and 0.026 at a specific sintering time of 3.0 under free sintering, with 1 MPa of pressure and 5 MPa of pressure, respectively.

The shrinkages of normalized thicknesses are shown in Figure 3a,b. The thicknesses of both layers decrease monotonically with time, with the rate of reduction gradually leveling off. The normalized thickness of the middle CGO_D layer decreases to about 0.94, 0.658, and 0.277 at a specific sintering time of 3 under free sintering, with 1 MPa of pressure and 5 MPa of pressure, respectively. The corresponding dimensional thicknesses are 0.028 mm, 0.0197, and 0.0083 mm, respectively. At the same time, the normalized thicknesses of the outside CGO_P layers decrease to 0.891, 0.594, and 0.221, respectively, and the corresponding dimensional thicknesses are 0.356 mm, 0.238 mm, and 0.0884 mm, respectively. The values for free sintering compare to the experimental data of 0.02 mm for the CGO_D layer and 0.28 mm for the CGO_P layers, respectively [28].

Figure 4 shows the kinetics of the in-plane dimension of the laminate. All three layers are assumed to deform together in the lamination model. The normalized dimension shrinks under free sintering, and the dimension reduces to about 0.824 at a specific sintering time of 3.0. The test result of Ni et al. [28] is about 0.7. The trends of kinetics of dimensional shrinkages and porosities are consistent with the experimental results of Ni et al. [28]. A number of factors may contribute to the difference between the theoretical predictions of this work and the experimental results of Ni et al. [28]. For example, grain growth may affect porosity and hence geometrical dimensions. Furthermore, the linear viscous sintering theory may not exactly describe the material behavior under sintering. For pressure-assisted sintering, the in-plane dimension monotonically increases with time under a 5 MPa pressure, and reaches 2.64 at a specific sintering time of 3. Under a pressure of 1 MPa, the normalized in-plane dimension initially decreases with time, reaches a minimum of 0.988 at a specific sintering time of 0.3, and then increases with time. The increase in the in-plane dimension under sinter forging is due to the uniaxial pressure in the thickness direction.

Figure 5a,b shows the normalized stresses developed in the CGO_P (layers 1 and 3) and CGO_D (layer 2) layers, respectively. Stresses developed in the individual layers because of a porosity mismatch, which results in differential free sintering strain rates. The CGO_P and CGO_D layers respond differently under free and pressure-assisted sintering. Under free sintering, the in-plane dimension shrinks, as shown in Figure 4. Hence, the outside CGO_P layers are subjected to tensile stress, which decreases with time. Compressive stress is developed in the middle CGO_D layer, and its magnitude also decreases with time. The magnitude of the stress in the CGO_D layer is significantly larger than that in the CGO_P layers, because the CGO_D layer is much thinner than the CGO_P layers as assumed. For sintering under the 1 MPa of uniaxial pressure, the individual layers respond similarly, although, generally, the magnitudes of the stresses become smaller. For sintering under the 5 MPa uniaxial pressure, however, the outside CGO_P layers are mostly subjected to compressive stress, and the middle CGO_D layer is mostly subjected to tensile stress.

A sensitivity analysis is performed to evaluate the dependence of the sintering kinetics on the thickness of the mid CGO_D layer. Figure 6 and Figure 7 show the evolutions of porosity and normalized thickness in the CGO_P (layers 1 and 3) and CGO_D (layer 2) layers, respectively, for a larger initial thickness of the mid layer of h₀₂ = 0.1 mm. Other parameters remain the same as those in Figure 2 and Figure 3. Comparing with the results in Figure 2 and Figure 3, it can be seen that this variation in the mid layer thickness has an insignificant influence on the porosity and normalized thickness evolutions.

6. Concluding Remarks

In this work, we present a continuum lamination model for pressure-assisted sintering of thin, multilayered, and porous structures based on a linear viscous constitutive theory of sintering and the classical lamination theory of plates. The lamination relations satisfied by the strain/curvature rates and force/moment resultants considering the applied out-of-plane pressure are derived. Application of the lamination model to sinter forging of a symmetric, tri-layer CGO_P/CGO_D/CGO_P structure under uniaxial pressure indicates that (i) both the porosities and thicknesses of the individual layers are significantly reduced by the applied pressure at a given sintering time, (ii) whereas the normalized in-plane dimension monotonically decreases with time under free sintering, the in-plane dimension may increase or decrease under sinter forging, depending on the magnitude of the applied pressure and sintering time, and (iii) the individual layers can exhibit complex stress characteristics under pressure-assisted sintering. The applied uniaxial pressure may produce tensile or compressive stress in a given individual layer, depending on the magnitude of the pressure. Finally, for the layered structure modelled here, the mid-layer thickness does not play a significant role in the densification.

This paper only considers sinter forging of multilayered structures under uniaxial pressures. In most pressure-assisted sintering, the in-plane movement of the material is normally constrained, which results in neither a uniaxial pressure nor isostatic pressure condition. The present model will be improved to include the effects of in-plane constraints in a future study. Overall, the lamination model presented in this paper is capable of predicting the effects of applied pressure on sintering kinetics, including evolutions of porosities and thicknesses of the individual layers, and the in-plane dimensions of the structure.