In order to find an exact geometrical solution, it is considered that the rigid body consists of an infinite number of material particles for which the distance between two points is maintained constant, regardless of the forces that act upon it. The elementary particles, characterized by elementary and infinitesimal mass, have a continuous distribution in the geometrical shape of the solid body.

When the density is constant inside the rigid structure, the solid rigid is homogeneous. When the integration limits applied on the geometrical outline are well defined, it results in a homogeneous body with a simple or regular shape. In this last case, the geometrical and mass integrals are applied.

An essential aspect in advanced dynamics is represented by mass properties [13] The mass and the position of the mass center is determined in relation with *O* 0*x* 0 *y* 0*z* 0 ≡ {0 } frame, as presented:

$$M = \int dm; \text{ where } dm = \{\rho\_{\mathcal{O}} \cdot d\sigma\}; \sigma = \{V; A; L\}, \tag{51}$$

$$
\overline{\rho}\_{\mathbb{C}}(t) = \frac{\int \overline{\rho}\_M(t) \cdot dm}{\int dm} = \frac{\int \overline{\rho}\_M(t) \cdot dm}{M},\tag{52}
$$

In the expression (51), ρσ defines the density of the material, where the symbol σ represents, by case, the volume (V), the area (A) or the length (L) of the rigid body subjected to study

The position of the mass center, with respect to the *O*0*x*0 *y*0*z*0 ≡ {0} reference frame, is expressed at first in a classical form, and then by using matrix exponentials as:

$$\begin{split} \overline{r}\_{\mathbb{C}}(t) &= \frac{\lceil \overline{r}\_{\mathbb{M}}(t) \cdot \text{du}}{\int \!\!\!/ \text{du}} = \frac{\lceil \overline{r}\_{\mathbb{M}}(t) \cdot \text{du}}{\text{M}} = \overline{r}\_{0}(t) + \overline{\rho}\_{\mathbb{C}}(t) = \overline{r}\_{0}(t) + \prescript{0}{}{\mathbb{S}}[\!\!R](t) \cdot \text{}^{\mathbb{S}}\overline{\rho}\_{\mathbb{C}} = \\ &= \overline{r}\_{0}(t) + \exp\left[\prod\_{\{\overline{\chi} = [u;v;w]\}} \exp[\overline{\chi} \times \delta\_{\chi}]\right] \cdot \text{}^{\mathbb{S}}\overline{\rho}\_{\mathbb{C}} = \overline{r}\_{0}(t) + \left\{\exp[\overline{u} \times a\_{u}] \cdot \exp[\overline{v} \times \beta\_{v}] \cdot \exp[\overline{w} \times \chi\_{v}]\right\} \cdot \text{}^{\mathbb{S}}\overline{\rho}\_{\mathbb{C}}. \end{split} \tag{53}$$

By applying the time derivative on (53), the linear velocity and acceleration of the mass center are determined:

$$\begin{split} \overline{v}\_{\mathcal{C}}(t) &= \dot{\overline{r}}\_{\mathcal{C}}(t) = \dot{\overline{r}}\_{0}(t) + \dot{\overline{\rho}}\_{\mathcal{C}}(t) = \dot{\overline{r}}\_{0}(t) + \frac{0}{\mathcal{S}}[\dot{R}](t) \cdot \frac{0}{\mathcal{S}}[\mathcal{R}]^{T}(t) \cdot \frac{0}{\mathcal{S}}[\mathcal{R}](t) \cdot {}^{\mathcal{S}}\overline{\rho}\_{\mathcal{C}} = \\ &= \overline{v}\_{0}(t) + \overline{w}(t) \times \overline{\rho}\_{\mathcal{C}}(t) = \frac{d}{dt} \Big[ \overline{r}\_{0}(t) + \left\{ \exp\Bigg[\prod\_{\{\boldsymbol{\gamma} = [w, \boldsymbol{\gamma}; w]\}} \exp[\overline{\boldsymbol{\chi} \times} \boldsymbol{\delta}\_{\boldsymbol{\lambda}}] \right] \right] \cdot {}^{\mathcal{S}}\overline{\rho}\_{\mathcal{C}} = \\ &= \frac{d}{dt} \Big[ \overline{r}\_{0}(t) + \{\exp[\overline{\boldsymbol{\mu}} \times \boldsymbol{\alpha}\_{\mathbb{I}}] \cdot \exp[\overline{\boldsymbol{\nu}} \times \boldsymbol{\beta}\_{\mathbb{I}}] \cdot \exp[\overline{\boldsymbol{\mu}} \times {}^{\mathcal{S}}\boldsymbol{\nu}\_{\mathbb{I}}] \cdot {}^{\mathcal{S}}\overline{\rho}\_{\mathcal{C}} \Big]. \end{split} \tag{54}$$

Using classical Formulations (52) and (53), the result is the expressions for the acceleration of mass center:

$$
\overline{a}\_{\mathbb{C}} = \dot{\overline{v}}\_{\mathbb{C}}(t) = \dot{\overline{v}}\_{0}(t) + \frac{d}{dt} [\overline{\omega}(t) \times \overline{\rho}\_{\mathbb{C}}(t)] = \overline{a}\_{0} + \overline{\varepsilon} \times \overline{\rho}\_{\mathbb{C}} + \overline{\omega} \times \overline{\omega} \times \overline{\rho}\_{\mathbb{C}} \tag{55}
$$

In dynamic modeling, the following expression is necessary:

$$\begin{array}{lcl}\frac{d}{dt}\Big(\overline{\rho}\_M(t)\Big) = \frac{d}{dt}\Big(\overline{\rho}\_\mathbb{C}(t) + \overline{\tau}'(t)\Big) = \frac{d}{dt}\Big(\prescript{0}{}{\mathbb{S}}[R](t) \cdot \left(\prescript{\mathbb{S}}{}{\overline{\rho}\_\mathbb{C}} + \prescript{\mathbb{S}}{}{\overline{\tau}'}\right)\Big) = \\\ = \overline{\omega}(t) \times \overline{\rho}\_M(t) = \overline{\omega}(t) \times \overline{\rho}\_\mathbb{C}(t) + \overline{\omega}(t) \times \overline{\tau}'(t). \end{array} \tag{56}$$

In the case of rigid bodies involved in rotation motion, the inertia property is highlighted by the mechanical moments of inertia [12–14]. According to Figure 2, the position of the elementary mass (*dm*) relative to the mass center is defined by means of the vector: *r* ∗ = *r* ∗ (*t*).

Considering (52), the following mass property is obtained:

$$\int \tilde{r}^\*(t) \cdot dm = \, ^0\_S[\mathbb{R}](t) \cdot \int \, ^S \tilde{r}^\* \cdot dm = 0. \tag{57}$$

Based on research [12,13], the inertial tensor and its variation law relative to the concurrent frames applied in the mass center: {*S*∗ } *and* {0∗ }*OR* ≡ {0 }*OR* is established, according to:

$$I\_S^\star = \int \left(\overline{r}^\star \times \right) \cdot \left(\overline{r}^\star \times \right)^T \cdot dm = \frac{0}{S} [R] \cdot \, ^S I\_S^\star \cdot \, ^0\_S [R]^T \, ^t \tag{58}$$

where *SI*∗ *S* is inertial tensor axial and centrifugal of the body (*S*) in with relation with {*S*∗ } frame, applied in the mass center (*C*), having the property: {*S*∗ }*OR* ≡ {*S*}*OR*:

$$^SI\_S^\* = \int \left(^S\vec{r}^\* \times \right) \cdot \left(^S\vec{r}^\* \times \right)^T \cdot dm,\tag{59}$$

Since the dynamic study refers to the absolute motion, the use of the following expressions of the inertial tensor are required:

$$I\_S' = {}^0\_S[R] \cdot {}^S I\_S \cdot {}^0\_S[R]^T = I\_{S\subset}' + I\_{S'}' \cdot \int \left( {}^S \overline{\rho}\_M \times \right) \cdot \left( {}^S \overline{\rho}\_M \times \right)^T \cdot dm = {}^S I\_S \text{ and } I\_{S\subset}' = M \cdot \left( \overline{\rho}\_\mathbb{C} \times \right) \cdot \left( \overline{\rho}\_\mathbb{C} \times \right)^T. \tag{60}$$

The matrix expression (60) represents the generalized variation law of the inertial tensor axial and centrifugal with respect to the frame {0 }. Further, in the expression (60), *I SC* is the inertia matrix axial and centrifugal of the mass center relative to {0 }. The inertial tensor axial and centrifugal is defined with respect to the absolute frame {0}*OR* = {0 }*OR*, as presented below:

$$I\_S = \int \left(\tilde{r}\_M \times \right) \cdot \left(\tilde{r}\_M \times \right)^T \cdot dm = M \cdot \left(\tilde{r}\_0 \times \right) \cdot \left(\tilde{r}\_0 \times \right)^T + I\_S' = I\_{SO} + I\_S' = I\_{SO} + I\_{S'}' + I\_{S'}' \tag{61}$$

In the Expression (61), *ISO* represents the inertia matrix containing axial and centrifugal inertia moments, which is determined with respect to the fixed system {0}.

#### **3. Higher Order Acceleration Energies**

The term, advanced notions, from analytical dynamics presented in this paper refers to the motion energies whose central functions are the higher order accelerations. They have been developing in any sudden and transitory motion of mechanical systems. By having Appell's function as the starting point [1,2], a function which was highlighted in 1899 and that is also known as the kinetic energy of accelerations, new mathematical formulations on the expressions for acceleration energies of the first, second, and third order have been developed [10,14]. In this section, they are presented only in an explicit form. From a physical point of view, the acceleration energy is expressed as - *J*/*s*2 . and represents the mechanical power developed by a mechanical system per time unit. By applying certain constraints, the acceleration energy can also be defined as the second order time derivative applied on the kinetic energy corresponding to the whole mechanical system. In the following, based on [10,14], the acceleration energies of the higher order are defined. They characterize the complex mechanical systems and include among others, the robots. The author established the acceleration energy corresponding to a rigid body involved in general motion, in a generalized form [10,14], also called the first order acceleration energy and defined with:

$$E\_A^{(1)} = \frac{1}{2} \cdot \int a\_M^2 \cdot dm = \frac{1}{2} \cdot \int \stackrel{?}{\overline{\nu}}\_M \cdot \stackrel{?}{\overline{\nu}}\_M \cdot dm = \frac{1}{2} \cdot \int \stackrel{?}{\overline{\nu}}\_M \cdot \stackrel{?}{\overline{\nu}}\_M \Big] \cdot dm, \quad \overline{a}\_M = \overline{a}\_0 + \varepsilon \times \overline{\rho}\_M + \overline{\omega} \times \overline{\omega} \times \overline{\rho}\_M \tag{62}$$

By developing Expression (71), the first three terms are obtained:

$$\frac{1}{2} \cdot \int \overline{a}\_0^T \cdot \overline{a}\_0 \cdot dm = \frac{1}{2} \cdot M \cdot \overline{a}\_0^T \cdot \overline{a}\_0 = \frac{1}{2} \cdot M \cdot a\_0^2 \tag{63}$$

$$\frac{1}{2} \cdot \int \overline{a}\_0^T \cdot \left(\overline{\boldsymbol{\varepsilon}} \times \overline{\boldsymbol{\rho}}\_M\right) \cdot dm = \frac{1}{2} \cdot \int \left(\overline{\boldsymbol{\varepsilon}} \times \overline{\boldsymbol{\rho}}\_M\right)^T \cdot \overline{a}\_0 \cdot dm = \frac{1}{2} \cdot M \cdot \overline{a}\_0 \cdot \left(\overline{\boldsymbol{\varepsilon}} \times \overline{\boldsymbol{\rho}}\_\mathbb{C}\right),\tag{64}$$

$$\frac{1}{2} \cdot \int \overline{a}\_0^T \cdot \left(\overline{\boldsymbol{w}} \times \overline{\boldsymbol{w}} \times \overline{\rho}\_M\right) \cdot dm = \frac{1}{2} \cdot \int \left(\overline{\boldsymbol{w}} \times \overline{\boldsymbol{w}} \times \overline{\rho}\_M\right)^T \cdot \overline{a}\_0 \cdot dm = \frac{1}{2} \cdot M \cdot \overline{a}\_0 \cdot \left(\overline{\boldsymbol{w}} \times \overline{\boldsymbol{w}} \times \overline{\rho}\_C\right). \tag{65}$$

*Symmetry* **2019**, *11*, 1077

The following three components define the resultant rotation. Among these, the first two contain the angular velocity and acceleration. The defining expressions are shown below:

$$\frac{1}{2} \cdot \int \left(\overline{\boldsymbol{\varepsilon}} \times \overline{\boldsymbol{\rho}}\_M\right)^T \cdot \left(\overline{\boldsymbol{\varepsilon}} \times \overline{\boldsymbol{\rho}}\_M\right) \cdot dm = \frac{1}{2} \cdot \overline{\boldsymbol{\varepsilon}}^T \cdot \left\{ \int \left(\overline{\boldsymbol{\rho}}\_M \times \right) \cdot \left(\overline{\boldsymbol{\rho}}\_M \times \right)^T \cdot dm \right\} \cdot \overline{\boldsymbol{\varepsilon}} = \frac{1}{2} \cdot \overline{\boldsymbol{\varepsilon}}^T \cdot \boldsymbol{l}\_S' \cdot \overline{\boldsymbol{\varepsilon}} \tag{66}$$

Performing some transformations and developing the terms from (76), the result is:

$$E\_A^{(1au)} = \frac{1}{2} \cdot \int \left(\overline{\boldsymbol{\varepsilon}} \times \overline{\boldsymbol{\rho}}\_M\right)^T \cdot \left(\overline{\boldsymbol{\omega}} \times \overline{\boldsymbol{\omega}} \times \overline{\boldsymbol{\rho}}\_M\right) \cdot dm = \frac{1}{2} \cdot \overline{\boldsymbol{\varepsilon}}^T \cdot \left(\overline{\boldsymbol{\omega}} \times\right) \cdot \left[\int \left(\overline{\boldsymbol{\rho}}\_M \times\right) \cdot \left(\overline{\boldsymbol{\rho}}\_M \times\right)^T \cdot dm\right] \cdot \overline{\boldsymbol{\omega}} = \frac{1}{2} \cdot \overline{\boldsymbol{\varepsilon}}^T \cdot \left(\overline{\boldsymbol{\omega}} \times \boldsymbol{l}\_S^{\prime} \cdot \overline{\boldsymbol{\omega}}\right). \tag{67}$$

The last component, corresponding to rotation motion, exclusively contains the angular velocity:

$$E\_A^{(1\omega^4)} = \frac{1}{2} \cdot \overline{\boldsymbol{\omega}}^T \cdot \left[ \overline{\boldsymbol{\omega}}^T \cdot \boldsymbol{l}\_S' \cdot \overline{\boldsymbol{\omega}} \right] \cdot \overline{\boldsymbol{\omega}} = \frac{1}{2} \cdot \overline{\boldsymbol{\omega}}^T \cdot \left[ \overline{\boldsymbol{\omega}}^T \cdot \operatorname{Trace} \left( \boldsymbol{l}\_{pS}' \right) \cdot \overline{\boldsymbol{\omega}} - \overline{\boldsymbol{\omega}}^T \cdot \boldsymbol{l}\_{pS}' \cdot \overline{\boldsymbol{\omega}} \right] \cdot \overline{\boldsymbol{\omega}}.\tag{68}$$

The acceleration energy of the first order can be written in the following final form:

$$\begin{array}{c} E\_A^{(1)} = \frac{1}{2} \cdot M \cdot \overline{a}\_0^T \cdot \overline{a}\_0 + M \cdot \overline{a}\_0^T \cdot \left(\overline{\varepsilon} \times \overline{\rho}\_{\text{C}}\right) + M \cdot \overline{a}\_0^T \cdot \left(\overline{\omega} \times \overline{\omega} \times \overline{\rho}\_{\text{C}}\right) + \frac{1}{2} \cdot \overline{\varepsilon}^T \cdot l\_{\text{S}}^{\prime} \cdot \overline{\varepsilon} + \\\ \qquad\quad + \overline{\varepsilon}^T \cdot \left(\overline{\omega} \times l\_{\text{S}}^{\prime} \cdot \overline{\omega}\right) + \frac{1}{2} \cdot \overline{\omega}^T \cdot \left[\overline{\omega}^T \cdot l\_{\text{S}}^{\prime} \cdot \overline{\omega}\right] \cdot \overline{\omega}. \end{array} \tag{69}$$

If the following conditions are met: *O* ≡ *C*, ρ*C* = 0, *and I S* ≡ *I* ∗ *S*, then the mobile reference frame is applied in the mass center of the rigid body and the position of the elementary particle (*dm*) with respect to this frame is known, as well as the inertial tensor axial and centrifugal *I S* ≡ *I* ∗ *S* with respect to the mass center. Thus, (69) defines the first order acceleration energy in a generalized form:

$$E\_A^{(1)} = \frac{1}{2} \cdot M \cdot \overline{a}\_{\mathbb{C}}^T \cdot \overline{a}\_{\mathbb{C}} + \frac{1}{2} \cdot \overline{\varepsilon}^T \cdot I\_{\mathbb{S}}^\* \cdot \overline{\varepsilon} + \overline{\varepsilon}^T \cdot \left(\overline{a} \times I\_{\mathbb{S}}^\* \cdot \overline{a}\right) + \frac{1}{2} \cdot \overline{\omega}^T \cdot \left[\overline{\omega}^T \cdot I\_{\mathbb{S}}^\* \cdot \overline{\omega}\right] \cdot \overline{\omega},\tag{70}$$

The expression (70) represents a reinterpretation of König's theorem for the acceleration energy of the first order. In the case of multibody systems (which also includes robots), the expression of acceleration energy of the first order (69) is determined according to:

$$\begin{split} \mathcal{E}\_{A}^{(1)} \Big[ \tilde{\mathcal{O}}(t); \stackrel{\scriptstyle}{\tilde{\mathcal{O}}}(t); \stackrel{\scriptstyle \mathcal{O}}{\tilde{\mathcal{O}}}(t) \right] &= \\ = (-1)^{\Delta\_{M}} \cdot \frac{1 - \Delta\_{M}}{1 + 3\Delta\_{M}} \sum\_{l=1}^{\tilde{\mathcal{U}}} \Big\{ \Big\{ \Big{v} \Big{1}\_{2} \cdot M\_{l} \cdot \big( {}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}}^{T} \big( {}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}} \big) + \Delta\_{M}^{2} \cdot \Big\{ \Big{v} \Big{1}\_{2} \cdot ({}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}}^{\*} \big( {}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}} \big{ {}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}} \big( {}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}} \big{ {}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}} \big( {}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}} \big{ {}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}} \big( {}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}} \big( {}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}} \big( {}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}} \big) \big{ v} \Big{ } {}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}} \big( {}^{(l)} \tilde{\mathcal{U}}\_{\mathcal{U}\_{l}} \big( {}^{($$

and *i* = 1 → *n* is the number of kinetic assemblies that compose the mechanical system.

Considering the advanced kinematics notions, the components of acceleration energy of the first order are:

*E*(1)*TR A* = 1 2 · '*n i*=1 *Mi* · *a*2 *Ci* = = 1 2 · '*n i*=1 *Mi*· ⎧ ⎪⎪⎪⎨ ⎪⎪⎪⎩ '*n j*=1 '*n p*=1 ⎡ ⎢⎢⎢⎢⎢⎢⎣ ∂2 (*m*) *rCi* ∂ (*m*) *qj* ·∂ (*m*) *qp* · .. *qj* · .. *qp* + 1 (*m*+1) 2 · ∂2 (*m*+1) *rCi* ∂ (*m*) *qj* ·∂ (*m*) *qp* · . *qj* · . *qp* + 1 *m*+1 · ∂ (*m*) *rCi* ∂ (*m*) *qj* · ∂ (*m*+1) *rCi* ∂ (*m*) *qp* · .. *qj* · . *qp* ⎤ ⎥⎥⎥⎥⎥⎥⎦ ⎫ ⎪⎪⎪⎬ ⎪⎪⎪⎭ (72) *E*(1)*ROT*ε *A* = 1 2 · (*n i*=1 ε*i T* · *I* ∗ *i* · ε*i* = ⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩ 1 2 · (*n i*=1 (*n j*=1 ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ ∂ (*m*) ψ*i* ∂ (*m*) *qj* · Δ*j* · .. *qj* + 1 *m* + 1 · ∂ (*m*+1) ψ*i* ∂ (*m*) *qj* · Δ*j* · . *qj* ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ *T* · *I* ∗ *i* · ε*i* ⎫ ⎪⎪⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎪⎪⎭ , (73)

*Symmetry* **2019**, *11*, 1077

$$I\_i^\* \cdot \overline{\varepsilon}\_i = I\_i^\* \cdot \sum\_{j=1}^n \left[ \frac{\partial \frac{(m)}{\partial \overline{\psi}\_i}}{\partial \overline{\vartheta}\_j^{(m)}} \cdot \Delta\_j \cdot \ddot{\overline{\eta}}\_j + \frac{1}{m+1} \cdot \frac{\partial \frac{(m+1)}{\partial \overline{\psi}\_i}}{\partial \overline{\vartheta}\_j^{(m)}} \cdot \Delta\_j \cdot \dot{\eta}\_j \right],\tag{74}$$

$$E\_A^{(1)ROTaw} = \sum\_{i=1}^n \overline{\varepsilon}\_i^T \cdot \left(\overline{a\_i} \times I\_i^\* \cdot \overline{\omega\_i}\right) = \sum\_{i=1}^n \sum\_{j=1}^n \left[\frac{\partial \overline{\overline{\psi}\_i}}{\partial \overline{\overline{q\_j}}} \cdot \Delta\_j \cdot \ddot{\overline{q}\_j} + \frac{1}{m+1} \cdot \frac{\partial \left.\overline{\overline{\psi}\_i}\right|\_j}{\partial \overline{q\_j}} \cdot \Delta\_j \cdot \dot{\overline{q}\_j}\right]^T \cdot E\_A^{(1)ROTaw},\tag{75}$$

$$E\_A^{(1)ROTaw} = \left(\overline{\boldsymbol{\alpha}}\_i \times I\_i^\* \cdot \overline{\boldsymbol{\alpha}}\_i\right) = \sum\_{i=1}^n \sum\_{j=1}^n \sum\_{p=1}^n \left[\frac{\partial \overline{\boldsymbol{\psi}}\_i^{(m)}}{\partial \boldsymbol{\dot{q}}\_j^{(m)}}\right] \times \left[I\_i^\* \cdot \frac{\partial \overline{\boldsymbol{\psi}}\_i^{(m)}}{\partial \boldsymbol{\dot{q}}\_p^{(m)}}\right] \cdot \Delta\_j \cdot \Delta\_p \cdot \boldsymbol{\dot{q}}\_j \cdot \boldsymbol{\dot{q}}\_p. \tag{76}$$

According research [10,14], the sudden motion of MBS (multibody system), the transient motion phases, as well as the mechanical systems subjected to the action of a system of external forces, with a time variation law [19,20], are characterized by linear and angular accelerations of the higher order. Thus, the acceleration energy of the second order was also developed.  According to [10,14], an explicit equation for the third order acceleration energy is developed:

*E*(3) *A* θ(*t*); . θ(*t*); .. θ(*t*); ... θ(*t*); .... θ (*t*) = (−1) Δ*M* · 1−Δ*M* 1+3·Δ*M* · '*n i*=1 1 2 · *Mi* · *i* ... *v T Ci* · *i* ... *v Ci* + +Δ2 *M* · '*n i*=1 " 1 2 · *i* ... ω*T i* · *i I*∗ *i* · *i* ... ω*i* + 3 · ω*T i* · /... ω*i* × *I* ∗ *pi* · .. ω*i* 0 +3 · . ω *T i* · /... ω*i* × *I* ∗ *pi* · . ω*i* 0 + 3 · .. ω *T i* · /... ω*i* × *I* ∗ *pi* · ω*i* 0 + +2 · ω*i* × ... ω*i T* · *I* ∗ *pi* · / . ω*i* × ω*i* 0 − 5 · . ω *T i* · ω*T i* · *I* ∗ *i* · ... ω*i* · ω*i* − ω*T i* · . ω *T i* · *I* ∗ *i* · ... ω*i* · ω*i*+ +ω*T i* · ω*T i* · *I* ∗ *pi* · ... ω*i* × ω*i* · ω*i* # (80)

An important aspect that has to be mentioned is that the component *E*(3) *A* θ(*t*); . θ(*t*); .. θ(*t*); ... θ 2 (*t*) of the acceleration energy of the third order is not included in the definition Equation (80). In the case of the multibody systems, the acceleration energies of the higher order are determined based on the following generalized expression:

*E*(*p*) *A* ⎡ ⎢⎢⎢⎢⎢⎢⎣ θ(*t*); . θ(*t*); ··· ; (*p*+1) θ(*t*) ⎤ ⎥⎥⎥⎥⎥⎥⎦ = 1 2 '*n i*=1 ⎧ ⎪⎪⎨ ⎪⎪⎩ , (*p*) *vi T* · (*p*) *vi* · *dm* ⎫ ⎪⎪⎬ ⎪⎪⎭ = 1 2 '*n i*=1 ⎧ ⎪⎪⎪⎨ ⎪⎪⎪⎩ , *Trace* ⎡ ⎢⎢⎢⎢⎢⎢⎣ (*p*+1) *ri* · (*p*+1) *r T i* ⎤ ⎥⎥⎥⎥⎥⎥⎦ · *dm* ⎫ ⎪⎪⎪⎬ ⎪⎪⎪⎭ = = 1 2 · '*n i*=1 *Trace* ⎧ ⎪⎪⎨ ⎪⎪⎩ (*p*+1) 0 *i* [*R*] · , *i r*∗ *i* · *i r*∗*T i* · *dm* +*i rCi* · *i rT Ci* · , *dm* · (*p*+1) 0 *i* [*R*] *T* ⎫ ⎪⎪⎬ ⎪⎪⎭ + 1 2 · '*n i*=1 *Trace* ⎡ ⎢⎢⎢⎢⎢⎢⎣ (*p*+1) *pi* · (*p*+1) *pT i* ⎤ ⎥⎥⎥⎥⎥⎥⎦ · , *dm* = = 1 2 · '*n i*=1 *Trace* ⎧ ⎪⎪⎨ ⎪⎪⎩ (*p*+1) 0 *i* [*R*] · *i I*∗ *pi* + *Mi* · *i rCi* · *i rT Ci* · (*p*+1) 0 *i* [*R*] *T* ⎫ ⎪⎪⎬ ⎪⎪⎭ + 1 2 · '*n i*=1 *Trace* ⎡ ⎢⎢⎢⎢⎢⎢⎣ (*p*+1) *pi* · (*p*+1) *pT i* ⎤ ⎥⎥⎥⎥⎥⎥⎦ · *Mi*, *where p* = 1, 2, 3, 4, ... (81)

$$where \, p \ge 1, \, k \ge 1, \, \{p; k\} = \{1; 2; 3; \ldots\}, \, (p) - order of \, acceleration, \, (k) - derivation \, or \, der \, (t - 1) \, \begin{bmatrix} 0 \\ \overline{\theta}(t); \dot{\overline{\theta}}(t); \cdots; \overline{\theta}(t) \end{bmatrix} = E\_A^{(p)} \left[ \overline{\theta}(t); \dot{\overline{\theta}}(t); \cdots; \overline{\theta}(t) \right] \tag{82}$$

Expression (81) includes the inertia tensor planar and centrifugal, *i I*∗ *pi* = , *i r*∗ *i* · *i r*∗*T i* · *dm*.

The expressions for the acceleration energies of the higher order, according to [13] are further included in the equations of advanced dynamics corresponding to the mechanical systems which are subjected to fast and sudden movements.

#### **4. Discussion on the Results**

This section was introduced exclusively for highlighting the existence of the accelerations and acceleration energies of the higher order. For this purpose, according to Figure 4, the rotation motion of the arm of the serial robot Fanuc LR Mate 100 iB, on the angular interval (0, π) was considered. The total time of the rotation motion was of 1.7 s.

**Figure 4.** Fanuc LR Mate 100 iB robot. The location of the transducer on the robot arm.

For highlighting the time variation law of the acceleration energies of the higher order, the polynomial interpolating functions of the fifth order have been applied in accordance with (48)–(51). This method requires, according to kinematical constraints (51), the establishment through measurements of the angular acceleration corresponding to robot arm, .. *q*3. As a result, in the first step, the time variation law corresponding to the tangential component of the acceleration for the characteristic point belonging to the robot arm has been determined at every 0.033 s. The experimental data was collected by using a mono-axial accelerometer. The robot arm (more precisely the third kinetic ensemble of the robot) was engaged in rotation motion at high speeds and by using the data collected by the accelerometer. The time variation law for the tangential component of the acceleration of the characteristic point was graphically represented. Considering the rotation motion, in the next step the time variation law of the angular acceleration .. *q*3(τ) was also represented. Both graphical representations are illustrated in Figure 5.

**Figure 5.** (**a**) Tangential accelerations variation law. (**b**) Generalized accelerations variation law. The above presented parameters were also determined on the same angular interval (0, π), by applying the sudden stop. The graphical representation is illustrated in Figure 6.

**Figure 6.** (**a**) The tangential first order accelerations variation law; (**b**) The generalized accelerations of first order variation law.

According to Figure 5, the maximum value for the linear acceleration exceeds 15 *m*/*s*2 . , which defines a sudden motion. The authors have considered, for this application, that the motion can be described by using the polynomial interpolating functions of the fifth order:

$$\stackrel{\text{e}}{\ddot{q}}\_{ji}(\tau) \;= \; \frac{\tau\_i - \tau}{t\_i} \cdot \stackrel{\text{e}}{\ddot{q}}\_{ji}(\tau\_{i-1}) \; + \; \frac{\tau - \tau\_{i-1}}{t\_i} \cdot \stackrel{\text{e}}{\ddot{q}}\_{ji}(\tau\_i) \; , \tag{83}$$

$$\ddot{q}\_{\ddot{j}\bar{i}}(\tau) = -\frac{\left(\tau\_{\bar{i}} - \tau\right)^{2}}{2 \cdot t\_{\bar{i}}} \cdot \dddot{\vec{q}}\_{\bar{j}\bar{i}-1} + \frac{\left(\tau - \tau\_{\bar{i}-1}\right)^{2}}{2 \cdot t\_{\bar{i}}} \cdot \dddot{\vec{q}}\_{\bar{j}\bar{i}} + a\_{\bar{j}\bar{1}},\tag{84}$$

$$\ddot{q}\_{\vec{\mu}}(\tau) = \frac{\left(\tau\_{\vec{i}} - \tau\right)^{3}}{6 \cdot t\_{\vec{i}}} \cdot \ddot{q}\_{\vec{j}\vec{i}-1} + \frac{\left(\tau - \tau\_{\vec{i}-1}\right)^{3}}{6 \cdot t\_{\vec{i}}} \cdot \dddot{q}\_{\vec{j}\vec{i}} + a\_{\vec{j}\vec{1}} \cdot \tau + a\_{\vec{j}\vec{2}}.\tag{85}$$

$$\dot{q}\_{\vec{j}\vec{i}}(\tau) = -\frac{\left(\tau\_{i} - \tau\right)^{4}}{24 \cdot t\_{i}} \cdot \dddot{q}\_{\vec{j}\vec{i}-1} + \frac{\left(\tau - \tau\_{i-1}\right)^{4}}{24 \cdot t\_{i}} \cdot \dddot{q}\_{\vec{j}\vec{i}} + a\_{\vec{j}\vec{1}} \cdot \frac{\tau^{2}}{2} + a\_{\vec{j}2} \cdot \tau + a\_{\vec{j}3} \tag{86}$$

$$q\_{\vec{\mu}}(\tau) = \frac{\left(\tau\_{i} - \tau\right)^{5}}{120 \cdot t\_{i}} \cdot \overset{\cdots}{q}\_{\vec{\mu}-1} + \frac{\left(\tau - \tau\_{i-1}\right)^{5}}{120 \cdot t\_{i}} \cdot \overset{\cdots}{q}\_{\vec{\mu}} + a\_{\vec{\mu}1} \cdot \frac{\tau^{3}}{6} + a\_{\vec{\mu}2} \cdot \frac{\tau^{2}}{2} + a\_{\vec{\mu}3} \cdot \tau + a\_{\vec{\mu}4}.\tag{87}$$

where *ajip*, *p* = 1 → 4, are the integration constants which are determined from the geometrical and kinematical constraints (50), with an important role in ensuring the continuity of the rotation motion on the angular interval (0, π), characterized by 51 interpolation segments. Every interpolation segment contains four unknowns, represented by the integration constants. It is obvious that the total number of the unknowns is 204. In order to determine the unknowns, according to (50), the numerical values for the angular acceleration (Figure 5b) were substituted in the functions (83)–(87) presented above. Considering the large number of results, more exactly 255 polynomial functions, only a few sequences are illustrated in Tables 1 and 2.

**Table 1.** The polynomial time functions for the generalized coordinates for the third rotation joint.

**Table 2.** The Polynomial Time Functions for Generalized Variables in Case of the Third Rotation Joint.

As a result, the numerical values for the angular acceleration (see Figure 5) were applied, finally obtaining all of the integration constants. In order to illustrate this step, some sequences containing the time variation laws of the functions described by Equations (83)–(87) are presented. These results are included in Tables 1 and 2. Similarly, the polynomial functions are obtained for all the 51 interpolation segments. As a result, the time variation law for the angular acceleration of the second and third order are also graphically represented according to Figure 7:

**Figure 7.** (**a**) The variation law of the angular accelerations of second order. (**b**) The variation law for the angular accelerations of third order.

Thus, by customizing Expressions (71), (77)–(80) for the considered robot structure, the following expressions characterizing the acceleration energies of the first, second and third order are:

$$E\_{A\bar{\rm ik}}^{(1)}(\tau) = \frac{1}{2} \cdot \left( M\_3 \cdot x\_{\gets3}^2 + M\_3 \cdot z\_{\gets3}^2 + {}^3I\_y \right) \cdot \left( \ddot{q}\_{3\bar{\rm ik}}^2(\tau) + \dot{q}\_{3\bar{\rm ik}}^4(\tau) \right) \tag{88}$$

$$E\_{A\&}^{(2)}(\mathbf{r}) = \frac{1}{2} \cdot \left(M\_3 \cdot \mathbf{x}\_{\subset3}^2 + M\_3 \cdot \mathbf{z}\_{\subset3}^2 + {}^3I\_9\right) \cdot \left[\overleftarrow{q}\_{3\text{ik}}^2(\mathbf{r}) - 2 \cdot \dot{\vec{q}}\_{3\text{ik}}^3(\mathbf{r}) \cdot \overleftarrow{\vec{q}}\_{3\text{ik}}(\mathbf{r}) + 9 \cdot \dot{\vec{q}}\_{3\text{ik}}^2(\mathbf{r}) \cdot \overrightarrow{\vec{q}}\_{3\text{ik}}^2(\mathbf{r}) + \dot{\vec{q}}\_{3\text{ik}}^6(\mathbf{r})\right] \tag{89}$$

$$E\_{A\bar{k}}^{(3)}(\tau) = \frac{1}{2} \cdot \left| \left( M\_3 \cdot \mathbf{x}\_{\text{C3}}^2 + M\_3 \cdot \mathbf{z}\_{\text{C3}}^2 + {}^3I\_y \right) \cdot \left[ \dot{q}\_{3\text{ik}}^8(\tau) - 8 \cdot \dot{q}\_{3\text{ik}}^5(\tau) \cdot \dddot{q}\_{3\text{ik}(\tau)} + 30 \cdot \dot{q}\_{3\text{ik}}^4(\tau) \cdot \ddot{q}\_{3\text{ik}}^2(\tau) - 4 \right] \right|^2 = \frac{1}{2} \cdot \left| \dot{q}\_{3\text{ik}}^2(\tau) \cdot \dddot{q}\_{3\text{ik}}^4(\tau) \cdot \dddot{q}\_{3\text{ik}}^2(\tau) - 12 \cdot \dot{q}\_{3\text{ik}}^2(\tau) \cdot \ddot{q}\_{3\text{ik}}(\tau) \cdot \dddot{q}\_{3\text{ik}}^2(\tau) + \dots \right| \cdot \tag{90}$$

The obtained results (see example from Tables 1 and 2) regarding the higher order polynomial functions have been included in the above expressions of the acceleration energies that characterize the fast rotation motion of the robot arm. The results are represented in Figures 8 and 9:

**Figure 8.** (**a**) The variation law for the kinetic energy. (**b**) The time variation law of the acceleration energy of the first order.

**Figure 9.** (**a**) The variation law for the acceleration energy of second order; (**b**) The time variation law of the acceleration energy of third order.

In the above graphical representations (Figures 8 and 9), the time variation laws for the kinetic energy are illustrated, as well as the acceleration energies of the first, second and third order for an application in which the serial Fanuc robot, presented in the Figure 4, is implemented.

#### **5. Conclusions**

The main objective of this paper was to highlight some new developments in the field of advanced dynamics of mechanical systems and higher order dynamic equations. These equations were developed by using matrix exponentials which proved to have indisputable advantages in the matrix study of any complex mechanical system. In this paper, new approaches have been proposed based on differential principles from analytical mechanics by using some important dynamics notions, regarding the acceleration energies of the first, second and third order. The study was extended to the equations of the higher order, which provided the possibility of applying the initial motion conditions in the positions, velocities and accelerations of the first and second order. In order to determine the time variation laws for the generalized variables and acceleration energies of the higher order, the time polynomial functions of the fifth order were applied. In Newtonian dynamics, the establishment of dynamic equations as second order differential equations are performed by considering the kinetic energy as the main function. In the case of fast and sudden movements, the presence of higher order accelerations leads to higher order differential equations which define the variation laws in relation to time in the case of the forces and the driving moments that impose the motion of the mechanical systems. As a result, the application of acceleration energies highlights, in an explicit and direct form, the presence of higher-order accelerations, unlike the case where the study was performed considering kinetic energy.

In accordance to inverse kinematics, also known as the control kinematics of the robots, the use of polynomial functions leads to the control of the mechanical motion, especially of transitory motions. These aspects were exemplified in the previous section. Therefore, it was noticed that on the time interval of 1.7 s, a number of 51 measurements of the tangential acceleration were performed at every 0.033 s. Regarding the polynomial functions, this aspect led to a very good approximation of the mechanical motion. At the same time, the polynomial functions were substituted in the expressions of the acceleration energies of the higher order. Therefore, their time variation laws were obtained as shown in graphical representations from the previous section.

#### **6. Contributions of the Authors**

The purpose of this paper was to present some original defining expressions for acceleration energies of the higher order. These expressions also include the accelerations of the higher order which define the fast and sudden motions of mechanical systems. The novelty of this paper was in determining in an explicit form of the expressions for the acceleration energies based on advanced kinematic notions and using polynomial interpolation functions of the higher order.

Therefore, the contributions of the main author are highlighted by a great number of expressions included in the following sections of this paper:

**Author Contributions:** Conceptualization, I.N.; methodology, I.N.; investigation and software A.-V.C; validation, I.N.; formal analysis, I.N; resources, I.N., A.-V.C.; data curation, I.N., A.-V.C.; writing—original draft preparation, A.-V.C; writing—review and editing, I.N., A.-V.C.; visualization, A.-V.C; supervision, I.N.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to thank Technical University of Cluj Napoca, Romania for providing the technical support in carrying out the experimental analysis.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **On the Partition of Energies for the Backward in Time Problem of Thermoelastic Materials with a Dipolar Structure**

#### **M. Marin 1,\*, S. Vlase 2, R. Ellahi 3,4 and M.M. Bhatti 5,6**

Received: 12 June 2019; Accepted: 26 June 2019; Published: 2 July 2019

**Abstract:** We first formulate the mixed backward in time problem in the context of thermoelasticity for dipolar materials. To prove the consistency of this mixed problem, our first main result is regarding the uniqueness of the solution for this problem. This is obtained based on some auxiliary results, namely, four integral identities. The second main result is regarding the temporal behavior of our thermoelastic body with a dipolar structure. This behavior is studied by means of some relations on a partition of various parts of the energy associated to the solution of the problem.

**Keywords:** backward in time problem; dipolar thermoelastic body; uniqueness of solution; Cesaro means; partition of energies

#### **1. Introduction**

In our study, we approach a thermoelastic body having a dipolar structure. This kind of structure falls within a more general theory, namely, the theory of bodies with microstructure. The first studies in this context were published by Eringen (see, for instance, references [1,2]). One may deduce the importance of the dipolar structure due to the large number of published studies dedicated to this topic, of which we can mention [3–7]. As such, our present work can be considered a continuation in this respect.

A continuation of the theories of microstructure, is a theory that takes into account the voids in the materials. It is considered that the initiators of this theory were Nunziato and Cowin, in their known paper [8]. After that, the number of studies within this topic has grown impressively. We want to enumerate some of these [9–16]: The first result for the backward in time problem belongs to Serrin, who approached this problem in the context of Navier–Stokes equations (see [17]). In the paper [17], we find some uniqueness in the results with regards to the forward in time problem. After that, the number of studies dedicated to the backward in time problem has increased considerably. Of particular importance are the works [18–27]. We have to point out that the results obtained by Ciarletta in [23], and Ciarletta and Chirita in [24] were improved by Quintanilla in [25]. In addition, Quintanilla approached the question of location in time for solutions to the backward in time problem, in the context of thermoelasticity of Green and Naghdi [26,27] and the theory of porous thermoelastic

bodies. The elastic porous bodies were also approached by Iovane and Passarella in [28]. The forward in time problem, in the context of theory for thermomicrostretch elastic solids, was approached by Passarella and Tibullo in [29]. It is worth noting that the idea of considering non-standard problems, in the context of the general theory of bodies having a dipolar structure, was inspired by Quintanilla and Straughan's work [30]. In [31,32] it is proved that the nonhomogenous temperature field has a profound influence on the nanobeam mechanics. Additionally, [33] is a recent contribution on stress-driven nonlocal modeling of thermoelastic nanostructures.

Here is the plane of our study. First of all, we summarize the main equations, the initial conditions, and the boundary data of the mixed problem. Then, we prove some estimates for the gradient of classical and dipolar displacements, and for the gradient of the function of voids. In the last part of our study we prove the main result, namely, the continuous dependence of solutions—with regards to the coefficients that couple the equations describing the dipolar deformation—with the equations that describe the behavior of voids. The description of the continuous dependence was possible due to the definition of an adequate measure.

#### **2. Basic Equations and Conditions**

In our paper, we approach a thermoelastic body having a dipolar structure. We will use an anisotropic body, which is situated in a regular domain *D*, included in the physical space *E*3, that is, the three-dimensional Euclidean space. Consider that the boundary of the domain is a piecewise smooth surface *∂D*. The closure of *D* is usual denoted by *D*¯ , *D*¯ = *D* ∪ *∂D*. An orthonormal system of references is introduced, and then tensors and vectors have components with Latin subscripts over 1, 2, 3. Typical conventions for summation over repeated indices and for derivation operations are implied. So, a subscript preceded by a comma is for a partial derivative with regards to corresponding spatial coordinate; while a superposed dot is for a derivative with regards to time variable. All the functions we use are assumed to be sufficiently regular as necessary. Additionally, if there is no possibility of confusion, then the dependence of function with regards to its spatial or time variables will be omitted. The evolution of the body with a dipolar structure will be described with the help of the following specific variables:

$$\min\_{i}(\mathbf{x},t), \ \phi\_{i\bar{j}}(\mathbf{x},t), \ \theta(\mathbf{x},t), \ (\mathbf{x},t) \in D \times [0, t\_0). \tag{1}$$

Here, we denoted by u*i* the components of the displacement vector field, by *φij* the components of the dipolar displacement tensor field, and by *θ* the absolute temperature.

Using the above variables u*i*(*x*, *t*), and *φij*(*x*, *t*) we will introduce the components of the tensors of strain, namely, *εij*, *κij*, and *χijk*, as follows:

$$
\mathfrak{L}\varepsilon\_{i\bar{j}} = \mathfrak{u}\_{\bar{j},i} + \mathfrak{u}\_{i,\bar{j}\prime}, \\
\kappa\_{i\bar{j}} = \mathfrak{u}\_{\bar{j},\bar{i}} - \phi\_{i\bar{j}\prime}, \\
\chi\_{i\bar{j}k} = \phi\_{\bar{i}\bar{j},k}.
\tag{2}
$$

All our considerations are made within a linear theory, therefore it is natural to consider that the Helmholtz's free energy is a quadratic form with regards to its independent constitutive variables. The Helmholtz's free energy in the reference configuration will be denoted by *W*. So, in accordance with the principle of conservation of energy, we develop in series the function *W* and we keep the terms only until the second order. Because the reference state was assumed to be free of loadings, we deduce that the Helmholtz's free energy per mass can be considered of the following form (see [30]):

$$\begin{array}{l} \mathcal{W} = \frac{1}{2} A\_{ijmn} \varepsilon\_{ij} \varepsilon\_{mn} + D\_{ijmn} \varepsilon\_{ij} \chi\_{mn} + F\_{ijmnr} \varepsilon\_{ij} \chi\_{mnr} + \frac{1}{2} B\_{ijmn} \kappa\_{ij} \chi\_{mn} \\ + G\_{ijmnr} \kappa\_{ij} \chi\_{mnr} + \frac{1}{2} C\_{ijkmnr} \chi\_{ijk} \chi\_{mnr} - a\_{ij} \varepsilon\_{ij} \theta - b\_{ij} \kappa\_{ij} \theta - c\_{ijk} \chi\_{ijk} \theta - \frac{1}{2} c \theta^2. \end{array} \tag{3}$$

We will use this form of free energy used in the entropy production inequality and deduce the motion equations. In addition, from the same inequality, the constitutive equations are obtained. These equations express the tensors of stress with the help of the tensors of deformation. We will denote the components of the stress measures by *τij*, *ηij*, and *μijk*. In this way, the constitutive equations establish a connection between the tensors *τij*, *ηij*, *μijk* and the tensors *εij*, *κij*, *χijk*.

We will use a procedure similar to that used by Green and Rivlin in [6], so that considering the Helmholtz' free energy (3) we deduce the next constitutive equations

$$\begin{aligned} \pi\_{\dot{l}\dot{l}} &= \frac{\partial \mathcal{W}}{\partial \varepsilon\_{\dot{i}}} = A\_{\dot{i}\dot{j}mn} \varepsilon\_{mn} + D\_{mn\dot{i}} \chi\_{mn} + F\_{mnrrj} \chi\_{mnr} - a\_{\dot{i}\dot{j}} \theta\_{\prime} \\ \eta\_{\dot{i}j} &= \frac{\partial \mathcal{W}}{\partial \kappa\_{\dot{i}j}} = D\_{\dot{i}jmn} \varepsilon\_{mn} + B\_{\dot{i}jmn} \kappa\_{mn} + G\_{\dot{i}jmnr} \chi\_{mnr} - b\_{\dot{i}j} \theta\_{\prime} \\ \mu\_{ijk} &= \frac{\partial \mathcal{W}}{\partial \chi\_{ijk}} = F\_{\dot{i}jkm} \varepsilon\_{mn} + G\_{mn\dot{i}\dot{j}} \kappa\_{mn} + G\_{\dot{i}jkmr} \chi\_{mnr} - c\_{\dot{i}j\dot{k}} \theta\_{\prime} \\ \eta &= -\frac{\partial \mathcal{W}}{\partial \theta} = a\_{\dot{i}\dot{j}} \varepsilon\_{\dot{i}j} + b\_{\dot{i}\dot{j}} \kappa\_{\dot{i}j} + c\_{\dot{i}j\dot{k}} \chi\_{ijk} + c \theta\_{\prime} \end{aligned} \tag{4}$$

which are satisfied in *D* × [0, *t*0). Here, we denoted by *η* the entropy per unit mass. For the vector of heat flux, having the components *qi* we have a classical constitutive relation, namely,

$$q\_{\dot{i}} = K\_{\dot{i}\dot{j}} \theta\_{\dot{i}\dot{\prime}} \tag{5}$$

where *Kij* is the thermal conductivity symmetric tensor.

Also, we can deduce the main equations that govern the thermoelasticity of bodies with a dipolar structure, namely (see [5,6]):

$$\begin{aligned} \left( \left( \mathbf{r}\_{i\dot{\jmath}} + \eta\_{i\dot{\jmath}} \right)\_{,\dot{\jmath}} + \rho f\_i = \rho \ddot{\mathbf{u}}\_{i\prime} \\ \mu\_{i\dot{\jmath}k,\dot{\imath}} + \eta\_{\dot{\jmath}k} + \rho g\_{\dot{\jmath}k} = I\_{kr} \ddot{\phi}\_{\dot{\jmath}\dot{\imath}} \end{aligned} \tag{6}$$

$$
\rho T \omicron \dot{\eta} = q\_{i,i} + \rho r.\tag{7}
$$

The signification of the notations that we introduced in preceding equations is as follows: *ρ*, the density of mass, which is a constant; *Iij*, the symmetric tensor of microinertia; *k*, the intrinsic inertia; *εij*, *κij*, *χijk*, the strain tensors; *τij*, *ηij*, *μijk*, the stress tensors; *fi*, the body forces; *gjk*, the dipolar charges; *Aijmn*, *Bijmn*, ..., *aij*, the functions that describe the properties of the material in terms of elasticity. Suppose the following symmetry relations take place:

$$\begin{array}{l} A\_{ijmn} = A\_{jimn} = A\_{mnij}, \; B\_{ijmn} = B\_{mnij}, \; a\_{ij} = a\_{ji}, \\ \mathbb{C}\_{ijkmnr} = \mathbb{C}\_{mnrijk}, \; F\_{ijkmn} = F\_{ijkmn}, \; D\_{ijmn} = D\_{ijnm}. \end{array} \tag{8}$$

Assuming that there are no supply terms and taking into account the constitutive Equations (4) and (5) and the kinematic Equation (2), Equations (6) and (7) become

$$\begin{aligned} \rho \ddot{\mathbf{u}}\_i &= \left[ \left( \mathbf{C}\_{ijmn} + \mathbf{G}\_{ijmn} \right) \mathbf{u}\_{n,m} + \left( \mathbf{G}\_{mnij} + \mathbf{B}\_{ijmn} \right) \left( \mathbf{u}\_{n,m} - \phi\_{mn} \right) + \\ &+ \left( \mathbf{F}\_{mnij} + \mathbf{D}\_{ijmnr} \right) \boldsymbol{\phi}\_{nr,m} - \left( a\_{ij} + b\_{ij} \right) \boldsymbol{\theta} \right]\_{,j}, \end{aligned} \tag{9a}$$

$$\begin{aligned} \mathbf{J}\_1 \ \ddot{\mathbf{g}}\_i &= \begin{bmatrix} \mathbf{F}\_{i1} & \mathbf{u}\_{n} \dots \mathbf{u}\_i \mathbf{D} & \cdots \left( \mathbf{u}\_{n} \dots \mathbf{u}\_{n} \dots \mathbf{f}\_{i1} & \boldsymbol{\theta}\_{n} \dots \mathbf{u}\_i \dots \boldsymbol{\theta} \end{bmatrix} \end{aligned} \tag{9a}$$

$$\begin{aligned} \mathbf{^L}\_{kr}\ddot{\boldsymbol{\Phi}}\_{jr} &= \left[ F\_{ijkmn}\mathbf{u}\_{n,m} + D\_{mmj\dot{k}} \left( \mathbf{u}\_{n,m} - \boldsymbol{\Phi}\_{mm} \right) + A\_{ij\dot{k}mnr} \boldsymbol{\Phi}\_{nr,m} - \mathbf{c}\_{ijk} \boldsymbol{\Phi} \right]\_{,i} + \\ &+ G\_{jkmn}\mathbf{u}\_{m,n} + B\_{jkmn} \left( \mathbf{u}\_{n,m} - \boldsymbol{\Phi}\_{mm} \right) + D\_{jkmmr} \boldsymbol{\Phi}\_{nr,m} - b\_{jk} \boldsymbol{\Phi}, \end{aligned} \tag{9b}$$

$$\mathcal{K}\_{i\dot{\jmath}} \left( \theta\_{,\dot{\jmath}} \right)\_{,\dot{\jmath}} = -T\_0 \left[ a\_{\dot{\imath}\dot{\jmath}} \dot{\mathbf{u}}\_{i,\dot{\jmath}} + b\_{\dot{\imath}\dot{\jmath}} \left( \dot{\mathbf{u}}\_{\dot{\jmath},\dot{\imath}} - \dot{\phi}\_{\dot{\imath}\dot{\jmath}} \right) + c\_{\dot{\imath}\dot{\jmath}k} \dot{\phi}\_{\dot{\imath}\dot{\jmath},k} + c\dot{\theta} \right]. \tag{9c}$$

From now, we will assume that the Equations (2), (4), and (9) will be satisfied on the interval (−∞, 0].

The outward unit normal to the surface *∂D* has the components *ni*. With the help of this normal we can define the surface traction's of components *ti*, the surface couple of components *μjk*, and the flux of heat, *q*. All of this makes sense in every point of regularity of the boundary *∂D* and has the following expressions

$$t\_i = \left(\pi\_{i\rangle} + \eta\_{i\rangle}\right) n\_{\rangle}, \; \mu\_{j\bar{k}} = \mu\_{ijk} n\_{i\prime}, \; q = q\_i n\_i. \tag{10}$$

In close relation to these surface tractions, we consider the following homogeneous boundary conditions:

$$\begin{aligned} \mathbf{u}\_i(\mathbf{x},t) &= 0, \ (\mathbf{x},t) \in \partial D\_\mathbf{u} \times (-\infty,0], \ t\_i = 0, \ (\mathbf{x},t) \in \partial D\_\mathbf{u}^\varepsilon \times (-\infty,0],\\ \phi\_{ij}(\mathbf{x},t) &= 0, \ (\mathbf{x},t) \in \partial D\_\Phi \times (-\infty,0], \ m\_{jk} = 0, \ (\mathbf{x},t) \in \partial D\_\Phi^\varepsilon \times (-\infty,0],\\ \theta(\mathbf{x},t) &= 0, \ (\mathbf{x},t) \in \partial D\_\theta \times (-\infty,0], \ q = 0, \ (\mathbf{x},t) \in \partial D\_\theta^\varepsilon \times (-\infty,0], \end{aligned} \tag{11}$$

where the surfaces *∂Du*, *∂Dφ*, *∂Dθ*, and its complements *∂Dc u*, *∂Dc φ*, *∂Dc θ* are subsurfaces of the border *∂D*, which are subject to the following restrictions:

$$\begin{aligned} \partial D\_{\mathfrak{u}} \cup \partial D\_{\mathfrak{u}}^{\varepsilon} &= \partial D\_{\mathfrak{\theta}} \cup \partial D\_{\mathfrak{\theta}}^{\varepsilon} = \partial D\_{\mathfrak{\theta}} \cup \partial D\_{\mathfrak{\theta}}^{\varepsilon} = \partial D\_{\mathfrak{\theta}}, \\ \partial D\_{\mathfrak{u}} \cap \partial D\_{\mathfrak{u}}^{\varepsilon} &= \partial D\_{\mathfrak{\theta}} \cap \partial D\_{\mathfrak{\theta}}^{\varepsilon} = \partial D\_{\mathfrak{\theta}} \cap \partial D\_{\mathfrak{\theta}}^{\varepsilon} = \mathcal{Q}. \end{aligned}$$

We still have to add the final restrictions. So, on the closed domain *D*¯ we have:

$$\begin{array}{l} \mathfrak{u}\_{i}(\mathbf{x},0) = \mathfrak{u}\_{i}^{0}(\mathbf{x}), \ \dot{\mathfrak{u}}\_{i}(\mathbf{x},0) = \mathfrak{u}\_{i}^{1}(\mathbf{x}), \ \theta(\mathbf{x},0) = \theta^{0}(\mathbf{x}),\\ \Phi\_{\dot{\mathbf{i}}j}(\mathbf{x},0) = \Phi\_{\dot{\mathbf{i}}j}^{0}(\mathbf{x}), \ \dot{\Phi}\_{\dot{\mathbf{i}}j}(\mathbf{x},0) = \Phi\_{\dot{\mathbf{i}}j}^{1}(\mathbf{x}), \end{array} \tag{12}$$

where u0 *i* (*x*), u1 *i* (*x*), *φ*0 *ij*(*x*), *φ*1 *ij*(*x*), and *θ*0(*x*) are continuous prescribed functions in all points where they are defined. Additionally, these functions are assumed be compatible with conditions (11) on the appropriate subsets of *∂D*.

Let us consider the internal energy density Ψ (see [30]), which has the following expression:

$$\Psi = \frac{1}{2} A\_{ijmn} \varepsilon\_{ij} \varepsilon\_{mn} + D\_{ijmn} \varepsilon\_{ij} \kappa\_{mn} + F\_{ijmnr} \varepsilon\_{ij} \chi\_{mnr} + \frac{1}{2} B\_{ijmnr} \kappa\_{ij} \chi\_{mnr} + \frac{1}{2} C\_{ijkmnr} \chi\_{ij} \chi\_{mnr}.$$

P is denoted the so-called boundary-final value problem, which consists of Equation (9), the boundary restrictions (11), and the final data (12).

To obtain the results we have proposed, we will have to impose some conditions on the functions we are dealing with.

So, if *Jm*(*x*) is the minimum eigenvalue of the inertia tensor *Iij*(*x*), then we need to assume that *Jm* and *ρ* are continuous functions and the constitutive coefficients are of class *C*1(*D*). We also assume that:

Based on hypothesis (b), we deduce that there exist two positive numbers, *Km* and *KM*, so that

$$\mathbb{K}\_m \theta\_{,i} \theta\_{,j} \le \mathbb{K}\_{i\uparrow} \theta\_{,i} \theta\_{,j} \le \mathbb{K}\_M \theta\_{,i} \theta\_{,j},\tag{14}$$

and, as a consequence of the hypothesis (c), we can find the positive constants *M*1 and *M*2 so that the next inequality is satisfied:

$$\frac{M\_1}{2} \left( \varepsilon\_{i\bar{j}} \varepsilon\_{i\bar{j}} + \kappa\_{i\bar{j}} \kappa\_{i\bar{j}} + \chi\_{i\bar{j}k} \chi\_{i\bar{j}k} \right) \le \Psi \le \frac{M\_2}{2} \left( \varepsilon\_{i\bar{j}} \varepsilon\_{i\bar{j}} + \kappa\_{i\bar{j}} \kappa\_{i\bar{j}} + \chi\_{i\bar{j}k} \chi\_{i\bar{j}k} \right). \tag{15}$$

These hypotheses are not considered as very restrictive, as they are commonly imposed in mechanics of continuous media.

It is not difficult to equate our boundary-final value problem P with a boundary-initial problem, denoted by P , by a convenient change of variables. In this regard, we set *h* (*t* ) = *h*(*t*), for *t* = −*t*. But, to simplify writing, we will give up the sign "prime" so that the P problem will be defined by the following conditions and equations:

$$\left(\mathbf{K}\_{\dot{l}\dot{l}}\left(\theta\_{,\dot{l}}\right)\_{,\dot{l}} = T\_0 \left[a\_{\dot{l}\dot{l}}\dot{\mathbf{u}}\_{\dot{l}\dot{l}} + b\_{\dot{l}\dot{l}}\left(\dot{\mathbf{u}}\_{\dot{l},\dot{l}} - \dot{\phi}\_{\dot{l}\dot{l}}\right) + c\_{\dot{l}\dot{l}\dot{k}}\dot{\phi}\_{\dot{l}\dot{l}}\dot{k} + c\dot{\theta}\right],\text{ in }D\times[0,\infty);\tag{16}$$

$$\begin{aligned} \mathbf{u}\_{i}(\mathbf{x},t) &= \mathbf{0}, \ (\mathbf{x},t) \in \partial D\_{\mathbf{u}} \times [0,\infty), \ t\_{i} = \mathbf{0}, \ (\mathbf{x},t) \in \partial D\_{\mathbf{u}}^{\mathbb{C}} \times [0,\infty), \\ \phi\_{ij}(\mathbf{x},t) &= \mathbf{0}, \ (\mathbf{x},t) \in \partial D\_{\theta} \times [0,\infty), \ m\_{jk} = \mathbf{0}, \ (\mathbf{x},t) \in \partial D\_{\theta}^{\mathbb{C}} \times [0,\infty), \\ \theta(\mathbf{x},t) &= \mathbf{0}, \ (\mathbf{x},t) \in \partial D\_{\theta} \times [0,\infty), \ q = \mathbf{0}, \ (\mathbf{x},t) \in \partial D\_{\theta}^{\mathbb{C}} \times [0,\infty). \end{aligned} \tag{17}$$

#### **3. Main Result**

We first establish some integral identities regarding a solution **u** = u*i*, *φij*, *θ* of the mixed problem P . These will be useful in obtaining the important results of our study.

**Proposition 1.** *If the ordered array* **u** = u*i*, *φij*, *θ satisfies the mixed problem* P *, then the following equality takes place:*

$$\begin{split} \int\_{B} \left[ \frac{1}{2} \left( \rho \dot{\mathbf{u}}\_{i}(t) \dot{\mathbf{u}}\_{i}(t) + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}(t) \dot{\boldsymbol{\phi}}\_{km}(t) \right) + \Psi(t) + \frac{1}{2} c \theta^{2}(t) \right] dV &= \\ = \int\_{B} \left[ \rho \dot{\mathbf{u}}\_{i}(0) \dot{\mathbf{u}}\_{i}(0) + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}(0) \dot{\boldsymbol{\phi}}\_{km}(0) + \Psi(0) + \frac{1}{2} c \theta^{2}(0) \right] dV + \\ &+ \int\_{0}^{t} \int\_{D} \frac{1}{T\_{0}} K\_{ij} \theta\_{,i}(\tau) \theta\_{,j}(\tau) dV d\tau, \,\forall t \in [0, \infty). \end{split} \tag{18}$$

**Proof.** Taking into account the kinematics compatibility relations (2) and the differential conditions of equilibrium (9a) and (9b), we obtain the following equality:

$$\begin{split} \frac{1}{2} \frac{\partial}{\partial t} \left( \rho \dot{\mathbf{u}}\_i(t) \dot{\mathbf{u}}\_i(t) + I\_{jk} \boldsymbol{\phi}\_{jm}(t) \boldsymbol{\phi}\_{km}(t) \right) &= \\ = \left[ \left( \boldsymbol{\pi}\_{ij} + \boldsymbol{\eta}\_{ij} \right) \dot{\mathbf{u}}\_j + \mu\_{ijk} \dot{\boldsymbol{\phi}}\_{jk} \right]\_{,i} - \left( \boldsymbol{\pi}\_{ij} \dot{\boldsymbol{\varepsilon}}\_{ij} + \boldsymbol{\eta}\_{ij} \dot{\mathbf{x}}\_{ij} + \mu\_{ijk} \dot{\boldsymbol{\chi}}\_{ijk} \right). \end{split} \tag{19}$$

Taking into account the constitutive Equation (4), the symmetry relations (8), and the expression of the internal energy density Ψ from (13), the last parentheses in the right-hand side of (19) becomes

$$
\left(\tau\_{i\bar{j}}\dot{\varepsilon}\_{l\bar{j}} + \eta\_{i\bar{j}}\dot{\kappa}\_{l\bar{j}} + \mu\_{i\bar{j}k}\dot{\chi}\_{l\bar{j}k}\right) = \frac{\partial}{\partial t}\left(\Psi + \frac{1}{2}c\theta^2\right) + \left(\frac{1}{T\_0}q\_{\bar{j}}\theta\right)\_{,\bar{j}} - \frac{1}{T\_0}K\_{\bar{j}l}\theta\_{,l}\theta\_{,\bar{l}}.\tag{20}
$$

We substitute Equation (20) into Equation (19), then the resulting equality is integrated on cylinder [0, *t*] × *D*. If we use the theorem of divergence and consider the conditions to the limit (17), we are led to equality (18), such that the proof of the proposition is finished.

In a similar way, one can demonstrate the identity that follows, as a complement to identity (18):

$$\begin{split} \int\_{D} \left[ \frac{1}{2} \left( \rho \dot{\mathbf{u}}\_{i}(t) \dot{\mathbf{u}}\_{i}(t) + I\_{jk} \dot{\boldsymbol{\theta}}\_{jm}(t) \dot{\boldsymbol{\theta}}\_{km}(t) \right) + \Psi(t) - \frac{1}{2} c \boldsymbol{\theta}^{2}(t) \right] dV &= \\ \quad + \int\_{D} \left[ \rho \dot{\mathbf{u}}\_{i}(0) \dot{\mathbf{u}}\_{i}(0) + I\_{jk} \dot{\boldsymbol{\theta}}\_{jm}(0) \dot{\boldsymbol{\theta}}\_{km}(0) + \Psi(0) - \frac{1}{2} c \boldsymbol{\theta}^{2}(0) \right] dV &- \\ \quad - \int\_{0}^{t} \int\_{D} \left\{ \dot{\mathbf{u}}\_{i}(\tau) \left[ \left( \dot{\mathbf{u}}\_{ij} + b\_{ji} \right) \boldsymbol{\theta}(\tau) \right]\_{,j} + \dot{\boldsymbol{\theta}}\_{ij}(\tau) \left[ c\_{ijk} \boldsymbol{\theta}(\tau) \right]\_{,k} - \\ \quad - b\_{ij} \dot{\boldsymbol{\theta}}\_{ij}(\tau) \boldsymbol{\theta}(\tau) + \frac{1}{T\_{0}} K\_{ij} \theta\_{,i}(\tau) \boldsymbol{\theta}\_{,j}(\tau) \right\} dV \tau\_{\prime}, \forall t \in [0, \infty). \end{split} \tag{21}$$

To simplify writing, we enter the notation

$$\begin{split} 2F(\mathbf{x},\mathbf{y}) &= A\_{ijmn}\varepsilon\_{ij}(\mathbf{x})\varepsilon\_{mn}(\mathbf{y}) + D\_{ijmn}\left[\varepsilon\_{ij}(\mathbf{x})\kappa\_{mn}(\mathbf{y}) + \varepsilon\_{ij}(\mathbf{y})\kappa\_{mn}(\mathbf{x})\right] + \\ &+ F\_{ijmnr}\left[\varepsilon\_{ij}(\mathbf{x})\chi\_{mnr}(\mathbf{y}) + \varepsilon\_{ij}(\mathbf{y})\chi\_{mnr}(\mathbf{x})\right] + B\_{ijmn}\kappa\_{ij}(\mathbf{x})\kappa\_{mn}(\mathbf{y}) + \\ &+ G\_{ijmnr}\left[\kappa\_{ij}(\mathbf{x})\chi\_{mnr}(\mathbf{y}) + \kappa\_{ij}(\mathbf{y})\chi\_{mnr}(\mathbf{x})\right] + C\_{ijkmr}\chi\_{ijk}(\mathbf{x})\chi\_{mnr}(\mathbf{y}). \end{split} \tag{22}$$

By using the symmetry relations (8), from (22) we deduce

$$F(\mathbf{x}, \mathbf{y}) = F(\mathbf{y}, \mathbf{x}). \tag{23}$$

By direct substitution in (22) and taking into account Equation (13), we also obtain

$$F(\mathbf{r}, \mathbf{r}) = \Psi(\mathbf{r}).\tag{24}$$

Now, we can prove a result similar to (18), but in the case of homogeneous initial conditions.

**Proposition 2.** *Consider a solution* u*i*, *φij*, *θ of the problem backward in time* P *, which corresponds to null initial data. Then, the next identity takes place*

$$\begin{split} \int\_{B} \left[ \frac{1}{2} \left( \rho \dot{\mathbf{u}}\_{i}(t) \dot{\mathbf{u}}\_{i}(t) + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}(t) \dot{\boldsymbol{\phi}}\_{km}(t) \right) - \frac{1}{2} c \theta^{2}(t) \right] dV &= \\ = \int\_{B} \left[ \frac{1}{2} A\_{ijmn} \varepsilon\_{ij}(t) \varepsilon\_{mn}(t) + D\_{ijmn} \varepsilon\_{ij}(t) \kappa\_{nm}(t) + \\ &+ F\_{ijmn} \varepsilon\_{ij}(t) \chi\_{mmr}(t) + \frac{1}{2} B\_{ijmn} \kappa\_{ij}(t) \kappa\_{mn}(t) + \\ &+ G\_{ijmnr} \kappa\_{ij}(t) \chi\_{mmr}(t) + \frac{1}{2} C\_{ijkmnr} \chi\_{ijk}(t) \chi\_{mmr}(t) \right] dV, \end{split} \tag{25}$$

*for all t* ∈ [0, ∞)*.*

**Proof.** By direct calculations, for a fixed *t* ∈ (0, ∞), we get the identity:

$$\begin{split} \frac{\partial}{\partial t} \left( \rho \ddot{\mathbf{u}}\_{i}(\tau) \ddot{\mathbf{u}}\_{i}(2t-\tau) + I\_{jk} \dot{\phi}\_{jm}(\tau) \dot{\phi}\_{km}(2t-\tau) - c\theta(\tau)\theta(2t-\tau) \right) &= \\ -\rho \ddot{\mathbf{u}}\_{i}(\tau) \ddot{\mathbf{u}}\_{i}(2t-\tau) + I\_{jk} \ddot{\phi}\_{jm}(\tau) \dot{\phi}\_{km}(2t-\tau) + c\theta(\tau) \dot{\theta}(2t-\tau) - \\ -\rho \ddot{\mathbf{u}}\_{i}(\tau) \ddot{\mathbf{u}}\_{i}(2t-\tau) + I\_{jk} \dot{\phi}\_{jm}(\tau) \ddot{\phi}\_{km}(2t-\tau) - c\dot{\theta}(\tau) \theta(2t-\tau). \end{split} \tag{26}$$

Taking into account the kinematic Equation (2), the constitutive Equation (4), the motion Equation (9), and the symmetry relations (8), we are led to the equality

$$\begin{split} \frac{\partial}{\partial t} \left( \rho \ddot{\mathbf{u}}\_{i}(\tau) \dot{\mathbf{u}}\_{i}(2t-\tau) + l\_{jk} \dot{\boldsymbol{\theta}}\_{jm}(\tau) \dot{\boldsymbol{\theta}}\_{km}(2t-\tau) - c\theta(\tau)\theta(2t-\tau) \right) &= \\ = \left[ \left( \tau\_{ij} + \eta\_{ij} \right) \left( \tau \right) \dot{\mathbf{u}}\_{i}(2t-\tau) - \left( \tau\_{ij} + \eta\_{ij} \right) \left( 2t-\tau \right) \dot{\boldsymbol{u}}\_{j}(\tau) + \\ \quad + \mu\_{ijk} \ddot{\boldsymbol{\theta}}\_{jk}(\tau) \dot{\boldsymbol{\theta}}\_{jk}(2t-\tau) - \mu\_{ijk} \ddot{\boldsymbol{\theta}}\_{jk}(2t-\tau) \dot{\boldsymbol{\theta}}\_{jk}(\tau) - \\ \quad - \frac{1}{T\_{0}} \theta(\tau) q\_{i}(2t-\tau) + \frac{1}{T\_{0}} \theta(2t-\tau) q\_{i}(\tau) \right]\_{,i} + F(\tau, 2t-\tau), \end{split} \tag{27}$$

the function *F*(., .) being defined in (22).

We just need to integrate this equality into [0, *t*] × *D*, to keep in mind that the initial data are null and to use the definition (13), so we get the desired equality (25) and the proof of proposition is complete.

The following two propositions are also useful in establishing the main outcomes of our study.

**Proposition 3.** *Consider that the ordered array* **u** = u*i*, *φij*, *θ satisfies the mixed problem* P *. Then, the following equality takes place*

$$\begin{split} \int\_{D} \left[ \rho \mathbf{u}\_{i}(t) \dot{\mathbf{u}}\_{i}(t) + I\_{jk} \boldsymbol{\theta}\_{jm}(t) \dot{\boldsymbol{\theta}}\_{km}(t) - \frac{1}{2T\_{0}} K\_{ij} \left( \int\_{0}^{t} \boldsymbol{\theta}(\tau) d\tau \right)\_{,i} \left( \int\_{0}^{t} \boldsymbol{\theta}(\tau) d\tau \right)\_{,j} \right] dV &= \\ = \int\_{D} \left[ \rho \mathbf{u}\_{i}(0) \dot{\mathbf{u}}\_{i}(0) + I\_{jk} \boldsymbol{\theta}\_{jm}(0) \dot{\boldsymbol{\theta}}\_{km}(0) \right] dV &+ \int\_{0}^{t} \int\_{D} \rho \eta(0) \boldsymbol{\theta}(\tau) dV d\tau + \\ &+ \int\_{0}^{t} \int\_{D} \left[ \rho \dot{\mathbf{u}}\_{i}(t) \dot{\mathbf{u}}\_{i}(t) + I\_{jk} \dot{\boldsymbol{\theta}}\_{jm}(t) \dot{\boldsymbol{\theta}}\_{km}(t) - 2\Psi(\tau) - c\theta^{2}(\tau) \right] dV d\tau, \end{split} \tag{28}$$

*for all t* ∈ [0, ∞)*.*

**Proof.** We will consider the constitutive relations (4), the geometric Equation (2), and the motion Equations (9a) and (9b), we deduce

$$\begin{split} \frac{\partial}{\partial t} \Big( \rho \mathbf{u}\_{i}(t) \mathbf{\dot{u}}\_{i}(t) + I\_{jk} \boldsymbol{\phi}\_{jm}(t) \dot{\boldsymbol{\phi}}\_{km}(t) \Big) &= \\ &= \rho \mathbf{\dot{u}}\_{i}(t) \mathbf{\dot{u}}\_{i}(t) + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}(t) \dot{\boldsymbol{\phi}}\_{km}(t) + \\ &\quad + \left[ \left( \mathbf{\dot{r}}\_{ij}(t) + \eta\_{ij}(t) \right) \mathbf{u}\_{j}(t) + \mu\_{ij\bar{k}}(t) \boldsymbol{\phi}\_{\bar{j}\bar{k}}(t) \right]\_{,j} - \\ &\quad - \left[ \left( \mathbf{\dot{r}}\_{ij}(t) + \eta\_{ij}(t) \right) \boldsymbol{\varepsilon}\_{ij}(t) + \eta\_{ij}(t) \kappa\_{ij}(t) + \mu\_{ijk}(t) \boldsymbol{\chi}\_{ijk}(t) \right]. \end{split} \tag{29}$$

Taking into account definition (13), the last parentheses from the right-hand side of (29) can be restated in the following form:

$$\begin{split} \pi\_{ij}(t)\varepsilon\_{ij}(t) + \eta\_{ij}(t)\kappa\_{ij}(t) + \mu\_{ijk}(t)\chi\_{ijk}(t) &= \\ &= \left(\frac{1}{2}c\theta^2(t) + \Psi(t)\right) + \left(\frac{1}{T\_0}\theta(t)\int\_0^t q\_i(\tau)d\tau\right)\_{,i} \\ &- \frac{1}{T\_0}\mathcal{K}\_{ij}\left(\int\_0^t \theta(\tau)d\tau\right)\_{,i}\left(\int\_0^t \theta(\tau)d\tau\right)\_{,i} - \rho\eta(0)\theta(t). \end{split} \tag{30}$$

Here, the expression of the entropy *η* is obtained by integrating its equation of evolution (19), with regards to time variable.

*Symmetry* **2019**, *11*, 863

We integrate the equality (30) on [0, *t*] × *D*, and take into account the data on the border (17) and the theorem of divergence. As such, we arrive at the proposed equality (29), and the demonstration of proposition is over.

**Proposition 4.** *Let us consider a solution* u*i*, *φij*, *θ of the mixed problem* P *. Then, the following equality takes place*

$$\begin{split} 2\int\_{D} \left[\rho \mathbf{u}\_{i}(t)\dot{\mathbf{u}}\_{i}(t) + I\_{jk}\dot{\boldsymbol{\rho}}\_{jm}(t)\dot{\boldsymbol{\rho}}\_{km}(t) - \frac{1}{2T\_{0}}K\_{ij}\left(\int\_{0}^{t} \theta(\tau)d\tau\right)\_{,j} \left(\int\_{0}^{t} \theta(\tau)d\tau\right)\_{,j} \right]dV \\ &= \int\_{D} \left[\rho \dot{\mathbf{u}}\_{i}(0)\dot{\mathbf{u}}\_{i}(2t) + I\_{jk}\dot{\boldsymbol{\rho}}\_{jm}(0)\dot{\boldsymbol{\rho}}\_{km}(2t)\right]dV + \\ &\quad + \int\_{D} \left[\rho \mathbf{u}\_{i}(0)\dot{\mathbf{u}}\_{i}(2t) + I\_{jk}\dot{\boldsymbol{\rho}}\_{jm}(0)\dot{\boldsymbol{\rho}}\_{km}(2t)\right]dV - \\ &\quad - \int\_{0}^{t} \int\_{D} \rho \eta(0) \left[\theta(t+\tau) - \theta(t-\tau)\right]dVd\tau, \end{split} \tag{31}$$

*for all t* ∈ [0, ∞)*.*

**Proof.** We will take into account the geometric relations (2) and the motion Equations (9a) and (9b), we deduce the following identity

$$\begin{split} \frac{\partial}{\partial t} \left[ \rho \left( \dot{\mathbf{u}}\_{i}(t+\tau) \mathbf{u}\_{i}(t-\tau) + \mathbf{u}\_{i}(t+\tau) \dot{\mathbf{u}}\_{i}(t-\tau) \right) \right] + \\ &+ \frac{\partial}{\partial t} \left[ I\_{jk} \left( \dot{\phi}\_{jm}(t+\tau) \phi\_{km}(t-\tau) + \phi\_{jm}(t+\tau) \dot{\phi}\_{km}(t-\tau) \right) \right] = \\ &= \left[ \left( \tau\_{ij} + \eta\_{ij} \right) \left( t+\tau \right) \mathbf{u}\_{i}(t-\tau) - \left( \tau\_{ij} + \eta\_{ij} \right) \left( t-\tau \right) \mathbf{u}\_{j}(t+\tau) + \\ &+ \mu\_{ijk}(t+\tau) \dot{\boldsymbol{\phi}}\_{jk}(t-\tau) - \mu\_{ijk}(t-\tau) \boldsymbol{\phi}\_{jk}(t+\tau) \right]\_{,i} - \\ &- \left[ \left( \tau\_{ij} + \eta\_{ij} \right) \left( t+\tau \right) \boldsymbol{\varepsilon}\_{ij}(t-\tau) - \left( \tau\_{ij} + \eta\_{ij} \right) \left( t-\tau \right) \boldsymbol{\varepsilon}\_{ij}(t+\tau) \right] - \\ &- \left[ \mu\_{ijk}(t+\tau) \boldsymbol{\chi}\_{ijk}(t-\tau) - \mu\_{ijk}(t-\tau) \boldsymbol{\chi}\_{jk}(t+\tau) \right]. \end{split} \tag{32}$$

On the other hand, by using the symmetry Equation (8) and the constitutive relations (4), the last two brackets receive the following form:

$$\begin{aligned} \left[ \left( \tau\_{ij} + \eta\_{ij} \right) \left( t + \tau \right) \varepsilon\_{ij} (t - \tau) - \left( \tau\_{ij} + \eta\_{ij} \right) \left( t - \tau \right) \varepsilon\_{ij} (t + \tau) \right] + \\ &+ \left[ \mu\_{ijk} (t + \tau) \chi\_{ijk} (t - \tau) - \mu\_{ijk} (t - \tau) \chi\_{jk} (t + \tau) \right] = \\ &= \frac{1}{T\_0} \left[ \theta(t + \tau) \int\_0^{t - \tau} q\_i(s) ds - \theta(t - \tau) \int\_0^{t + \tau} q\_i(s) ds \right]\_j - \\ &- \frac{1}{T\_0} K\_{ij} \left[ \left( \int\_0^{t + \tau} \theta(s) ds \right)\_j \left( \int\_0^{t - \tau} \theta(s) ds \right)\_j - \left( \int\_0^{t + \tau} \theta(s) ds \right)\_j \left( \int\_0^{t - \tau} \theta(s) ds \right)\_j \right] - \\ &- \rho \eta \left( 0 \right) \left[ \theta(t + \tau) - \theta(t - \tau) \right]. \end{aligned} (33)$$

Here, the expression of the entropy *η* is obtained by integrating its equation of evolution (19), with regards to time variable.

Let us integrate the equality (33) on [0, *t*] × *D*, and take into account the border data (17) and the theorem of divergence. As such, we arrive at the proposed equality (31) and the demonstration of the proposition is over.

If we combine the results from Equations (25) and (31), then we obtain a new useful equality

$$\begin{split} 2\int\_{B} \left[ \left( \rho \dot{\mathbf{u}}\_{i}(t) \dot{\mathbf{u}}\_{i}(t) + I\_{j\mathbb{R}} \dot{\boldsymbol{\phi}}\_{jm}(t) \dot{\boldsymbol{\phi}}\_{km}(t) \right) - c\theta^{2}(t) \right] dV &= \\ = -2 \int\_{D} \left[ \rho \mathbf{u}\_{i}(0) \dot{\mathbf{u}}\_{i}(0) + I\_{j\mathbb{R}} \boldsymbol{\Phi}\_{jm}(0) \dot{\boldsymbol{\phi}}\_{km}(0) \right] dV + \\ &+ \int\_{D} \left[ \rho \dot{\mathbf{u}}\_{i}(0) \dot{\mathbf{u}}\_{i}(2t) + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}(0) \dot{\boldsymbol{\phi}}\_{km}(2t) \right] dV + \\ &+ \int\_{D} \left[ \rho \mathbf{u}\_{i}(0) \dot{\mathbf{u}}\_{i}(2t) + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}(0) \dot{\boldsymbol{\phi}}\_{km}(2t) \right] dV + \\ &- \int\_{0}^{t} \int\_{D} \rho \eta(0) \left[ 2\theta(\tau) + \theta(t + \tau) - \theta(t - \tau) \right] dV d\tau, \end{split} \tag{34}$$

for all *t* ∈ [0, ∞).

Based on the previously demonstrated integral identities, we are able to address the main results of our study. First, we established a result of uniqueness for the solution of the backward in time problem. As a consequence, we approach the question of localization of the solutions of the backward in time problem.

**Theorem 1.** *At most, an ordered array* **u** = u*i*, *φij*, *θ can satisfy the equations and conditions of the backward problem* P *.*

**Proof.** As usual, we will assume, by absurdum, that the problem would admit two solutions. The difference of the two solutions is also a solution, because the problem P is a linear one. Suffice it to show that this difference is null. For this we have to show that the problem P , for which the boundary and initial data are null, admits the null solution. It is clear that for the difference of two solutions, the boundary and initial conditions become homogeneous.

To simplify writing, we introduce the function *M* as a measure of the solution, defined by

$$M(t) = \int\_{D} \left[ \frac{\varepsilon}{2} \left( \rho \dot{\mathbf{u}}\_{i}(t) \dot{\mathbf{u}}\_{i}(t) + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}(t) \dot{\boldsymbol{\phi}}\_{km}(t) \right) + (\varepsilon + 2) \Psi(t) + \frac{\varepsilon}{2} c \theta^{2}(t) \right]. \tag{35}$$

Here, *ε* is a small positive number.

Based on assumptions (a), (b), and (c), it can be deduced that the function *M* is positive. Because the initial data are zero, the identity (18) received the simpler form:

$$\begin{cases} \int\_{B} \left[ \frac{1}{2} \left( \rho \dot{\mathbf{u}}\_{i}(t) \dot{\mathbf{u}}\_{i}(t) + I\_{jk} \dot{\boldsymbol{\theta}}\_{jm}(t) \dot{\boldsymbol{\theta}}\_{km}(t) \right) + \Psi(t) + \frac{1}{2} c \theta^{2}(t) \right] dV = \\ \int\_{0}^{t} \int\_{D} \frac{1}{\tau\_{0}} K\_{ij} \theta\_{,i}(\tau) \theta\_{,j}(\tau) dV d\tau, \; \forall t \in [0, \infty). \end{cases} \tag{36}$$

Analogously, the identity (21) becomes

$$\begin{split} \int\_{B} \left[ \frac{1}{2} \left( \rho \dot{\mathbf{u}}\_{i}(t) \dot{\mathbf{u}}\_{i}(t) + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}(t) \dot{\boldsymbol{\phi}}\_{km}(t) \right) + \Psi(t) - \frac{1}{2} \epsilon \theta^{2}(t) \right] dV &= \\ = - \int\_{0}^{t} \int\_{D} \left\{ \dot{\mathbf{u}}\_{i}(\tau) \left[ \left( a\_{ij} + b\_{ji} \right) \theta(\tau) \right]\_{,j} + \dot{\theta}\_{ij}(\tau) \left[ c\_{ij} \mathbf{k} \theta(\tau) \right]\_{,k} - \\ - b\_{ij} \dot{\theta}\_{ij}(\tau) \theta(\tau) + \frac{1}{T\_{0}} K\_{ij} \theta\_{,i}(\tau) \theta\_{,j}(\tau) \right\} dV d\tau, \ \forall t \in [0, \infty). \end{split} \tag{37}$$

If we use Equations (36) and (37), then the function *M* from (35) becomes

$$\begin{split} M(t) &= -\int\_{0}^{t} \int\_{D} \left\{ 2\dot{\mathfrak{u}}\_{i}(\boldsymbol{\tau}) \left[ \left( a\_{ij} + b\_{ji} \right) \boldsymbol{\theta}(\boldsymbol{\tau}) \right]\_{,j} + \dot{\phi}\_{ij}(\boldsymbol{\tau}) \left[ c\_{ijk} \boldsymbol{\theta}(\boldsymbol{\tau}) \right]\_{,k} - \\ &- b\_{ij} \dot{\phi}\_{ij}(\boldsymbol{\tau}) \boldsymbol{\theta}(\boldsymbol{\tau}) + \frac{1-\varepsilon}{T\_{0}} K\_{ij} \boldsymbol{\theta}\_{j}(\boldsymbol{\tau}) \boldsymbol{\theta}\_{,j}(\boldsymbol{\tau}) \right\} dV d\boldsymbol{\tau}, \; \forall t \in [0, \infty). \end{split} \tag{38}$$

*Symmetry* **2019**, *11*, 863

By direct derivation with regards to the variable *t* in (37), we are led to the following equality:

$$\begin{cases} \frac{dM(t)}{dt} = -2\int\_{D} \left\{ \dot{\mathbf{u}}\_{i}(\boldsymbol{\tau}) \left[ \left( a\_{ij} + b\_{ji} \right) \theta(\boldsymbol{\tau}) \right]\_{,j} + \dot{\theta}\_{ij}(\boldsymbol{\tau}) \left[ c\_{ijk} \theta(\boldsymbol{\tau}) \right]\_{,k} - \\\ -b\_{ij} \phi\_{ij}(\boldsymbol{\tau}) \theta(\boldsymbol{\tau}) + \frac{1-\varepsilon}{2T\_{0}} K\_{ij} \theta\_{,i}(\boldsymbol{\tau}) \theta\_{,j}(\boldsymbol{\tau}) \right\} dV, \ \forall t \in [0, \infty). \end{cases} \tag{39}$$

With the help of Schwarz' inequality and by using the arithmetic–geometric mean inequality, from (39) we can deduce the following inequality:

$$\frac{dM(t)}{dt} \le \mathcal{C}\_1 \int\_D \left[\rho \dot{\mathbf{u}}\_i(t)\dot{\mathbf{u}}\_i(t) + I\_{jk}\dot{\boldsymbol{\phi}}\_{jm}(t)\dot{\boldsymbol{\phi}}\_{km}(t) + c\theta^2(t)\right]dV + \\ + \frac{\delta - 1 + \varepsilon}{T\_0} \int\_D K\_{ij}\theta\_{,i}(\tau)\theta\_{,j}(\tau)dV, \ \forall t \in [0, \infty). \tag{40}$$

Now, we take into account that the internal energy density Ψ is a positive definite quadratic form, according to the hypothesis (c), using the definition (35) of the function *M* and choose *δ* ≤ 1 − *ε*. Then, from (40) we obtain the inequality:

$$\frac{dM(t)}{dt} \le \frac{\mathcal{C}\_1}{\varepsilon} \int\_D \varepsilon \left[ \rho \dot{\mathbf{u}}\_i(t) \dot{\mathbf{u}}\_i(t) + I\_{jk} \dot{\phi}\_{jm}(t) \dot{\phi}\_{km}(t) + c\theta^2(t) \right] dV \le \frac{\mathcal{C}\_1}{\varepsilon} M,\tag{41}$$

and it is then clear that a solution to this inequality meets the next inequality:

$$0 \le M(t) \le M(0)e^{\mathcal{C}\_1/\varepsilon}.\tag{42}$$

We recall that for the difference of the two supposed solutions, the initial conditions are homogeneous, then we have *M*(0) = 0, so from (42) we get

$$M(t) = 0, \,\forall t \in [0, \infty)$$

and this together with the assumptions leads to the conclusion that our problem has only the solution

$$\mathfrak{u}\_i(t) = 0, \; \mathfrak{\phi}\_{ij}(t) = 0, \; \mathfrak{\theta}(t) = 0, \; \forall t \in [0, \infty), i$$

and the proof of Theorem 1 is concluded.

Our final result is dedicated to the partition of various energies associated with the solution of the backward in time problem P∗. We recall that this problem consists of the equations of motion (9), the constitutive relations (4), the kinematic Equation (2), the initial data (12), and the boundary restrictions in their homogeneous form (11).

First, using the procedure outlined at the end of Section 2, we transform the boundary-final value problem P∗ into the boundary-initial value problem P . In this way, in what follows, we will be able to make the considerations only on the problem P . Let us denote by T the set of of those thermoelastodynamic processes defined in the cylinder (−∞, 0] × *D*, which satisfy the restriction

$$\int\_{D} \frac{1}{T\_0} \mathbb{K}\_{\vec{i}\vec{j}} \theta\_{,i}(\tau) \theta\_{,\vec{j}}(\tau) dV \le \mathbb{C}\_1. \tag{43}$$

for all *t* ∈ [0, ∞). Here, *C*1 is a given positive constant.

*Symmetry* **2019**, *11*, 863

We will introduce the known Cesaro means necessary to evaluate the various types of energies that can be attached to a solution to the problem P . So, if u*i*, *φij*, *θ* is a solution of the mixed problem, then the Cesaro means are:

$$\begin{split} K(t) &= \frac{1}{t} \int\_{0}^{t} \int\_{D} \left[ \rho \dot{\mathbf{u}}\_{i}(\tau) \dot{\mathbf{u}}\_{i}(\tau) + I\_{jk} \dot{\phi}\_{jm}(\tau) \dot{\phi}\_{jm}(\tau) \right] dV d\tau; \\ S(t) &= \frac{1}{t} \int\_{0}^{t} \int\_{D} \left[ \frac{1}{2} A\_{ijmn} \varepsilon\_{ij}(\tau) \varepsilon\_{mn}(\tau) + D\_{ijmn} \varepsilon\_{ij}(\tau) \chi\_{mn}(\tau) + F\_{ijmn} \varepsilon\_{ij}(\tau) \chi\_{mmr}(\tau) + \\ &\quad + \frac{1}{2} B\_{ijmn} \kappa\_{ij}(\tau) \chi\_{mn}(\tau) + G\_{ijmn} \kappa\_{ij}(\tau) \chi\_{mmr}(\tau) + \frac{1}{2} C\_{ijkmn} \chi\_{ijk}(\tau) \chi\_{mmr}(\tau) \right] dV d\tau; \\ R(t) &= \frac{1}{t} \int\_{0}^{t} \int\_{D} \frac{1}{2c} c \,\theta^{2}(\tau) dV d\tau; \\ D(t) &= \frac{1}{t} \int\_{0}^{t} \int\_{0}^{\tau} \int\_{D} \frac{1}{T\_{0}} K\_{ij} \theta\_{,j}(s) \theta\_{,j}(s) dV ds \, d\tau. \end{split} \tag{44}$$

In the particular case when meas(*∂Du*)= 0 and meas(*∂Dφ*)= 0, it can be determined a rigid displacement, a rigid dipolar displacement, and a null temperature, which satisfy the equations of motion (9), the constitutive relations (4), the kinematic relations (2), and verify the homogeneous boundary conditions (17). In this case, the initial data can be decomposed as follows:

$$\begin{aligned} \mathbf{u}\_i^0 &= \mathbf{u}\_i' + V\_i^0 \; \prime \; \dot{\mathbf{u}}\_i^0 = \dot{\mathbf{u}}\_i' + V\_i^0 \; \prime \\ \boldsymbol{\Phi}\_{ij}^0 &= \boldsymbol{\Phi}\_{ij}' + \Psi\_{ij\prime}^0 \; \dot{\boldsymbol{\Phi}}\_{ij}^0 = \dot{\boldsymbol{\Phi}}\_{ij}' + \Psi\_{ij}^0 \; \prime \end{aligned} \tag{45}$$

The rigid displacements u *i* , ˙u *i* , *φ ij*, *φ*˙ *ij* can be computed using the functions

$$\mathcal{W}\_{\vec{l}}(\omega) = \int\_{D} \rho \omega\_{\vec{l}} Dv,\ \mathcal{W}\_{\vec{j}\vec{k}}(\psi) = \int\_{D} I\_{mn} \psi\_{nj} \psi\_{nk} dV.$$

such that we have

$$\begin{aligned} \mathcal{W}\_i(V\_i^0) &= 0, \; \mathcal{W}\_i(\dot{V}\_i^0) = 0, \\ \mathcal{W}\_{jk}(\Psi\_{ij}^0) &= 0, \; \mathcal{W}\_{jk}(\Psi\_{ij}^0) = 0. \end{aligned} \tag{46}$$

Together with known notation *C*1(*D*), we will use the notation *Wn*(*D*) for a Sobolev space defined on the domain D, and **W***n*(*D*)=[*Wn*(*D*)]3.

Other new notations:

$$\begin{aligned} \mathcal{C}^1(D) &= \left\{ (\mathbf{u}\_i, \boldsymbol{\phi}\_{\hat{\boldsymbol{\eta}}}) \in \mathbb{C}^1(D)^3 \times \mathbb{C}^1(D)^9 : \mathbf{u}\_i = 0 \text{ on } \partial D\_{\mathbf{u}}, \ \boldsymbol{\phi}\_{\hat{\boldsymbol{\eta}}} = 0 \text{ on } \partial D\_{\boldsymbol{\theta}}; \\ \text{if } \text{meas}(\partial D\_{\boldsymbol{u}}) &= 0 \text{ and } \text{meas}(\partial D\_{\boldsymbol{\phi}}) = 0, \text{ then } \boldsymbol{W}\_i(\mathbf{u}\_i) = 0, \ \boldsymbol{W}\_{jk}(\boldsymbol{\phi}\_{\hat{\boldsymbol{\eta}}}^0) = 0 \right\}; \\ \tilde{\mathcal{C}}^1(D) &= \left\{ \theta \in \mathcal{C}^1(D) : \theta = 0 \text{ on } \partial D\_{\boldsymbol{\theta}} \right\}; \\ \mathcal{W}\_1(D) &\text{ the completion of } \mathcal{C}^1(D); \\ \tilde{\mathcal{W}}\_1(D) &\text{ the completion of } \mathcal{C}^1(D), \end{aligned}$$

the completion is by means of the original norm of the respective Sobolev space. Based on definition (13), and taking into account the hypothesis (15), we can obtain the following inequality (see [34]):

$$\int\_{D} \Psi(\mathbf{u})dV \ge \frac{M\_1}{2} \int\_{D} \left[\rho \mathbf{u}\_i \mathbf{u}\_i + I\_{jk} \phi\_{jm} \phi\_{km}\right] dV. \tag{47}$$

for any **u** = u*i*, *φij* ∈ W1(*D*). *M*1 > 0 is a convenient chosen constant.

On the other hand, if we take into account the hypothesis (14), we can obtain the next inequality of Poincare' type:

$$\int\_{D} \mathcal{K}\_{i\bar{j}} \theta\_{,i} \theta\_{,\bar{j}} dV \geq M\_2 \int\_{D} \theta^2 dV,\tag{48}$$

for any *θ* ∈ *W*˜ 1(*D*). *M*2 > 0 is a convenient chosen constant.

Let us consider a solution u*i*, *φij*, *θ* to the problem P in the particular case when meas(*∂Du*)= 0 and meas(*∂Dφ*)= 0. We can represent this solution in the following form:

$$\begin{cases} \mathbf{u}\_{i}(t, \mathbf{x}) = u\_{i}(t, \mathbf{x}) + u\_{i}^{\prime}(t, \mathbf{x}) + t \dot{u}\_{i}^{\prime}(t, \mathbf{x}), \; \boldsymbol{\theta}(t, \mathbf{x}) = \boldsymbol{\theta}(t, \mathbf{x}),\\ \boldsymbol{\phi}\_{\overline{i}\overline{j}}(t, \mathbf{x}) = \boldsymbol{\psi}\_{\overline{i}\overline{j}}(t, \mathbf{x}) + \boldsymbol{\phi}\_{\overline{i}\overline{j}}^{\prime}(t, \mathbf{x}) + t \dot{\boldsymbol{\phi}}\_{\overline{i}\overline{j}}^{\prime}(t, \mathbf{x}), \; \boldsymbol{\theta}(t, \mathbf{x}) \in [0, \infty) \times D,\end{cases} \tag{49}$$

in which *ui*, *ψij*, *ϑ* ∈ W1(*D*) × *W*˜ 1(*D*). In addition, the ordered array *ui*, *ψij*, *ϑ* satisfies the problem P , which corresponds to the following initial data:

$$\begin{aligned} \mu\_i(0, \mathbf{x}) &= V\_i(\mathbf{x}), \; \dot{\imath}\_i(0, \mathbf{x}) = \dot{V}\_i(\mathbf{x}), \; \theta(t, \mathbf{x}) = \theta(\mathbf{x}),\\ \Phi\_{\dot{\imath}\jmath}(0, \mathbf{x}) &= \Psi\_{\dot{\imath}\jmath}(\mathbf{x}), \; \dot{\phi}\_{\dot{\imath}\jmath}(0, \mathbf{x}) = \Psi\_{\dot{\imath}\jmath}(\mathbf{x}), \; \forall \mathbf{x} \in D. \end{aligned}$$

We now have everything prepared to address the problem of the equipartition of the various types of energies associated with the solution to the problem P .

**Theorem 2.** *Consider a solution* u*i*, *φij*, *θ of the backward in time problem* P *. If the initial data satisfy the following conditions:*

$$\begin{aligned} \mathbf{u} &= (\mathbf{u}\_i) \in \mathbf{W}\_1(D), \; \dot{\mathbf{u}} = (\dot{\mathbf{u}}\_i) \in \mathbf{W}\_0(D), \\ \boldsymbol{\Phi} &= (\dot{\boldsymbol{\phi}}\_{\dot{\boldsymbol{\eta}}}) \in \mathbf{W}\_1(D), \; \dot{\boldsymbol{\Phi}} = (\dot{\boldsymbol{\phi}}\_{\dot{\boldsymbol{\eta}}}) \in \mathbf{W}\_0(D), \; \boldsymbol{\theta} \in \mathbf{W}\_0(D), \end{aligned}$$

*then the following three statements are true:*

*i) The thermal component of energy, R, vanishes as t* → ∞*:*

$$\lim\_{t \to \infty} R(t) = 0;\tag{50}$$

*ii) if meas(∂Du)*= 0 *and meas(∂Dφ)*= 0*, then we have*

$$\lim\_{t \to \infty} K(t) = \lim\_{t \to \infty} S(t) + \frac{1}{2} \int\_{D} \left[ \rho \dot{\mathbf{u}}\_{i}^{\prime} \dot{\mathbf{u}}\_{i}^{\prime} + I\_{jk} \dot{\phi}\_{jm}^{\prime} \dot{\phi}\_{km}^{\prime} \right] dV\_{\cdot} \tag{51}$$

$$\begin{split} \lim\_{t \to \infty} D(t) &= 2 \lim\_{t \to \infty} K(t) - \frac{1}{2} \int\_{D} \left[ \rho \dot{\mathbf{u}}\_{i}^{\prime} \dot{\mathbf{u}}\_{i}^{\prime} + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}^{\prime} \dot{\boldsymbol{\phi}}\_{km}^{\prime} \right] dV - E(0) = \\ &= 2 \lim\_{t \to \infty} S(t) + \frac{1}{2} \int\_{D} \left[ \rho \dot{\mathbf{u}}\_{i}^{\prime} \dot{\mathbf{u}}\_{i}^{\prime} + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}^{\prime} \dot{\boldsymbol{\phi}}\_{km}^{\prime} \right] dV - E(0); \end{split} \tag{52}$$

*iii) if meas(∂Du)*= 0 *or meas(∂Dφ)*= 0*, then we have*

$$\begin{aligned} \lim\_{t \to \infty} \mathcal{K}(t) &= \lim\_{t \to \infty} \mathcal{S}(t), \\ \lim\_{t \to \infty} D(t) &= 2 \lim\_{t \to \infty} \mathcal{K}(t) - E(0) = \\ &= 2 \lim\_{t \to \infty} \mathcal{S}(t) - E(0), \end{aligned} \tag{53}$$

*where, to simplify writing in (52) and (53), we used the notation:*

$$E(t) = \int\_{D} \left[ \frac{1}{2} \left( \rho \dot{\mathbf{u}}\_{i}(t) \dot{\mathbf{u}}\_{i}(t) + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}(t) \dot{\boldsymbol{\phi}}\_{km}(t) \right) + \Psi(t) + \frac{1}{2} c \theta^{2}(t) \right] dV,\tag{54}$$

*Symmetry* **2019**, *11*, 863

*the internal energy density* Ψ *being defined in (13).*

**Proof.** Taking into account the relations (18), (54), and (44) we immediately deduce the following equality:

$$K(t) + \mathcal{S}(t) + R(t) = D(t) + E(0), \ \forall t \in (0, \infty). \tag{55}$$

Now, we consider the restrictions (43) and (48), the definitions (44), and the notation (54) in order to deduce the following estimation:

$$R(t) \le \frac{T\_0 C\_1}{2M\_2 t} \max\_{\mathbf{x} \in D} \{ c(\mathbf{x}) \}, \ \forall t \in (0, \infty). \tag{56}$$

Clearly, since *c*(*x*) is a continuous function, we deduce that max *x*∈*D* {*c*(*x*)} is bounded, such that after we pass to the limit in (56) for *t* → ∞, we obtain the result (50).

If we consider the identity (34) and use the notations (44), we are led to the the following identity:

$$\begin{split} &K(t) - S(t) - R(t) = \\ &= -\frac{1}{2t} \int\_{D} \left[ \rho \mathbf{u}\_{i}(0) \dot{\mathbf{u}}\_{i}(0) + I\_{jk} \Phi\_{jm}(0) \dot{\Phi}\_{km}(0) \right] dV + \\ &\quad + \frac{1}{4t} \int\_{D} \left[ \rho \dot{\mathbf{u}}\_{i}(0) \mathbf{u}\_{i}(2t) + I\_{jk} \dot{\Phi}\_{jm}(0) \Phi\_{km}(2t) \right] dV + \\ &\quad + \frac{1}{4t} \int\_{D} \left[ \rho \mathbf{u}\_{i}(0) \dot{\mathbf{u}}\_{i}(2t) + I\_{jk} \Phi\_{jm}(0) \dot{\Phi}\_{km}(2t) \right] dV + \\ &\quad - \frac{1}{4t} \int\_{0}^{t} \int\_{D} \rho \eta(0) \left[ 2\theta(\tau) + \theta(t+\tau) - \theta(t-\tau) \right] dV d\tau, \end{split} \tag{57}$$

for all *t* ∈ [0, ∞).

On the other hand, taking into account the hypothesis (15), the inequalities (43) and (48), and the identity (18) we obtain the estimates

$$\int\_{D} \rho \dot{\mathbf{u}}\_{i}(t) \dot{\mathbf{u}}\_{i}(t) dV \le 2(\mathbb{C}\_{1} + E(0)),\\\int\_{D} I\_{jk} \dot{\Phi}\_{jm}(t) \dot{\Phi}\_{km}(t) dV \le 2(\mathbb{C}\_{1} + E(0)),\\\int\_{D} 2\Psi(t) dV \le 2(\mathbb{C}\_{1} + E(0)),\tag{58}$$

where *c*0 is from hypothesis (a), *C*1 is from (43), and *E* from (54).

In (57) we use the Schwarz's inequality and taking into account the estimates (58), we pass the limit—for *t* → ∞—so we get the equality

$$\lim\_{t \to \infty} K(t) = \lim\_{t \to \infty} S(t) + \lim\_{t \to \infty} \frac{1}{4t} \int\_{D} \left[ \rho \dot{\mathbf{u}}\_i(0) \mathbf{u}\_i(2t) + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}(0) \dot{\boldsymbol{\phi}}\_{km}(2t) \right] dV. \tag{59}$$

In this way, it will be easy to demonstrate the relation (53), if we show that the integral from the right-hand side of the identity (59) is bounded.

For this aim, we will use the fact that meas(*∂Du*)= 0 or meas(*∂Dφ*)= 0 and *ui*, *ψij* ∈ W1(*D*). Furthermore, we consider the relations (47), (54), and (18) in order to get the following estimates:

$$\begin{aligned} M\_1 \int\_D \rho \mathbf{u}\_i(t) \mathbf{u}\_i(t) dV &\le \int\_D 2\Psi(t) dV \le 2(\mathcal{C}\_1 + E(0)),\\ M\_1 \int\_D I\_{jk} \phi\_{jm}(t) \phi\_{km}(t) dV &\le \int\_D 2\Psi(t) dV \le 2(\mathcal{C}\_1 + E(0)), \end{aligned}$$

such that it is easy, after we apply the Schwarz's inequality, to deduce

$$\lim\_{t \to \infty} \frac{1}{4t} \int\_{D} \left[ \rho \dot{\mathbf{u}}\_{i}(0) \mathbf{u}\_{i}(2t) + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}(0) \boldsymbol{\phi}\_{km}(2t) \right] dV = 0. \tag{60}$$

From (59) and (60) we obtain the first result from (53), and after we consider the equality (55), the second equality from (53) is proven.

Finally, we will prove the equalities (51) and (52). In this regard, we take into account that meas(*∂Du*)= 0 and meas(*∂Dφ*)= 0 such that with the help of the decompositions (45) and (49), and conditions (46), we are led to the equality

$$\begin{split} \frac{1}{4t} \int\_{D} \left[ \rho \dot{\mathbf{u}}\_{i}(0) \mathbf{u}\_{i}(2t) + I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}(0) \boldsymbol{\phi}\_{km}(2t) \right] dV &= \\ = \frac{1}{4t} \int\_{D} \rho \dot{\mathbf{u}}\_{i}^{\prime} \dot{\mathbf{u}}\_{i}^{\prime}(2t) dV + \frac{1}{4t} \int\_{D} \rho \left( \dot{\mathbf{u}}\_{i}^{\prime} + \dot{V}\_{i}^{0} \right) \boldsymbol{u}\_{i}(2t) dV + \\ &+ \frac{1}{2} \int\_{D} \rho \dot{\mathbf{u}}\_{i}^{\prime} \dot{\mathbf{u}}\_{i}^{\prime} dV + \frac{1}{4t} \int\_{D} I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}^{\prime} \boldsymbol{\phi}\_{km}^{\prime} dV + \\ &+ \frac{1}{4t} \int\_{D} I\_{jk} \left( \dot{\mathbf{s}}\_{jm}^{\prime} + \dot{\mathbf{v}}\_{jm}^{0} \right) \boldsymbol{\psi}\_{jk}(2t) dV + \frac{1}{2} \int\_{D} I\_{jk} \dot{\boldsymbol{\phi}}\_{jm}^{\prime} \dot{\boldsymbol{\phi}}\_{km}^{\prime} dV. \end{split} \tag{61}$$

Let us observe that relations (47), (54), and (18) involve the estimates

$$\begin{aligned} M\_1 \int\_D \rho u\_i(t) u\_i(t) dV &\le 2(C\_1 + E(0)), \\ M\_1 \int\_D I\_{jk} \psi\_{jm}(t) \psi\_{km}(t) dV &\le 2(C\_1 + E(0)), \end{aligned}$$

such that from (61) we are led to the equality

$$\begin{array}{ll}\displaystyle\lim\_{t\to\infty}\frac{1}{4t}\int\_{D}\left[\rho\dot{\mathbf{u}}\_{i}(0)\mathbf{u}\_{i}(2t)+I\_{jk}\dot{\boldsymbol{\phi}}\_{jm}(0)\boldsymbol{\phi}\_{km}(2t)\right]dV=\\ =\frac{1}{2}\int\_{D}\rho\dot{\mathbf{u}}\_{i}^{\prime}\dot{\mathbf{u}}\_{i}^{\prime}dV+\frac{1}{2}\int\_{D}I\_{jk}\dot{\boldsymbol{\phi}}\_{jm}^{\prime}\dot{\boldsymbol{\phi}}\_{km}^{\prime}dV.\end{array} \tag{62}$$

Now, we substitute (62) into Equation (59) and then obtain the relation (51). Lastly, we consider the relations (50), (51), and (55) in order to obtain (52). With this, the proof of the theorem is completed.

#### **4. Conclusions**

We first formulate the mixed backward in time problem in the context of thermoelasticity for dipolar materials. To prove the consistency of this mixed problem, our first main result is regarding the uniqueness of the solution for this problem. This is obtained based on some auxiliary results, namely, four integral identities. The second main result is regarding the temporal behavior of our thermoelastic body with dipolar structure. This behavior is studied by means of some relations on partition of various parts of the energy associated to the solution of the problem. After we introduce the Cesaro means for all parts of the total energy, we can evaluate the asymptotic partition of these parts. We have to say that the kinetic energy and potential energy become asymptotically equal when the variable time tends to infinity.

**Author Contributions:** All four authors conceived the framework and structured the whole manuscript, checked the results and completed the revision of the paper. The authors have equally contributed to the elaboration of this manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to express their gratitude to the anonymous reviewers for their constructive comments and suggestions, which led to an improved form of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Two-Dimensional Finite Element in General Plane Motion Used in the Analysis of Multi-Body Systems**

#### **Eliza Chircan 1, Maria-Lumini¸ta Scutaru 1,\* and Cătălin Iulian Pruncu 2,3**

Received: 5 June 2019; Accepted: 27 June 2019; Published: 1 July 2019

**Abstract:** A description of the motion equation of a single two-dimensional finite element used to model a multi-body system with elastic elements is made in the article. To establish them, the Lagrange's equations are used. Obtaining the dynamic response of a system with deformable components has become important for technical applications in recent decades. These engineering applications are characterized by high applied loads and high acceleration and velocities. A study of such mechanical systems leads to the identification of different mechanical phenomena (due to high deformations, resonance phenomena, and stability). Coriolis effects and relative motions significantly modify the motion equations and, implicitly, the dynamic response. These effects are highlighted in this paper for plane motion.

**Keywords:** multi-body system; finite element method (FEM); linear elastic elements; Lagrange's equations; two-dimensional finite element; plane motion

#### **1. Introduction**

Finite Element Analysis (FEA) of elastic bodies with overall rigid motion began in the 1970s. A number of papers analyzed different cases and types of finite elements used in this type of analysis [1–3]. One-dimensional and three-dimensional finite elements have been developed with plane motion or three-dimensional motion [4,5]. The field has remained interesting for researchers so far, with many papers on different aspects of the issue being published. The new field of multi-body systems (MBS), which has developed in recent decades, has reached the study of mechanical systems with (some) linear/nonlinear elastic elements. Why is this type of analysis necessary? If the velocities, accelerations, or loads that a constitutive element of a multi-body system supports are high, the elasticity of the elements may cause instability phenomena. The classical hypothesis of the rigid elements, commonly used in the dynamic analysis and synthesis of multi-body systems, no longer corresponds to the reality. The phenomena of resonance and loss of stability represent classic forms of manifestation of the elasticity properties. The classical approach involves the use of continuous mechanical models. Unfortunately, this type of analysis leads to a nonlinear system of differential equations that cannot be solved analytically in common engineering applications. The most convenient approach is FEA, a time-validated and frequently used method for analyzing mechanical systems with elastic elements. However, classical models of static or quasi-static systems do not correspond to the analysis of multi-body systems, which have rigid motion. The relatively large amount of motion between elements and Coriolis effects lead to the existence of additional terms in motion equations that cannot be neglected. These terms are usually strongly nonlinear. As a consequence, proper modeling of the problem is required [6,7]. The many aspects of a one-dimensional problem are described in [8]. In the following research, we develop a method for a two-dimensional finite element with plane motion. The motion equations for this type of element are established.

#### **2. Two-Dimensional Finite Element**

#### *2.1. Two-Dimensional Model*

Let us consider a two-dimensional finite element with a body with the general rigid plane motion of an MBS. To determine the governing equations for this element, the Lagrange's equations are used. The main phase in this type of approach is to obtain the Lagrangian for the one two-dimensional finite element. The terms necessary to build the Lagrangian are the kinetic energy of the two-dimensional finite element, the internal energy, and the work done by the concentrated and distributed loads. The shape function (determined in every case by the type of finite element chosen and the hypothesis made) will finally determine the form of the matrix coefficients and, as a consequence, the differential system of the governing motion equations. A basic hypothesis in this study is that the deformations of the elements of MBS are small enough not to influence the plane rigid motion of the system. Both the problem of the plane multi-body rigid motion of the system and velocities and the acceleration field distribution for each element can be solved using the classical method (see [9–11]).

The finite element chosen is related to the local reference frame Ox1x2 (participating with the finite element in the general plane rigid motion) and to the global reference frame O'X1X2. It is considered that the vector velocity and acceleration of the origin O related to the global reference system are known from a previous dynamical analysis. All the bodies are considered to be rigid as is the angular velocity and angular acceleration of the local reference system. In this way, the fields of velocity and acceleration are known for each point of each constitutive element of the MBS [12–18].

For one single two-dimensional finite element, the generalized independent coordinates can be the nodal displacements, depending on the shape function and hypothesis used to express the displacement of a point in the finite element. The final number of independent coordinates will depend on this hypothesis and shape function.

Consider the displacement {δ(*u, v*)} of an arbitrary point M of the domain of the finite element. If we use the shape function *Nij* and the vector of the nodal displacements δ*e,j*, the displacements in the local coordinate system are

$$u = \delta\_1 = N\_{1\rangle} \delta\_{\mathfrak{e}, \mathfrak{j}}; \quad v = \delta\_2 = N\_{2\rangle} \delta\_{\mathfrak{e}, \mathfrak{j}}; \quad j = 1, p. \tag{1}$$

or

$$
\delta\_{\mathbf{i}} = N\_{\mathbf{i}\mathbf{j}} \delta\_{\mathbf{t}, \mathbf{j}} \quad \mathbf{i} = \mathbf{1}, \mathbf{2}; \quad \mathbf{j} = \overline{\mathbf{1}, p}, \tag{2}
$$

where the nodal displacement vector of the finite element is labeled e, and δ*e,j* depends on the independent coordinate chosen to define the element. The number of independent coordinates is *p*. For a triangular finite element we have

$$
\delta\_\varepsilon^T = [\boldsymbol{\mu}\_1 \ \boldsymbol{\upsilon}\_1 \ \boldsymbol{\mu}\_2 \ \boldsymbol{\upsilon}\_2 \ \boldsymbol{\mu}\_3 \ \boldsymbol{\upsilon}\_3].\tag{3}
$$

The two lines of the shape function matrix **N** correspond to the displacements *u* and *v* and are named N(u) = N(1) and N(v) = N(2):

$$N = \begin{bmatrix} N\_{(u)} \\ N\_{(v)} \end{bmatrix} = \begin{bmatrix} N\_{(1)} \\ N\_{(2)} \end{bmatrix} = \begin{bmatrix} N\_{\bar{i}\bar{j}} \end{bmatrix}, \quad i = 1, 2; \quad j = \overline{1, p}. \tag{4}$$

The displacement of the point M(*x*1,*x*2), which is arbitrarily chosen, becomes, after deformation, M'(*x'*1,*x'*2) and can be expressed through its components:

$$\mathbf{x}'\_1 = \mathbf{x}\_1 + \boldsymbol{\mu} = \mathbf{x}\_1 + \boldsymbol{\delta}\_1;\ \mathbf{x}'\_2 = \mathbf{x}\_2 + \boldsymbol{\upsilon} = \mathbf{x}\_2 + \boldsymbol{\delta}\_2,\tag{5}$$

or with respect to the global reference system:

$$\begin{aligned} X\_1' &= X\_{10} + r\_{1i} \mathbf{x}\_i' = X\_{1o} + r\_{1i} \mathbf{x}\_i + r\_{1i} \delta\_i = X\_{1o} + r\_{1i} \mathbf{x}\_i + r\_{1i} N\_{i\bar{j}} \delta\_{\bar{e},\bar{j}\prime} \\ X\_2' &= X\_{20} + r\_{2i} \mathbf{x}\_i' = X\_{2o} + r\_{2i} \mathbf{x}\_i + r\_{2\bar{i}} \delta\_{\bar{i}} = X\_{2o} + r\_{2i} \mathbf{x}\_i + r\_{2\bar{i}} N\_{\bar{i}\bar{j}} \delta\_{\bar{e},\bar{j}} \delta\_{\bar{e},\bar{j}} \end{aligned} \tag{6}$$

or

$$X\_k' = X\_{k\nu} + r\_{k\bar{\imath}} x\_{\bar{\imath}} + r\_{k\bar{\imath}} N\_{\bar{\imath}\bar{\jmath}} \delta\_{\varepsilon\_{\bar{\imath}}\bar{\jmath}} \quad k, \bar{\imath} = 1, 2, \quad j = \overline{1, p}. \tag{7}$$

#### *2.2. Lagrangian of an Element*

Using (6), the velocity of point M' can be obtained after differentiation with respect to the time:

$$\dot{X}\_k = \dot{X}\_{k\flat} + \dot{r}\_{\text{ki}} \mathbf{x}\_{\text{i}} + \dot{r}\_{\text{ki}} \mathbf{N}\_{\text{ij}} \delta\_{\mathfrak{e},\text{j}} + r\_{\text{ki}} \mathbf{N}\_{\text{ij}} \dot{\delta}\_{\mathfrak{e},\text{j}}, \quad k, \mathfrak{i} = 1, 2, \quad \mathfrak{j} = \overline{1, p}. \tag{8}$$

The kinetic energy for one single element due to this velocity is

$$\begin{aligned} E\_{\rm tI} &= \frac{1}{2} \int \rho A \dot{\mathbf{X}}\_{k}^{\prime} \dot{\mathbf{X}}\_{k}^{\prime} dA \\ &= \frac{1}{2} \int \rho t \Big( \dot{\mathbf{X}}\_{k\nu} + \dot{r}\_{\rm li} \mathbf{x}\_{i} + \dot{r}\_{\rm li} \mathbf{N}\_{\rm ij} \delta\_{\varepsilon,j} + r\_{\rm li} \mathbf{N}\_{\rm ij} \dot{\delta}\_{\varepsilon,i} \Big) \Big( \dot{\mathbf{X}}\_{k\nu} + \dot{r}\_{\rm li} \mathbf{x}\_{\rm I} + \dot{r}\_{\rm I} \mathbf{N}\_{\rm Im} \delta\_{\varepsilon,\mathbf{m}} + r\_{\rm iI} \mathbf{N}\_{\rm Im} \dot{\delta}\_{\varepsilon,\mathbf{m}} \Big) dA \end{aligned} \tag{9}$$

where *t* is the thickness of the element, and ρ is the mass density.

In plane motion, the rotation of the local system of coordinates is expressed by the plane rotation tensor *rij*. The columns of this tensor define the positions of the unit vectors of the local reference frame Oxyz. For plane rotation, these coefficients can be expressed as

$$r\_{11} = r\_{22} = \cos\theta; \; r\_{12} = r\_{21} = -\sin\theta. \tag{10}$$

The ortho-normality condition of these unit vectors leads to

$$r\_{\rm ij}r\_{k\rm j} = r\_{\rm jk}r\_{\rm ji} = \delta\_{\rm ij}, \text{ ijk} = 1,2 \tag{11}$$

where δ*ij* is the Kronecker delta (derived from the general three-dimensional transformation). By differentiation, the following equation will result:

$$
\dot{r}\_{i\bar{j}}r\_{k\bar{j}} + r\_{i\bar{j}}\dot{r}\_{k\bar{j}} = 0, \quad i,k = 1,2. \tag{12}
$$

Denote

$$
\omega\_{i\bar{k}} = \dot{r}\_{i\bar{j}} r\_{k\bar{j}\prime}.\tag{13}
$$

the skew-symmetric tensor angular velocity. The relation (12) becomes

$$
\omega\_{ik} + \omega\_{ki} = 0, \quad i,k = 1,2. \tag{14}
$$

To this corresponds the angular velocity vector components defined by

$$
\omega\_3 = \omega\_{21} = -\omega\_{12} = \omega,\ \omega\_1 = \omega\_{32} = -\omega\_{23} = 0,\ \omega\_2 = \omega\_{13} = -\omega\_{31} = 0. \tag{15}
$$

The angular velocity vector and the angular acceleration vector have, in plane motion, the same components in both the local reference system and the global reference system.

The angular acceleration skew symmetric operator it is defined as

$$
\varepsilon\_{\rm ik} = \dot{\omega}\_{\rm ik} = \bar{r}\_{\rm ij} r\_{k\rm j} + \dot{r}\_{\rm ij} \dot{r}\_{k\rm j}, \,\mathrm{ijk} = \mathbf{1}, \mathbf{2}. \tag{16}
$$

The angular acceleration vector is defined by

$$
\varepsilon\_3 = \varepsilon\_{21} = -\varepsilon\_{12} = \dot{\omega} = \varepsilon. \tag{17}
$$

The other components are null. We shall have

$$
\varepsilon\_{i\mathbf{k}} = \omega\_{i\mathbf{k}} = \ddot{r}\_{\mathbf{i}\mathbf{j}} r\_{\mathbf{k}\mathbf{j}} + \dot{r}\_{\mathbf{i}\mathbf{j}} \dot{r}\_{\mathbf{k}\mathbf{j}} = \ddot{r}\_{\mathbf{i}\mathbf{j}} r\_{\mathbf{k}\mathbf{j}} + \dot{r}\_{\mathbf{i}\mathbf{j}} r\_{\mathbf{j}} r\_{\mathbf{k}\mathbf{l}} \dot{r}\_{\mathbf{k}\mathbf{m}} = \ddot{r}\_{\mathbf{i}\mathbf{j}} r\_{\mathbf{k}\mathbf{j}} - \omega\_{\mathbf{i}\mathbf{l}} \omega\_{\mathbf{l}\mathbf{k}} \tag{18}
$$

from where ..

$$
\ddot{r}\_{\text{ij}} r\_{k\text{j}} = \varepsilon\_{\text{ik}} + \omega\_{\text{il}} \omega\_{\text{lk}}, \text{ ij}
\text{l} = 1, 2. \tag{19}
$$

The internal energy stored is obtained via the relation

$$E\_p = \frac{1}{2} \int\_0^L \sigma\_{ij} \varepsilon\_{ij} dV \quad i, j = 1, 2. \tag{20}$$

From (20), it is possible to obtain the stiffness matrix depending on the hypothesis and type of finite element (and consequently, the shape functions) chosen. Here, we present a case as an example; any other case can be treated in the same manner. The relationship between stress and strain is

$$
\sigma\_{\rm ij} = D\_{\rm ik} \varepsilon\_{\rm kj}. \tag{21}
$$

If we choose, for example, a rectangular finite element with sides *a* and *b*, a displacement field of the form

$$u = a\_1 + a\_2 \mathbf{x} + a\_3 \mathbf{x} \mathbf{y} + a\_4 \mathbf{y} \mathbf{v} = \beta\_1 + \beta\_2 \mathbf{x} + \beta\_3 \mathbf{x} \mathbf{y} + \beta\_4 \mathbf{y} \tag{22}$$

can be used to obtain the stiffness tensor.

The boundary conditions (*x* = 0, *y* = 0 = > *u* = *u*1, *v* = *v*1; *x* = *a*, *y* = 0 = > *u* = *u*2, *v* = *v*2; *x* = *a*, *y* = *b* = > *u* = *u*3, *v* = *v*3; *x* = 0, *y* = *b* = > *u* = *u*4, *v* = *v*4;) lead to the following fields for the displacements:

$$\begin{aligned} u &= (1 - \underline{\xi})(1 - \eta)u\_1 + \underline{\xi}(1 - \eta)u\_2 + \underline{\xi}\eta u\_3 + \eta(1 - \underline{\xi})u\_4 \\ v &= (1 - \underline{\xi})(1 - \eta)v\_1 + \underline{\xi}(1 - \eta)v\_2 + \underline{\xi}\eta v\_3 + \eta(1 - \underline{\xi})v\_4. \end{aligned} \tag{23}$$

Applying the classical laws of linear elasticity [2], the relation between the strain and nodal displacement are, in this case,

$$\begin{array}{c} \varepsilon\_{11} = -\frac{1-\eta}{a}u\_{1} + \frac{1-\eta}{a}u\_{2} + \frac{\eta}{a}u\_{3} - \frac{\eta}{a}u\_{4},\\ \varepsilon\_{22} = -\frac{1-\xi}{b}v\_{1} + -\frac{\xi}{b}v\_{2} + \frac{\xi}{b}v\_{3} + \frac{1-\xi}{b}v\_{4},\\ \varepsilon\_{12} = \varepsilon\_{21} = -\frac{1-\xi}{b}u\_{1} - \frac{1-\eta}{a}v\_{1} - \frac{\xi}{b}u\_{2} + \frac{1-\eta}{a}v\_{2} + \frac{\xi}{b}u\_{3} + \frac{\eta}{a}v\_{3} + \frac{1-\xi}{b}u\_{4} - \frac{\eta}{a}v\_{4} \end{array} \tag{24}$$

Using (20), (21), and (24), after some elementary calculus, it is possible to obtain the internal energy of the form (see [2])

$$E\_p = \frac{1}{2} \delta\_i k\_{\mathbf{e},ij} \delta\_j \quad i, j = \overline{1, p} \tag{25}$$

The external work of the concentrated load is

$$\mathcal{W} = q\_{\mathfrak{e}\_{\mathfrak{e}\_{\mathfrak{e}}}} \delta\_{\mathfrak{e}, \mathfrak{i}}, \quad \mathfrak{i} = 1, 2, \quad \mathfrak{j} = \overline{1, p}. \tag{26}$$

Here, *qej* is the vector of the concentrated loads acting in nodes.

*Symmetry* **2019**, *11*, 848

The external work of the distributed loads is

$$\mathcal{W} = \bigcup\_{0}^{L} (p\_i \delta\_i) dA = \left( \int\_{0}^{L} p\_i \mathcal{N}\_{ij} dA \right) \delta\_{\varepsilon, j} = q\_{\varepsilon, j}^\* \delta\_{\varepsilon, j}, \quad i = 1, 2, \quad j = \overline{1, p}. \tag{27}$$

where the notation

$$q\_{\mathfrak{c},j}^{\*} = \left(\int\_{0}^{L} p\_{i} \mathbf{N}\_{i\bar{j}} d\mathbf{x}\right) \delta\_{\mathfrak{c},j\circ} \quad i = 1,2, \quad j = \overline{1,p} \tag{28}$$

is used. We have all the data to build the Lagrangian [8,19]

$$L = E\_{\mathfrak{c}} - E\_p + \mathcal{W}^d + \mathcal{W}^c,\tag{29}$$

or, taking into account the relations (9), (24), (25), and (26),

$$L = \frac{1}{2} \int \rho t \left( \dot{\mathbf{X}}\_{ko} + \dot{r}\_{ki} \mathbf{x}\_{l} + \dot{r}\_{li} \mathbf{N}\_{lj} \delta\_{e,j} + r\_{li} \mathbf{N}\_{lj} \dot{\delta}\_{e,j} \right) \left( \dot{\mathbf{X}}\_{ko} + \dot{r}\_{kl} \mathbf{x}\_{l} + \dot{r}\_{kl} \mathbf{N}\_{lm} \delta\_{e,m} + r\_{li} \mathbf{N}\_{lm} \dot{\delta}\_{e,m} \right) dA \tag{30}$$

$$- \frac{1}{2} \delta\_{e,i} \mathbf{k}\_{e,ij} \delta\_{e,j} + q\_{e\mathbf{I}\_{e\uparrow}}^{\*} \delta\_{e,j} + q\_{e\mathbf{J}} \delta\_{e,i}.$$

#### **3. Motion Equations**

**Theorem**: *The motion equations written in the local coordinate system for two-dimensional finite element take the form*

$$m\_{\varepsilon,\bar{\imath}\rangle}\ddot{\boldsymbol{\delta}}\_{\varepsilon,\bar{\jmath}} + 2\boldsymbol{\varepsilon}\_{\varepsilon,\bar{\jmath}}^{\mu}\dot{\boldsymbol{\delta}}\_{\varepsilon,\bar{\jmath}} + \left(k\_{\varepsilon,\bar{\imath}\rangle} + k\_{\varepsilon,\bar{\jmath}}^{\varepsilon} + k\_{\varepsilon,\bar{\jmath}}^{\mu\bar{\jmath}}\right)\boldsymbol{\delta}\_{\varepsilon,\bar{\jmath}} = q\_{\varepsilon,i} + q\_{\varepsilon,i}^{\*} - q\_{\varepsilon,i}^{\varepsilon} - q\_{\varepsilon,i}^{\mu\bar{\jmath}} - m\_{\varepsilon,i}^{o}\ddot{\boldsymbol{x}}\_{\bar{\jmath}\boldsymbol{\varepsilon}}\tag{31}$$

*where*

$$m\_{\varepsilon;ij} = \int\_0^L \rho t N\_{\rm ki} N\_{\rm kj} dA, \quad i, j = \overline{1, p}, \quad k = 1, 2\\\varepsilon\_{\varepsilon;ij}^{\omega} = 2 \int\_0^L N\_{\rm k} \omega\_{\rm km} N\_{\rm mj} \rho t dA; \tag{32}$$

$$k\_{\varepsilon,ij}^{\varepsilon} = \bigcup\_{0}^{L} \mathbf{N}\_{kl} \varepsilon\_{km} \mathbf{N}\_{ml} \rho \mathbf{t} dA; \quad k\_{\varepsilon,ij}^{\mu 2} = \bigcup\_{0}^{L} \mathbf{N}\_{kl} \omega\_{km} \omega\_{ml} \mathbf{N}\_{lj} \rho \mathbf{t} dA; \quad m\_{\varepsilon,ij}^{\mu} = \int\_{0}^{L} \rho A \mathbf{N}\_{jl} d\mathbf{x}; \mathbf{j}, j = \overline{1, p}, \quad k, m, l = 1, 2. \tag{33}$$

**Proof**: We apply the Lagrange's equations in the form

$$
\frac{d}{dt} \frac{\partial L}{\partial \dot{\delta}\_{\varepsilon,i}} - \frac{\partial L}{\partial \delta\_{\varepsilon,i}} = 0. \tag{34}
$$

We obtain the motion equations (an extended presentation can be found in Appendix A):

$$\begin{aligned} \left[\int\_{0}^{L} \rho t \mathbf{N}\_{\bar{\text{R}}} \mathbf{N}\_{\bar{\text{R}}} dA \right] \ddot{\boldsymbol{\delta}}\_{c,j} &+ 2 \left[\int\_{0}^{L} \mathbf{N}\_{\bar{\text{R}}} \boldsymbol{\omega}\_{\text{km}} \mathbf{N}\_{\text{m}/\bar{\text{p}}} \rho t dA \right] \ddot{\boldsymbol{\delta}}\_{c,j} + \begin{bmatrix} \boldsymbol{k}\_{c,ij} + \int\_{0}^{L} \mathbf{N}\_{\bar{\text{R}}} \boldsymbol{\varepsilon}\_{\text{km}} \mathbf{N}\_{\text{m}/\bar{\text{p}}} \rho t dA + \\ \int\_{0}^{L} & \mathbf{L} \\ + \int\_{\bar{\text{R}}} \mathbf{N}\_{\bar{\text{R}}} \boldsymbol{\omega}\_{\text{km}} \boldsymbol{\omega}\_{\text{ml}} \mathbf{N}\_{\bar{\text{I}}} \rho t dA \end{bmatrix} \boldsymbol{\delta}\_{c,j} \\ &= q\_{c,j} + q\_{c,i}^{\*} - q\_{c,i}^{\*} - q\_{c,i}^{\mu^{2}} - \left(\int\_{0}^{L} \rho A \mathbf{N}\_{\bar{\text{J}}} \mathbf{x} \right) \ddot{\mathbf{x}}\_{\bar{\text{J}}}. \end{aligned} \tag{35}$$

Using the notation mentioned above, it is possible to obtain (30).

#### **4. Conclusions**

In engineering practice, there are frequent cases in which the elements of a device or a machine are required to operate with high forces or have high accelerations and speeds. In these cases, the elasticity of the elements can be manifested in such a way as to negatively influence the working process of the assembly. An analysis of the MBS motion in such conditions is required. In this work, we determined the motion equations for a two-dimensional finite element with plane motion. This case often occurs in practice. Most of the technical systems have elements operating in planar motion. In this case, a system of differential equations with matrix coefficients that can be strongly nonlinear is obtained. The mass tensor and the classical rigidity tensor are symmetrical matrices, while the Coriolis effects tensor is skew-symmetric. There are also terms that change the classical stiffness of the element and additional inertial terms. If we use FEA to obtain the motion equations for one single element in plane motion, we obtain a system of second-order differential equations, as presented in (30). Generally, this system is strongly nonlinear, with the matrix coefficients of the system depending on time and on the position of the system. For a usual engineering application, some aspects of these equations can cause difficulties. We highlight some properties of this system:

The common method used to solve this system is to linearize this system, considering the tensor coefficients as being constant for very short time intervals (rigid motion freezing). In this way, it is possible to obtain a system of differential equations with constant coefficients. To solve this, usual and well-known methods can be used. There can be singularities due to the inertia term affecting the stiffness matrix [20].

**Author Contributions:** Conceptualization, M.L.S.; methodology, E.C., S.M.L. and C.P.; investigation and software M.L.S.; validation, E.C., M.L.S.; formal analysis, M.L.S.; resources, M.L.S., C.P.; data curation, M.L.S., C.P.; writing—original draft preparation, M.L.S.; writing—review and editing, M.L.S.; visualization, M.L.S.; supervision, C.P.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to thank to the Transylvania University of Brasov and the Imperial College of London for providing support to help us to successfully complete the work.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

The derivatives of the Lagrangian:

$$\begin{split} \frac{\partial L}{\partial \dot{\delta}\_{c,i}} &= \int\_0^L \rho \dot{X}\_{\text{kr}} r\_{\text{kl}} N\_{\text{lm}} t dA + \int\_0^L \rho A \dot{r}\_{\text{kl}} \ge r\_{\text{kl}} N\_{\text{lm}} t dA + \int\_0^L \rho \dot{r}\_{\text{kl}} N\_{\text{lj}} \delta\_{c,j} r\_{\text{kl}} N\_{\text{lm}} t dA \\ &+ \int\_0^L \rho r\_{\text{kl}} N\_{\text{l}j} \dot{\delta}\_{c,j} r\_{\text{kl}} N\_{\text{lm}} t dA \\ -\frac{\partial L}{\partial \dot{\delta}\_{c,i}} &= -\dot{X}\_{\text{kr}} \dot{r}\_{\text{kl}} \Big{{} \Big{{}}^L \mathcal{N}\_{\text{lm}} \rho t dA - \dot{r}\_{\text{kl}} \dot{r}\_{\text{kl}} \Big{{}}^L \mathcal{N}\_{\text{l}i} N\_{\text{lm}} \rho t dA \\ &- \dot{r}\_{\text{kl}} r\_{\text{kl}} \Big{{}}^L \Big{{}} N\_{\text{l}j} N\_{\text{lm}} \rho t dA \Big{{} \Big{{}} \mathcal{N}\_{\text{l}c} \dot{r}\_{\text{kl}} \Big{{}}^L \mathcal{N}\_{\text{l}c} \dot{r}\_{\text{kl}} \Big{{}} \delta\_{c,j} \end{split} \tag{A1}$$

$$\begin{aligned} \text{Applying } \frac{d}{dt} \Big( \frac{\rho L}{\partial \dot{\phi}\_{ij}} \Big) - \frac{\rho L}{\partial \dot{\phi}\_{ij}} &= 0 \text{ obtains} \\\\ \left( \int\_0^L N\_{\text{Im}} \rho t dA \right)\_{\text{Fl}} \bar{X}\_{\text{Ib}} + \bar{r}\_{\text{li}} r\_{\text{li}} \Big( \int\_0^L \bar{x} N\_{\text{Im}} \rho t dA \Big) + \bar{r}\_{\text{li}} r\_{\text{li}} \Big( \int\_0^L N\_{\text{ij}} N\_{\text{Im}} \rho t dA \Big) \delta\_{c,j} + 2 \bar{r}\_{\text{li}} r\_{\text{li}} \Big( \int\_0^L N\_{\text{ij}} N\_{\text{ij}} \rho t dA \Big) \delta\_{c,j} \\ &+ r\_{\text{li}} r\_{\text{li}} \Big( \int\_0^L N\_{\text{li}} N\_{\text{Im}} \rho t dA \Big) \bar{\delta}\_{c,j} - (q\_{c,i} + q\_{c,i}^\*) + k\_{c,i} \delta\_{c,j} = 0; \\ \left( \int\_0^L \rho t N\_{\text{ik}} N\_{\text{jk}} d\Lambda \right) \bar{\delta}\_{c,j} &+ 2\omega \Big( \int\_0^L \rho t N\_{\text{2}} N\_{\text{1}} \rho A d - \int\_0^L \rho t N\_{\text{1i}} N\_{\text{2}} \rho A d \Big) \delta\_{c,j} + \\ &\left[ k\_{c,j} + \varepsilon \Big( \int\_0^L \rho t N\_{\text{2}} N\_{\text{1i}} M\_{\text{1}} - \int\_0^L \rho t N\_{\text{1i}} N\_{\text{1i}} M\_{\text{1}} + \int\_0^L \rho t N\_{\text{2}} N\_{\text{2}} \rho A d \$$

or

$$m\_{\varepsilon;\bar{j}}\ddot{\delta}\_{\varepsilon,\bar{j}} + 2\omega \Big( m\_{\mathrm{i}\bar{j}}^{(21)} - m\_{\mathrm{i}\bar{j}}^{(12)} \Big) \dot{\delta}\_{\varepsilon,\bar{j}} + \left[ k\_{\varepsilon;\bar{j}} + \varepsilon \Big( m\_{\mathrm{i}\bar{j}}^{(21)} - m\_{\mathrm{i}\bar{j}}^{(12)} \Big) - \omega^{2} \Big( m\_{\mathrm{i}\bar{j}}^{(11)} + m\_{\mathrm{i}\bar{j}}^{(22)} \Big) \right] \delta\_{\varepsilon,\bar{j}} = \tag{A3}$$

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Inference about the Ratio of the Coe**ffi**cients of Variation of Two Independent Symmetric or Asymmetric Populations**

#### **Zhang Yue 1 and Dumitru Baleanu 2,3,\***

Received: 23 May 2019; Accepted: 20 June 2019; Published: 21 June 2019

**Abstract:** Coefficient of variation (CV) is a simple but useful statistical tool to make comparisons about the independent populations in many research areas. In this study, firstly, we proposed the asymptotic distribution for the ratio of the CVs of two separate symmetric or asymmetric populations. Then, we derived the asymptotic confidence interval and test statistic for hypothesis testing about the ratio of the CVs of these populations. Finally, the performance of the introduced approach was studied through simulation study.

**Keywords:** coefficient of variation; ratio; symmetric and asymmetric distributions; test of hypothesis

#### **1. Introduction**

Based on the literature, to describe a dataset (random variable), three main characteristics containing central tendencies, dispersion tendencies and shape tendencies, are used. A central tendency (or measure of central tendency) is a central or typical value for a random variable that describes the way in which the random variable is clustered around a central value. It may also be called a center or location of the distribution of the random variable. The most common measures of central tendency are the mean, the median and the mode. Measures of dispersion like the range, variance and standard deviation tell us about the spread of the values of a random variable. It may also be called a scale of the distribution of the random variable. The shape tendencies such as skewness and kurtosis describe the distribution (or pattern) of the random variable.

The division of the standard deviation to the mean of population, *CV* = σ μ , is called as coefficient of variation (CV) which is an applicable statistic to evaluate the relative variability. This free dimension parameter can be widely used as an index of reliability or variability in many applied sciences such as agriculture, biology, engineering, finance, medicine, and many others [1–3]. Since it is often necessary to relate the standard deviation to the level of the measurements, the CV is a widely used measure of dispersion. The CVs are often calculated on samples from several independent populations, and questions about how to compare them naturally arise, especially when the distributions of the populations are skewed. In real world applications, the researchers may intend to compare the CVs of two separate populations to understand the structure of the data. ANOVA and Levene tests can be used to investigate the equality of CVs of populations in case the means or variances of the populations are equal. It is obvious that in many situations two populations with different means and variances may have an equal CV. For the normal case, the problems of interval estimating the CV or comparison of two or several CVs have been well addressed in the literature. Due to possible small differences of two small CVs and no strong interpretation, the ratio of CVs is more accurate

than the difference of CVs. Bennett [4] proposed a likelihood ratio test. Doornbos and Dijkstra [5] and Hedges and Olkin [6] presented two tests based on the non-central t test. A modification of Bennett's method was provided by Shafer and Sullivan [7]. Wald tests have been introduced by [8–10]. Based on Renyi's divergence, Pardo and Pardo [11] proposed a new method. Nairy and Rao [12] applied the likelihood ratio, score test and Wald test to check that the inverse CVs are equal. Verrill and Johnson [13] applied one-step Newton estimators to establish a likelihood ratio test. Jafari and Kazemi [14] developed a parametric bootstrap (PB) approach. Some statisticians improved these tests for symmetric distributions [15–23]. The problem of comparing two or more CVs arises in many practical situations [24–26]. Nam and Kwon [25] developed approximate interval estimation of the ratio of two CVs for lognormal distributions by using the Wald-type, Fieller-type, log methods, and the method of variance estimates recovery (MOVER). Wong and Jiang [26] proposed a simulated Bartlett corrected likelihood ratio approach to obtain inference concerning the ratio of two CVs for lognormal distribution.

In applications, it is usually assumed that the data follows symmetric distributions. For this reason, most previous works have focused on the comparison of CVs in symmetric distributions, especially normal distributions. In this paper, we propose a method to compare the CVs of two separate symmetric or asymmetric populations. Firstly, we propose the asymptotic distribution for the ratio of the CVs. Then, we derive the asymptotic confidence interval and test statistic for hypothesis testing about the ratio of the CVs. Finally, the performance of the introduced approach is studied through simulation study. The introduced approach seems to have many advantages. First, it is powerful. Second, it is not too computational. Third, it can be applied to compare the CVs of two separate symmetric or asymmetric populations. We apply a methodology similar to that which has been used in [27–33]. The comparison between the parameters of two datasets or models has been considered in several works [34–40]

#### **2. Asymptotic Results**

Assume that *X* and *Y* are uncorrelated variables with non-zero means μ*X* and μ*Y*, and the finite *i th* central moments:

$$\mu\_{iX} = E(X - \mu \chi)^l, \quad \mu\_{iY} = E(Y - \mu\_Y)^l, \quad i \in \{2, 3, 4\},$$

respectively. Also assume two samples *X*1, ... , *Xm*, and*Y*1, ... , *Yn*, distributed from *X* and*Y*, respectively. From the motivation given in the introduction, the parameter:

$$\gamma = \frac{CV\_Y}{CV\_X}\_{\prime}$$

is interesting to inference, where *CVY* and *CVX* are the CVs corresponding to *Y* and *X*, respectively.

Assume *miX* = 1 *m* '*m k*=1 (*xk* − *x*) *i* , *miY* = 1 *n* '*n k*=1 (*yk* − *y*) *i* , *i* ∈ {2, 3, 4}. *CVX* and *CVY*, are consistently estimated [41] by *CV*ˆ *X* = √*m*2*X X* and *CV*ˆ *Y* = √*m*2*Y Y* , respectively. So, it is obvious that

$$\mathfrak{P} = \frac{\mathfrak{C}V\_Y}{\mathfrak{C}V\_X}$$

can reasonably estimate the parameter γ. For simplicity, let *m* = *n*. When *n* = *m*, let *n*∗ = *min*(*m*, *n*) instead of *m* and *n* in the following discussions.

**Lemma 1.** *If the above assumptions are satisfied, then:*

$$\sqrt{n}(\heartsuit V\_X - \heartsuit V\_X) \xrightarrow{L} N(0, \delta\_X^2)\_{\prime} \quad \text{as} \quad n \to \infty, \tau$$

*Symmetry* **2019**, *11*, 824

*where:*

$$
\delta\_X^2 = \left[ \frac{\mu\_{4X} - \mu\_{2X}^2}{4\mu\_X^2 \mu\_{2X}} - \frac{\mu\_{3X}}{\mu\_X^3} + \frac{\mu\_{2X}^2}{\mu\_{4X}} \right].
$$

*is the asymptotic variance.*

**Proof.** The outline of proof can be found in [41].

The next theorem corresponds to the asymptotic distribution of γˆ. This theorem will be applied to construct the confidence interval and perform hypothesis testing for the parameter γ.

**Theorem 1.** *If the previous assumptions are satisfied, then:*

$$\sqrt{n}(\wp - \wp) \xrightarrow{L} \mathcal{N}(0, \lambda^2), \quad \text{as} \ n \to \infty,$$

*where:*

$$
\lambda^2 = \frac{1}{CV\_X^2} (\nu^2 \delta\_X^2 + \delta\_Y^2),
$$

*and:*

$$
\delta\_Y^2 = \left[ \frac{\mu\_{4Y} - \mu\_{2Y}^2}{4\mu\_Y^2 \mu\_{2Y}} - \frac{\mu\_{3Y}}{\mu\_Y^3} + \frac{\mu\_{2Y}^2}{\mu\_{4Y}} \right].
$$

**Proof.** By using Lemma 1, we have:

$$\sqrt{n}(\operatorname{\mathbf{C}}V\_X - \operatorname{\mathbf{C}}V\_X) \xrightarrow{L} \operatorname{N}(\operatorname{0}, \delta\_X^2)\_{\prime} \quad \text{as} \ n \to \infty,$$

and:

$$\sqrt{n}(\mathring{C}V\_Y - CV\_Y) \xrightarrow{L} N(0, \delta\_Y^2)\_{\prime} \quad \text{as} \ n \to \infty.$$

Slutsky's Theorem gives:

$$
\sqrt{n}\left[\begin{pmatrix} \hat{C}V\_X\\ \hat{C}V\_Y \end{pmatrix} - \begin{pmatrix} CV\_X\\ CV\_Y \end{pmatrix} \right] \xrightarrow{L} N\left(0, \begin{pmatrix} \delta\_X^2 & 0\\ 0 & \delta\_Y^2 \end{pmatrix} \right).
$$

for independent samples [41].

Now define *f* : *R*2 → *R* as *f*(*x*, *y*) = *y x* . Then we have:

$$
\nabla f(\mathbf{x}, \mathbf{y}) = \left(-\frac{\mathcal{Y}}{\mathbf{x}^2}, \frac{1}{\mathbf{x}}\right).
$$

where ∇*f*(*x*, *y*) is the gradient function. Consequently, we have ∇*f*(*CVX*,*CVY*) '(∇*f*(*CVX*,*CVY*))*T* = λ2. Because of continuity of ∇*f* in the neighbourhood of (*CVX*,*CVY*), by using Cramer's Rule:

$$\sqrt{n}\Big(f\Big(\hat{\mathcal{C}}V\_{\mathcal{X}}, \hat{\mathcal{C}}V\_{Y}\Big) - f\Big(\mathcal{C}V\_{\mathcal{X}}, \mathbb{C}V\_{Y}\Big)\Big) = \sqrt{n}(\circ - \circ) \stackrel{L}{\to} \mathsf{N}\Big(0, \lambda^{2}\Big), \ n \to \infty,$$

the proof ends.

Thus, the asymptotic distribution can be constructed as:

$$T\_n = \sqrt{n} \Big( \frac{\mathfrak{P} - \mathfrak{I}'}{\lambda} \Big) \xrightarrow{L} N(0, 1), \quad \text{as} \ n \to \infty. \tag{1}$$

#### *2.1. Constructing the Confidence Interval*

As can be seen, the parameter λ depends on *CVX*, δ2 *X*, δ2 *Y* and γ which are unknown parameters in practice. The result of the next theorem can be applied to construct the confidence interval and to perform the hypothesis testing for the parameter γ.

**Theorem 2.** *If the previous assumptions are satisfied, then:*

$$T\_n^\* = \sqrt{n} \Big(\frac{\hat{\mathcal{I}} - \mathcal{V}}{\hat{\mathcal{A}}}\Big) \xrightarrow{L} N(0, 1), \quad \text{as } \to \infty,\tag{2}$$

*where:*

$$
\lambda^2 = \frac{1}{\widehat{\circ}V\_X^2} (\widehat{\circ}^2 \delta\_X^2 + \delta\_Y^2),
$$

$$
\delta\_X^2 = \left[ \frac{m\_{4X} - m\_{2X}^2}{4\overline{X}^2 m\_{2X}} - \frac{m\_{3X}}{\overline{X}^3} + \frac{m\_{2X}^2}{m\_{4X}} \right],
$$

$$
\delta\_X = \delta\_X \qquad \gamma = 1
$$

*and:*

$$
\delta\_Y^2 = \left[ \frac{m\_{4Y} - m\_{2Y}^2}{4\overline{Y}^2 m\_{2Y}} - \frac{m\_{3Y}}{\overline{Y}^3} + \frac{m\_{2Y}^2}{m\_{4Y}} \right].
$$

**Proof.** From the *Weak Law of Large Numbers*, it is known that:

$$
\overline{X} \xrightarrow{p} \mu\_{X\prime} \quad \overline{Y} \xrightarrow{p} \mu\_{Y\prime} \quad m\_{iX} \xrightarrow{p} \mu\_{iX\prime} \quad m\_{iY} \xrightarrow{p} \mu\_{iY\prime} \quad i \in \{2,3,4\},
$$

as *n* → ∞.

Consequently, by applying Slutsky's Theorem, we have λˆ *p* → λ, as *n* → ∞. Appliying Theorem 1 the proof is completed.

Now, *T*∗ *n* is a pivotal quantity for γ. In the following, this pivotal quantity is used to construct asymptotic confidence interval for γ.

$$\left(\hat{\gamma} - \frac{\lambda}{\sqrt{n}} Z\_{\frac{\varphi}{\xi}}, \hat{\gamma} + \frac{\lambda}{\sqrt{n}} Z\_{\frac{\varphi}{\xi}}\right). \tag{3}$$

#### *2.2. Hypothesis Testing*

In real word applications, researchers are interested in testing about the parameter γ. For example, the null hypothesis *H*0 : γ = 1 means that the CVs of two populations are equal. To perform the hypothesis test *H*0 : γ = γ0, the test statistic:

$$T\_0 = \sqrt{n} \left(\frac{\hat{\mathcal{Y}} - \mathcal{Y}\_0}{\lambda^\*}\right) \tag{4}$$

is generally applied, such that:

$$
\lambda^{\*2} = \frac{1}{\hat{C}V\_X^2} (\nu\_0^2 \delta\_X^2 + \delta\_Y^2).
$$

If the null hypothesis *H*0 : γ = γ0 is satisfied, then the asymptotic distribution of *T*0 is standard normal.

#### *2.3. Normal Populations*

Naturally, many phenomena follow normal distribution. This distribution is very important in natural and social sciences. Many researchers focused on the comparison between the CVs of two independent normal distributions. Nairy and Rao [12] reviewed and studied several methods such as likelihood ratio test, score test and Wald test that could be used to compare the CVs of two independent normal distributions. If the parent distributions *X* and *Y* are normal, then:

$$
\mu\_{\mathfrak{P}X} = \mu\_{\mathfrak{P}Y} = 0, \quad \mu\_{4X} = 3\mu\_{2X'}^2 \quad \mu\_{4Y} = 3\mu\_{2Y'}^2.
$$

Consequently, for normal distributions, δ2 *X* and δˆ2 *X* can be rewritten as:

$$
\delta\_X^2 = \frac{\mu\_{2X}\mu\_X^2 + 2\mu\_{2X}^2}{2\mu\_X^4},
$$

and:

$$
\hat{\delta}\_X^2 = \frac{m\_{2X}\overline{X}^2 + 2m\_{2X}^2}{2\overline{X}^4},
$$

respectively.

#### **3. Simulation Study**

In this section, the accuracy of the given theoretical results is studied and analyzed by different simulated datasets. For the populations *X* and *Y*, we respectively simulated different samples from symmetric distribution (normal) and asymmetric distributions (gamma and beta) with different CV values, (*CVX*,*CVY*) ∈ - (1, 1),(1, 2),(2, 3),(2, 5) , which are equivalent to γ ∈ {1, 2, 1.5, 2.5}. Figures 1–3 show the plots of probability density function (PDF) for the considered distributions.

**Figure 1.** Probability density function (PDF) of normal (μ, σ2) distribution with different coefficient of variation (CV) values. Black: μ = 1, σ = 1, CV = 1; red: μ = 1, σ = 2, CV = 2; green: μ = 1, σ = 3, CV = 3; blue: μ = 1, σ = 5, CV = 5.

**Figure 2.** PDF of gamma (α, λ) distribution with different CV values. Black: α = 1, λ = 0.001, CV = 1; red: α = 0.25, λ = 0.001, CV = 2; green: α = 0.11, λ = 0.001, CV = 3; blue: α = 0.04, λ = 0.001, CV = 5.

**Figure 3.** PDF of beta (α, β) distribution with different CV values. Black: α = 0.94, β = 30.39, CV = 1; red: α = 0.21, β = 6.87, CV = 2; green: α = 0.08, β = 2.51, CV = 3; blue: α = 0.009, β = 0.285, CV = 5.

The simulations are accomplished after 1000 runs using the *R 3.3.2* software (*R Development Core Team*, 2017) on a PC (Processor: Intel(R) CoreTM(2) Duo CPU T7100 @ 1.80GHz 1.80GHz, RAM: 2.00GB, System Type: 32-bit).

To check the accuracy of Equations (3) and (4), we estimated the coverage probability,

$$\text{CP} = \frac{\text{number of runs that Equation (3) contained true } \underline{\text{y}}}{1000}$$

for each parameter setting. We also computed the value of the test statistic in Equation (4), for each run. Then we considered the Shapiro–Wilk's normality test and the Q–Q plots to verify normality assumption for the proposed test statistic. Table 1 summarizes the CP values for different parameter settings.

**Table 1.** The CP values for different parameter settings.

As Table 1 indicates, the CP are very close to the considered level (1 − α = 0.95), especially when sample size was increased, and consequently the proposed method controlled the type I error. In other words, about 95% of simulated confidence intervals contained true γ and consequently it can be accepted that Equation (3) is asymptotically confidence interval for γ. The values of CPU times (in seconds) for different parameter settings given in Table 2, verify that this approach is not too time consuming. Furthermore, Figure 4 and Table 3 illustrate the Q-Q plots and the p-values of Shapiro–Wilk's test, respectively, to study the normality of the introduced test statistic.

First column: