Numerical Aspects of a Continuum Sintering Model Formulated in the Standard Dissipative Framework

: Robust and computationally efﬁcient numeric algorithms are required to simulate the sintering process of complex ceramic components by means of the ﬁnite element method. This work focuses on a thermodynamically consistent sintering model capturing the effects of both, viscosity and elasticity, within the standard dissipative framework. In particular, the temporal integration of the model by means of several implicit ﬁrst and second order accurate one step time integration methods is discussed. It is shown in terms of numerical experiments on the material point level that the ﬁrst order schemes exhibit poor performance when compared to second order schemes. Further numerical experiments indicate that the results translate directly to ﬁnite element simulations.


Introduction
The sintering of ceramics is accompanied by substantial volume reductions and, in some cases, large shape distortions. In order to predict such phenomena and optimize the production process of ceramic components, finite element simulations are in use for at least three decades, see, e.g., Zipse [1], Kraft and Riedel [2] and the references cited therein. These simulation methods have been developed with the primary goal to predict the shrinkage and warpage of ceramic components during sintering, and to adjust the dimensions of the green body accordingly. Though still not widely used in an industrial context, this approach has the potential to reduce the need for the usual time-consuming and expensive iterative experimental trial and error approach in process optimization.
Recently, there is growing interest in the modeling of sintering processes due to the rapid progress of additive as well as multi-material manufacturing methods. In this context, undesired shape distortions as well as the generation of residual stresses (i) due to the presence of inhomogeneous temperature fields at the viscous-elastic transition during cooling and (ii) as a result of the mismatch between material properties in multi-material situations become important. Given the complexity of these phenomena and the geometry of the manufactured components, the computational performance of sintering models becomes relevant; and it is therefore desirable to implement robust and computationally efficient numerical schemes. In this context, it is remarked that the present work focuses on the robustness and efficiency of individual steps of temporal integration schemes. Substantial further reductions of computational cost are possible by additionally utilizing classical adaptive schemes like adaptive time stepping and hp-adaptive finite element methods or model order reduction approaches. An example for the use of the latter in the context of sintering simulations can be found in Sarbandi et al. [3]. A discussion of such adaptive schemes is, however, beyond the scope of this contribution.
A thermodynamically consistent large strain/large deformation sintering model capturing both, the viscous effects near the sintering temperature and the elastic effects at temperatures below has been proposed by Stark and Neumeister [4]. In view of the need for numerically efficient schemes, this model has been formulated within the standard dissipative framework, which facilitates the numerical solution by utilizing the variational structure of the problem [5]. Furthermore, the use of a Lagrangian approach allows to circumvent some of the conceptual difficulties associated with the usual rate-type formulations of sintering models and their implementation into finite element codes. This concerns in particular the consistent incorporation of elasticity into the model, which is non-trivial if rate-type constitutive equations are used as a starting point [6] as well as the difficulty involved in choosing appropriate work-conjugate strain measures and objective stress rates [7]. A drawback of the Lagrangian approach is that the constitutive equations are in a form, which is somewhat less intuitive to interpret than rate-type equations.
The focus of this work is to address some numerical aspects of the model of Stark and Neumeister [4]. In particular, after a brief description of the model and the space-time continuous equations, methods for spatial and temporal discretization are discussed, with a focus on several generic first and second order accurate one step time integration schemes for standard dissipative continua proposed in [5]. It is shown by means of a first set of numerical experiments at the material point level that the theoretical rates of convergence of the temporal discretization are indeed achieved in practice. Furthermore it is found that the second order accurate schemes perform substantially better than the first order schemes at about the same numerical cost. A second set of numerical experiments involving the co-sintering of two different materials indicates that these results translate to finite element simulations.

Sintering Model
Let Ω ⊂ R 3 be the reference configuration of a three-dimensional material body, with X ∈ Ω denoting a material point. Furthermore, let T ⊂ R be the time interval under consideration, with t ∈ T denoting an instant of time. In particular, T = [t 0 , t e ) is assumed, with t 0 and t e being the initial and the final time, respectively. Based on the above definitions, the thermodynamic state variables u : : are introduced. Here, u is the mechanical displacement field, which describes how the placements of material points in the reference configuration are mapped onto their placements in the current configuration. In particular, x(X, t) = X + u(X, t), with x(X, t) being the placement of the material point X at time t. U i is the (symmetric) inelastic stretch tensor, which describes the viscous deformation of the body during sintering, and G represents the average grain size of the ceramic material. These quantities are supplemented by the initial conditions with the constant G ref being the initial grain size. Furthermore, u may be subject to Dirichlet type conditions depending on the problem studied, while it is assumed that no such conditions are imposed on U i and G.
In addition to the variables introduced above, the temperature history in the body needs to be taken into account for the description of sintering processes. In this regard, it is for simplicity assumed that the temperature field is homogeneous in the body and known a priori. As a consequence, the temperature field T = T(X, t) enters the below equations as a parameter.
The sintering model is formulated in the standard dissipative framework [8][9][10][11] and thus based on a Helmholtz free energy function(al), a dissipation function(al) and a power function(al). The Helmholtz free energy is assumed to be of the form where dV denotes the volume element and F = ∇x the deformation gradient (Due to the principle of material frame-indifference, ψ is required to depend on F only through the right Cauchy-Green deformation tensor C = F · F. However, as this work focuses on numerics rather than physics, this aspect is not made explicit in the equations. For a detailed discussion of physical aspects including material frame-indifference and material symmetry, the reader is referred to [4]). The particular form for the Helmholtz free energy density ψ used in this work is described in detail in [4] and briefly summarized in Appendix A.
For the dissipation functional, the generic relation is assumed, where( ) denotes the time derivative. For the local dissipation functionδ, the form proposed in [4] is used, see also Appendix A for a brief summary. The power functional P[u,U i ,Ġ; t] accounts for external loads. It is assumed to be linear in the ratesu,U i ,Ġ and is not further specified at this point. It is noted that it is possible to additionally include a parametric dependency of P on the thermodynamic state variables if needed.
Using the Helmholtz free energy functional, the dissipation functional and the power functional, the space time continuous problem is Find the unknowns (u, H) such that u(X, t 0 ) = 0, and H(X, t 0 ) = H 0 (X), and for all (δu, δḢ) and all t ∈ (T \ t 0 ).
In these equations and in the following, U i and G are jointly denoted by H in order to shorten the equations. Furthermore, the notation is used here and henceforth that the first variation of a functional I[x 1 , . . . , x n ; y 1 , . . . , y n ] is δI[x 1 , . . . , x n , δx 1 , . . . , δx n ; y 1 , . . . , y n ]. It is remarked in this context that δP does not depend onu andḢ since P is linear in the rates of the thermodynamic state variables.
The incremental variational principle (1) implies the mechanical equilibrium equations and a set of local constitutive relations, which are discussed in detail in [4] and summarized in Appendix A. The latter discussion does, however, not involve direct contributions of the power functional to the evolution equations for H as such contributions are deemed unphysical. The only reason for including a dependency of P onḢ into the power functional is to enable for later application of the manufactured solution approach, which is used to test the convergence behavior of the time integration methods by comparison with an analytical solution.

Spatial and Temporal Discretization
The finite element method is used for the spatial discretization of the fields u and H. In particular, u is discretized with standard Lagrange elements of degree k, while discontinuous Lagrange elements of degree k + 1 are used for all components of H. By choosing the support points of the latter elements to coincide with the quadrature points, the degrees of freedom associated with H become local to the quadrature points. Firstly, this leads to a substantial reduction of entries in the sparsity pattern of the finite element matrix and, secondly, it allows for eliminating H at the quadrature point level just as in common plasticity formulations as will be discussed below.
With regard to temporal discretization, the schemes described in [5] are slightly adapted to the current situation. For self-containedness of the present work, these will be briefly described in the following, although the reader is referred to [5] for details.
It is assumed throughout that the time interval [t 0 , t e ] is divided by the discrete time points t n = t 0 + n∆t, where n = 0, 1, . . . N and ∆t = (t e − t 0 )/N. The numerically approximated values of u and H at time t n are denoted by u n and H n , respectively; and the initial values at time t 0 of these quantities are given by u(X, t 0 ) = u 0 (X) = 0 and H(X, t 0 ) = H 0 (X).
The method represents a first order accurate time integrator provided sufficient regularity of the problem. An advantage of this formulation is that the variational structure of the space-time continuous problem is preserved, which makes the method particularly attractive from the mathematical point of view. This concerns, e.g., the fact that the resulting finite element systems are symmetric.

α-Family
The "α-family" introduces a real parameter α ∈ (0, 1], where α → 0 corresponds to a Forward Euler type scheme, α = 1/2 to a Crank-Nicolson-type scheme, and α = 1 to a Backward Euler-type scheme (The α-family closely resembles the well-known θ-family for the solution of ordinary differential equations; and in certain special cases the methods of the α-family are equivalent to those of the θ-family. However, this is not the case in general. Therefore, the term α-family is preferred to distinguish it from the standard θ-family.). Moreover, the following quantities are defined: Then, the time-discrete version of (1) reads Find the unknowns (u n+1 , H n+1 ) such that for all (δu n+1 , δH n+1 ) and all n ∈ {0, 1, . . . N − 1}.
For a sufficiently regular problem, the method is second order accurate for α = 1/2 and first order accurate else. However, the method is non-convergent for α < 1/2 for the problem considered here. Furthermore, the finite element systems are unsymmetrical due to the dependency of∆ on H α n .
is solved. Based on the resulting predicted values, the corrector problem is Find the unknowns (u n+1 , H n+1 ) such that for all (δu n+1 , δH n+1 ) and all n ∈ {0, 1, . . . , For the problem considered here, the method provides with the same theoretical rate of convergence in time as the α-family, while retaining the symmetric, variational structure, so that the finite element systems are symmetric.

Numerical Solution Algorithm and Implementation
In principle, a standard Newton-Raphson method can be used to solve (2)-(5). However, without modifications, this method exhibits poor convergence in some situations for the problem studied here. The reason is that ψ contains a penalty type contribution of the parameter, see also Appendix A. Physically, this term ensures that no further densification is possible as soon as all pores of the material all closed. Numerically, the term has the effect that the step, which can be taken during a Newton-Raphson iteration, is often limited by a few quadrature points for which the relative density ρ rel = ρ rel 0 /J i approaches 1; and, hence, a large number of iterations is necessary until the termination criterion is met or the algorithm does not converge at all. There are several remedies to this situation. One possibility is to adapt classical interior point algorithms with a logarithmic barrier function and drive µ globally to the desired small value in several steps. However, this substantially increases the number of assemblies and solutions of the global finite element system. Therefore, the approach used in this work utilizes the fact that the degrees of freedom related to H are local to each quadrature point as indicated above and can, therefore, be eliminated at the quadrature point level by using a second local Newton-Raphson cycle. This approach is essentially the same as in classical plasticity and involves the local computation of consistent tangents. This procedure does of course exhibit the same problem type at the quadrature point level. However, due to the small number of unknowns, the local Newton-Raphson cycle is still well-behaved for reasonably small values of µ.
The actual implementation is based on the libraries GalerkinTools and IncrementalFE, which leverage the capabilities of the open source finite element code deal.II [12,13] for the case of standard dissipative formulations, see also [14]. In particular, the library GalerkinTools provides with a "high-level" algorithm, which implements the local elimination procedure just described in a generic way without reference to the actual relations used for the functionals. With these libraries, the only problem-specific part of the implementation are the functions ψ(F, H; T) andδ(Ḣ; H, T) together with their gradients and Hessians; and the implementation of the gradients and Hessians has been carefully checked by comparison with results obtained from numerical differentiation. For further details, the reader is referred to the documentation of GalerkinTools and In-crementalFE and, in particular, to the documentation of the C++ class ScalarFunctionalLo-calElimination of GalerkinTools (The source codes for the libraries GalerkinTools and Incre-mentalFE are available together with corresponding documentations from the repositories https://github.com/sebastian-stark/GalerkinTools, accessed on 16 May 2023, and https://github.com/sebastian-stark/IncrementalFE, accessed on 16 May 2023, respectively).

Normalization and Initial Conditions
For the numerical experiments discussed below, the reference values T * = 1500 K, G * = 0.1 µm, t * = 1 h, σ * = 10 MPa and L * = 1 mm are used for normalization of the equations. With this choice, the normalized quantities arẽ Concerning the required initial condition for the grain size, G ref = 0.17 µm is chosen.

Numerical Experiments
In this section, two examples are considered, which demonstrate the performance of the time integration schemes. While the first example corresponds to free sintering of a homogeneous block of material, the second example addresses the more complex situation of co-sintering of a disc consisting of two different materials.

General Approach
The same general approach is used to discuss both examples. In particular, the temperature history is homogeneously prescribed according to which represents a linear ramp up between the temperatures T min and T max for a duration of t e /2 followed by a linear ramp down between T max and T min for the same duration.
In particular, T min = 293.15 K, T max = 1623.15 K, and t e = 5 h are chosen. As a starting point, numerical experiments with the choice P = 0 for the power functional have been conducted for both examples in order to determine a well-converged "reference solution". For the first example, this concerns only the temporal discretization. In contrast, for the second example, both the temporal as well as the spatial consideration need to be considered. With regard to the temporal discretization, it has empirically been found that the α-family with α = 1/2 and 1536 equally spaced time increments for the entire time interval between t = 0 and t = t e results in highly accurate numerical solutions for both examples. The spatial discretization used to obtain a reference solution for the second example will be discussed later.
Based on the reference solutions, manufactured analytical solutions to a time-continuous (but, in the case of the second example, spatially discrete) problem have been obtained for both examples by the procedure described in the following. For this purpose, let q n,i denote the finite element solution for the i-th finite element degree of freedom of the reference calculation at time instant t n . Then, the data points (t n , q n,i ) have been used to obtain a time-continuous function q a i (t) for each finite element degree of freedom by cubic spline interpolation (for the actual implementation of the spline interpolation, the cubic spline interpolation library of Kluge [15] has been used). By combining the resulting timecontinuous functions with the finite element interpolation, the analytical solution u a (X, t), H a (X, t) is obtained. The corresponding power functional consistent with this analytical solution follows from (1) and reads P[u,Ḣ; t] = δΨ[u a (X, t), H a (X, t),u,Ḣ; t] + δ∆[Ḣ a (X, t),Ḣ; H a (X, t), t].
This power functional has then been used for all subsequent calculations. It is remarked that this manufactured solution approach has the advantage to allow for a comparison of the numerical results with an analytical solution while ensuring that an example is considered, which is representative of a realistic physical situation. In particular, the power P associated with the analytical solution remains very close to zero throughout the entire time interval under consideration since the reference solution is well converged for the choice P = 0, provided that there are no errors in the implementation. Such errors would, however, most likely become evident in subsequent convergence studies using the manufactured solution.
Based on the manufactured solutions, convergence studies have been performed for both examples using the power functional according to Equation (6). The starting point for these calculations is to take a single time increment for the time range 0 ≤ t < 1/3t e . Subsequently, the time range 1/3t e < t ≤ 2/3t e is split into 128 time increments, followed by a single time increment again between 2/3t e < t ≤ t e . This approach is suitable because the elastic behavior is dominant below temperatures of approximately T min + 2/3(T max − T min ), such that long time steps are permissible. In contrast, smaller time steps are needed in the viscous range at higher temperatures due to rapid sintering. It is remarked in this context that the number of time increments required to obtain reasonable results in the latter range could be substantially reduced by adaptive time stepping. This would, however, complicate the discussion and has, therefore, not been implemented. Based on the computation with 130 time increments in total, the number of time increments has been increased to N t = 130 · 2 m t by uniform refinement, where m t is the refinement cycle. In particular, calculations with m t = 0, 1, . . . , m t,max have been performed, with m t,max = 10 for the variationally consistent method as well as the α-family and the modified α-family with α = 1; and with m t,max = 8 for the α-family and the modified α-family with α = 1/2. With further temporal refinement, no further increase of the accuracy of the solution has been observed with the methods with α = 1/2, which is most likely caused by the limited numerical accuracy related to the use of double precision numbers. In order to actually assess the convergence behavior in time by comparison with the manufactured analytical solution, the error e (m t ) is, for each m t = 0, 1, . . . , m t,max , evaluated according to withF a being the deformation gradient corresponding toũ a , the subscript N referring to the solution after the last time increment (i.e., at t = t N = t e ), and the bracketed superscript referring to the number of refinements in time.

Numerical Experiment 1
In the first numerical experiment, the free sintering behavior (i.e., without applied tractions and displacement constraints other than those minimally needed to remove rigid body motions) of a homogeneous body is considered. For this situation, the shape of the body is not of interest and, therefore, the situation can be reduced to a single material point. In practice, this has been implemented using a single finite element having unit volume. Figure 1 shows the solution obtained with the α-family and α = 1/2 with m t = 3 refinements in time. It is clearly visible, that the sintering happens rapidly in the time range 2 <t < 2.5, and this sintering is accompanied by rapid grain growth. Furthermore, it is remarked that the differences in the total volume change as described byJ = det(F) and the inelastic volume change as described byJ i = det(Ũ i ) are attributable to thermal expansion.

Results
The temporal convergence behavior is discussed in terms of the results shown in Figure 2. In particular, Figure 2a shows the errors e (m t ) for the different cases considered. It can be seen that all results are consistent with the expected rates of convergence, i.e., the variationally consistent method and the other methods with α = 1 exhibit a rate of convergence of k t = 1, while the methods with α = 1/2 are associated with a rate of convergence of k t = 2. While the variationally consistent method is associated with the highest accuracy at a given time step size among the first order methods, the α-family produces slightly more accurate results than the modified α-family for the second order methods with α = 1/2. Figure 2b is used to further discuss the absolute error levels. For this purpose, the transient grain size solution is compared for the time interval 2.5 ≤t ≤ 3 for selected cases. In particular, the α-family is considered for α = 1 and α = 1/2; and for both cases the solutions obtained for m t = 0 and m t = 3 are compared. Hardly a difference can be noticed between the solutions for m t = 0 and m t = 3 for α = 1/2. I.e., for α = 1/2, the solution can be considered converged for m t = 0 already. This contrasts with the results obtained for α = 1. In this case, substantial error is involved for m t = 0, and even for m t = 3 there is still a slight difference to the corresponding solution obtained with α = 1/2, which may be taken as the converged reference solution. Taking into account that the computation time needed for a single time increment is virtually the same for all values of α, it becomes clear that the α-family performs substantially better with α = 1/2 than with α = 1. In practice, this translates to savings of computational time of more than one order of magnitude for a pre-defined error level if α = 1/2 is used. For the case of the modified α-family with α = 1/2, the advantage is nominally slightly less because the corrector step requires additional Newton-Raphson iterations. However, depending on the linear solver, this disadvantage can be compensated for by utilizing the fact that the linear systems associated with the modified α-family are, in contrast to the α-family, symmetric. In addition, positive-definiteness of the system can usually be expected, although this property cannot be guaranteed in the large deformation context. Concerning the first order methods, the variationally consistent method and the modified α-family with α = 1 have similar advantages over the α-family with α = 1 in that the linear systems are symmetric and usually positive definite.

Numerical Experiment 2
In the second numerical experiment, the co-sintering of two layers with different material properties is considered. To simplify matters, an axisymmetric situation is assumed as illustrated in Figure 3. The geometry is a circular cylinder with radius R a = 20 mm and total thickness d = 1 mm. The spatial region occupied by the cylinder is described by 0 < R ≤ R a and −d/2 ≤ Z ≤ d/2 with regard to the radial coordinate R and the axial coordinate Z. The computational domain Ω is split into the two subdomains Ω 1 and Ω 2 , with Ω 1 being associated with Z ≤ 0 and Ω 2 with Z > 0. In Ω 2 , the same constitutive relations are used as in the first numerical experiment, while the viscosity tensor η is scaled by a factor of 2 in Ω 1 compared to the values given in Appendix A (i.e., the quantity η 0 introduced in Appendix A is multiplied by 2). This causes substantially slower sintering in Ω 1 , thus leading to bending of the cylinder towards the positive Z-direction. In order to eliminate the axial rigid body translation, which is still possible in the axisymmetric setting, the point R = 0, Z = 0 is fixed in the axial direction. Furthermore, the problem described in Sections 2 and 3 for the fully three-dimensional case is appropriately modified to account for the axisymmetric situation. As a consequence, two-dimensional finite elements can be used for spatial discretization, which facilitates the convergence studies by reducing computational times and, therefore, allows for a higher number of refinements of the discretization. For the calculations shown below, k = 2 is generally chosen. I.e., quadratic Lagrange finite elements are used to discretize u. These perform much better than their linear counterparts for the bending situation considered here.

Finite Element Grid
In order to obtain a high quality finite element grid, several grid sizes have been attempted for the reference calculation. Based on a comparison of the solutions, the grid shown in Figure 4 has finally been used to determine the manufactured solution and for all subsequent calculations. The figure also includes the outlines of the deformed shape of the disc at the end of the simulation for this grid as well as a grid obtained after three uniform global grid refinement steps. It can be seen that the deformed shapes are indistinguishable and, therefore, the unrefined grid is considered sufficiently accurate.

Results
In the same way as for the numerical experiment 1, the errors e (m t ) are shown in Figure 5. A comparison of the latter figure with Figure 2a shows that the results are very similar apart from the absolute error level and the fact that now all first order methods exhibit virtually the same convergence behavior. Consequently, the same conclusions can be drawn for the component level as for the material point level. Variationally consistent α-family, α = 1 modified α-family, α = 1 α-family, α = 1/2 modified α-family, α = 1/2 Figure 5. Errors e (m t ) for different methods for numerical experiment 2.

Concluding Remarks
In this work, a previously published model for the sintering of ceramic materials has been investigated with regard to temporal discretization. Several first and second order accurate one-step schemes have been considered at the material point level and at the component level. It has been shown by numerical experiments that the theoretically expected rates of convergence are achieved in practice as well. The most important conclusion from this work is that the second order accurate methods discussed in this work are superior to the first order accurate methods. As a consequence, substantial savings of computational time are possible with these schemes; and their use is recommended whenever possible. Aspects not discussed in this work are the application of the methods to situations with thermomechanical coupling, spatial discretization, and adaptive time stepping. All of these are of substantial practical importance and need to be addressed in the future.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A. Constitutive Equations of the Sintering Model
In the following, the particular forms for the Helmholtz free energy density ψ and the local dissipation functionδ as well as the resulting local constitutive equations are briefly summarized. It is remarked that the equations stated below are representative of a particular alumina ceramic material; and details of the constitutive equations are discussed in [4].
Appendix A.1. Helmholtz Free Energy Density ψ The Helmholtz free energy density is assumed to be split according to where C = F · F is the right Cauchy-Green deformation tensor, J i = det(U i ), ψ el represents an elastic contribution, and ψ if represents a microstructural contribution, which constitutes the main driving force for sintering.
The form of the elastic contribution is motivated by an assumed multiplicative split of the deformation gradient F into a thermoelastic part F te and a rotation-free inelastic part U i according to F = F te · U i . F te may be further split into a rotation R te and a thermoelastic stretch U te according to F te = R te · U te , where R te may represent large rotations. In contrast, the thermal expansion strains and the elastic strains are, in the context of sintering, usually "small", such that U te ≈ I, with I being the unit tensor. Hence, the elastic contribution to the Helmholtz free energy density can be based on the usual small strain expression. In particular, where C el (J i ; T) is the elasticity tensor, ε t (J i ; T) the thermal expansion tensor, and E te = 1 2 F te · F te − I = 1 2 U i −1 · C · U i −1 − I is the thermoelastic Lagrange-Green strain tensor (which is small compared to unity). The elasticity tensor and the thermal expansion tensor are assumed to be isotropic. I.e., where Y = Y(J i ; T) is Young's modulus, ν = ν(J i ; T) is Poisson's ratio, ε t = ε t (J i ; T) describes the isotropic thermal expansion, and the tensor product A = a sym b corresponds to the coordinate form A klmn = (a km b ln + a lm b kn )/2.
With regard to the microstructural contribution, the relation The particular functional forms chosen for Y, ν, ε t and ψ if,0 read In these equations, 1 + e −19(ρ rel −0.757) − 1 1 + e −4.617 , and ρ rel = ρ rel (J i ) = ρ rel,0 J i is the relative density, with ρ rel,0 = 0.584 being the initial relative density. It is noted that the logarithmic penalty term −µ ln 1 − ρ rel in ψ if,0 ensures that ρ rel < 1. The penalty parameter µ > 0 must be chosen as small as possible in order to ensure that the term affects the constitutive behavior only for values of ρ rel very close to 1, and as large as necessary in order to avoid pathological numerical behavior. For the value of µ chosen here, the impact of the penalty term on the sintering stress remains below 1% for ρ rel 0.99 [4].
being the inelastic stretching tensor, η the viscosity tensor, and ξ a "micro dissipation coefficient". The viscosity tensor is assumed to be isotropic. In particular, thatψ(E te , J i , G; T) = ψ(F, U i , G; T). It is pointed out that the sintering stress also involves contributions related to the change of the elastic material behavior as the relative density of the material evolves, which naturally result from the inherent thermodynamic consistency of the approach. However, these are typically insignificant compared to the contributions related to the changes of interfacial energy density. It is evident that (A2) has the form of the usual relation between stress and strain rate applicable for purely viscous materials. Furthermore, the evolution Equation (A3) for the grain size suggests a cubic grain growth behavior at fixed temperature and relative density, i.e., G 3 − G * 3 = K(t − t * ), where G and G * denote the grain size at time t and t * , respectively, and K is a constant depending on temperature and relative density.