%% 
%% Copyright 2019-2021 Elsevier Ltd
%% 
%% This file is part of the 'CAS Bundle'.
%% --------------------------------------
%% 
%% It may be distributed under the conditions of the LaTeX Project Public
%% License, either version 1.2 of this license or (at your option) any
%% later version.  The latest version of this license is in
%%    http://www.latex-project.org/lppl.txt
%% and version 1.2 or later is part of all distributions of LaTeX
%% version 1999/12/01 or later.
%% 
%% The list of all files belonging to the 'CAS Bundle' is
%% given in the file `manifest.txt'.
%% 
%% Template article for cas-sc documentclass for 
%% single column output.

\documentclass[a4paper,fleqn]{cas-sc}

% If the frontmatter runs over more than one page
% use the longmktitle option.

%\documentclass[a4paper,fleqn,longmktitle]{cas-sc}

\usepackage[numbers]{natbib}
%\usepackage[authoryear]{natbib}
%\usepackage[authoryear,longnamesfirst]{natbib}

%%%Author macros
\def\tsc#1{\csdef{#1}{\textsc{\lowercase{#1}}\xspace}}
\tsc{WGM}
\tsc{QE}
%%%

% Uncomment and use as if needed
%\newtheorem{theorem}{Theorem}
%\newtheorem{lemma}[theorem]{Lemma}
%\newdefinition{rmk}{Remark}
%\newproof{pf}{Proof}
%\newproof{pot}{Proof of Theorem \ref{thm}}

\begin{document}
\let\WriteBookmarks\relax
\def\floatpagepagefraction{1}
\def\textpagefraction{.001}

% Short title
\shorttitle{Model-assisted guided-wave based approach for debonding assessment in honeycomb sandwich structures}

% Short author
\shortauthors{P. Fiborek, P. Kudela}  

% Main title of the paper
\title [mode = title]{Model-assisted guided-wave based approach for debonding assessment in honeycomb sandwich structures}

% Title footnote mark
% eg: \tnotemark[1]
%\tnotemark[1] 

% Title footnote 1.
% eg: \tnotetext[1]{Title footnote text}
%\tnotetext[1]{<footnote text>} 

% First author
%
% Options: Use if required
% eg: \author[1,3]{Author Name}[type=editor,
%       style=chinese,
%       auid=000,
%       bioid=1,
%       prefix=Sir,
%       orcid=0000-0000-0000-0000,
%       facebook=<facebook id>,
%       twitter=<twitter id>,
%       linkedin=<linkedin id>,
%       gplus=<gplus id>]

\author[1]{Piotr Fiborek}[type=editor,
orcid=0000-0002-5030-3312]

% Corresponding author indication
\cormark[1]

% Footnote of the first author
%\fnmark[1]

% Email id of the first author
\ead{pfiborek@imp.gda.pl}

% URL of the first author
\ead[url]{www.imp.gda.pl}

% Credit authorship
% eg: \credit{Conceptualization of this study, Methodology, Software}
%\credit{<Credit authorship details>}

% Address/affiliation
\affiliation[1]{organization={Institute of Fluid Flow Machinery, Polish Academy of Sciences},
            addressline={Fishera 14 st.}, 
            city={Gdańsk},
%          citysep={}, % Uncomment if no comma needed between city and postcode
            postcode={80-231},
            country={Poland}}

\author[1]{Pawel Kudela}[orcid=0000-0002-5130-6443]

% For a title note without a number/mark
%\nonumnote{}

\begin{abstract}[S U M M A R Y]
	This study aims to assess the debonding of the carbon skin from the honeycomb core in the sandwich panel based on guided wave propagation.
	The effect of damage on wave propagation in the slab was determined using numerical simulations.
	The two most common honeycomb core models found in the literature are a homogenized, orthogonal continuous medium and a real geometry core modelled by finite elements, usually using commercial software.
	In a simplified model, the results are in good agreement only for low frequencies.
	In this paper, the panel has been modelled by the time-domain spectral element method taking into account the real geometry of the honeycomb core.
	
	The presented model was compared with the homogenized model of the 
	The result of the parametric study is a function of damage influence on the amplitude of propagating waves.
\end{abstract}
\begin{keywords}
Honeycomb, Debonding, Computational modelling, Spectral element methods, 
\end{keywords}

\maketitle

\section{Introduction}
\label{sec:intro}
Honeycomb Sandwich Composite~(HSC) is a multi-layered structure composed of the mid-core with the geometry of honeycomb sandwiched between thin skins.
They are widely used in the aerospace, marine and automotive industries due to their high strength-to-weight ratio, high energy absorption capability and effective acoustic insulation.
However, these complex structures are exposed to various types of damage that are not found in metal alloy materials, e.g. hidden disbonds between the skin and the core, delamination of the skin plates, or the core impact damage.
The damage can occur during a manufacturing process, storage or in-service life; thus, advanced methods for online damage detection are required.

The Guided Wave (GW) based method is one of the structural health monitoring (SHM) techniques with very high potential for damage detection and localization in HSC. \cite{mustapha2011assessment, sikdar2016guided, sikdar2016ultrasonic,radzienski2016assessment, yu2019core}.
Among numerous techniques developed for damage detection and localization, the most popular are pitch-catch \cite{sikdar2017structural}, pulse-echo \cite{}, phase array \cite{lu2006crack, ostachowicz2008elastic}, and time-reversal mirror \cite{fink1992time, eremin2016analytically}.


mustapha2006debonding
mustapha2011assessment,
mustapha2014leaky,
baid2008detection
baid2015dispersion,
sikdar2016guided,
sikdar2016ultrasonic,
radzienski2016assessment,
sikdar2017structural,
yu2019core,
wang2020improved


qi2008ultrasonic,
hosseini2013numerical


The most common numerical modeling of the phenomenon of GW in HSC found in the literature is a calculation of the effective material properties of the honeycomb structure by the homogenization process~\cite{shi1995derivation, qi2008ultrasonic, mustapha2014leaky, baid2015dispersion, sikdar2016guided}.
The homogenized model significantly simplifies the grid formulation and time consuming but does not accurately represent the wave propagation phenomenon in such material. 
A model with a real geometry of the hexagonal cell was presented by Ruzzenne et al. \cite{ruzzene2003wave}.
The finite element model and the theory of periodic structures was used to evaluate the dynamic behaviour of the honeycomb and cellular structures.
Recently, the simulations of the wave propagation in the HSC have been conducted with the commercially available finite element code~\cite{song2009guided, hosseini2013numerical, tian2015wavenumber, zhao2018wave}.
However, using the finite element method (FEM) to model GW is inefficient as it requires a significant amount of memory and is time-consuming.

The computational efficiency of the FEM in case of GW modelling in the HSC can be improved by using the time-domain spectral element method (SEM).
The SEM was originally used for the numerical solution of the fluid--flow in a channel by Patera \cite{patera1984spectral} but has also been successfully developed for elastic wave propagation~\cite{ostachowicz2011guided}.
Kudela proposed a model of the GW in HSC by the parallel implementation of the SEM \cite{kudela2016parallel}.
However, this model had a large number (1.5 million) of degrees-of-freedom (DOF), because cells of the core and skin plate were modelled by the three-dimensional (3D) spectral elements.
Although, the simulation was limited to only one skin plate and a small dimension of the HSC (\(179 \times 159 \) mm).

Above mentioned drawbacks were motivation to propose a new model of the HSC. In the present paper, the skin plates, adhesive layers and each wall of the hexagonal core were modelled by two-dimensional (2D) spectral elements.
However, 2D elements have nodes only in a mid-plane, therefore, there is no direct connection between the two adjacent structures.
This connection was implemented by interface elements based on Lagrange multipliers \cite{ashwin2014formulation, fiborek20192d}.
Additionally, the signal was generated and recorded with piezoelectric transducers (PZT), while Kudela excited the waves with an external force applied at the point of the panel \cite{kudela2016parallel}.
A non-matching interface between the transducers and the skin was used to avoid too complex mesh, likewise to the interfaces developed for the FEM \cite{flemisch2000elasto, flemisch2012non}. To the best of authors’
knowledge, it has not been implemented for the SEM.

The parametric study conducted in the paper leads to the determination of a model-assisted damage identification function (MADIF), which determines the influence of the size of the composite defect on wave propagation.
In this case, the defect is assumed to be a disbond between the skin and the core.

\section{The time-domain spectral element method formulation}
\label{sec:time_SEM}
\subsection{The spectral element method}
\label{sec:sem}
The general concept of the SEM is based on the idea of the FEM.
The similarity of both methods lies in the fact that the modelled domain is divided into non-overlapping finite elements, and external forces and arbitrary boundary conditions are imposed in the particular nodes.
The main difference between those methods is a choice of the shape function \( N=N(\xi )\), which is interpolated by Lagrange polynomial that passes through the element nodes. The nodes are localized on the endpoint of an interval, \(\xi\in[-1,1]\) and the roots of the first derivative of Legendre polynomial of degree \(p-1\):

\begin{eqnarray}
	(1-\xi^2)P'_{p-1}(\xi)=0
	\label{eq:nodes}
\end{eqnarray}

The approximation of an integral over the elements is achieved according to Gauss-Lobatto-Legendre (GLL) rule at points coinciding with the element nodes, 
and the weights \(w=w(\xi)\) calculated as:
\begin{eqnarray}
	{w(\xi)} = \frac{2}{p(p-1)(P_{p-1}(\xi))^2}
	\label{eq:weights}
\end{eqnarray}
This approach guarantees a diagonal mass matrix.
The shape functions and the weights for 2D or 3D elements are obtained by the Kronecker product of vectors of individual axes, denoted by \(\otimes\) as follows:
\begin{eqnarray}
	N(\xi,\eta) = N(\xi)\otimes N(\eta), & N(\xi,\eta,\zeta) = N(\xi)\otimes N(\eta)\otimes N(\zeta) \nonumber\\
	w(\xi,\eta) = w(\xi)\otimes w(\eta), & w(\xi,\eta,\zeta) = w(\xi)\otimes w(\eta)\otimes w(\zeta) 
	\label{eq:3Dshape_weights}
\end{eqnarray}
\subsection{2D spectral element}
\label{sec:2D_SEM}

According to the first-order shear deformation theory~\cite{reissner1945effect, mindlin1951influence} the vector of displacements is expressed as:
\begin{eqnarray}
	\left \{ \begin{array}{c}
		\textbf{u}^e(\xi,\eta,\zeta) \\
		\textbf{v}^e(\xi,\eta,\zeta) \\
		\textbf{w}^e(\xi,\eta,\zeta)
	\end{array} \right\} = 
	\left \{ \begin{array}{c}
		\textbf{u}_0^e(\xi,\eta) + z\boldsymbol{\varphi}_x^e(\xi,\eta)\\
		\textbf{v}_0^e(\xi,\eta) + z\boldsymbol{\varphi}_y^e(\xi,\eta)\\
		\textbf{w}_0^e(\xi,\eta) \\
	\end{array} \right\}
\end{eqnarray}
where \(\textbf{u}_0^e\), \(\textbf{v}_0^e\) and \(\textbf{w}_0^e\) are nodal displacements, \(\boldsymbol{\varphi}_x^e\), \(\boldsymbol{\varphi}_y^e\) are the rotations of the normal to the mid-plane with respect to axes \textit{x} and \textit{y}, respectively.
\begin{eqnarray}
	\left \{\begin{array}{c}
		\textbf{u}_0^e(\xi,\eta) \\
		\textbf{v}_0^e(\xi,\eta) \\
		\textbf{w}_0^e(\xi,\eta) \\
		\boldsymbol{\varphi}_x^e(\xi,\eta) \\
		\boldsymbol{\varphi}_y^e(\xi,\eta)
	\end{array} \right\} 
	 = \textbf{N}^e(\xi,\eta)\widehat{\textbf{d}}^e
	 = \sum_{n=1}^q\sum_{m=1}^p\textbf{N}_m^e(\xi)\textbf{N}_n^e(\eta)
	\left [ \begin{array}{c}
		\widehat{\textbf{u}}_0^e \\
		\widehat{\textbf{v}}_0^e \\
		\widehat{\textbf{w}}_0^e \\
		\widehat{\boldsymbol{\varphi}}_x^e \\
		\widehat{\boldsymbol{\varphi}}_y^e
	\end{array} \right]
\end{eqnarray}
The nodal bending strain-displacement relations are given in the form:
\begin{eqnarray}
	\boldsymbol{\epsilon}_b^e =
	\textbf{B}_b^e\widehat{\textbf{d}}^e = 
	\left [
	\begin{array}{ccccc}
		\frac{\partial N^e}{\partial x} & 0 & 0 & 0 & 0\\
		0 & \frac{\partial N^e}{\partial y} & 0 & 0 & 0\\
		\frac{\partial N^e}{\partial y} & \frac{\partial N^e}{\partial x} & 0 & 0 & 0\\
		0 & 0 & 0 & -\frac{\partial N^e}{\partial x} & 0\\
		0 & 0 & 0 & 0 & -\frac{\partial N^e}{\partial y}\\
		0 & 0 & 0 & -\frac{\partial N^e}{\partial y} & -\frac{\partial N^e}{\partial x}
	\end{array} \right]
	\left \{ \begin{array}{c}
		\textbf{u}_0^e \\
		\textbf{v}_0^e \\
		\textbf{w}_0^e \\
		\boldsymbol{\varphi}_x^e \\
		\boldsymbol{\varphi}_y^e
	\end{array} \right\}
\end{eqnarray}
The nodal shear strain-displacement relations are given in the form:
\begin{eqnarray}
	\boldsymbol{\epsilon}_s^e =
	\textbf{B}_s^e\widehat{\textbf{d}}^e = 
	\left [
	\begin{array}{ccccc}
		0 & 0 & \frac{\partial N^e}{\partial y} & -1 & 0\\
		0 & 0 & \frac{\partial N^e}{\partial y} & 0 & -1
	\end{array} \right]
	\left \{ \begin{array}{c}
		\textbf{u}_0^e \\
		\textbf{v}_0^e \\
		\textbf{w}_0^e \\
		\boldsymbol{\varphi}_x^e \\
		\boldsymbol{\varphi}_y^e
	\end{array} \right\}
\end{eqnarray}
\subsection{3D model of the PZT transducer}
\label{sec:3D_SEM}
The displacements vector of the PZT transducer is composed of three translational displacements, and it is defined as:
\begin{eqnarray}
	\left \{ \begin{array}{c}
		\textbf{u}^e(\xi,\eta,\zeta) \\
		\textbf{v}^e(\xi,\eta,\zeta) \\
		\textbf{w}^e(\xi,\eta,\zeta)
	\end{array} \right\}
	= \textbf{N}^e(\xi,\eta, \zeta)\widehat{\textbf{d}}^e\nonumber\\
	= \sum_{l=1}^r\sum_{n=1}^q\sum_{m=1}^p\textbf{N}_m^e(\xi)\textbf{N}_n^e(\eta)\textbf{N}_l^e(\zeta)
	\left \{ \begin{array}{c}
		\widehat{\textbf{u}}^e(\xi_m,\eta_n,\zeta_l) \\
		\widehat{\textbf{v}}^e(\xi_m,\eta_n,\zeta_l) \\
		\widehat{\textbf{w}}^e(\xi_m,\eta_n,\zeta_l)
	\end{array} \right\}
	\label{eq:3D_displ}
\end{eqnarray}
where \(\widehat{\textbf{u}}^e\), \(\widehat{\textbf{v}}^e\) and 
\(\widehat{\textbf{w}}^e\) are displacements of the element nodes in \(\xi,\eta\) and \(\zeta\) direction.

The nodal strain-displacement relations are given as \cite{kudela20093d}:
\begin{eqnarray}
	\boldsymbol{\epsilon}^e=\textbf{B}_{d}^e\widehat{\textbf{d}}^e=
	\left [
	\begin{array}{ccc}
		\frac{\partial N^e}{\partial x} & 0 & 0\\
		0 & \frac{\partial N^e}{\partial y} & 0\\
		0 & 0 & \frac{\partial N^e}{\partial z}\\
		0 & \frac{\partial N^e}{\partial z} & \frac{\partial N^e}{\partial y}\\
		\frac{\partial N^e}{\partial z} & 0 & \frac{\partial N^e}{\partial x}\\
		\frac{\partial N^e}{\partial y} & \frac{\partial N^e}{\partial x} & 0
	\end{array} \right]
	\left \{ \begin{array}{c}
		\textbf{u}^e \\
		\textbf{v}^e \\
		\textbf{w}^e
	\end{array} \right\}
\end{eqnarray}

The electromechanical coupling is governed by the linear constitutive equation of piezoelectric material according to~\cite{giurgiutiumicromechatronics, rekatsinas2017cubic}, and it is defined as:
\begin{eqnarray}
	\left [ 
	\begin {array}{c}
	\boldsymbol{\sigma}\\
	\textbf{D}
\end{array}\right ]=
\left [ 
\begin{array}{cc}
	\textbf{c}^E & -\textbf{e}^T \\
	\textbf{e} & \epsilon^S 
\end{array} \right ]
\left[ 
\begin{array}{c}
	\textbf{S}\\
	\textbf{E} 
\end{array} \right ]
\end{eqnarray}
where \(\textbf{c}^E\) is the stiffness coefficients measured at zero electric field, \textbf{e} is the piezoelectric coupling tensor and \(\boldsymbol{\epsilon}^S\) electric permittivity, \textbf{E} and \textbf{D} are the electric field and electric displacement measured at zero strain.
The superscript T denotes transpose matrix.
The eElectric field is defined as:
\begin{eqnarray}
\textbf{E}^e=-\textbf{B}_\phi^e \widehat{\phi}^e = \left[ \begin{array}{c}
	\frac{\partial N^e}{\partial \xi}\\
	\frac{\partial N^e}{\partial \eta}\\
	\frac{\partial N^e}{\partial \zeta}
\end{array} \right] \widehat{\phi}^e
\end{eqnarray}
\subsection{Elementary governing equations of motion}
\label{sec:motion}
The classical equations of motion \(\textbf{M}\ddot{\textbf{d}} + \textbf{C}\dot{\textbf{d}} + \textbf{K}\textbf{d} = \textbf{F}\) known from FEM are complemented by piezoelectric and interface coupling. Thus, the governing equations are defined as:
\begin{eqnarray}
	\textbf{M}_{dd} \widehat{\ddot{\textbf{d}}} + \textbf{C}_{dd} \widehat{\dot{\textbf{d}}} + \textbf{K}_{dd} \widehat{\textbf{d}} + \textbf{K}_{d\phi} \widehat{\boldsymbol{\phi}} = \textbf{F} - \textbf{G}^T \boldsymbol{\lambda},
	\label{eq:motion}
\end{eqnarray}
\begin{eqnarray}
	\textbf{K}_{\phi d}\widehat{\textbf{d}} + \textbf{K}_{\phi \phi} \widehat{\boldsymbol{\phi}} = \textbf{Q}
	\label{eq:piezocoupling}
\end{eqnarray} 

where \(\textbf{M}_{dd}\), \(\textbf{C}_{dd}\), \(\textbf{K}_{dd}\) are structural mass, damping and stiffness matrices, respectively; \(\textbf{K}_{\phi d}={\textbf{K}_{d\phi}}^T\) are piezoelectric coupling matrices; \(\textbf{K}_{\phi \phi}\) is the dielectric permittivity matrix, \(\widehat{\textbf{d}}\) is the vector of unknown nodal displacements, and \(\widehat{\boldsymbol{\phi}}\) is the electric potential vector; \(\textbf{F}\) the nodal external force vector, \(\textbf{Q}\) is the nodal charge vector, \(\boldsymbol{\lambda}\) is the Lagrange multipliers vector, and \(\textbf{G}\) is the interface coupling matrix; \((\dot{\ })=\frac{\partial}{\partial t}\).
The formulae of the matrices are provided in App.~\ref{app:matrices}.
The coupling is realised by imposing the traction forces, represented by a vector of Lagrange multipliers. 
\subsection{Transformation of the core elements}
\label{sec:transformation}
All core elements are rotated relative to both skins, so it is necessary to transform degrees of freedom from the local core's coordinate system to the global system.
For this purpose, an additional sixth global dof has been incorporated, i.e. $\boldsymbol{\varphi}_z^e$ rotation with respect to the z-axis .
Firstly, a part of the mass matrix accounted for rotary inertia (see App. \ref{app:matrices}) is globalised as follows:
\begin{eqnarray}
	\textbf{M}_R^g=\textbf{a}\,\textbf{M}_R^l\,\textbf{a}^{-1}
	\label{eq:inertia}
\end{eqnarray}
where \(\textbf{a}\) is general rotation matrix obtained from the multiplication of basic rotation matrices, and for the regular hexagonal core a is equal to:
\begin{eqnarray}
	\textbf{a}=\left [ 
	\begin{array}{ccc}
		m & -n & 0\\
		0 & 0 & -1\\
		n & m & 0\\
	\end{array}
	\right ]
	\label{eq:rotation}
\end{eqnarray}
where \(m=\cos(\alpha),\:n=\sin(\alpha)\), for \(\alpha\in{[0^\circ,\,60^\circ,\,120^\circ,180^\circ,\,240^\circ,\,300^\circ]}\) depending on which wall of the cell it is applied to.

The stiffness--strain relationship in local and global coordinate systems is given by:
\begin{eqnarray}
	\begin{array}{ccc}
		\left [
		\begin{array}{c}
			\sigma^l_{11}\\
			\sigma^l_{22}\\ 
			\sigma^l_{33}\\ 
			\sigma^l_{23}\\
			\sigma^l_{13}\\
			\sigma^l_{12}\\
		\end{array}
		\right ]=
		\textbf{c}\,\left [
		\begin{array}{c}
			\epsilon^l_{11}\\
			\epsilon^l_{22}\\ 
			\epsilon^l_{33}\\
			\gamma^l_{23}\\
			\gamma^l_{13}\\
			\gamma^l_{12}\\
		\end{array}
		\right ], & and &
		\left [
		\begin{array}{c}
			\sigma^g_{11}\\
			\sigma^g_{22}\\ 
			\sigma^g_{33}\\ 
			\sigma^g_{23}\\
			\sigma^g_{13}\\
			\sigma^g_{12}\\
		\end{array}
		\right ]=
		\bar{\textbf{c}}\,\left [
		\begin{array}{c}
			\epsilon^g_{11}\\
			\epsilon^g_{22}\\ 
			\epsilon^g_{33}\\
			\gamma^g_{23}\\
			\gamma^g_{13}\\
			\gamma^g_{12}\\
		\end{array}
		\right ]
	\end{array}
	\label{eq:stress_global}
\end{eqnarray}
The global form of the stiffness matrix \(K\), involves the determination of stiffness tensor \(\bar{c}\) from Eq.~\ref{eq:stress_global}.
The transformation of the stress tensor is expressed as:
\begin{eqnarray}
	\left [ 
	\begin{array}{ccc}
		\sigma^g_{11} & \sigma^g_{12} & \sigma^g_{13}\\
		\sigma^g_{21} & \sigma^g_{22} & \sigma^g_{23}\\
		\sigma^g_{31} & \sigma^g_{32} & \sigma^g_{33}\\
	\end{array}
	\right ]
	=
	\textbf{a}\,
	\left [ 
	\begin{array}{ccc}
		\sigma^l_{11} & \sigma^l_{12} & \sigma^l_{13}\\
		\sigma^l_{21} & \sigma^l_{22} & \sigma^l_{23}\\
		\sigma^l_{31} & \sigma^l_{32} & \sigma^l_{33}\\
	\end{array}
	\right ]
	\,\textbf{a}^T
	\label{eq:sigma_tensor}
\end{eqnarray}
After the multiplication of matrices and using the symmetry of the stress tensor Eq.~\ref{eq:sigma_tensor} can be rewritten in a simplified form:

\begin{eqnarray}
	\left [
	\begin{array}{c}
		\sigma^g_{11}\\
		\sigma^g_{22}\\ 
		\sigma^g_{33}\\ 
		\sigma^g_{23}\\
		\sigma^g_{13}\\
		\sigma^g_{12}\\
	\end{array}
	\right ]=
	\textbf{T}\,\left [
	\begin{array}{c}
		\sigma^l_{11}\\
		\sigma^l_{22}\\ 
		\sigma^l_{33}\\
		\sigma^l_{23}\\
		\sigma^l_{13}\\
		\sigma^l_{12}\\
	\end{array}
	\right ]
	\label{eq:stress}
\end{eqnarray}

Analogously, a strain relationship also is given in the form:
\begin{eqnarray}
	\begin{array}{ccc}
		\left [
		\begin{array}{c}
			\epsilon^g_{11}\\
			\epsilon^g_{22}\\ 
			\epsilon^g_{33}\\ 
			\epsilon^g_{23}\\
			\epsilon^g_{13}\\
			\epsilon^g_{12}\\
		\end{array}
		\right ]=
		\textbf{T}\,\left [
		\begin{array}{c}
			\epsilon^l_{11}\\
			\epsilon^l_{22}\\ 
			\epsilon^l_{33}\\
			\epsilon^l_{23}\\
			\epsilon^l_{13}\\
			\epsilon^l_{12}\\
		\end{array}
		\right ], & and & \left [
		\begin{array}{c}
			\epsilon_{11}\\
			\epsilon_{22}\\ 
			\epsilon_{33}\\ 
			\gamma_{23}\\
			\gamma_{13}\\
			\gamma_{12}\\
		\end{array}
		\right ]=
		\textbf{R}\,\left [
		\begin{array}{c}
			\epsilon_{11}\\
			\epsilon_{22}\\ 
			\epsilon_{33}\\
			\epsilon_{23}\\
			\epsilon_{13}\\
			\epsilon_{12}\\
		\end{array}
		\right ]
	\end{array}
	\label{eq:strain}
\end{eqnarray}
where \(\textbf{R}\) is the Reuter's matrix, defined as:
\begin{eqnarray}
	\textbf{R} = \left [
	\begin{array}{cccccc}
		1 & 0 & 0 & 0 & 0 & 0\\
		0 & 1 & 0 & 0 & 0 & 0\\
		0 & 0  & 1 & 0 & 0 & 0\\
		0 & 0 & 0 & 2 & 0 & 0\\
		0 & 0 & 0 & 0 & 2 & 0\\
		0 & 0 & 0 & 0 & 0 & 2
	\end{array}
	\right ]
	\label{eq:reuters}
\end{eqnarray}
Using the Eqs.~\ref{eq:rotation} through \ref{eq:reuters}, the relationship of stress--strain in the global frame can be expressed as:
\begin{eqnarray}
	\left [
	\begin{array}{c}
		\sigma^g_{11}\\
		\sigma^g_{22}\\ 
		\sigma^g_{33}\\ 
		\sigma^g_{23}\\
		\sigma^g_{13}\\
		\sigma^g_{12}\\
	\end{array}
	\right ]=
	\textbf{T}\,\textbf{c}\,\textbf{R}\,\textbf{T}^{-1}\,\textbf{R}^{-1}\,
	\left [
	\begin{array}{c}
		\epsilon^g_{11}\\
		\epsilon^g_{22}\\ 
		\epsilon^g_{33}\\
		\gamma^g_{23}\\
		\gamma^g_{13}\\
		\gamma^g_{12}\\
	\end{array}
	\right ]
	\label{eq:stress-strain}
\end{eqnarray}
Therefore, transformed stiffness tensor for the core's elements is equal:
\begin{eqnarray}
	\bar{\textbf{c}}=\textbf{T}\,\textbf{c}\,\textbf{R}\,\textbf{T}^{-1}\,\textbf{R}^{-1}
	\label{eq:c_global}
\end{eqnarray}


%% Loading bibliography style file
%\bibliographystyle{model1-num-names}
\bibliographystyle{model1-num-names}

% Loading bibliography database
\bibliography{model_hc}

\end{document}

