homogenization.tex

\chapter{Homogenization and Justification}

Mechanical models of highly periodic structures within the asymptotic homogenization theory, which defines a multi-scale method for determining the effective moduli of periodic media, have been studied profoundly for example by \textit{Bakhvalov, Panasenko} \cite{bakhvalov1989homogenisation:}. The approach exploits a separation of scales within the composite material, deriving a leading order homogenized equations governing the effective macroscopic behavior.

Given the assumptions for the multiscale elastodynamic model of cortical bone, the objective in this section will be to rigorously give a justification of the existence of such multiscale PDE system motivated from \cite{altenbach2018generalized}, \cite{Parnell2008} then show the schematic procedure to obtain the so-called macroscopic mechanical behavior of bone which will be used for the numerical simulations and finally give a convergence result of the solution between the multiscale model and the homogenized one. 

First, let us define the adequate spaces for our problem. Given the domain $\Omega \subset \mathbb{R}^d$, it will be denoted by $\mathbf{L}^2(\Omega)$ the $d$-dimensional vector functions, being each component on $L^2(\Omega)$. Similarly, it is defined the space $\mathbf{H}^1(\Omega)$ of d-dimensional vector functions, where as before, each component belongs to $H^1(\Omega)$ defined as usual in the literature \cite{evans2010partial}. 
In particular, denoting the vector trace operator by $\gamma: \mathbf{H}^1(\Omega) \rightarrow \mathbf{L}^2(\partial \Omega)$, we can define the space:
\begin{equation*}
    \mathbf{H}^1(\Omega, \Gamma_D) = \big \{ v \in \mathbf{H}^1(\Omega) \, \vert \, \gamma (v) \vert_{\Gamma_D} = \mathbf{0} \big \}
\end{equation*}

Moreover, by applying the \textit{Rellich-Kondrachov} theorem, the embedding from $\mathbf{H}^1(\Omega, \Gamma_D)$ into $\mathbf{L}^2(\Omega)$ is compact \footnote{From a more abstract point of view, it is possible to define $\kappa: \mathbf{L}^2( \Omega) \longrightarrow \mathbf{H}^{-1/2}(\Gamma_D)$ being the trace operator, i.e., $\kappa (u) = u \vert_{\Gamma_D}$. Then define an adequate space of square integrable vector functions with homogeneous Dirichlet condition on $\Gamma_D$ by $\mathbf{L}^2(\Omega, \Gamma_D) := \kappa^{-1}\big( \{ \mathbf{0}\}\big)$. Note that the operator $\kappa$ defined is linear and continuous in the corresponding spaces, moreover, since $\{\mathbf{0}\}$ is a close subset, then $\mathbf{L}^2(\Omega, \Gamma_D)$ is a close subspace of $\mathbf{L}^2(\Omega)$.
In particular, $\mathbf{L}^2(\Omega, \Gamma_D)$ is separable, so that there exist a Hilbertian base associated.}.

Given $\epsilon > 0$, let us fix a porosity $p \in (0,1)$ and consider $(C_{ijkl})_{ijkl}:=\mathbf{C}(p) \in \text{ lin}\big(\textbf{Sym}^{d\times d})$, the space of linear, continuous operators on the $d\times d$ symmetric matrices space, being uniformly elliptic and bounded elasticity tensors. Also, denote by $\rho^{\epsilon}(\mathbf{x}) = \rho \big( \frac{\mathbf{x}}{\epsilon}\big)$ an uniformly bounded density, $C_{ijkl}^{\epsilon}(\mathbf{x}) = C_{ijkl}(\frac{\mathbf{x}}{\epsilon})$, $\mathbf{F}(t, \mathbf{x})$ a force independent of $\epsilon$ and consider for fixed parameter $T > 0$ the following evolution PDE problem modelling the displacements:
\begin{equation}
    \label{MainPDE}
    \left \{
    \begin{array}{cc}
        \rho^{\epsilon} \partial_{tt} u^{\epsilon} - \nabla\cdot \sigma^{\epsilon}(u^{\epsilon})= \mathbf{0} & \text{ in } (0,T) \times \Omega \\
        \sigma^{\epsilon}_{ij}(u^{\epsilon}) = C_{ijkl}^{\epsilon} \mathbf{e}_{kl}(u^{\epsilon}) & \text{ in } (0,T)\times \Omega \\
        u^{\epsilon} = \mathbf{0} & \text{ on } (0,T) \times \Gamma_D \\
        \sigma^{\epsilon}_{ij} n_j = \mathbf{F} & \text{ on } (0,T) \times \Gamma_N \\
        \partial_t u^{\epsilon} = u^{\epsilon} = \mathbf{0} & \text{ on } \{t=0\} \times \Omega
    \end{array}
    \right.
\end{equation}
\begin{prop}
Assuming $0 < \rho_0 \leq \rho\big( \frac{\mathbf{x}}{\epsilon} \big) < + \infty$ and $\mathbf{F} \in L^{\infty}(0,T;\mathbf{L}^2(\Omega,\Gamma_N))$, there exist a unique solution $u^{\epsilon}$ to (\ref{MainPDE}) for each $\epsilon > 0$, such that
\begin{equation*}
    u^{\epsilon} \in \mathcal{C}^0(0,T;\mathbf{H}^1(\Omega,\Gamma_D)) \cap \mathcal{C}^1(0,T;\mathbf{L}^2(\Omega))
\end{equation*}
\end{prop}

\begin{rem}
Inspired from \cite{raviart1983introduction}, a constructive proof is given based in the spectral decomposition of the elastic operator on $\mathbf{H}^1(\Omega, \Gamma_D)$ and solutions of the reduced ODE's associated to each eigenvalue\footnote{The formulation of the elastodynamic problem with Dirichlet boundary conditions had been studied in literature, in particular on the multiscale works from \cite{panasenko2005multi-scale}, \cite{bakhvalov1989homogenisation:}. On the other hand, the formal study of Neumann boundary conditions has been done mainly on \cite{oleinik1992mathematical}.}. % of the elastic operator.
\end{rem}

In the following, it is denoted the dot product $(u,v)_{\Omega} := \int_{\Omega} \rho^{\epsilon} u\cdot v \, d\mathbf{x}$ for $u, v \in \mathbf{L}^2(\Omega)$, which is equivalent to the usual one of $\mathbf{L}^2(\Omega)$ for each $\epsilon$ by the assumptions of boundness on the density. On the other hand, $(\cdot\,, \cdot)_{\Gamma_N}$ denotes the usual dot product of $\mathbf{L}^2(\Gamma_N)$ vector functions.

Under such notation, problem (\ref{MainPDE}) can be rewritten in variational form satisfied in distributions over $(0,T)$ by: 
\begin{equation}
    \label{MainTimePDE}
    \begin{array}{cc}
        \text{Find } u^{\epsilon} \in \mathcal{C}^0 (0,T;\mathbf{H}^1(\Omega,\Gamma_D)) \cap \mathcal{C}^1(0,T;\mathbf{L}^2(\Omega)) & \text{ s.t. }\\
        \partial_{tt} (u^{\epsilon}(t),v)_{\Omega} + \mathcal{I}_{C}(u^{\epsilon}(t),v) = (\mathbf{F}(t),v)_{\Gamma_N}&  \forall v \in \mathbf{H}^1(\Omega,\Gamma_D) \\
        \partial_{t} u^{\epsilon}(0) = u^{\epsilon}(0) = \mathbf{0} & \\
    \end{array}
\end{equation}
being $\mathcal{I}_{C}(u,v) := \int_{\Omega} C_{ijkl}^{\epsilon}\mathbf{e}_{kl}(u^{\epsilon}(t)) \partial_{x_j} v_i$.
Let us prove the existence and uniqueness of solution for (\ref{MainTimePDE}).
\begin{proof}
\begin{enumerate}
    \item To this end, it will be considered approximate solutions to (\ref{MainTimePDE}). It is defined the subspace $V_m$ generated by the first $m \in \mathbb{N}$ eigenvectors $\{w_1, \dots, w_m \}$ being $(w_i)_{i \in \mathbb{N}} \subset \mathcal{C}^{\infty}(\Omega)$ a Hilbertian base of $\mathbf{H}^1(\Omega, \Gamma_D)$ were the regularity follows from bootstrap and by the Sobolev embedding.
    Then, let us consider the problem in time, defined by:
    \begin{equation}
        \label{ApproxTimePDE}
        \begin{array}{cc}
            \text{Find } u^{\epsilon}_m: t \in (0,T) \longrightarrow u_m(t) \in V_m & \text{ s.t. } \\
            \partial_{tt}(u_m^{\epsilon}(t),v)_{\Omega} + \mathcal{I}_{C}(u^{\epsilon}_m(t),v) = (\mathbf{F}(t),v)_{\Gamma_N} & \forall v \in \mathbf{H}^1(\Omega,\Gamma_D) \\
            \partial_{t} u^{\epsilon}_m(0) = u^{\epsilon}_m(0) = \mathbf{0} & \\
        \end{array}
    \end{equation}
    \begin{rem}
    Note that the linear operator $b(t)(v) := (\mathbf{F}(t),v)_{\Gamma_N}$ is well-defined and moreover continuous on $\mathbf{H}^1(\Omega,\Gamma_D)$ for each $t \in (0,T)$.
    This follows from the trace theorem applied on $\mathbf{H}^1(\Omega, \Gamma_D)$ since $\forall t \geq 0$, we have the bounds:
    \begin{align*}
        \vert b(t)(v) \vert & \leq \Vert \mathbf{F}(t) \mathbb{I}_{\Gamma_N} \Vert_{\mathbf{L}^2(\partial \Omega)} \Vert \gamma (v) \Vert_{\mathbf{L}^2(\Gamma_N)} \\
        & \leq \Vert \mathbf{F}\Vert_{\mathbf{L}^{\infty}(0,T;\mathbf{L}^2(\Gamma_N))} \Vert \gamma(v) \Vert_{\mathbf{L}^2(\partial \Omega)} \\
        & \leq \Vert \mathbf{F} \Vert_{\mathbf{L}^{\infty}(0,T;\mathbf{L}^2(\Gamma_N))} \Vert v \Vert_{\mathbf{H}^1(\Omega)}\\
        & = \Vert \mathbf{F}\Vert_{\mathbf{L}^{\infty}(0,T;\mathbf{L}^2(\Gamma_N))} \Vert v \Vert_{\mathbf{H}^1(\Omega, \Gamma_D)}
    \end{align*}
    \end{rem}
    Defining $u^{\epsilon}_m(t) = \sum_{i=1}^m \alpha_i^{\epsilon}(t) w_i$ with $\alpha_i^{\epsilon} (t) = (u^{\epsilon}_m(t),w_i)_{\Omega}$ it follows that the functions $\alpha_i^{\epsilon}$ are solutions to the ODE system:
    \begin{equation}
        \label{AlphaODE}
        \begin{array}{cc}
            \partial_{tt} \alpha_i^{\epsilon}(t) + \lambda_i \alpha_i^{\epsilon}(t) = (\mathbf{F}(t),w_i)_{\Gamma_N}& \forall t \in (0,T) \\
            \partial_t \alpha_i^{\epsilon}(0) = \alpha_i^{\epsilon}(0) = 0 & 
        \end{array}
    \end{equation}
    being $(\lambda_i)_{i \geq 1}$ a positive, increasing sequence of eigenvalues associated to the decomposition of the operator $\mathcal{I}_C(\cdot, \cdot)$.\\
    Let us note by applying the variation of parameters technique, that the solution to (\ref{AlphaODE}) is given in the form:
    \begin{equation}
        \label{AlphaODEsol}
        \alpha_i^{\epsilon} (t) = \frac{1}{\sqrt{\lambda_i}} \int\limits_0^t sin(\sqrt{\lambda_i} (t-s)) (\mathbf{F}(t),w_i)_{\Gamma_N} \, ds
    \end{equation}
    so that, defining the matrix for $\omega \in \mathbb{R}$ by:
    \begin{equation*}
        Q(\omega) =
        \begin{bmatrix}
        cos(\omega) & sin(\omega) \\
        -sin(\omega) & cos(\omega)
        \end{bmatrix}
    \end{equation*}
    it is obtained the following relation for the solution and their derivative in the form:
    \begin{equation}
        \label{MatrixODEsol}
        \begin{bmatrix}
        \sqrt{\lambda_i} \alpha_i^{\epsilon}(t) \\
        \partial_{t} \alpha_i^{\epsilon}(t) 
        \end{bmatrix}
        = \int \limits_0^t Q\big(\sqrt{\lambda_i}(t-s)\big)
        \begin{bmatrix}
        0 \\
        (\mathbf{F}(t),w_i)_{\Gamma_N}
        \end{bmatrix}
    \end{equation}
    
    \item Next, the sequence $(u^{\epsilon}_m)_{m \in \mathbb{N}}$ is of Cauchy type on the spaces $\mathcal{C}^0(0,T;\mathbf{H}^1(\Omega, \Gamma_D))$ and $\mathcal{C}^1(0,T; \mathbf{L}^2(\Omega))$.\\
    Let $m,p$ be two integers such that $p > m \geq 1$ then from (\ref{AlphaODEsol}) it follows that
    \begin{equation*}
        \mathcal{I}_C(u_p^{\epsilon}(t)- u_m^{\epsilon}(t),u_p^{\epsilon}(t)- u_m^{\epsilon}(t)) + \vert \partial_t (u_p^{\epsilon}(t)- u_m^{\epsilon}(t))_{\Omega} \vert^2 = \sum_{i=m+1}^p \big( \lambda_i \vert \alpha_i^{\epsilon}(t) \vert^2 + \vert \partial_t \alpha_i^{\epsilon}(t) \vert^2 \big) 
    \end{equation*}
    and since $Q(\omega)$ is an orthogonal matrix, from (\ref{MatrixODEsol}) it follows that
    \begin{equation}
        \label{AlphaBound}
        \big( \lambda_i \vert \alpha_i^{\epsilon}(t) \vert^2 + \vert \partial_t \alpha_i^{\epsilon}(t) \vert^2 \big)^{1/2} \leq \int \limits_0^t \vert (\mathbf{F}(t),w_i)_{\Gamma_N} \vert \, ds 
    \end{equation}
    so that, by using the Cauchy-Schwatz inequality, it can be obtained
    \begin{align*}
        \lambda_i \vert \alpha_i^{\epsilon}(t) \vert^2 + \vert \partial_t \alpha_i^{\epsilon}(t) \vert^2 &\leq 2 \big( \int_0^t \vert (\mathbf{F}(t),w_i)_{\Gamma_N} \vert \, ds \big)^2 \\
        & \leq  2 t \int_0^t \vert (\mathbf{F}(t),w_i)_{\Gamma_N} \vert^2 \, ds
    \end{align*}
    from which it can be deduced the bound
    \begin{equation*}
        \mathcal{I}_C \big(u_p^{\epsilon}(t) - u_m^{\epsilon}(t),u_p^{\epsilon}(t) - u_m^{\epsilon}(t) \big) + \vert \partial_t (u_p^{\epsilon}(t) - u_m^{\epsilon}(t)) \vert^2 \leq 2 \sum_{i=m+1}^p t \int_0^t \vert (\mathbf{F}(s),w_i)_{\Gamma_N} \vert^2 \, ds
    \end{equation*}
    Now, since $\mathbf{F} \in L^{\infty}(0,T;\mathbf{L}^2( \Gamma_N))$ it follows
    \begin{equation*}
        \underset{p,m \longrightarrow + \infty}{\text{lim}} \sum_{i=m+1}^p \big \{ T \int_0^T \vert (\mathbf{F}(s),w_i)_{\Gamma_N}\vert^2 \, ds \big \} = 0
    \end{equation*}
    and by using the uniform ellipticity of the tensor $\mathbf{C}$, it can be concluded that $(u_m^{\epsilon}(t))_{m \in \mathbb{N}}$ is a Cauchy sequence on the spaces $\mathcal{C}^0(0,T; \mathbf{H}^1(\Omega, \Gamma_D))$ and $\mathcal{C}^1(0,T; \mathbf{L}^2(\Omega, \Gamma_D))$.
    
    
    \item Since the above spaces are complete, there exist $u^{\epsilon}(t)$ limit as $m \longrightarrow +\infty$ of $u^{\epsilon}_m(t)$ belonging to the spaces $\mathcal{C}^0(0,T; \mathbf{H}^1(\Omega, \Gamma_D))$ and $\mathcal{C}^1(0,T; \mathbf{L}^2(\Omega))$.\\
    Let us see that $u^{\epsilon}$ effectively solves the problem (\ref{MainPDE}). To do this, let $m \geq 1$ and consider $\psi \in \mathcal{D}(0,T)$, since $u_m(t)$ solves (\ref{AlphaODE}) it follows:
    \begin{equation*}
        \int_0^T (u_m^{\epsilon}(t),v)_{\Omega} \partial_{tt}\psi(t) + \int_0^T \mathcal{I}_C (u_m^{\epsilon}(t),v) \psi(t)  = \int_0^T (\mathbf{F}(t),v)_{\Gamma_N} \psi(t) 
    \end{equation*}
    so that, in the limit as $m \longrightarrow 0$, the solution $u^{\epsilon}(t)$ solves the problem
    \begin{equation*}
        \int_0^T (u^{\epsilon}(t),v)_{\Omega} \partial_{tt}\psi(t) + \int_0^T \mathcal{I}_C (u^{\epsilon}(t),v) \psi(t) = \int_0^T (\mathbf{F}(t),v)_{\Gamma_N}\psi(t) 
    \end{equation*}
    with $\partial_t u^{\epsilon}(0) = u^{\epsilon}(0) = 0$ and the desired results follows. Moreover, continuity with respect to boundary data is obtained, since from (\ref{AlphaBound}) the solution $u^{\epsilon}(t)$ in the limit satisfies the inequality 
    \begin{equation*}
        \big( \mathcal{I}_C(u^{\epsilon}(t), u^{\epsilon}(t))+ \vert \partial_t u^{\epsilon}(t) \vert^2 \big)^{1/2} \leq \int_0^T \Vert \mathbf{F}(t) \Vert_{\mathbf{L}^2(\Gamma_N)} \, dt \leq T\, \Vert \mathbf{F}\Vert_{L^{\infty}(0,T;\mathbf{L}^2(\Gamma_N)}
    \end{equation*}
    and by the uniform boundness of the tensor $C$, the following bound is obtained
    \begin{equation*}
        \Vert u^{\epsilon} \Vert_{L^{\infty}(0,T;\mathbf{H}^1(\Omega, \Gamma_D))}  + \Vert \partial_t u^{\epsilon} \Vert_{L^{\infty}(0,T;\mathbf{L}^2(\Omega))} \lesssim T \, \Vert \mathbf{F}\Vert_{L^{\infty}(0,T; \mathbf{L}^2(\Gamma_N))}
    \end{equation*}
    Note in particular that the above bound is independent of $\epsilon > 0$ since there is no dependency on the right-hand side of the inequality, from which follows well-possedness for each multiscale elastic problem with parameters $\epsilon > 0$.
\end{enumerate}

\end{proof}

%%%%%
%%%%% THE ABOVE ADDED 21/10/2018
%%%%%

Now, the idea is to obtain a result of time regularity for the solution of the multiscale elastodynamic model. The motivation comes from \textit{Panasenko} works \cite{panasenko2005multi-scale}, who obtained such type of results in a general viscoelastic case with full \textit{Dirichlet} boundary conditions.

In this case, a bootstrap type method is used to obtain the time-regularity for the mixed boundary problem, with the bounds derived from spectral theory.

\begin{prop}
\label{BootstrapingProp}
Let $u^{\epsilon} \in \mathcal{C}^{0}(0,T; \mathbf{H}^1(\Omega, \Gamma_D))$ denote the unique solution of the mixed boundary multiscale problem and suppose also the regularity for the parameters and surface force given by $C_{ijkl} \in L^{\infty}(\mathbf{Y})$ and $\mathbf{F} \in \mathcal{C}^p(0,T; \mathbf{L}^{2}(\Gamma_N))$, then for each $p \geq 1$ we have that:
\begin{equation*}
    u^{\epsilon}, \partial_t u^{\epsilon} \in \mathcal{C}^p(0,T; \mathbf{L}^2(\Omega)).
\end{equation*}
\end{prop}
\begin{rem}
In particular, from the proof it's also possible to obtain that
\begin{equation*}
    u^{\epsilon}, \partial_t u^{\epsilon} \in \mathcal{C}^p(0,T; \mathbf{H}^1(\Omega, \Gamma_D)).
\end{equation*}
\end{rem}

\begin{proof}
Let us develop the main idea of the proof:\\
\begin{enumerate}
    \item Using the spectral theory for the existence of multiscale problem, it is obtained a solution $u^{\epsilon}$ having regularity:
    \begin{equation*}
        u^{\epsilon} \in \mathcal{C}^0(0,T;\mathbf{H}^1(\Omega, \Gamma_D)) \cap \mathcal{C}^1(0,T;\mathbf{L}^2(\Omega))
    \end{equation*}
    thus the results being valid for initial cases $p = 0, 1$.
    \item Let us take $\vert h \vert \ll 1$ such that $t+h \in (0,T)$ and define the difference with parameter $\epsilon > 0$ fixed as:
    \begin{equation*}
        u_h^{\epsilon}(t, \mathbf{x}, \frac{\mathbf{x}}{\epsilon}) := u^{\epsilon} (t+h, \mathbf{x}, \frac{\mathbf{x}}{\epsilon}) - u^{\epsilon}(t, \mathbf{x}, \frac{\mathbf{x}}{\epsilon})
    \end{equation*}
    Now, defining the difference of surface forces as
    \begin{equation*}
        \mathbf{F}_h(t):= \mathbf{F}(t+h) - \mathbf{F}(t)
    \end{equation*}
    it follows that the problem satisfied for functions $u_h^{\epsilon}$ is given by:
    \begin{equation*}
        (P_{\epsilon}) \left \{
        \begin{array}{cc}
            \rho^{\epsilon} \partial_{tt} u_h^{\epsilon} - \nabla \cdot \sigma^{\epsilon}( u_h^{\epsilon} ) = \mathbf{0}  &  \text{ in } (0,T)\times \Omega \\
            \sigma^{\epsilon}_{ij} (u^{\epsilon}_h) = C^{\epsilon}_{ijkl} \mathbf{e}_{kl}(u^{\epsilon}_h) & \text{ in } (0,T) \times \Omega \\
            C^{\epsilon}_{ijkl}\mathbf{e}_{kl}( u_h^{\epsilon})n_j = \mathbf{F}_h & \text{ on } (0,T)\times \Gamma_N \\
            u^{\epsilon}_h = \mathbf{0} & \text{ on }(0,T)\times \Gamma_D
        \end{array}
        \right .
    \end{equation*}
    with initial condition at rest (i.e. $\partial_t u_h^{\epsilon} = u_h^{\epsilon} = \mathbf{0}$ on $\{t=0\} \times \Omega$).
    
    \item Using the assumption of regularity for the surface force $\mathbf{F}$, it can be obtained from spectral theory the bound:
    \begin{equation*}
        \Vert u_h^{\epsilon} \Vert_{\mathbf{H}^1(\Omega, \Gamma_D)} + \Vert \partial_t u_h^{\epsilon}\Vert_{\mathbf{L}^2(\Omega)} \lesssim \int \limits_t^{t+h} \Vert \mathbf{F}(s) \Vert_{\mathbf{L}^2(\Gamma_N)} \, ds
    \end{equation*}
    so, from the continuity it is obtained that as $h \rightarrow 0$, the terms 
    \begin{equation*}
        \Vert u_h^{\epsilon} \Vert_{\mathbf{H}^1(\Omega, \Gamma_D)}, \Vert \partial_t u_h^{\epsilon} \Vert_{\mathbf{L}^2 (\Omega)} \rightarrow 0
    \end{equation*}
    then $u^{\epsilon} \in \mathcal{C}^2(0,T; \mathbf{L}^2(\Omega))$.
    
    
    \item Now, let us observe that we can take a time derivative to the full problem $(P_{\epsilon})$, and obtain a similar multiscale problem. Applying again the spectral theory, it is possible to obtain a solution $v_h^{\epsilon} = \partial_t u_h^{\epsilon}$, satisfying same bounds as before. In such a way, it follows that:
    \begin{equation*}
        \Vert \partial_{tt} u_h^{\epsilon} \Vert_{\mathbf{H}^1(\Omega, \Gamma_D)} + \Vert \partial_t u_h^{\epsilon} \Vert_{\mathbf{L}^2 (\Omega)} \lesssim \int_t^{t+h} \Vert \partial_t \mathbf{F}(t) \Vert_{\mathbf{L}^2(\Gamma_N)}
    \end{equation*}
    so that, as $h \rightarrow 0$ we conclude the regularity
    \begin{equation*}
        \partial_t u_h^{\epsilon} \in \mathcal{C}^0(0,T;\mathbf{H}^1(\Omega, \Gamma_D)) \cap \mathcal{C}^1(0,T;\mathbf{L}^2(\Omega)) 
    \end{equation*}
    then
    \begin{equation*}
        u_h^{\epsilon} \in \mathcal{C}^3(0,T; \mathbf{L}^2(\Omega))
    \end{equation*}
    and by a bootstrap type argument, the results follows for each $p \geq 1$.
\end{enumerate}

\end{proof}


%%%%%
%%%%% THE ABOVE ADDED 21/10/2018
%%%%%

\section{Homogenization Procedure}

 In this section, it is described the procedure to obtain effective (or macroscopic) equations derived from the multiscale model (\ref{MainPDE}), by means of asymptotic approximation of the displacement $u^{\epsilon}(\mathbf{x},t)$ proposed by the two-scale homogenization theory inspired from \cite{altenbach2018generalized}. This framework enables us to obtain an effective PDE model governing the overall macroscopic mechanical behavior that incorporates the highly oscillatory microstructure variation by using the cell problems. It provides moreover an algorithmic procedure suitable for numerical implementation\footnote{It keeps the mechanical behavior explicitly stated in all the derivation of the homogenized model, thus keeping its physical interpretation in the context of solid mechanics.}.

To give jargon regarding homogenization literature, the variable $\mathbf{x}$ is defined as slow or global coordinate, while $\mathbf{y}$ denote the fast or local variable, related by $\mathbf{y} = \epsilon^{-1}\mathbf{x}$ as used before. In particular, the $\epsilon$ parameter represent the highly-oscillatory periodicity assumed on the material components described by the cell structure.

\subsection{Two-Scale Asymptotic Homogenization}
It is assumed that displacement solutions to problem (\ref{MainPDE}) can be expressed as expansion at different $\epsilon$ scales. Thus, the solution is found by an asymptotic expression at each time $t \in \mathbb{R}_+$ in the form:
\begin{equation}
    \label{AsymptoticExpansion}
    u^{\epsilon}(t, \mathbf{x}) = \sum_{a=0}^{\infty} \epsilon^a u^{(a)}(t, \mathbf{x},\mathbf{y}) 
\end{equation}
where the vector functions $u^{(a)}(\mathbf{x}, \cdot)$ are assumed Y-periodic, spanning the full domain $\Omega$ for each $a>0$. Moreover, throughout the developments it is assumed regularity for each $t \in (0,T)$ in the form:
\begin{equation*}
    u^{(a)}(t, \mathbf{x},\mathbf{y}) \in \mathbf{H}^1\big(\Omega; \, \mathbf{H}^1(\mathbf{Y})\big) \quad \forall a \in \mathbb{N}
\end{equation*}
where $\mathbf{H}^1(\mathbf{Y})$ denotes the space of d-dimensional vector function on $H^1(\mathbf{Y})$.
In what follows, it will be imposed several restrictions for the proposed solution (\ref{AsymptoticExpansion}) to effectively satisfy the problem (\ref{MainPDE}). In particular, the regularity assumption will naturally arise from restriction imposed for each $u^{(a)}$. In all that follows, it is dropped the time dependency in the functions for easiness of exposition.

Since there is an explicit relation between slow and fast variables, it is considered expressions of strain rate dependent on each one. Let us define for $\Phi \in \mathcal{C}^{\infty}(\Omega \times \mathbf{Y})$ the strain on each variable by:
\begin{equation*}
    \mathbf{e}_{kl,\alpha} (\Phi(\alpha)) = \frac{1}{2}(\partial_{\alpha_l} \Phi_k (\alpha) + \partial_{x_k} \Phi_l (\alpha)) \quad \forall \alpha \in \{\mathbf{x}, \mathbf{y}\}
\end{equation*}
thus, after applying it over the strain tensor $\mathbf{e}(\cdot)$ it follows the relation at components $k,l$ given in the form:
\begin{equation}
    \label{Multiscale-Strain}
    \mathbf{e}_{kl} ( u^{(a)}(\mathbf{x}, \frac{\mathbf{x}}{\epsilon})) = \mathbf{e}_{kl,x}( u^{(a)} (\mathbf{x},\mathbf{y})) + \epsilon^{-1} \mathbf{e}_{kl,y} (u^{(a)}(\mathbf{x},\mathbf{y})).
\end{equation}
%%%%%%%%%%%%
In the following, it is used a separation of scales to obtain the so-called homogenized coefficients, approximating the overall macroscopic (effective) behavior by applying the asymptotic solution, i.e., seeking
\begin{equation*}
    \rho^{\epsilon}(\mathbf{x}) u(\mathbf{x})^{\epsilon} + P^{(\epsilon)}(u^{\epsilon}(\mathbf{x})) = \mathcal{O}(\epsilon^{\gamma}) \quad \gamma \geq 0
\end{equation*}
where the elasticity operator $P^{(\epsilon)}$ is expected to satisfy under some corrections
\begin{equation}
    \label{P-Operator}
    P^{(\epsilon)}(u^{\epsilon}(\mathbf{x}) = - div \, \big[\mathbf{C}^{hom}(\mathbf{x}): \mathbf{e}(u^{(0)}(\mathbf{x}))\big] + \mathcal{O}(\epsilon)
\end{equation}
where $\mathbf{C}^{hom}(\mathbf{x})$ is some four-order tensor to be found and the right-hand side $\mathcal{O}(\epsilon^{\gamma})$ denotes the error expected from such an approximation. 

To handle derivatives from (\ref{P-Operator}), let us define the operator 
\begin{equation*}
    L_{\alpha \beta} (\cdot) = - \partial_{\alpha_j} \big( C_{ijkl} (\mathbf{y}) \mathbf{e}_{kl, \beta}(\cdot) \big), \quad \alpha, \beta \in \{ \mathbf{x},\mathbf{y} \}
\end{equation*}
And recall that by the chain rule $\forall \, \Phi(\mathbf{x},\mathbf{y})$ enough regular vector functions, it follows:
\begin{equation*}
    \partial_{x_j} (\Phi (\mathbf{x}, \frac{\mathbf{x}}{\epsilon})) = \big \{ \partial_{x_j} \Phi (\mathbf{x}, \mathbf{y}) + \frac{1}{\epsilon} \partial_{y_j} \Phi(\mathbf{x},\mathbf{y}) \big \}_{\mathbf{y}= \epsilon^{-1}\mathbf{x}}.
\end{equation*}

Now, using (\ref{AsymptoticExpansion}), the idea is to regroup in powers of $\epsilon$ thus impose constrains over each term $\epsilon^{(a)}$ to satisfy the main formulation (\ref{MainPDE}). To this end, note that from (\ref{P-Operator}) it follows:
\begin{equation*}
    P^{(\epsilon)}(u(\mathbf{x},\epsilon)) = A + B + C + \mathcal{O}(\epsilon)
\end{equation*}
where 
\begin{equation*}
    \begin{array}{cc}
        A &= P^{(\epsilon)}(u^{(0)}(\mathbf{x},\mathbf{y})) \\
        B &= \epsilon P^{(\epsilon)}(u^{(1)}(\mathbf{x},\mathbf{y})) \\
        C &= \epsilon^2 P^{(\epsilon)}(u^{(2)}(\mathbf{x},\mathbf{y})) \\
    \end{array}
\end{equation*}
Expanding every term and using (\ref{Multiscale-Strain}), it follows:
\begin{align*}
    A &= - \partial_{x_j} \big( C_{ijkl}\mathbf{e}(u^{(0)}) \big) \\
    &=- \partial_{x_j} \big( C_{ijkl} \mathbf{e}_{kl,x} (u^{(0)}) + \frac{1}{\epsilon}C_{ijkl}\mathbf{e}_{kl,y}(u^{(0)}) \big)\\
    &= - L_{xx}u^{(0)} - \frac{1}{\epsilon} L_{yx}u^{(0)} - \frac{1}{\epsilon} L_{xy}u^{(0)} - \frac{1}{\epsilon^2}L_{yy}u^{(0)}
\end{align*}
Similarly, for the other terms it can be obtained:
\begin{align*}
    B &= -\epsilon L_{xx} u^{(1)} - L_{yx}u^{(1)} - L_{xy} u^{(1)} - \frac{1}{\epsilon} L_{yy}u^{(1)} \\
    C &= -\epsilon^2 L_{xx} u^{(2)} - \epsilon L_{yx}u^{(2)} - \epsilon L_{xy} u^{(2)} - L_{yy}u^{(2)} 
\end{align*}
So, the following necessary conditions for powers of $\epsilon$ are derived:
\begin{equation}
    \label{Epsilon-Separation}
    \begin{array}{ccc}
        \epsilon^{-2} \longrightarrow & L_{yy} u^{(0)} &= \mathbf{0} \\
        \epsilon^{-1} \longrightarrow & L_{xy}u^{(0)} + L_{yx} u^{(0)} + L_{yy} u^{(1)} &= \mathbf{0} \\
        \epsilon^{0} \longrightarrow & L_{xx} u^{(0)} + L_{xy} u^{(1)} + L_{yx} u^{(1)} + L_{yy} u^{(2)} + \rho^{\epsilon} u^{(0)} &= \mathbf{0}
    \end{array}
\end{equation}

Let us consider the asymptotic expansion (\ref{AsymptoticExpansion}) as an approximation for the exact solution of the original problem (\ref{MainPDE}) where their boundary conditions had been replaced to the behavior at $\mathcal{O}(1)$.
\begin{equation*}
    \begin{array}{cc}   
        u^{(0)}(\mathbf{x},\mathbf{y}) = \mathbf{0} & \forall (\mathbf{x}, \mathbf{y}) \in \Gamma_D\times \mathbf{Y}\\
        \big(\mathbf{C}(\mathbf{y}): \mathbf{e} (u^{(0)}(\mathbf{x}, \mathbf{y}) \big) \cdot n = \hat{\mathbf{F}}(\mathbf{x}) & \forall (\mathbf{x},\mathbf{y}) \in \Gamma_N \times \mathbf{Y}
    \end{array}
\end{equation*}
and the remaining terms for each $a \in \mathbb{N}$ in the expansion are assigned by
\begin{equation*}
    \begin{array}{cc}
        u^{(a)}(\mathbf{x},\mathbf{y}) = \mathbf{0} & \forall (\mathbf{x}, \mathbf{y}) \in \Gamma_D \times \mathbf{Y} \\
        \big( \mathbf{C}(\mathbf{y}): \mathbf{e} (u^{(a)}(\mathbf{x},\mathbf{y}) \big) \cdot n = \mathbf{0} & \forall  (\mathbf{x}, \mathbf{y}) \in \Gamma_N\times \mathbf{Y} 
    \end{array}
\end{equation*}

Let us recall a classical results for elliptic problems. It relates the existence for problems in (\ref{Epsilon-Separation}) by some compatibility conditions.
\begin{lem}
\label{ExistenceLemma}
Let $f(\cdot)$ be a square integrable function over $\mathbf{Y}$. Consider the problem:
\begin{equation*}
    L_{yy} \Phi(\mathbf{y}) = f(\mathbf{y}) \text{ in } \mathbf{Y}
\end{equation*}
where $\Phi$ is $\mathbf{Y}$-periodic function. Then it holds:
\begin{enumerate}
    \item There exist a $\mathbf{Y}$-periodic solution $\Phi$ iff $\langle f \rangle = \mathbf{0}$
    \item If a $Y$-periodic solution $\Phi$ exists, then it's unique up to a constant vector $\mathbf{c} \in \mathbb{R}^d$.
\end{enumerate}
\end{lem}

\begin{rem}
It is being used as notation
\begin{equation*}
    \langle f \rangle := \frac{1}{\vert Y \vert} \int_{\mathbf{Y}} f(\mathbf{y}) \, d\mathbf{y}
\end{equation*}
where $\vert Y \vert$ denotes the measure of the set $\mathbf{Y}$.
\end{rem}

\subsection{Contribution at Second Order}
For the contribution at order $\mathcal{O}(\epsilon^{-2})$ recall that problem (\ref{MainPDE}) with their boundary conditions states:
\begin{equation}
    \label{Order-2VectorPDE}
    \left \{
    \begin{array}{cc}
        L_{yy} u^{(0)}( \mathbf{x},\mathbf{y}) = \mathbf{0} & \text{ in } \Omega \times \mathbf{Y}\\
        u^{(0)} (\mathbf{x},\mathbf{y}) = \mathbf{0} & \text{ in } \Gamma_D \times \mathbf{Y} \\
        \big( \mathbf{C}(\mathbf{y}) :\mathbf{e}(u^{(0)}(\mathbf{x},\mathbf{y}) \big) \cdot n = \mathbf{F}(t, \mathbf{x}) & \text{ in } \Gamma_N \times \mathbf{Y} \\
    \end{array}
    \right .
\end{equation}
By the Lemma (\ref{ExistenceLemma}) it can be deduced that $u(\mathbf{x},\mathbf{y})$ is solution of (\ref{Order-2VectorPDE}) iff it is constant with respect to the $\mathbf{y}$-variable. It implies then:
\begin{equation}
    \label{IndepencyofY}
    u^{(0)}(t, \mathbf{x},\mathbf{y}) = v(t, \mathbf{x})
\end{equation}
i.e. being independent with respect to the microstructure (or fast variable) $\mathbf{Y}$. In particular, the contribution of the boundary conditions is associated to such vector function $v(t, \mathbf{x})$, at each time $t \in (0,T)$ in the form:
\begin{equation*}
    \begin{array}{cc}
        v(\mathbf{x}) = \mathbf{0} & \text{ on } \Gamma_D\\
        \big(\mathbf{C}(\mathbf{y}):\mathbf{e}(v(\mathbf{x})) \big) \cdot n = \mathbf{F}(\mathbf{x}) & \text{ on } \Gamma_N \times \mathbf{Y}
    \end{array}
\end{equation*}

\subsection{Contributions at First Order}
For the contribution at order $\mathcal{O}(\epsilon^{-1})$ recall by (\ref{IndepencyofY}) it follows $u^{(0)}(\mathbf{x},\mathbf{y}) = v(\mathbf{x})$ then it can deduced by definition:
\begin{equation*}
    L_{xy} v(\mathbf{x}) = \mathbf{0}
\end{equation*}
So, the problem formulation in (\ref{Epsilon-Separation}) is reduced to
\begin{equation}
    \label{Order-1VectorPDE}
    \begin{array}{cc}
        L_{yy} u^{(1)}(\mathbf{x},\mathbf{y}) = - L_{yx} v(\mathbf{x}) & \text{ in } \Omega \times \mathbf{Y}
    \end{array}
\end{equation}
Now, using the Lemma (\ref{ExistenceLemma}) on (\ref{Order-1VectorPDE}) taking into account (\ref{IndepencyofY}), the $Y$-periodicity of $\mathbf{C}(\mathbf{y})$ and the diverge theorem, it follows:
\begin{equation*}
    \big\langle - L_{yx} u^{(0)}(\mathbf{x}, \cdot) \big\rangle = \mathbf{0}
\end{equation*}
thus, the existence of solution for problem \ref{Order-1VectorPDE} is satisfied.

Now, by separation of variables and the second condition of lemma (\ref{ExistenceLemma}), a general solution to the system of equations (\ref{Order-1VectorPDE}) can be given by:
\begin{equation}
    \label{Order-1Ansatz}
    u^{(1)}(\mathbf{x},\mathbf{y}) = \mathbf{N}^{rs}(\mathbf{y}) \mathbf{e}_{rs,x}(v(\mathbf{x})) 
\end{equation}
where $N^{rs} \in \mathbf{H}^1 (\mathbf{Y})$ the so-called cell solution.

Then replacing (\ref{IndepencyofY}), (\ref{Order-1Ansatz}) into (\ref{Order-1VectorPDE}) results on the so-called cell problems. Explicitly, note that replacing the terms described before, it follows:
\begin{align*}
    &\,L_{yy} \big( \mathbf{N}^{rs} (\mathbf{y}) \mathbf{e}_{rs,x} (v(\mathbf{x})) \big) + L_{yx}(v(\mathbf{x}) ) = \mathbf{0} \\
    \implies& -\partial_{y_j} \big( C_{ijkl}\mathbf{e}_{kl,y}(\mathbf{N}^{rs}(\mathbf{y}) ) \mathbf{e}_{rs,x}(v(\mathbf{x})) = \partial_{y_j} \big( C_{ijrs}\mathbf{e}_{rs,x}(v(\mathbf{x})) \big) \\
    \implies& - \partial_{y_j} \big( C_{ijkl} \mathbf{e}_{kl,y} (\mathbf{N}^{rs}(\mathbf{y})) \big) = \partial_{y_j} (C_{ijrs})
\end{align*}
and applying lemma (\ref{ExistenceLemma}), it follows that $\mathbf{N}^{rs}$ is $Y$-periodic where it has been added the boundary conditions taken from (\ref{Order-1VectorPDE}), obtaining:
\begin{equation*}
    \left \{
    \begin{array}{cc}
         - \partial_{y_j} \big( C_{ijkl} \mathbf{e}_{kl,y} (\mathbf{N}^{rs}(t,\mathbf{y})) \big) = \partial_{y_j} (C_{ijrs}) & \text{ in } (0,T)\times \mathbf{Y} \\
        \big( \mathbf{C}(\mathbf{y}) : \mathbf{e}(\mathbf{N}^{rs}(t,\mathbf{y})) \big) \cdot n = \mathbf{0} & \text{ on } (0,T)\times \partial \mathbf{Y}\\
        \mathbf{N}^{rs} (0, \mathbf{y}) = \mathbf{0} &  \text{ in } \mathbf{Y}
    \end{array}
    \right.
\end{equation*}

the above deduced by construction using the ansatz (\ref{Order-1Ansatz}) define a solution to the PDE system (\ref{Order-1VectorPDE}).

\subsection{Contribution at the zero order}
For the contribution at order $\mathcal{O}(\epsilon^0)$ problem (\ref{MainPDE}) after regrouping the terms it follows:
\begin{equation}
    \label{Order-0VectorPDE}
    \begin{array}{cc}
        L_{yy} u^{(2)} = \rho^{\epsilon}(\mathbf{x})\partial_{tt} v(\mathbf{x}) - \tilde{P}(u^{(1)},u^{(2)}) (\mathbf{x},\mathbf{y}) & \text{ in } \Omega \times \mathbf{Y}
    \end{array}
\end{equation}
where it has been denoted 
\begin{equation*}
    \tilde{P}(u^{(0)}, u^{(1)}) (\mathbf{x},\mathbf{y}) :=  L_{xx} u^{(0)} + L_{xy} u^{(1)} + L_{yx} u^{(1)}
\end{equation*}
then by Lemma \ref{ExistenceLemma}, it follows existence of a $Y$-periodic solution of (\ref{Order-0VectorPDE}) iff 
\begin{equation}
    \label{Order-0ExistenceCond}
    \big \langle \rho(\mathbf{y}) \partial_{tt} v(\mathbf{x}) - L_{xx} u^{(1)} (\mathbf{x},\mathbf{y}) - L_{xy} u^{(1)}(\mathbf{x},\mathbf{y}) - L_{yx} u^{(1)}(\mathbf{x},\mathbf{y}) \big \rangle = \mathbf{0}
\end{equation}
From \ref{Order-0ExistenceCond} the homogenized equation is obtained, which can be deduced explicitly in the form:
\begin{align*}
    & \quad \langle \rho(\mathbf{y}) \rangle \partial_{tt} v(\mathbf{x}) + \langle L_{xx} u^{(0)}(\mathbf{x},\mathbf{y}) +L_{xy} u^{(1)}(\mathbf{x},\mathbf{y}) \rangle = \mathbf{0} \\
    \Leftrightarrow & \quad \langle \rho(\mathbf{y}) \rangle \partial_{tt} v(\mathbf{x}) + \langle L_{xx} u^{(0)}(\mathbf{x}, \mathbf{y}) + L_{xy}\big( \mathbf{N}^{rs}(\mathbf{y})\mathbf{e}_{rs,x}(v(\mathbf{x})) \big) = \mathbf{0} \\
    \Leftrightarrow & \quad \langle \rho\rangle \partial_{tt} v_i(\mathbf{x}) - \partial_{x_j}\big \langle \big(C_{ijrs}^{hom}(\mathbf{N}^{rs})\big) \mathbf{e}_{rs,x}(\hat{v}(\mathbf{x})) \big \rangle  = \mathbf{0} \quad \forall i \in \{1,\dots, d\}\\
    \Leftrightarrow & \quad \rho^{0} \partial_{tt} v_i(\mathbf{x}) - C_{ijkl}^{hom} \partial_{x_j} \mathbf{e}_{rs,x} (v(\mathbf{x})) = \mathbf{0} \quad \forall i \in \{1,\dots,d\}
\end{align*}
where it is denoted $\rho^{0} = \langle \rho \rangle$ and $C_{ijrs}^{hom}$ the so-called homogenized elastic coefficients defined by 
\begin{equation*}
    C_{ijrs}^{hom} = \big \langle C_{ijrs}(\mathbf{y}) + C_{ijkl}\mathbf{e}_{kl,y}\big(\mathbf{N}^{rs}(\mathbf{y})\big) \big \rangle 
\end{equation*}
With the above definition and the homogenized equation obtained from (\ref{Order-0VectorPDE}) it can be obtained the effective behavior of the system defined by:
\begin{equation}
    \label{HomogenizedPDE}
    \left \{
    \begin{array}{cc}
        \rho^{0} \partial_{tt} v - \nabla \cdot \sigma^{0} (v) = \mathbf{0} & \text{ in } (0,T) \times \Omega \\
        \sigma^{0}_{ij}(v) = C^{hom}_{ijkl}\mathbf{e}_{kl,x}(v) & \text{ in } (0,T)\times\Omega \\
        v = \mathbf{0} & \text{ on } (0,T)\times\Gamma_D \\
        \sigma^{0}(v) \cdot n = \mathbf{F} & \text{ on }(0,T)\times \Gamma_N
    \end{array}
    \right .
\end{equation}
In particular, the PDE problem (\ref{HomogenizedPDE}) is well-posed and the mechanical behavior maintains the linear elasticity property from the material.

Finally, it will be useful in the next chapters to recall that the homogenized elastic operator defines a bilinear form on $\mathbf{H}^1(\Omega, \Gamma_D)$ by:
\begin{equation*}
    a(u,v) := \int \limits_{\Omega} \sigma^{hom}_{ij}(u(\mathbf{x}) \partial_{x_j} v_i \, dx \quad u,v \in \mathbf{H}^1(\Omega, \Gamma_D)
\end{equation*}
moreover, since the homogenized coefficients are bounded and uniformly elliptic it can be applied the classical theorem of spectral decomposition \cite{raviart1983introduction}.
\begin{prop}
\label{EigenValuesProp}
Let $V \subset H$ Hilbert spaces, such that $V$ is dense and continuously embedded in $H$. Suppose the canonical injection from $V$ on $H$ is compact, and the bilineal form $a(\cdot, \cdot)$ is symmetric, V-elliptic. Then there exist an increasing sequence which tend to $+ \infty$ of eigenvalues
\begin{equation*}
    0 < \lambda_1 < \lambda_2  \leq \dots \leq \lambda_m \leq \cdots 
\end{equation*}
and a hilbertian orthonormal base of $H$ given by eigenvector $w_m$ such that:
\begin{equation*}
    \forall v \in V, \quad a(w_,, v) = \lambda_m (w_m, v), \quad \forall m = 1, 2, \dots
\end{equation*}
\end{prop}



\section{Justification of the homogenization}
Having obtained the homogenized model, in this section it is developed the justification procedure of convergence for the solution in the space $\mathbf{H}^1(\Omega, \Gamma_D)$. It corresponds to an extension of the work done by \cite{oleinik1992mathematical} in the static case. To handle easily the derivatives, it will be reformulated the problem (\ref{MainPDE}) to a canonical form described below\footnote{The choice of a canonical formulation for the mixed boundary elastodynamic model is based on straightforward usage of space derivatives and estimations to assure the justification of the asymptotic solution, in the sense that the two-scale approximate solution $u^{\epsilon}(t, \mathbf{x})$ converge to the real (experimental solution) $u(t, \mathbf{x})$ in some space with enough regularity.}.

Over the bounded smooth domain $\Omega \subset \mathbb{R}^d$ it is considered the multiscale problem of second order in time:
\begin{equation}
    \label{MainMultiPDE}
    \left \{
    \begin{array}{cc}
        \mathcal{L}_{\epsilon}(u^{\epsilon}) := \rho\big( \frac{\mathbf{x}}{\epsilon} \big) \partial_{tt} u^{\epsilon} - \partial_{x_h} \big( A^{hk}\big( \frac{\mathbf{x}}{\epsilon} \big) \partial_{x_k} u^{\epsilon} \big)  = \mathbf{0} & \text{ in } (0,T)\times \Omega  \\
        \sigma^{\epsilon}(u^{\epsilon})\cdot n := A^{hk}\big( \frac{\mathbf{x}}{\epsilon} \big) \partial_{x_k} u^{\epsilon} n_k  = \mathbf{F}(t,\mathbf{x}) & \text{ on } (0,T) \times \Gamma_N \\
        u^{\epsilon} =  \mathbf{0} & \text{ on } (0,T) \times \Gamma_D \\
        \partial_t u^{\epsilon} = u^{\epsilon} = \mathbf{0} & \text{ on } \{ t=0 \} \times \Omega 
    \end{array}
    \right.
\end{equation}
being $n$ the unit outward normal to $\partial \Omega$, the disjoint decomposition $\partial \Omega = \Gamma_D \dot \cup \Gamma_N$ and assumed enough regularity for the boundaries (at least \textit{Lipschitz}). It is also assumed that the coefficients matrices $A^{hk}(\mathbf{y})$ are smooth functions, $1$-periodic that satisfy conditions of uniformly ellipticity and boundness in $\mathbf{y}$.\\
Recall also from (\ref{HomogenizedPDE}) the homogenized mixed boundary problem, rewritten in canonical form is given by:
\begin{equation}
    \label{HomMultiPDE}
    \left \{
    \begin{array}{ccc}
        \mathcal{L}_0 (u^0) := \rho^0 \partial_{tt} u^0 - \partial_{x_h}\big( A^{hk}_{hom} \partial_{x_k} u^0 \big) = \mathbf{0} & \text{ in } (0,T)\times \Omega \\
        \sigma^0(u^0) \cdot n := A^{hk}_{hom} \partial_{x_k}u^0 n_h = \mathbf{F}(t, \mathbf{x}) & \text{ on } (0,T) \times \Gamma_N \\
        u^0 = \mathbf{0} & \text{ on } (0,T) \times \Gamma_D \\
        \partial_t u^0 = u^0 = \mathbf{0} & \text{ on } \{ t=0 \} \times \Omega
    \end{array}
    \right .
\end{equation}
where the matrices $A^{hk}_{hom}$ (denoting the homogenized material coefficients) are defined by the formulas:
\begin{equation}
    \label{HomCoeffs}
    A^{pq}_{hom} = \int \limits_{\mathbf{Y}} A^{pq}(\mathbf{y}) + A^{pj} (\mathbf{y}) \partial_{y_j} \mathbf{N}^q (\mathbf{y}) \, d\mathbf{y}
\end{equation}
being the matrices $\mathbf{N}^q(\mathbf{y})$ the so-called cell solutions of the following boundary value problem:
\begin{equation}
    \label{CellProblems}
    \left \{
    \begin{array}{cc}
        \partial_{y_k} \big( A^{kj}(\mathbf{y}) \partial_{y_j} \mathbf{N}^q \big) = -\partial_{y_k} \big( A^{kq}(\mathbf{y}) \big) & \text{ in } \mathbf{Y} \\
        \mathbf{N}^q(\mathbf{y}) \quad 1\text{-periodic in } \mathbf{y}, & \int_{\mathbf{Y}}  \mathbf{N}^q (\mathbf{y}) \, d\mathbf{y} = \mathbf{0}
    \end{array}
    \right .
\end{equation}
where $\mathbf{Y} = (0,1)^d$ is the characteristic microstructure and $d \in \mathbb{N}^*$ denotes the dimension, usually $2$ or $3$.\\

\begin{rem}
 Let us note from (\ref{CellProblems}) that by applying  coercivity and boundness conditions over the coefficients $A^{kj}(\mathbf{y})$, it follows continuity and boundness conditions on (\ref{HomCoeffs}). Thus the homogenized problem (\ref{HomMultiPDE}) is of elastic type, and moreover it can be applied the same proof as in the case of prop. \ref{BootstrapingProp} obtaining by bootstrapping method:
 \begin{equation*}
     u^0, \, \partial_{t} u^0 \in \mathcal{C}^p (0,T; \mathbf{H}^1(\Omega, \Gamma_D))
 \end{equation*}
\end{rem}

It is assumed an approximate solution to the problem (\ref{MainMultiPDE}) in the form:
\begin{equation}
    \label{Asymptotic}
    \tilde{u}(t,\mathbf{x}) = u^0 (t,\mathbf{x}) + \epsilon \mathbf{N}^{s} \big(\frac{\mathbf{x}}{\epsilon} \big) \partial_{x_s} u^0(t,x)
\end{equation}
where $u^0$ is the solution to problem (\ref{HomMultiPDE}) and $\mathbf{N}^{s} (\mathbf{y})$ solutions to the so-called cell problems (\ref{CellProblems}).

%The function $\varphi(\mathbf{x})$ denotes the truncation function satisfying the conditions 
%\begin{enumerate}
%    \item $\varphi \in \mathcal{C}^{\infty}(\Omega)$ with $\vert \nabla \varphi\vert \lesssim \epsilon^{-1}$,
%    \item $\varphi = 0$ in $\Gamma_D$ and $\varphi = 1$ outside an $\epsilon$-neighborhood of $\Gamma_D$.
%\end{enumerate}
%Such type of functions exists, in particular, we can consider the function $\varphi(\mathbf{x}) := w(\epsilon^{-1} \rho(\mathbf{x},\partial \Omega))$ being

%\begin{equation*}
%    w(t) = 
%    \left \{
%    \begin{array}{cc}
%        t & \text{ if } 0 \leq t \leq 1 \\
%        1 & \text{ if } t > 1
%    \end{array}
%    \right .
%\end{equation*}

%\begin{rem}
%The truncation function $\varphi$ enters in the expression of $\tilde{u}$ since the matrices $\mathbf{N}^s \big(\frac{\mathbf{x}}{\epsilon}\big)$ are in general not defined in a neighborhood of $\Gamma_D$ (the other boundaries $\Gamma_N$ are associated with \textit{Neumann} type boundary conditions).
%\end{rem}

Now, let us enunciate the main result, where the sketch of proof will be presented in the next section.
\begin{theo}
Assuming that $\mathbf{F} \in L^{\infty}(0,T;\mathbf{L}^{2}(\Gamma_N))$ and also $\mathbf{N}^q \in L^{\infty}(0,T; \mathbf{H}^1(\Omega))$ for each $q \in \{1,\dots, N\}$. Then the solutions $u^{\epsilon}$ and $u^0$ of problems (\ref{MainMultiPDE}) and (\ref{HomMultiPDE}) respectively, satisfy the following inequality
\begin{equation*}
    \label{MainInequality}
    \Vert u^{\epsilon} - u^0 - \epsilon \mathbf{N}^s \big(\frac{\mathbf{x}}{\epsilon} \big) \partial_{x_s} u^0 \Vert_{L^{\infty}(0,T; \mathbf{H}^1(\Omega, \Gamma_D))} \lesssim \epsilon^{1/2} \Vert u^0 \Vert_{L^{\infty}(0,T; \mathbf{H}^{3}(\Omega))}
\end{equation*}
\end{theo}

\section{Sketch of the Proof}
The result is constructed using the work done on \cite{oleinik1992mathematical}, extending it from the elastostatic case to the elastodynamic model with mixed boundary conditions by controlling the behavior of each time-dependent that appears in the deduction. In this sense, the proof corresponds to an extension of the work done by \textit{Oleinik}.

The idea will be to study at fixed time $t \in (0,T)$ the elastodynamic operator $\mathcal{L}_{\epsilon}$.
Let us apply the operator $\mathcal{L}_{\epsilon}$ to a vector valued function $u^{\epsilon}-\tilde{u}$, where $\tilde{u}$ is defined by (\ref{Asymptotic}), then:
\begin{align*}
    \mathcal{L}_{\epsilon} (u^{\epsilon}-\tilde{u}) =&\, \rho^{\epsilon} \partial_{tt} (u^{\epsilon}-\tilde{u}) - \partial_{x_h} \big( A^{hk}\partial_{x_k} u^{\epsilon} \big) + \partial_{x_h} \big( A^{hk}\partial_{x_k} (u^0 + \epsilon \mathbf{N}^s \partial_{x_s}u^0 )\big) \\
    \overset{(*)}{=}& - (\rho^{\epsilon}-\rho^0) \partial_{tt} u^0 - \epsilon \rho^{\epsilon} \partial_{tt}\big(\mathbf{N}^s \partial_{x_s}u^0 \big)  \\
    & - \partial_{x_h} \big[ \big( A_{hom}^{hk} - A^{hk} - \epsilon A^{hj}\partial_{x_j} \mathbf{N}^k \big) \partial_{x_k} u^0 \big] + \epsilon \partial_{x_h} \big( A^{hk} \mathbf{N}^s \partial_{x_k x_s} u^0 \big) 
\end{align*}
where in $(*)$ it was rewritten, relabeled and separated term with $\epsilon$ order.
Taking into account the equation in (\ref{CellProblems}) for the $\mathbf{N}^s$ it follows:
\begin{align*}
    \mathcal{L}_{\epsilon} (u^{\epsilon} - \tilde{u})  :=& - (\rho^{\epsilon}-\rho^0) \partial_{tt} u^0 - \epsilon \rho^{\epsilon} \partial_{tt}\big(\mathbf{N}^s \partial_{x_s}u^0 \big) + \epsilon A^{hk}\mathbf{N}^s \partial_{x_h x_k x_s}u^0 \\
    & - \big[ A^{pq}_{hom} - A^{pq} - A^{pj} \partial_{y_j} \mathbf{N}^q - \partial_{y_h}(A^{hp} \mathbf{N}^q) \big] \partial_{x_p x_q} u^0 
    \end{align*}

From the above expression, let us define the matrices $\mathbf{N}^{hk}(\mathbf{y})$ as weak solutions of the following boundary value problem
\begin{equation}
    \label{SecondCellProblem}
    \left \{
    \begin{array}{cc}
        \partial_{y_k} \big( A^{kj}\partial_{y_j} \mathbf{N}^{pq}\big) = -\partial_{y_k} \big( A^{kp} \mathbf{N}^q\big) - A^{pj}\partial_{y_j} \mathbf{N}^q - A^{pq} + A^{pq}_{hom} & \text{ in } \mathbf{Y}\\
        \sigma(N^{pq})\cdot n := A^{kl}\partial_{x_k}\mathbf{N}^{pq}n_l = - n_k A^{pk}\mathbf{N}^q & \text{ on } \partial \mathbf{Y} \\
        \mathbf{N}^{hk}(\mathbf{y}) \text{ $\mathbf{Y}$-periodic}, \quad  \int_{\mathbf{Y}} N^{pq} d \mathbf{y} = \mathbf{0} &
    \end{array}
    \right .
\end{equation}

So, using (\ref{SecondCellProblem}) it can be obtained:
\begin{align*}
    \mathcal{L}_{\epsilon} (u^{\epsilon} -\tilde{u})  = & - (\rho^{\epsilon}-\rho^0) \partial_{tt} u^0 - \epsilon \rho^{\epsilon} \partial_{tt}\big(\mathbf{N}^s \partial_{x_s}u^0 \big)  \\
    & + \epsilon A^{hk} \mathbf{N}^s \partial_{x_h x_k x_s} u^0 + \epsilon A^{kj} \partial_{y_j} \mathbf{N}^{pq} \partial_{x_p x_q x_k} u^0  \\
    & - \epsilon \partial_{x_k} \big[ A^{kj} \partial_{y_j} \mathbf{N}^{pq} \partial_{x_p x_q} u^0\big] 
\end{align*}

Thus, it follows:
\begin{equation}
    \mathcal{L}_{\epsilon} (u^{\epsilon}-\tilde{u}) = G_{tt}^0 + \epsilon G_{tt}^1 + \epsilon F_0 + \epsilon \partial_{x_k} F_k
\end{equation}
where each term is defined by:
\begin{equation}
    \label{Variables}
    \begin{aligned}
        G_{tt}^0 = & -(\rho^{\epsilon} - \rho^0) \partial_{tt}u^0 \\
        G_{tt}^1 = & -\rho^{\epsilon} \partial_{tt}(\mathbf{N}^s \partial_{x_s} u^0) \\
        F_0 = & \, A^{hk}N^s \partial_{x_h x_k x_s} u^0  + A^{kj}\partial_{y_j} \mathbf{N}^{pq} \partial_{x_p x_q x_k} u^0 \\
        F_k = & - A^{kj} \partial_{y_j} \mathbf{N}^{pq} \partial_{x_p x_q} u^0 
    \end{aligned}
\end{equation}

\begin{rem}
Let us recall the two-scale asymptotic solution proof, by the hypothesis of continuity for the $\mathcal{O}(1)$ term in the form:
\begin{equation*}
    v \in \mathcal{C}^2(0,T; \mathbf{H}^1(\Omega, \Gamma_D))
\end{equation*}
It follows that terms $G^0_{tt}, G_{tt}^1$ are bounded.
\end{rem}

Let us now obtain a expression for the \textit{Neumann} condition at $\Gamma_N$. To this end, applying $\sigma_{\epsilon}$ to the difference $u^{\epsilon} - \tilde{u}$ it follows that:
\begin{align*}
    \sigma_{\epsilon} (u^{\epsilon} - \tilde{u}) & =  \sigma_{\epsilon}(u^{\epsilon}) - \sigma_{\epsilon} (\tilde{u}) \\
    & = F_h(t) n_h - A^{hk} \partial_{x_k} \big( u^0 + \epsilon N^s \partial_{x_s}u^0 \big)n_h \\
    & = \big(A^{hk}_{hom} - A^{hk} \big) \partial_{x_k} u^0 n_h - A^{hl} \partial_{y_l} N^s \partial_{x_s} u^0 n_h - \epsilon A^{hk} N^s \partial_{x_k x_s} u^0 n_h
\end{align*}

Taking into account the obtained expression, let us define:
\begin{equation*}
    \alpha^{is}(\mathbf{y}) := A^{is}_{hom} - A^{is}(\mathbf{y}) - A^{ij}(\mathbf{y}) \partial_{y_j} \mathbf{N}^s(\mathbf{y}), \quad i,s \in \{1,\dots, d\}
\end{equation*}
That allows us to obtain the expression at the \textit{Neumann} boundary in the form:
\begin{equation}
    \label{NeumannExp}
    \begin{aligned}
    \sigma_{\epsilon} (u^{\epsilon}-\tilde{u}) &= \alpha^{hs} \partial_{x_s} u^0 n_h - \epsilon A^{hk}N^s \partial_{x_k x_s} u^0 n_h \\
    & = I_0 + \epsilon I_1
    \end{aligned}
\end{equation}
where $I_0 = \alpha^{hs} \partial_{x_s} u^0 n_h$ and $I_1= -A^{hk}N^s \partial_{x_k x_s} u^0 n_h$. 
The function $\alpha^{hs}$ has already been studied in the similar context by \textit{Oleinik} \cite{oleinik1992mathematical}, who shown a continuity property: For each $v \in \mathbf{H}^1(\Omega)$ 
\begin{equation}
    \label{OleinikLemma2.2}
    \left \vert \int_{\partial \Omega} \alpha^{ik}(\frac{\mathbf{x}}{\epsilon}) v_k n_i \,ds \right \vert \lesssim\epsilon^{1/2} \Vert \nabla v \Vert_{\mathbf{L}^2(\Omega)}
\end{equation}
which applies to our case of functions $v \in \mathbf{H}^1(\Omega, \Gamma_D)$, whereas the $I_2$ term can be shown to satisfy similarly the bound:
\begin{equation}
    \label{I1-bound}
    \left \vert \int_{\partial \Omega} I_2 \right \vert \lesssim \epsilon^{1/2} \Vert \nabla u_0 \Vert_{\mathbf{L}^2(\Omega)}
\end{equation}
The intuition behind  (\ref{I1-bound}) relies in using that $\partial \Omega \subset \partial B_{\epsilon}(\Omega)$ where $B_{\epsilon}(\Omega) = \{ \mathbf{x} \in \Omega \, \vert \, \rho(\mathbf{x}, \Omega) \leq \epsilon \}$ from which the $\epsilon^{1/2}$ bound follows by applying \textit{Lemma 2.2} from \cite{oleinik1992mathematical} with the uniform bound on $\mathbf{Y}$ of $A^{hk}$ and $N^s$.


On the other hand, at the boundary $\Gamma_D$ for the difference $u^{\epsilon} - \tilde{u}$ it follows:
\begin{equation*}
    u^{\epsilon} - \tilde{u} = - \epsilon N^s \partial_{x_s}u^0 \equiv \psi_{\epsilon}
\end{equation*} 
Such kind of function satisfies moreover
\begin{equation}
    \label{BoundDirichlet}
    \Vert \psi_{\epsilon} \Vert_{\mathbf{H}^{1/2}(\Gamma_D)} \lesssim \epsilon^{1/2} \Vert u^0 \Vert_{\mathbf{H}^2(\Omega)}
\end{equation} 
To prove (\ref{BoundDirichlet}), it suffices to find a function $\Psi_{\epsilon} \in H^1(\Omega)$ such that $\Psi_{\epsilon} \vert_{\Gamma_D} = \psi_{\epsilon} \in \mathbf{H}^{1/2}(\Gamma_D)$ with $\Vert \Psi_{\epsilon} \Vert_{\mathbf{H}^1(\Omega)} \lesssim \epsilon^{1/2} \Vert u^0 \Vert_{\mathbf{H}^2(\Omega)}$.

Let us define such $\Psi_{\epsilon}(\mathbf{x})$ using a cutoff function. To this end, let $\varphi_{\epsilon} \in \mathcal{C}^{\infty}(\Omega)$ be such that:
\begin{equation*}
    \left \{
    \begin{aligned}
    \varphi_{\epsilon}(\mathbf{x}) = 1 & \text{ if } \rho(\mathbf{x}, \Gamma_D) \leq \epsilon \\
    \varphi_{\epsilon}(\mathbf{x}) = 0 & \text{ if } \rho(\mathbf{x}, \Gamma_D) \geq 2\epsilon 
    \end{aligned}
    \right.
\end{equation*}
with $0 \leq \varphi_{\epsilon}(\mathbf{x}) \leq 1$ and $\vert \nabla \varphi_{\epsilon} \vert \lesssim \epsilon^{-1}$.
Then, by construction it follows that $\Psi_{\epsilon} := \epsilon \varphi_{\epsilon} N^s \partial_{x_s}u^0 \in \mathbf{H}^1(\Omega)$.
Thus, from $\varphi_{\epsilon}$ and the uniformly boundness in $\epsilon$ of $N^s, \partial_{y_j} N^s$ for each $s,j \in \{1,\dots, d\}$ it follows:
\begin{equation*}
    \Vert \Psi_{\epsilon} \Vert_{H^1(\Omega)} \leq \Vert u^0 \Vert_{H^1(B_{\epsilon}(\Gamma_D))} + \epsilon \Vert u^0 \Vert_{\mathbf{H}^2(\Omega)}
\end{equation*}
being $B_{\epsilon}(\Gamma_D) = \{\mathbf{x}\in \Omega \, :\, \rho(\mathbf{x},\Gamma_D) \leq 2\epsilon \}$. From which it can be deduced the desired result after applying the  result (Lemma 1.5 at \cite{oleinik1992mathematical}). %$\Vert u^0 \Vert_{\mathbf{H}^1(\Omega^{\epsilon})} \lesssim \epsilon^{1/2} \Vert u^0 \Vert_{\mathbf{H}^2(\Omega)}$.

From the developments done, a mixed boundary evolution problem for the variable $u^{\epsilon} - \tilde{u}$ is obtained given by:
\begin{equation}
    \label{MixedFormulationDiff}
    \left \{
    \begin{aligned}
        \mathcal{L}_{\epsilon}(u^{\epsilon} - \tilde{u}) = G_{tt}^0 + \epsilon G_tt^1 + \epsilon \partial_{x_k} F_k + \epsilon F_0   & \text{ in } (0,T) \times \Omega \\
        u^{\epsilon}-\tilde{u} = \psi_{\epsilon} & \text{ on } (0,T)\times \Gamma_D\\
        \sigma_{\epsilon}(u^{\epsilon}-\tilde{u}) = I_0 + \epsilon I_1 & \text{ on } (0,T)\times\Gamma_N
    \end{aligned}
    \right .
\end{equation}
where each term $G_{tt}^0, G_{tt}^1, F_k^0$ are uniformly bounded in $\epsilon$ and time, $I_0, I_1$ bounded in (\ref{OleinikLemma2.2}), (\ref{I1-bound}) respectively
and $\psi_{\epsilon}$ bounded in (\ref{BoundDirichlet}). Moreover, the operator $\mathcal{L}_{\epsilon}(\cdot)$ is of linear elastic type.
It follows that (\ref{MixedFormulationDiff}) defines a well-posed evolution problem, which can be solved as in the case of the previous multiscale problems by using the spectral decomposition method \cite{raviart1983introduction}. Thus, the existence of solution is guaranteed and the continuity condition with respect to the initial data implies in particular the bound
\begin{equation*}
    \Vert u^{\epsilon} - u^0 - \epsilon N^s \partial_{x_s}u^0 \Vert_{L^{\infty}(0,T; \mathbf{H}^1(\Omega))} \lesssim \epsilon^{1/2} \Vert u^0 \Vert_{L^{\infty}(0,T; \mathbf{H}^3(\Omega)}
\end{equation*}


\begin{rem}
Let us note from the above bound that for each $t \in (0,T)$ it can be deduced:
\begin{equation}
    \label{TimeBound}
    \Vert u^{\epsilon} - u^0 - \epsilon \varphi \mathbf{N}^s\big( \frac{\mathbf{x}}{\epsilon} \big) \partial_{x_s} u^0 \Vert_{\mathbf{H}^1(\Omega, \Gamma_N)} \lesssim \epsilon^{1/2} \Vert \mathbf{F}\Vert_{L^{\infty}(0,T;\mathbf{L}^{2}(\Gamma_N))}
\end{equation}
\end{rem}

\section{Energy Estimations}
The convergence results shown before contains information of the elastic operator behavior as $\epsilon$ varies, thus the convergence of its energy operator. The results is inspired from the energy result of \textit{Oleinik} \cite{oleinik1992mathematical}, extending it to a time-dependent model.

Let us define the energy operators $E_{\epsilon}$ and $E_{0}$ associated to the multiscale and homogenized problems defined in the form:
\begin{equation*}
    E_{\epsilon}(u^{\epsilon}) := \int_0^T \int_{\Omega} \frac{1}{2}u_i^{\epsilon}u_i^{\epsilon} + \int_0^T \int_{\Omega} \big( \partial_{x_j} u^{\epsilon}, A^{jk} \big( \frac{\mathbf{x}}{\epsilon}\big) \partial_{x_k} u^{\epsilon} \big)
\end{equation*}
and analogously:
\begin{equation*}
    E_0(u^0) := \int_0^T \int_{\Omega} \frac{1}{2}u_i^{0}u_i^{0} + \int_0^T \int_{\Omega} \big( \partial_{x_j} u^0, A^{jk}_{hom} \partial_{x_k} u^0 \big)
\end{equation*}
Thus, the following can be obtained
\begin{lem}
Assuming the hypothesis of the above theorem, we have the convergence of the energy operators, i.e., 
\begin{equation*}
    \vert E_{\epsilon} (u^{\epsilon}) - E_0 (u^0) \vert \lesssim \epsilon^{1/2} \Vert \mathbf{F} \Vert_{L^{\infty}(0,T; \mathbf{L}^{2}(\Gamma_N))} 
\end{equation*}
\end{lem}
\begin{proof}
Let us fix $t \in (0,T)$, it follows from the result (\ref{TimeBound}) the bounds for $u^{\epsilon}$ and $\partial_{x_i} u^{\epsilon}$ in the form:
\begin{equation*}
    u^{\epsilon} = u^0 + \epsilon \mathbf{N}^s\big( \frac{\mathbf{x}}{\epsilon} \big) \partial_{x_s} u^0 + r^{\epsilon}(\mathbf{x},t) 
\end{equation*}
where the remaining term is bounded for
\begin{equation*}
    \Vert r^{\epsilon}(\cdot, t) \Vert_{\mathbf{H}^1(\Omega; \Gamma_D)} \lesssim \epsilon^{1/2} \Vert \mathbf{F}\Vert_{L^{\infty}(0,T;\mathbf{L}^{2}(\Gamma_N))}
\end{equation*}
and for each $i \in \{1,2,3\}$:
\begin{equation*}
    \partial_{x_i} u^{\epsilon} = \partial_{x_i} u^0 + \epsilon \partial_{y_i} \mathbf{N}^s \big( \frac{\mathbf{x}}{\epsilon} \big) \partial_{x_s} u^0 + q^{\epsilon}_i(\mathbf{x},t)
\end{equation*}
where now, the remaining term satisfies 
\begin{equation*}
    \Vert q^{\epsilon}_i(\cdot, t) \Vert_{\mathbf{H}^1(\Omega, \Gamma_D)} \lesssim \epsilon^{1/2} \Vert \mathbf{F}\Vert_{L^{\infty}(0,T; \mathbf{L}^{2}(\Gamma_N))}
\end{equation*}
So, using the definition of the energy operator it follows that:
\begin{align*}
    E_{\epsilon}(u^{\epsilon}) &= \int_0^T \int_{\Omega} \frac{1}{2}(u^0 + \epsilon \mathbf{N}^s \big( \frac{\mathbf{x}}{\epsilon} \big) \partial_{x_s} u^0, u^0 + \epsilon \mathbf{N}^s \big( \frac{\mathbf{x}}{\epsilon} \big) \partial_{x_s} (u^0) \big) + R(r^{\epsilon}) \\
    & + \int_0^T \int_{\Omega} \big( \partial_{x_j} u^0 + \epsilon \partial_{x_j} \mathbf{N}^s \big( \frac{\mathbf{x}}{\epsilon} \big) \partial_{x_s} u^0, A^{jk} \big(\frac{\mathbf{x}}{\epsilon} \big) \big[\partial_{x_k} u^0 + \epsilon \partial_{x_k} \mathbf{N}^s \big( \frac{\mathbf{x}}{\epsilon}\big) \partial_{x_s} u^0 \big] \big) + Q(q^{\epsilon})
\end{align*}
where the remaining integral terms can be bounded by
\begin{equation*}
    \vert R(r^{\epsilon}) \vert \lesssim \epsilon^{1/2} \Vert \mathbf{F}\Vert_{L^{\infty}(0,T;\mathbf{L}^{2}(\Gamma_N))} 
\end{equation*}
and moreover
\begin{equation*}
    \vert Q(q^{\epsilon} \vert \lesssim \epsilon^{1/2} \Vert \mathbf{F}\Vert_{L^{\infty}(0,T;\mathbf{L}^{2}(\Gamma_N))}
\end{equation*}
since the functions $u^0, \partial_{x_k}u^0, \mathbf{N}^s, \partial_{x_k} \mathbf{N}^s$ are bounded.

Let us consider a base for the $d\times d$ matrices $(m_{ij})_{ij}$, using the notation $m_s = (m_{s1}, \dots, m_{sd})$. Thus, consider for each $s, t \in \{1,2,3\}$ the matrices 
\begin{equation*}
    H^{st}(\mathbf{y}) = \partial_{y_i}\big(\mathbf{N}^s + y_s m_s\big) A^{ij}(\mathbf{y}) \partial_{y_j} \big(\mathbf{N}^t +  y_t m_t\big) - A^{st}_{hom}
\end{equation*}
It follows then by the homogenized coefficient definition that $\langle H^{st} \rangle_{\mathbf{Y}} = \mathbf{0}$.
Using the definition of $H^{st}$ over the multiscale energy expression, it can be deduced that
\begin{equation*}
    E_{\epsilon}(u^{\epsilon}) - E_{0}(u^{0}) = \int_0^T \int_{\Omega} \frac{1}{2}(u^0, u^0) + \partial_{x_s}u^0 H^{st}(\frac{\mathbf{x}}{\epsilon}) \partial_{t}u^0  + \tilde{R}(r^{\epsilon}) + Q(q^{\epsilon})
\end{equation*}
being $\vert \tilde{R}(r^{\epsilon})\vert \lesssim \epsilon^{1/2}\Vert \mathbf{F} \Vert_{L^{\infty}(0,T;\mathbf{L}^2(\Gamma_N))}$, obtained using the same technique as in \textit{Theorem 1.3} \cite{oleinik1992mathematical}.
\end{proof}

\section{Continuity on the Porosity}
The definition of homogenized coefficients contains an explicit dependency on the porosity level since it defines the characteristic microstructure. In this section, it will be studied the continuity of the homogenized coefficients with respect to the porosity level addressing the well-posedness of the inverse problem formulation. Moreover, differentiability is also assessed\footnote{A similar study has been done for the multiscale \textit{Calderón} problem by with high-oscillation conductivity coefficients, oriented to the inverse problem identification of homogenized coefficients \textit{A. Abudlle and A. di Blasio} \cite{Abdulle:221051}} by means of the cell problem definition. 

More explicitly, it is given estimations for such effective coefficients with a continuity property under reasonable hypothesis on the porosity\footnote{Such reasonable hypothesis are oriented in a continuity condition on the definition of the multiscale coefficients $C_{ijkl}(\mathbf{y},p)$, for example of \textit{Lipschitz} type on $p$.}. To this end, let us recall that our formulation for the elasticity coefficients at a cell point $\mathbf{y}\in \mathbf{Y}$ is given by:
\begin{equation*}
    C_{ijkl}(\mathbf{y}, p) = C_{ijkl}^{m} \mathbb{I}_{\mathbf{y} \in \mathbf{Y}_m(p)} + C_{ijkl}^{f} \mathbb{I}_{\mathbf{y} \in \mathbf{Y}_f(p)}
\end{equation*}
being $\mathbf{Y}_m(p)$ and $\mathbf{Y}_f(p)$ the matrix and fluid-like parts of the cell, where the explicit dependency on the porosity fraction $p \in (0,1)$ implies \textit{a posteriori} also the dependence of $p$ on the homogenized coefficient.

\subsection{Definition and Estimates}
In this section is defined the space in which the highly oscillatory coefficients $C_{ijkl}(\cdot)$ belongs dependent on the porosity parameter.
\begin{defn}
Let us define the set of operators on the feasible porosity interval $(0,1)$, associated to a cell-domain $\mathbf{Y}$ by:
\begin{equation}
    \label{T-definition}
    \mathbf{T}((0,1) \times \mathbf{Y}) := \big \{ \mathbf{C}(\cdot, \cdot) \in L^{\infty}((0,1); \, \text{lin}(\textbf{Sym}^d))\, : \mathbf{C}(p, \mathbf{y})^* = \mathbf{C}(p, \mathbf{y}) \, \forall \mathbf{y} \in \mathbf{Y}\big \}
\end{equation}
where it is used the notation of $lin(\mathbf{Sym}^d)$ as the space of linear, continuous operators from the $d\times d$ symmetric matrices space with $\mathbf{L}^2(\mathbf{Y})$ entries into itself, being $d=2,3$ and $*$ denoting the adjoint. 
\end{defn}
\begin{rem}
Note that, the linear elastic second-order tensors $\mathbf{C} = (C_{ijkl})_{ijkl}$ belong to $\mathbf{T}((0,1)\times \mathbf{Y})$. In particular, the symmetry of the tensors expressed in \textit{Voigt} engineering notation described as $C_{IJ}=C_{JI}$ is given by the self-adjoint property, by identifying the set $lin(\mathbf{Sym}^d)$ with the space of arrays indexed in the form $\mathbf{C}(p, \cdot) = (C_{ijkl} (p, \cdot) )_{ijkl}$.
Over such considerations, the set (\ref{T-definition}) becomes:
\begin{equation*}
    \mathbf{T}((0,1)\times \mathbf{Y}) = \big \{ C_{ijkl}(\cdot, \mathbf{y}) \in L^{\infty}((0,1); \mathbb{R})  \,:\, C_{ijkl}(\mathbf{y}) = C_{ijlk}(\mathbf{y}) = C_{klij}(\mathbf{y}) \in L^2(\mathbf{Y})\big \}
\end{equation*}
where it is assumed implicitly the dependency on $p$ of each coefficient.
\end{rem}
With the above it can be defined the space of bounded and coercive tensors associated to parameters $(\alpha,\beta) \in \mathbb{R}^2_+$ given as a subset of $\mathbf{T}((0,1)\times \mathbf{Y})$ with conditions specified below.
Let $(C_{ijkl})_{ijkl} \in \mathbf{T}((0,1)\times \mathbf{Y})$, consider $\alpha, \beta$ two positive parameters and define the following properties:
\begin{enumerate}
    \item[(H1)] If $\xi \in \mathbf{R}^{d\times d} \setminus \{\mathbf{0}_{d \times d}\}$ then $\alpha \xi_{ij} \, \xi_{ij} \leq C_{ijkl}(p,\mathbf{y})  \xi_{kl}\,\xi_{ij}\, \text{a.e.} \mathbf{y} \in \mathbf{Y}$ for each $p \in (0,1)$.
    \item[(H2)] If $\xi \in \mathbf{R}^{d\times d} \setminus \{\mathbf{0}_{d \times d}\}$ then $C_{ijkl}(p,\mathbf{y}) \xi_{kl} \xi_{ij} \leq \beta  \xi_{ij} \xi_{ij} \, \text{ a.e. }\mathbf{y} \in \mathbf{Y}$ for each $p \in (0,1)$.
\end{enumerate}
The above enables to consider the suitable space of linear elastic tensors, related to the generalized \textit{Hookes} law.
\begin{defn}
A tensor $\mathbf{C} = (C_{ijkl})_{ijkl} \in \mathcal{T}(\alpha, \beta, (0,1)\times \mathbf{Y})$
if $\mathbf{C}$ is a feasible tensor, i.e., belongs to $\mathbf{T}((0,1)\times \mathbf{Y})$ and satisfies the properties $(H1)-(H2)$.
\end{defn}
\begin{rem}
In the two-scale elastodynamic model formulation, the linear elastic tensor is expressed on the macroscopic (or slow) variable by
\begin{equation*}
    \mathbf{C}^{\epsilon}(p,\mathbf{x}) =\mathbf{C}(p,\frac{\mathbf{x}}{\epsilon})
\end{equation*} being $\mathbf{x}\in \Omega$, in such a way that $\mathbf{C} \in \mathcal{T}(\alpha, \beta, (0,1)\times \mathbf{Y})$. Recall also that, the $\epsilon$ parameter indicates the relation of micro-to-macro scale (\ref{TimeDif-Scheme}).
\end{rem}
Following the two-scale homogenization procedure studied in the sections before, the homogenized problem governing the macroscopic mechanical behavior is of elastic type with effective coefficients defined at porosity $p \in (0,1)$ by:
\begin{equation}
    \label{eq:homogenized-coeff}
    C^{hom}_{ijrs}(p) = \frac{1}{\vert \mathbf{Y}\vert} \int\limits_{\mathbf{Y}} C_{ijrs}(p,\mathbf{y}) \, dy + \frac{1}{\vert \mathbf{Y}\vert} \int\limits_{\mathbf{Y}} C_{ijkl}(p,\mathbf{y}) \mathbf{e}_{kl,y}( N^{rs}(p))\,dy
\end{equation}
for each $ i,j,r,s \in \{1,2,3\}$, where the vector functions $(N^{rs})_{rs}$ are defined as the unique $\mathbf{Y}$-periodic solution to the cell problems, satisfying for each $i \in \{1,2,3\}$, the equality:
\begin{equation}
    \label{eq:cell}
    -\partial_{y_j} \big[ C_{ijkl}(p,\mathbf{y}) \mathbf{e}_{kl}(N^{rs}(p,\mathbf{y}) \big] = \partial_{y_j} \big[ C_{ijrs}(p,\mathbf{y}) \big] 
\end{equation}
with the normalization condition $\langle N^{rs} \rangle_{\mathbf{Y}} = \mathbf{0}$. Such a condition is necessary to obtain a well-posed elliptic problem on the space $\mathbf{H}^1(\mathbf{Y})$.
Explicitly, vector solutions $N^{rs}(\cdot,p)$ to the cell problem are found on   
\begin{equation*}
    \mathbf{H}^1_{ 0} (\mathbf{Y}) := \big \{ N(p) \in  \mathbf{H}^1(\mathbf{Y}) \, : \, \langle N(p) \rangle_{\mathbf{Y}}=\mathbf{0} \, \big \}
\end{equation*}

\subsubsection{Variational Formulations}
Solutions to problem (\ref{eq:cell}) are found via variational formulation of such problem for each $r,s \in \{1,2,3\}$. In all that follows, it is assumed a condition of radial symmetry on the coefficients\footnote{Such kind of condition, which is naturally presented in periodic and compound media is assumed throughout this study, moreover the numerical results developed are obtained over a formulation of bone satisfying that property since the model is defined as compounds material with the inclusion being of different kind where the surrounding is a known fixed material. It is important to observe that under such assumptions the term associated to $\partial \mathbf{Y}$ vanishes obtaining the expression observed on (\ref{eq:cell_variational}).} $C_{ijkl}(\mathbf{y})$ which applied to (\ref{eq:cell}) gives us the desired equivalent problem:
\begin{equation}
    \label{eq:cell_variational}
    \left \{
    \begin{array}{cc}
        \text{Find a function } N^{rs} \in \mathbf{H}^1_{0}(\mathbf{Y}) \text{ such that } & \quad \\
        \int\limits_{\mathbf{Y}} C_{ijkl}(p,\mathbf{y})\mathbf{e}_{kl,y} (N^{rs}) \, \partial_{y_j}(v_i)\,d\mathbf{y} = -\int\limits_{\mathbf{Y}} C_{ijrs}(p, \mathbf{y}) \partial_{y_j}(v_i)\,d\mathbf{y} &  \forall v \in \mathbf{H}^1_{0}(\mathbf{Y}).
    \end{array}
    \right.
\end{equation}
Note the variational formulation (\ref{eq:cell_variational}) can be rewritten by the symmetry of the operator $\mathbf{C}(\cdot,p)$ in the form:
\begin{equation}
    \label{eq:cell_symmetric}
    \left \{
    \begin{array}{cc}
        \text{Find a function } N^{rs}(p) \in \mathbf{H}^1_{0}(\mathbf{Y}) \text{ such that } & \quad\\
        \int\limits_{\mathbf{Y}} C_{ijkl}(p,\mathbf{y})\mathbf{e}_{kl,y} (N^{rs})\, \mathbf{e}_{kl,y}(v)\,d\mathbf{y} = -\int\limits_{\mathbf{Y}} C_{ijrs}(p,\mathbf{y}) \mathbf{e}_{ij,y}(v)\,d\mathbf{y} &\forall v \in \mathbf{H}^1_{0}(\mathbf{Y}).
    \end{array}
    \right.
\end{equation}
Let us now define the following useful notation for the contraction of the last two indexes:
\begin{defn}
Let us consider $\mathbf{a} = (a_{ij})_{ij}$ and $\mathbf{b} = (b_{ij})_{ij}$ elements in $\mathcal{M}_{d\times d}(\mathbb{R})$, and consider $\mathbf{A}=(A_{ijkl})_{ijkl}$ a element of $\mathbf{T}([0,1]\times \mathbf{Y})$ identified with their multidimensional array. It is defined the contraction of the (last) two indices by
\begin{equation*}
    \mathbf{a}:\mathbf{b} := a_{ij}b_{ij}
\end{equation*}
and naturally extending the definition to three elements in the form:
\begin{equation*}
    \mathbf{A}:\mathbf{a}:\mathbf{b} := A_{ijkl}a_{kl}b_{ij}
\end{equation*}
\end{defn}
\begin{rem}
In particular, since $\mathbf{A} \in \mathbf{T}((0,1)\times \mathbf{Y})$ then it can be deduced the symmetry conditions $A_{ijkl}=A_{klij}=A_{ijlk}$, which allows to obtain a conmutative property between $\mathbf{a},\mathbf{b}$ two matrices:
\begin{equation}
    \mathbf{A}: \mathbf{a}:\mathbf{b}= \mathbf{A}:\mathbf{b}:\mathbf{a}
\end{equation}
\end{rem}

%%%%%%%%%%%%%%%%%%%%%%%%%5

\subsubsection{Symmetries}
As before, let us denote the base of $\mathbf{Sym}^d$ by the set $(m_{ij})_{ij}$ for each $i,j \in \{1,\dots d\}$, i.e., the set of $d \times d$ matrices with $1$ at the entry $(i,j)$ and $0$ otherwise. Observe that the homogenized coefficient (\ref{eq:homogenized-coeff}) (where it is assumed w.l.o.g $\vert \mathbf{Y} \vert = 1$) can be rewritten with the previously defined notation in the form:
\begin{equation}
    \label{eq:effective_coef}
    \begin{array}{ccc}
        C_{ijrs}^{hom}(p) &=& \int\limits_{\mathbf{Y}} \mathbf{C}(\mathbf{y},p):m_{rs}:m_{ij} + \int\limits_{\mathbf{Y}} \mathbf{C}(\mathbf{y},p):\mathbf{e}(N^{rs}):m_{ij} \, d\mathbf{y} \\
         &=&\int\limits_{\mathbf{Y}} \mathbf{C}(\mathbf{y},p):\big(m_{rs} + \mathbf{e}(N^{rs})\big):m_{ij} \, d\mathbf{y}
    \end{array}
\end{equation}

Moreover, expression (\ref{eq:cell_symmetric}) can be rewritten in the form:
\begin{equation}
    \label{eq:cell_v}
    \int\limits_{\mathbf{y}}\mathbf{C}(\mathbf{y},p):\mathbf{e}(v): \mathbf{e}(N^{rs}) \, d\mathbf{y} = - \int\limits_{\mathbf{y}} \mathbf{C}(\mathbf{y},p):\mathbf{e}(v):m_{rs}\, d\mathbf{y} \quad \forall v \in \mathbf{H}^1_{0}(\mathbf{Y}) \, d\mathbf{y}
\end{equation}

So, by using the above equality (\ref{eq:cell_v}) with $v = N^{ij}$ in the expression (\ref{eq:effective_coef}), it is obtained the homogenized coefficient in the form:
\begin{equation}
    \label{eq:effective_coef_symmetric}
    C_{ijkl}^{hom}(p) = \int\limits_{\mathbf{Y}} \mathbf{C}(\mathbf{y},p):(m_{rs}+\mathbf{e}(N^{rs})):(m_{ij}+\mathbf{e}(N^{ij})) \, d\mathbf{y}
\end{equation}
It follows then from (\ref{eq:effective_coef_symmetric}) the full symmetric properties of the elements $C^{hom}_{ijkl}$, since by (\ref{eq:cell_symmetric}) follows the that $N^{rs}(p,\mathbf{y}) = N^{sr}(p, \mathbf{y})$ thus taking into account the symmetry of $\mathbf{C}(p)$ it can be deduced:
\begin{equation*}
    C^{hom}_{ijkl} = C^{hom}_{klij} = C^{hom}_{ijlk} \quad \forall i,j,k,l \in \{1,2,3\}
\end{equation*}
which implies in particular the usage of engineering \textit{Voigt} notation for the homogenized equations, obtaining a standard elasticity tensor.

In all that follows, it will be denoted by $M(\mathbf{Y})$ the space of $d \times d$ matrices $L^2(\mathbf{Y})$ entries with norm $\Vert \cdot \Vert_{M(\mathbf{Y})}$ defined for $A = (A_{kl})_{kl} \in M(\mathbf{Y})$ by
\begin{equation}
    \Vert A \Vert_{M(\mathbf{Y})} = \big(\sum_{k,l=1}^d \Vert A_{kl} \Vert^2_{L^2(\mathbf{Y})} \big)^{1/2}
\end{equation}
Let us recall the coercivity condition over $\mathbf{C}(p, \cdot)$, applying a duality argument over the space $M(\mathbf{Y})$ it follows that:
\begin{equation}
    \label{eq:Estimate-Nrs}
    \Vert \mathbf{e}(N^{rs}(p))\Vert_{M(\mathbf{Y})} \leq \frac{\beta}{\alpha}
\end{equation}
which will be useful in the estimation in the following subsections. Observe with the above definition, the elastic tensors belong to $lin(M(\mathbf{Y}))$.
%%%%%%%%%%%%%%
\subsection{About the Continuity}

To study the continuity condition on the effective coefficients $C^{hom}_{ijkl}(\cdot)$, let $p, q \in (0,1)$ be two porosity levels. It will be first shown a useful estimation of $N^{rs}(p)-N^{rs}(q)$ for each $r,s \in \{1,2,3\}$. 
Using (\ref{eq:cell_symmetric}) rewritten in tensor notation, the following equality can be obtained that relates the cell solutions to their coefficients associated to the porosity level in the form:
\begin{equation}
    \label{eq:diff_Nrs}
    \begin{aligned}
        \int\limits_{\mathbf{Y}} \mathbf{C}(\mathbf{y},p):\mathbf{e}\big( N^{rs}(p)-N^{rs}(q) \big) : \mathbf{e}(v) & =  -\int\limits_{\mathbf{Y}} \mathbf{C}(p)-\mathbf{C}(q) : m_{rs}:\mathbf{e}(v) \\
        & \quad - \int\limits_{\mathbf{Y}}\mathbf{C}(p)-\mathbf{C}(q):\mathbf{e}(N^{rs}(q)):\mathbf{e}(v).
    \end{aligned}
\end{equation}
So that, using (\ref{eq:diff_Nrs}) it follows the representation $\forall v \in \mathbf{H}^1_{0}(\mathbf{Y})$
\begin{equation}
    \label{eq:diff_Nrs_full}
    \int\limits_{\mathbf{Y}} \mathbf{C}(p):\mathbf{e}\big(N^{rs}(p)-N^{rs}(q)\big):\mathbf{e}(v) = -\int\limits_{\mathbf{Y}} \mathbf{C}(p)-\mathbf{C}(q):m_{rs}+\mathbf{e}(N^{rs}(q)):\mathbf{e}(v).
\end{equation}
Applying then a duality argument on (\ref{eq:diff_Nrs_full}) since $M(\mathbf{Y})^* = M(\mathbf{Y})$ it is deduced the first estimate for the difference between the cell solutions, expressed by:
\begin{equation*}
    \Vert \mathbf{C}(p):\mathbf{e}\big(N^{rs}(p)-N^{rs}(q)\big) \Vert_{M(\mathbf{Y})} \leq \Vert \mathbf{C}(p)-\mathbf{C}(q):m_{rs}+\mathbf{e}(N^{rs}(q))\Vert_{M(\mathbf{Y})}
\end{equation*}
from which after taking H\"{o}lder inequality and recalling the coercivity of $\mathbf{C}(p)$ it is possible to obtain:
\begin{equation}
    \label{eq:ineq_diff_N}
    \Vert \mathbf{e}\big(N^{rs}(p)-N^{rs}(q)\big) \Vert_{M(\mathbf{Y})} \leq \alpha^{-1}\Vert \mathbf{C}(p)-\mathbf{C}(q) \Vert_{lin (M(\mathbf{Y}))} \Vert m_{rs}+\mathbf{e}(N^{rs}(q)) \Vert_{M(\mathbf{Y})}.
\end{equation}

Let us then estimate the difference between the homogenized coefficients associated to the porosities $p, q$. By applying the definition and rearranging terms it follows:
\begin{equation*}
    \begin{aligned}
        C^{hom}_{ijrs}(p) - C^{hom}_{ijrs}(q) = & \int_{\mathbf{Y}} \mathbf{C}(p):m_{rs} + \mathbf{e}(N^{rs}(p)):m_{ij} - \int_{\mathbf{Y}} \mathbf{C}(q):m_{rs} + \mathbf{e}(N^{rs}(q)):m_{ij}\\
        = & \int_{\mathbf{Y}} \mathbf{C}(p) - \mathbf{C}(q):m_{rs} + \mathbf{e}(N^{rs}(p)):m_{ij} \\
        & + \int_{\mathbf{Y}} \mathbf{C}(q):\mathbf{e}(N^{rs}(p)) - \mathbf{e}(N^{rs}(q)):m_{ij} 
    \end{aligned}
\end{equation*}
and by applying the H\"{o}lder inequality and using (\ref{eq:ineq_diff_N}) it follows:
\begin{equation}
    \label{eq:diff-HomCoeffs}
    \begin{aligned}
        \vert C^{hom}_{ijrs}(p) - C^{hom}_{ijrs}(q)\vert \leq & \, \Vert \mathbf{C}(p)-\mathbf{C}(q) \Vert_{lin(M(\mathbf{Y}))} \Vert m_{rs}+e(N^{rs}(p)) \Vert_{M(\mathbf{Y})} \\
         & +  \frac{\beta}{\alpha} \Vert \mathbf{C}(p)-\mathbf{C}(q) \Vert_{lin(M(\mathbf{Y}))} \Vert m_{rs}+\mathbf{e}(N^{rs}(q)) \Vert_{M(\mathbf{Y})} 
    \end{aligned}
\end{equation}

Finally, applying the estimation (\ref{eq:Estimate-Nrs}) over (\ref{eq:diff-HomCoeffs}) it follows the continuity property for the coefficients:
\begin{equation}
    \label{ContinuityPropHom}
    \vert C^{hom}_{ijrs}(p) - C^{hom}_{ijrs}(q)\vert \leq \Vert \mathbf{C}(p) - \mathbf{C}(q) \Vert_{lin(M(\mathbf{Y}))} (1+\frac{\beta}{\alpha})\frac{\beta}{\alpha}
\end{equation}

\subsection{About the Derivative}
Our interest now is to obtain some estimates for the derivatives with respect to the porosity $p$. By using the expression (\ref{eq:effective_coef_symmetric}) for fixed $i,j,r,s \in \{1,2,3\}$ and $\overline{p} \in (0,1)$ it can be deduced:
\begin{equation}
    \label{eq:DerivHomCoeff}
    \partial_{p}  C_{ijkl}^{hom} (\overline{p}) = \int\limits_{\mathbf{Y}} \partial_{p}\big( \mathbf{C}(\mathbf{y},\overline{p}) \big)\, : \,m_{rs}+\mathbf{e}(N^{rs}(\overline{p}))\, :\,m_{ij}+\mathbf{e}(N^{ij}(\overline{p})) \, d\mathbf{y} 
\end{equation}
which results by distributing each term, rules of differentiation and (\ref{eq:cell_v}).

It follows then from (\ref{eq:DerivHomCoeff}) after applying holder inequality that:
\begin{equation*}
    \vert \partial_p C^{hom}_{ijrs}(\overline{p}) \vert \leq \Vert \mathbf{C}(\overline{p}) \Vert_{lin(M(\mathbf{Y}))} \Vert m_{rs} + \mathbf{e}(N^{rs}(\overline{p})) \Vert_{M(\mathbf{Y})} \Vert m_{ij} + \mathbf{e}(N^{ij}(\overline{p})) \Vert_{M(\mathbf{Y})}
\end{equation*}

which after applying the estimate (\ref{eq:Estimate-Nrs}) the following inequality for the derivative can be deduced:
\begin{equation}
    \label{EstimateDerivHomCoeff}
    \vert \partial_p C^{hom}_{ijrs}(p) \vert \leq (1+\frac{\beta}{\alpha})^2 \Vert \partial_p \mathbf{C}\Vert_{L^{\infty}((0,1); lin(M(\mathbf{Y}))}
\end{equation}

The above estimates enable us to obtain a bound for the difference of the homogenized elastic tensors.

Now, let $\overline{p}, \overline{q}$ be two porosity level on $(0,1)$.
Applying the equality (\ref{eq:DerivHomCoeff}), rearranging terms and using directly H\"{o}lder inequality, it follows:
\begin{align*}
    \vert \partial_p C_{ijrs}(\overline{p}) - & \partial_p C_{ijrs}(\overline{q}) \vert \leq  \quad \\
    &\quad \Vert \partial_p \mathbf{C}(\overline{p}) - \partial_p \mathbf{C}(\overline{q}) \Vert_{lin(M(\mathbf{Y}))} \Vert m_{rs} + \mathbf{e}(N^{rs}(\overline{p})) \Vert_{M(\mathbf{Y})} \Vert m_{ji} + \mathbf{e}(N^{ij}(\overline{p}))\Vert_{M(\mathbf{Y})} \\
    & + \Vert \partial_p \mathbf{C}(\overline{q}) \Vert_{lin(M(\mathbf{Y}))} \Vert \mathbf{e}(N^{rs}(\overline{p})) -\mathbf{e}(N^{rs}(\overline{q})) \Vert_{M(\mathbf{Y})} \Vert m_{ij} + \mathbf{e}(N^{ij}(\overline{q}) \Vert_{M(\mathbf{Y})}  \\
    & + \Vert \partial_p \mathbf{C}(\overline{q}) \Vert_{lin(M(\mathbf{Y}))} \Vert m_{rs} + \mathbf{e}(N^{rs}(\overline{p})) \Vert_{M(\mathbf{Y})} \Vert \mathbf{e}(N^{ij}(\overline{p})) - \mathbf{e}(N^{ij}(\overline{q})) \Vert_{M(\mathbf{Y})}
\end{align*}
From which can be deduced taking into consideration the bounds for the cell solutions (\ref{eq:Estimate-Nrs}), (\ref{eq:ineq_diff_N}) that:
\begin{equation}
    \label{EstimateDiffDeriv}
    \begin{aligned}
        \vert \partial_p C_{ijrs}(\overline{p}) - \partial_p C_{ijrs}(\overline{q}) \vert & \,  \leq (1+ \frac{\beta}{\alpha})^2\Vert \partial_p \mathbf{C}(\overline{p}) - \partial_p \mathbf{C}(\overline{q}) \Vert_{lin(M(\mathbf{Y}))} \\
        & + 2\frac{\beta}{\alpha}(1+\frac{\beta}{\alpha})^2 \Vert \partial_p \mathbf{C}(\overline{q}) \Vert_{lin(M(\mathbf{Y}))} \Vert \mathbf{C}(\overline{p}) - \mathbf{C}(\overline{q}) \Vert_{lin(M(\mathbf{Y}))}
    \end{aligned}
\end{equation}

As conclusion, it is obtained the following properties:
\begin{prop}
Under the assumptions of regularity for the multiscale elastic coefficients $\mathbf{C} = (C_{ijkl})_{ijkl} \in \mathcal{T}(\alpha, \beta, (0,1)\times \mathbf{Y})$ of type:
\begin{enumerate}
    \item[HI] A Lipschitz continuity assumption for the elastic coefficients in the form:
    \begin{equation*}
        \Vert \mathbf{C}(\overline{p}) - \mathbf{C}(\overline{q}) \Vert_{lin(M(\mathbf{Y}))} \leq L \vert \overline{p} - \overline{q} \vert^{l} 
    \end{equation*}
    for some $L, l > 0$ positive constants, and each $\overline{p},\overline{q} \in (0,1)$.
\end{enumerate}
Then it follows for each $i,j,r,s \in \{1,2,3\}$ and $\overline{p},\overline{q}$ porosity levels:
\begin{enumerate}
    \item[PI] A continuity property for the homogenized coefficients:
    \begin{equation*}
        \vert C_{ijrs}^{hom}(\overline{p})-C^{hom}_{ijrs}(\overline{q}) \vert \leq L \vert \overline{p} - \overline{q} \vert^{l} \frac{\beta}{\alpha} (1+\frac{\beta}{\alpha})
    \end{equation*}
\end{enumerate}

Moreover, if it is satisfied:
\begin{enumerate}
    \item[HII] Uniformly boundness over the porosity range for the derivative, i.e., $\Vert \partial_p \mathbf{C}(\overline{q}) \Vert_{lin(M(\mathbf{Y}))} < \infty$ for each $\overline{q} \in (0,1)$.
    \item[HIII] Lipschitz continuity condition for the derivative, in the form:
    \begin{equation*}
        \Vert \partial_p \mathbf{C}(\overline{p})- \partial_p \mathbf{C}(\overline{q}) \Vert_{lin(M(\mathbf{Y}))} \leq \overline{L} \vert \overline{p}-\overline{q}\vert^{\overline{l}}
    \end{equation*}
    for $\overline{L}, \overline{l} > 0$ constants and each $\overline{p},\overline{q} \in (0,1)$.
\end{enumerate}
Then it follows that each $i,j,r,s \in \{1,2,3\}$ and $\overline{p},\overline{q}$ porosity levels:
\begin{enumerate}
    \item[PII] A continuity property for the derivatives of hom. coeffients:
    \begin{equation*}
        \vert \partial_p C_{ijrs}^{hom}(\overline{p})-\partial_p C_{ijrs}(\overline{q}) \vert \leq C(\alpha, \beta) \vert \overline{p}-\overline{q} \vert^{\overline{l}} + \tilde{C}(\alpha, \beta) \vert \overline{p}-\overline{q} \vert^{l}
    \end{equation*}
    being $C(\alpha,\beta), \tilde{C}(\alpha,\beta) > 0$ constants.
\end{enumerate}
\end{prop}