%~Mouliné par MaN_auto v.0.32.4 2024-11-13 10:29:07
\documentclass[CRMECA, Unicode, XML, screen, biblatex,published]{cedram}

\addbibresource{Linn.bib}


\usepackage{changes}
\usepackage{subcaption}


%\newcommand{\bbR}{\mathbb{R}}
%\newcommand{\bbC}{\mathbb{C}}
%\newcommand{\bbZ}{\mathbb{Z}}
\newcommand{\rme}{\mathrm{e}}
\newcommand{\rmm}{\mathrm{mm}}
\newcommand{\rth}{\mathrm{th}}
\newcommand{\rint}{\mathrm{int}}
\newcommand{\rmd}{\mathrm{d}}
\newcommand{\rmI}{\mathrm{I}}
\newcommand{\rmc}{\mathrm{c}}
\newcommand{\rmh}{\mathrm{h}}
\newcommand{\rmT}{\mathrm{T}}
\newcommand{\rmK}{\mathrm{K}}
\newcommand{\rmp}{\mathrm{p}}
\newcommand{\rext}{\mathrm{ext}}
\newcommand{\MPa}{\mathrm{MPa}}
\newcommand{\bvarepsilon}{\boldsymbol{\varepsilon}}
\newcommand{\bsigma}{\boldsymbol{\sigma}}
\newcommand{\bv}{\boldsymbol{v}}
\newcommand{\bb}{\boldsymbol{b}}
\newcommand{\bA}{\boldsymbol{A}}
\newcommand{\bphi}{\boldsymbol{\phi}}
\newcommand{\bB}{\boldsymbol{B}}

\makeatletter
\renewcommand\thesubfigure{(\alph{subfigure})}
\DeclareCaptionLabelFormat{opening}{#2}\subcaptionsetup[figure]{labelformat=opening}
\makeatother

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

%\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\makeatletter
\def\TITREspecial{\relax}\def\cdr@specialtitle@english{Advances in finite fracture mechanics -- A tribute to Dominique Leguillon's scientific achievements}
%\def\cdr@specialtitle@french{}
\makeatother{}


\title{Closed-form analysis of a thermally loaded single-layer system on a rigid foundation}

\alttitle{Analyse de forme fermée d'un système monocouche thermiquement chargé sur une fondation rigide}


\author{\firstname{Debora} \lastname{Linn}\IsCorresp}
\address{Technical University of Darmstadt, Institute of Applied Dynamics, Otto-Berndt-Str.~2, 64287 Darmstadt, Germany}
\email[D. Linn]{linn@ad.tu-darmstadt.de}

\author{\firstname{Wilfried} \lastname{Becker}}
\address[1]{Technical University of Darmstadt, Institute of Applied Dynamics, Otto-Berndt-Str.~2, 64287 Darmstadt, Germany}
\email[W. Becker]{becker@fsm.tu-darmstadt.de}



\keywords{Analytical model, Stress prediction, Interlaminar stresses, Transversal cracks, Crack pattern}

\altkeywords{Modèle analytique, prédiction des contraintes, contraintes interlaminaires, fissures transversales, schéma de fissuration}

\begin{abstract}
A closed-form analytical model of a single layer with a linear elastic material behavior on a rigid foundation subjected to thermal loading is investigated. The closed-form analytical model is based on a higher-order displacement approach that takes the singularity exponent into account. Two applications are considered in this work. First, the interlaminar stresses are analyzed at the interface between the substrate and the material layer. These stresses are an indicator of the formation of interlaminar cracks, which cause the individual layer to peel off. Based on the interlaminar stresses, the closed-form analytical approach, which considers the singularity exponent in the displacement approach, is compared to a model with a second order displacement approach and a FEM model. In the second part, the development of transversal cracks is considered within the framework of Finite Fracture Mechanics, using a coupled stress and energy criterion. Transverse cracks often occur in thin brittle layers such as ceramic coatings and form in periodic patterns. In this work, a representative unit cell is considered -- i.e. the material layer between two cracks. This unit cell is used to determine the cooling temperature at which transverse cracks develop. Furthermore, the resulting distance between two cracks can be determined when larger cooling temperatures are applied.
\end{abstract}

\begin{altabstract}
Un modèle analytique de forme fermée d'une couche unique avec un comportement de matériau élastique linéaire sur une fondation rigide soumise à une charge thermique est étudié. Le modèle analytique de forme fermée est basé sur une approche de déplacement d'ordre supérieur qui prend en compte l'exposant de singularité. Deux applications sont considérées dans ce travail. Premièrement, les contraintes interlaminaires sont analysées à l'interface entre le substrat et la couche de matériau. Ces contraintes sont un indicateur de la formation de fissures interlaminaires, qui entraînent le décollement de la couche individuelle. Sur la base des contraintes interlaminaires, l'approche analytique de forme fermée, qui prend en compte l'exposant de singularité dans l'approche de déplacement, est comparée à un modèle avec une approche de déplacement du second ordre et à un modèle FEM. Dans la deuxième partie, le développement de fissures transversales est étudié dans le cadre de la mécanique de la rupture finie, en utilisant un critère couplé de contrainte et d'énergie. Les fissures transversales se produisent souvent dans des couches fragiles minces telles que les revêtements céramiques et se forment selon des schémas périodiques. Dans ce travail, une cellule unitaire représentative est considérée -- c'est-à-dire la couche de matériau entre deux fissures. Cette cellule unitaire est utilisée pour déterminer la température de refroidissement à laquelle les fissures transversales se développent. En outre, la distance résultante entre deux fissures peut être déterminée lorsque des températures de refroidissement plus élevées sont appliquées.
\end{altabstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\dateposted{2025-01-08}
\begin{document}
\maketitle

\section{Introduction}\label{sec: introduction}

In the field of engineering, there are many applications in which an adhesive bonding is used to join two components. There is great interest in understanding the mechanics behind such adhesive joints in order to design them better and find parameters to improve joint strength. In civil engineering, for example, concrete structures to which a patch of finite size is bonded are analyzed. If concrete structures are damaged due to age or environmental impact, patches can be glued onto the concrete surface to repair or generally reinforce the structure. These patches typically have a high strength and are made of carbon fiber reinforced plastic (CFRP), for example. Frhaan et al. give an overview of the reinforcement of concrete structures with CFRP patches in~\cite{Frhaan_2021}. Various bonding methods of concrete beam structures are examined by Martinelli et al. in~\cite{Martinelli_2023,Martinelli_2023b}. On the one hand, the CFRP reinforcement bars are bonded directly to the concrete surface and on the other hand, grooves are inserted into the concrete before the CFRP bars are applied, providing a form of anchoring. The detachment behavior of these reinforcement patches is investigated experimentally and numerically. Steel reinforcement patches are also conceivable. Awassa et al. tested shafts reinforced with steel rods from the underground core network in a 4-point bending test in~\cite{Awassa_2023}. The typical mode of failure that appears in such a configuration of patch and substrate is the delamination of the patch. Due to the geometry and the different material properties, very high stress concentrations occur locally in the area of the free edges. Calculations of these stresses are essential in order to be able to design such configurations more precisely with regard to their safety.

In the field of structural mechanics, many models for the analysis of bonded joints are presented. A good overview and comparison is given by da Silva et al. in~\cite{daSilva_2009a,daSilva_2009b}. The modeling developed in this paper is similar to that of Methfessel and Becker, which is described in~\cite{Methfessel_2022,Methfessel_2023}. There, a closed-form analytical model is presented, which is based on the generalized model of Bigwood and Crocombe~\cite{Bigwood_1989} and has been extended. In view of a simple and pragmatic model of Methfessel and Becker (\cite{Methfessel_2022,Methfessel_2023}) only the overlap area is considered, in which the internal forces and moments and displacements are applied at the edges, depending on the load case. First-order displacement approaches are used for modeling the adherends. Similar to the work of Ojalvo and Eidinoff in~\cite{Ojalvo_1978}, the displacements in the adhesive layer are composed of the displacements of the adherends, but, unlike in~\cite{Ojalvo_1978}, are extended by higher-order terms. In the current work a higher order displacement approach using the singularity exponent of the system is considered.

The second part of this work focuses on the formation of transverse cracks in thin layers.Thin layers of a solid, brittle material are often applied to protect a component surface. For example, brittle coatings made of a ceramic material on a flexible polymer are used for electronic devices. Another example are automobile catalytic converters, to which $\gamma$-alumina layers are applied as catalytic support. A typical failure mode that occurs in such coatings is the formation of transverse cracks that appear in crack patterns. In some cases, such cracks in thickness direction lead to a loss of stiffness or even delamination of the layer. All in all, this leads to a loss of functionality of the coating. In~\cite{Leguillon_2014}, Leguillon et al. modeled a $\gamma$-alumina layer on a rigid substrate and analyzed the formation of transverse cracks using a coupled stress and energy criterion~\cite{Leguillon_2002}. A coupled stress and energy criterion is also used in~\cite{Leguillon_2013,Leguillon_2017,
Leguillon_2018} by Leguillon et al. to predict crack initiation in thickness direction. There are some more publications that are studying the formation of transversal cracks like~\cite{Bahr_2010}, where Bahr et al. investigated the crack growth or scaling behavior during shrinkage processes, such as cooling or drying, using a fracture mechanics bifurcation analysis with two plausible scaling assumptions. Another investigation that is worth mentioning is~\cite{Jenkins_2009}. There, Jenkins considered a material layer that solidifies starting from a liquid. The solidification process causes the layer to shrink until cracks form in the layer due to excessive internal stresses. The distance between the cracks is then determined using energy minimization. In~\cite{Shao_2010}, Shao et al. carried out experimental tests with ceramic plates and ceramic slabs in which crack patterns form under thermal shock (quenching).


\section{Generic model}\label{sec: model}


\subsection{Mechanical situation}\label{sec: mech model}

Before focusing on a special material and application, this section will begin by considering a generic model where a thin layer of thickness $t$ is perfectly attached to a substrate and subjected to an external thermal loading $\Delta T$. To keep the model simple there are some assumptions made. First, linear elastic, isotropic material behavior of the layer with Young's modulus $E$, Poisson's ratio $\nu$ and thermal expansion coefficient $\alpha_\rth$ is supposed. Second, the substrate is idealized as a rigid body. The model is shown in Figure~\ref{fig:modell_mech}. Due to the fact that geometry and loading do not change in $y$-direction plane-strain behavior is assumed. Thus the displacement $u_y=v$ in $y$-direction vanishes and it is possible to reduce the system to a two-dimensional model with length~$l$. Later, two applications are analyzed where specific parameters are chosen (see Section~\ref{sec:application}).

\begin{figure}[!htbp]
\centering
\includegraphics[width=1\textwidth]{1_Mechanical_Model.pdf}
\vspace*{-20pt}
\caption{Mechanical model of the thin layer on a rigid substrate}\label{fig:modell_mech}
\end{figure}


\subsection{Singularity exponent}

Before the analytical model is derived, a pre-study is performed. Looking more closely at the situation at hand, it can be seen that stress singularities occur at the free edges at the layer substrate interface. In order to understand these singularities more precisely a pre-study is carried out using the method of complex potentials. This method in essence was formulated by Kolosow~\cite{Kolosow_1909}, and was later expanded and widely described by Muschelišvili~\cite{Muschelisvili_1971}. With this method of complex potentials it is possible to determine the singularity exponent. The basic idea is to represent the real field quantities by means of two complex potentials $\Phi$ and $\Psi$ depending on the complex variable $\zeta = x + i z$. The decisive benefit then is the automatic fulfillment of the equilibrium conditions, kinematics and Hooke's law when the stresses and displacements are represented in the following manner (Kolosow's equations):
\begin{align}
\sigma_\varphi + i \tau_{r\varphi} &= \Phi'(\zeta) + \overline{\Phi'(\zeta)} + \zeta \Phi''(\zeta)+ \frac{\zeta}{\overline{\zeta}} \Psi'(\zeta) \ ,\label{eq:kolosov1}
\\
u_r + i u_{r\varphi} &= \frac{\rme^{-i\varphi}}{2 G} \left(\kappa \Phi(\zeta) - \zeta \overline{\Phi'(\zeta)} - \overline{\Psi(\zeta)}\right), \quad \text{with: } \quad \kappa=3-4\nu \ \text{(plane strain)} \ .\label{eq:kolosov2}
\end{align}
The first line contains the relation for the stresses $\sigma_\varphi$ and $\tau_{r\varphi}$, and the second line contains the relation for the displacements. $\zeta$ is the complex coordinate $\zeta=r\rme^{i\varphi}$ (see also Figure~\ref{fig:modell_mech}) and $\Phi(\zeta)$ and $\Psi(\zeta)$ are the complex potentials for which in the present case the following representations are chosen.
\begin{equation}
\Phi(\zeta) = \left(a_1 + i a_2 \right) \zeta^\lambda \ ; \quad \Psi(\zeta) = \left(b_1 + i b_2\right) \zeta^\lambda
\end{equation}
$a_1$, $a_2$, $b_1$ and $b_2$ are real constants and $\lambda$ is the singularity exponent which is to be determined. When inserting the chosen representations into the equations~\eqref{eq:kolosov1} and~\eqref{eq:kolosov2} from the given boundary conditions an algebraic equation system is obtained. Namely for the free edge on the left ($\varphi=\pi/2$ no traction applied) the stresses must be zero and for the bottom edge ($\varphi=0$ rigid substrate) the displacements must be zero. Evaluating these boundary conditions and separating real and imaginary parts leads to a homogeneous system of equations. For a non-trivial solution of this system the determinant of the coefficient matrix has to be equal zero. This results in the characteristic polynomial:
\begin{equation}
-4 \lambda^2 + 2 \kappa \cos(\pi \lambda) + \kappa + 1 = 0 \ .
\end{equation}
From this characteristic polynomial the singularity exponent $\lambda$ can be obtained. For the special case of a $\gamma$-alumina layer with a Poisson's ratio of $\nu=0.2$ for instance a singularity exponent of $\lambda=0.7811$ can be calculated. This means that the displacements locally (for very small radii $r$) behave like $r^{0.7811}$ and the strains and stresses behave like $r^{-0.2789}$. So there are relatively strong stress singularities, which should be taken into account in some way in further investigations.


\subsection{Approximate displacement approach}

For the given mechanical situation now an approximate displacement approach is suggested. Goal is to get an approximate closed-form analytical description of the deformation, strain and stress field. Therefore the displacements $u$ in horizontal and $w$ in vertical direction within the layer are described with a second-order approach extended with an additional term. In this additional term the singularity exponent of the free edges is taken into account. This additional term is used to describe the displacements and stresses more precisely, especially at the ends of the layer, as will be shown later in a comparison of the interlaminar stresses calculated using simpler approaches. With this motivation the horizontal and vertical displacements in the layer are represented as follows
\begin{align}
u(x,z) &= u_1(x) z + u_2(x) z^2 + u_3(x) z^{\lambda} \ ,\label{eq: displacement_approach_u}
\\
w(x,z) &= w_1(x) z + w_2(x) z^2 + w_3(x) z^{\lambda} \ .\label{eq: displacement_approach_w}
\end{align}
Herein $u_1(x)$, $u_2(x)$, $u_3(x)$, $w_1(x)$, $w_2(x)$ and $w_3(x)$ are unknown functions of only the coordinate~$x$ and have to be determined in such a way that representations~\eqref{eq: displacement_approach_u} and~\eqref{eq: displacement_approach_w} are a good approximation of the real deformations. Next step is to determine the unknown displacement functions. For a simpler description, the arguments $x$ and $z$ of the displacement functions are omitted in the following. To determine the unknown displacement functions dependent on $x$ the minimum total energy principle is used
\begin{equation}
\Pi = \Pi_\rint + \Pi_\rext = \min \ ,
\end{equation}
which means that the total energy, consisting of the internal potential $\Pi_{\rint}$ and the external potential~$\Pi_{\rext}$, is required to become minimal. As we only consider thermal and no mechanical loads the external potential is equal zero and the potential energy only consists of the internal energy (which means the total thermally induced elastic strain energy):
\begin{equation}
\begin{split}
\Pi &= \frac{1}{2} \int_{V} \bsigma : \bvarepsilon \ \rmd V \\
&= \frac{1}{2} \int_{0}^{l} \int_{0}^{h} \left[\sigma_{x} \left(\varepsilon_x - \alpha_\rth \Delta T \right) + \sigma_{y} \left(- \alpha_\rth \Delta T\right) + \sigma_{z} \left(\varepsilon_z - \alpha_\rth \Delta T\right) + \tau_{xz} \gamma_{xz}\right] \rmd z \rmd x \ .
\end{split}
\end{equation}
This minimum total energy principle postulates that a system is in equilibrium when the potential energy $\Pi$ becomes minimal.
\begin{equation}
\delta \Pi = 0 \ .
\end{equation}
Using the kinematic relations
\begin{align}
\varepsilon_{x} &= \frac{\partial u}{\partial x} = {u_1}' z + {u_2}' z^2 + {u_3}' z^{\lambda} \ , \\
\varepsilon_{z} &= \frac{\partial w}{\partial z} = w_1 + 2 w_2 z + \lambda w_3 z^{\lambda-1} \ , \\
\gamma_{xz} &= \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} = u_1 + 2 u_2 z + \lambda u_3 z^{\lambda-1} + {w_1}' z + {w_2}' z^2 + {w_3}' z^{\lambda} \ , \\
\varepsilon_{y} &= \gamma_{xy} = \gamma_{yz} = 0 \ ,
\end{align}
Hooke's law
\begin{equation}
\left(
\begin{array}{cccc} \sigma_{x} \\ \sigma_{y} \\ \sigma_{z} \\ \tau_{xz} \\
\end{array} \right) = \frac{\mathit{E}}{2(1+\nu)(1-2 \nu)} \left(
\begin{array}{cccc}
1-\nu & \nu & \nu & 0 \\
\nu & 1-\nu & \nu & 0 \\
\nu & \nu & 1-\nu & 0 \\
0 & 0 & 0 & \frac{1-2\nu}{2}
\end{array}\right)\left(
\begin{array}{rrrc} \varepsilon_{x} - \alpha_\rth \Delta T \\ - \alpha_\rth \Delta T \\ \varepsilon_{z} - \alpha_\rth \Delta T \\ \gamma_{xz} \
\end{array}\right) \ ,
\end{equation}
and applying the variation for all displacement variables $u_1$, $u_2$, $u_3$, $w_1$, $w_2$, $w_3$, we get a non-homogeneous differential equation system of second order (see Appendix~\ref{ap: DGL}). Next step is solving this differential equation system. Therefore we have to reduce it to a system of first order of the following kind:
\begin{equation}
\bA \bphi + \bB \boldsymbol{\dot{\phi \ }} = \bb \ .
\end{equation}
where the vector $\bphi$ contains all displacement functions and their derivatives which are all dependent on $x$.
\begin{equation}
\bphi = \left(u_1, u_2, u_3, w_1, w_2, w_3, u_1', u_2', u_3', w_1', w_2', w_3' \right)^\rmT
\end{equation}
As the system is non-homogeneous, the general solution consists of a homogeneous part $\bphi_\rmh$ and a particular part $\bphi_\rmp$.
\begin{equation}
\bphi = \bphi_\rmh + \bphi_\rmp \ .
\end{equation}
While the particular solution can be obtained using the method of undetermined coefficients,
\begin{equation}
\bphi_\rmp = \bA^{-1} \bb \ ,
\end{equation}
the homogeneous solution is derived from solving an eigenvalue problem. For the homogeneous solution the equation system is transformed to this form,
\begin{equation}
\boldsymbol{\dot{\phi}}_\rmh = -\bB^{-1} \bA \bphi_\rmh \ ,
\end{equation}
where $-\bB^{-1}\bA$ is the new system matrix. By inserting the general exponential representation~$\bv \rme^{\mu x}$ into the system of differential equations an eigenvalue system is obtained. Consequently, the 12 eigenvalues $\mu_{1\dots12}$ can be determined by setting the coefficient determinant equal to zero. Thus eventually with the corresponding eigenvectors $\bv_{1\dots12}$ the general solution can be written as follows
\begin{equation}
\bphi = C_1 \rme^{\mu_1 x} \bv_1 + C_2 \rme^{\mu_2 x} \bv_2 +\dots + C_{12} \rme^{\mu_{12} x} \bv_{12} + \bA^{-1} \bb \ .
\end{equation}
$C_{1\dots12}$ describe the 12 still unknown constants, which in the next step are determined by using the 12 boundary conditions. The following six equations result from the variation approach and are evaluated once at the left edge ($x=0$) and once at the right edge ($x=l$) and in total form a linear system of equations with 12 equations for the 12 unknown constants $C_{1\dots12}$.
\begin{multline}
3\left(1-\nu\right) \left(\frac{2}{3} h^2 u_1^\prime + \frac{1}{2} h^3 u_2^\prime + \frac{2} {\lambda+2} h^{\lambda+1} u_3\prime \right) + 3 h \nu w_1 + 4 h^2 \nu w_2 + \frac{6 \lambda}{\lambda+1} h^{\lambda} \nu w_3 \\
= 3 \left(1+\nu\right) \Delta T \alpha_\rth h 
\end{multline}
\begin{multline}
3 \left(1-\nu\right) \left(\frac{1}{2} h^2 u_1^\prime + \frac{2}{5} h^3 u_2^\prime + \frac{2}{\lambda+3} h^{\lambda+1} u_3^\prime\right) + 2 h \nu w_1 + 3 h^2 \nu w_2 + \frac{6 \lambda }{\lambda+2} h^{\lambda} \nu w_3 \\
= 2 \left(1+\nu\right) \Delta T \alpha_\rth h
\end{multline}
\begin{multline}
\left(1-\nu\right) \left(\frac{2 h^{2}}{\lambda+2} u_1^\prime + \frac{2 h^{3}}{\lambda+3} u_2^\prime + \frac{2 h^{\lambda+1}}{2 \lambda+1} u_3^\prime\right) + \frac{2 h \nu}{\lambda+1} w_1 + \frac{4 h^{2} \nu}{\lambda+2} w_2 + h^{\lambda} \nu w_3 \\
= \frac{2}{\lambda+1}\left(1+\nu\right)\Delta T \alpha_\rth h
\end{multline}
\begin{align}
h u_1 + \frac{4}{3} h^2 u_2 + \frac{2 \lambda}{\lambda+1} h^{\lambda} u_3 + \frac{2}{3} h^2 w_1^\prime + \frac{1}{2} h^3 w_2^\prime + \frac{2}{\lambda+2} h^{\lambda+1} w_3^\prime &= 0 \\
\frac{2}{3} h u_1 + h^2 u_2+\frac{2\lambda}{\lambda+2} h^{\lambda} u_3 + \frac{1}{2} h^2 w_1^\prime + \frac{2}{5} h^3 w_2^\prime + \frac{2}{\lambda+3} h^{\lambda+1} w_3^\prime &= 0 \\
\frac{2}{\lambda+1} h u_1 + \frac{4}{\lambda+2} h^{2} u_2 + h^{\lambda} u_3 + \frac{2}{\lambda+2} h^{2} w_1^\prime + \frac{2}{\lambda+3} h^{3} w_2^\prime + \frac{2}{2\lambda+1} h^{\lambda+1} w_3^\prime &= 0
\end{align}
Solving this system of equations provides the unknown constants, which can be inserted into the displacement functions. The displacement field is thus determined and the stresses are also available via the kinematics and the constitutive law.


\subsection{Finite Element Model}

To compare and validate the analytical solution, a numerical model is created using the finite element method in Simulia Abaqus. The model is discretized in a way that at the lower left region, where the stress singularity is expected, the mesh is the finest (see Figure~\ref{fig:modell_num}). In this area, square elements with an edge length of~$2*10^{-4}\ \rmm$ are chosen, for comparison the layer thickness is~$0.1\ \rmm$. Starting from this area, the mesh becomes coarser towards the other edges. The model uses quadratic basis functions and has about $410,000$ degrees of freedom.


\begin{figure}[!htbp]
\centering
\includegraphics[width=1\textwidth]{2_Finite_Element_Model.pdf}
\caption{Finite Element Model}\label{fig:modell_num}
\end{figure}


\section{Applications}\label{sec:application}


Based on the model for a single layer of material on a rigid substrate under a temperature load, two cases are considered in the following. Firstly, the interlaminar stresses along the interface between the layer and the substrate are analyzed. These stresses are of particular interest as they may lead to the formation of interlaminar cracks and consequently delamination of the layer. In the second part, the development of transverse cracks is considered. In both cases, a thin $\gamma$-alumina layer with Young's modulus~$E=250,000\ \MPa$, Poisson's ratio~$\nu=0.2$ and thermal expansion coefficient~$\alpha_\rth=6*10^{-6}\ \rmK^{-1}$ is considered. The layer is modeled with a thickness of~$t=0.1\ \rmm$ and length~$l=5\ \rmm$. It is perfectly bonded to the rigid base and subjected to an external thermal loading of~$\Delta T=-200\ \rmK$. The model is shown in Figure~\ref{fig:modell_mech} and all parameters are summarized in Table~\ref{tab:parameter}.

\begin{table}
\centering
\caption{Model parameters of the $\gamma$-alumina layer on the substrate}\label{tab:parameter}
\begin{tabular}{lcc}
\hline Parameter & Value & Unit \\
\hline Layer thickness $t$ & $0.1$ & $\rmm$ \\
Length $l$ & $5.0$ & $\rmm$ \\
Elastic modulus $E$ & $250,000$ & $\MPa$ \\
Poisson's ration $\nu$ & $0.2$ & - \\
Thermal expansion coefficient $\alpha_\rth$ & $6*10^{-6}$ & $\rmK^{-1}$ \\
Thermal loading $\Delta T$ & $-200$ & $\rmK$ \\
Singularity exponent $\lambda$ & $0.7811$ & $-$
\end{tabular}
\end{table}


\subsection{Analysis of interlaminar stresses}\label{sec: results}

\begin{figure}[!htbp]
\centering
\includegraphics[width=1\textwidth]{3_Interlaminar_Stresses.pdf}
\vspace*{-30pt}
\caption{Interlaminar stresses $\sigma_z$ and $\tau_{xz}$ along the interface between layer and substrate}\label{fig:interl_stresses}
\end{figure}

The interlaminar stresses along the interface between layer and substrate are now evaluated. Figure~\ref{fig:interl_stresses} on the right side shows the interlaminar peel stress $\sigma_{z}$ and the shear stress $\tau_{xz}$ up to the center of the model, and on the left side the region of the edge where the stress singularity is expected in detail. Now the stress distribution resulting from the analytical model when the singularity exponent is considered can be compared with the numerical model and the analytical model that uses only a second order polynomial approach for the displacements. It can be seen that the expected stress singularity at the free edge is approximated better with the new displacement approach than with the second order approach. Further it shows a better agreement to the calculation with the FEM. Especially the interlaminar peel stress $\sigma_{z}$ shows an almost perfect agreement between the new displacement approach and the FEM.


\subsection{Formation of transversal cracks}\label{sec: transversal cracks}

Beside the delamination of the layer there is another failure mode that can be mentioned. This is the formation of transverse cracks in thin, brittle layers like the $\gamma$-alumina layer. For this purpose, a thin~$\gamma$-alumina layer is first considered as infinitely extended in the longitudinal direction and is loaded with a temperature load~$\Delta T$ (see Figure~\ref{fig:4a}). Now, failure is considered to be the formation of two cracks in the thickness direction, which have a characteristic distance $l_c$ between them (see Figure~\ref{fig:4b}). If just the section of layer between two cracks with the distance~$l_\rmc$ is considered, the generic initial situation is given again (see Figure~\ref{fig:modell_mech}).


\begin{figure}
\centering
%\includegraphics[width=1\textwidth]{4_Model_Transversal_Cracks.pdf}
\subcaptionbox{\label{fig:4a}}{\includegraphics[width=0.45\textwidth]{4_Model_Transversal_Cracks-a.png}}\hfill
\subcaptionbox{\label{fig:4b}}{\includegraphics[width=0.45\textwidth]{4_Model_Transversal_Cracks-b.png}}
\vspace*{-10pt}
\caption{(a): Single layer with infinite length without cracks. (b): Model of the formation of transversal cracks in pattern}\label{fig:transv_cracks_modell}
\end{figure}

Namely, a layer with linear elastic properties on a rigid substrate with length~$l_\rmc$ and thickness~$t$. Consequently, the displacement field is determined using the previously described approach and the stress field within the layer can be derived using kinematic relations and the constitutive law. The failure prediction can be performed with Finite Fracture Mechanics using a coupled stress and energy criterion. Here, a quadratic stress criterion
\begin{equation}
f(\sigma_x,\sigma_{\rmc}) = \left(\frac{\sigma_{x}}{\sigma_\rmc}\right)^2 \geq 1 \ ,
\end{equation}
is used, where the stress~$\sigma_x$ is used at the center of the infinitely extended layer. The present structural situation is modeled in plane elasticity (plane strain), so the stress~$\sigma_x$ can be expressed as follows:
\begin{equation}
\sigma_x = \frac{E}{ (1 - \nu)^2 } \alpha_\rth \Delta T \ .
\end{equation}
According to this term, the stress and consequently the stress criterion depend only on the temperature change~$\Delta T$. Consequently, the evaluation of the stress criterion results in a constant failure temperature change, which is independent of the layer thickness.

For the energy criterion the following linear criterion is used:
\begin{equation}
g(\overline{G}_\rmI, G_\rmc) = \frac{\overline{G}_\rmI}{G_\rmc} \geq 1 \quad \text{with:} \quad \overline{G}_\rmI = \frac{1}{t} \int_{0}^{t} \sigma_x u_x \ \rmd z \ .
\end{equation}
Herein the quantity $\overline{G}_\rmI$ is the so-called incremental energy release rate and the quantity $G_\rmc$ is the fracture toughness of the given material. The coupled criterion postulates the formation of a crack through the total layer thickness~$t$ if both criteria for this length are fulfilled.

The coupled criterion is now evaluated, so that the critical temperature change $\Delta T_\rmc$ can be calculated for different layer thicknesses. This situation has also been investigated by Leguillon~\cite{Leguillon_2014}, so that a comparison of the results can be drawn (see Figure~\ref{fig:transv_cracks} (left)).


\begin{figure}[!htbp]
\centering
\includegraphics[width=1\textwidth]{5_Plot_Transversal_Cracks.pdf}
\vspace*{-20pt}
\caption{Critical temperature reduction (left) and fragment-length (right) for different layer-thicknesses}\label{fig:transv_cracks}
\end{figure} 

For layer thicknesses below~$0.1\ \rmm$, the energy criterion is decisive, that means the critical cooling increases with thinner layer thickness. For layer thicknesses larger than~$0.1\ \rmm$, the stress criterion is governing failure and cracks develop at a cooling of~$179\ \rmK$. For layer thicknesses of exactly~$0.1\ \rmm$, both criteria return a critical load of~$-179\ \rmK$. When comparing the results to the investigations of Leguillon, extremely good agreements can be seen. Both curves match each other almost perfectly, even the transition to a constant failure temperature occurs with both at a layer thickness of~$0.1\ \rmm$.

In the following, the layer is loaded with temperature loads that are~$k$ times larger than the critical temperature load~$\Delta T_\rmc$ so that the resulting fragment length~$l_\rmc$, which is the characteristic distance between two cracks is determined. These fragment lengths are plotted in Figure~\ref{fig:transv_cracks} (right) for different layer thicknesses. A comparison of the curves for different temperature loads shows that smaller fragment lengths occur with higher temperature loads (i.e. with a higher factor~$k$). Qualitatively, the curves are the same for different temperature loads: For layer thicknesses up to~$0.1\ \rmm$, the fragment length increases with increasing layer thickness and at layer thicknesses larger than~$0.1\ \rmm$, the fragment length is decreasing with increasing layer thickness.


\section{Conclusion}\label{sec: conclusion}

In this work, a single layer with elastic material behavior on a rigid substrate has been modeled under temperature loading. An analytical model has been derived to describe the displacement and stress field in an approximate closed-form. Thereby the singularity exponent is taken into account in the displacement approach. One advantage over calculations with the FEM is that parameter studies can be carried out much more quickly and easily. In order to evaluate the accuracy of the closed-form analytical model, a comparison of the interlaminar stresses calculated with different displacement approaches and FEM has been carried out. It has been shown that the results are in better agreement with the approach that takes the singularity exponent into account than with a simple second-order approach. Another important point that should be mentioned is that very high stresses occur at the free edge. These high interlaminar stresses are a cause for the development of interlaminar cracks that lead to a delamination of the layer. An analysis of the development of such interlaminar cracks is a necessary future work and could be performed for example by Finite Fracture Mechanics. The second application investigated in this work is the prediction of transverse cracks, i.e. cracks in the thickness direction that form in a pattern. In this analysis, the coupled criterion has been used to predict the failure temperature change and the fragment length. Although a very simple model is used, the results show good agreement with Leguillon's findings.



\section*{Conflicts of interest}

The authors declare no competing financial interest.

\section*{Dedication}

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

\appendix

\section{Differential equation system}\label{ap: DGL}

\begin{multline}
\left(2\nu-1\right)\left(hu_1+h^2u_2+h^\lambda u_3\right) + \left(2\nu-\frac{1}{2}\right) h^ 2 w_1^\prime + \left(2\nu-\frac{1}{3}\right) h^3 w_2^\prime + \left(2\nu-\frac{1}{\lambda+1}\right) h^{\lambda+1} w_3^\prime \\
-\left(\nu-1\right)\left(\frac{2}{3}h^3u_1^{\prime\prime}+\frac{1}{2}h^4 u_2^{\prime\prime} + \frac{2}{\lambda+2}h^{\lambda+2} u_3^{\prime\prime}\right) = 0\end{multline}
\begin{multline}
\left(2\nu-1\right)\left(h^2u_1+\frac{4}{3}h^3u_2+\frac{2\lambda}{\lambda+1}{h^{\lambda+1}u}_3\right) + 2 \left(\nu-\frac{1}{3}\right) h^3 w_1^\prime+ 2\left(\nu-\frac{1}{4}\right) h^4 w_2^\prime + 2 \left(\nu-\frac{1}{\lambda+2}\right) h^{\lambda+2} w_3^\prime \\
 -\left(\nu-1\right)\left(\frac{1}{2} h^4 u_1^{\prime\prime} + \frac{2}{5} h^5 u_2^{\prime\prime} + \frac{2}{\lambda+3} h^{\lambda+3} u_3^{\prime\prime}\right) = 0
\end{multline}
\begin{multline}
\left(2\nu-1\right)\left(h^\lambda u_1 + \frac{2\lambda}{\lambda+1}h^{\lambda+1} u_2 + \frac{\lambda^2}{2\lambda-1} h^{2\lambda-1} u_3\right) + \left(2\nu-\frac{\lambda}{\lambda+1}\right) h^{\lambda+1} w_1^\prime + \left(2\nu-\frac{\lambda}{\lambda+2}\right) h^{\lambda+2} w_2^\prime \\
+ \left(2\nu-\frac{1}{2}\right) h^{2\lambda} w_3^\prime - 2 \left(\nu-1\right) \left(\frac{1}{\lambda+2} h^{\lambda+2} u_1^{\prime\prime} + \frac{1}{\lambda+3} h^{\lambda+3} u_2^{\prime\prime} + \frac{1}{2\lambda+1} h^{2\lambda+1}u_3^{\prime\prime}\right)=0
\end{multline}
\begin{multline}
 2\left(\nu+1\right)\left(hw_1+h^2w_2+h^\lambda w_3\right) + \left(2\nu-\frac{1}{2}\right) h^2 u_1^\prime+ 2 \left(\nu-\frac{1}{3}\right) h^3 u_2^\prime + \left(2\nu-\frac{1}{\lambda+1}\right) h^{\lambda+1} u_3^\prime\\
+ \left(2\nu-1\right) \left(\frac{1}{3} h^3 w_1^{\prime\prime} + \frac{1}{4} h^4 w_2^{\prime\prime} + \frac{1}{\lambda+2} h^{\lambda+2} w_3^{\prime\prime}\right) = 2\left(\nu+1\right) \Delta T \alpha_\rth
\end{multline}
\begin{multline}
-2\left(\nu-1\right) \left(h^2w_1+\frac{4}{3} h^3 w_2 + \frac{2\lambda}{\lambda+1} h^{\lambda+1} w_3\right) + \left(2\nu-\frac{1}{3}\right) h^3 u_1^\prime + \left(2\nu-\frac{1}{2}\right) h^4 u_2^\prime + \left(2\nu-\frac{1}{\lambda+2}\right) h^{\lambda+2} u_3^\prime \\
+ \left(2\nu-1\right) \left(\frac{1}{4} h^4 w_1^{\prime\prime} + \frac{1}{5} h^5 w_2^{\prime\prime} + \frac{1}{\lambda+3} h^{\lambda+3} w_3^{\prime\prime}\right) = -2\left(\nu-1\right) h^2 \Delta T \alpha_\rth
\end{multline}
\begin{multline}
-2 \left(\nu-1\right)\left(h^\lambda w_1 + \frac{2\lambda}{\lambda+1} h^{\lambda+1} w_2 + \frac{\lambda^2}{2\lambda-1} h^{2\lambda-1} w_3\right) + \left(2\nu-\frac{1}{\lambda+1}\right) h^{\lambda+1} u_1^\prime + 2 \left(\nu-\frac{1}{\lambda+2}\right) h^{\lambda+2} u_2^\prime \\
+ \left(2\nu-\frac{1}{2}\right) h^{2\lambda} u_3^\prime + 2\left(\nu-\frac{1}{2}\right) \left(\frac{1}{\lambda+2} h^{\lambda+2} w_1^{\prime\prime} + \frac{1}{\lambda+3} h^{\lambda+3} w_2^{\prime\prime} + \frac{1}{2\lambda+ 1} h^{2\lambda+1} w_3^{\prime\prime}\right)\\
= 2\left(\nu+1\right) h^\lambda\Delta T \alpha_\rth
\end{multline}



%\nocite{*}


\section*{Declaration of interests}
The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\printbibliography
\end{document}
