Gibbs–Appell Equations in Finite Element Analysis of Mechanical Systems with Elements Having Micro-Structure and Voids

Abstract

In this paper, the authors propose the application of the Gibbs–Appell equations to obtain the equations of motion in the case of a mechanical system that has elements with a micro-polar structure, containing voids. Voids can appear as a result of the processing or manufacturing of the parts, or can be intentionally introduced. This research involves a model of the considered solid material containing voids. To determine the dynamic behavior of such a system, the Gibbs–Appell (GA) method is used to obtain the evolution equations, as an alternative to Lagrange’s classical description. The proposed method can be applied to any mechanical system consisting of materials with a micro-polar structure and voids. The study of such systems is interesting because the literature shows that even a reduce number of small voids can produce significant variations in physical behavior. The proposed method requires a smaller number of mathematical operations. To apply this method, the acceleration energy is calculated, which is then used to derive the equations. The method comes with advantages in the application to multibody systems having the mentioned properties and, in particular, in the study of robots and manipulators. Using the GA method, it is necessary to do a fewer differentiation operations than applying the Lagrange’s equations. This leads to a reduced amount of computation for obtaining the evolution equations.

1. Introduction

Recent years have led to an impressive development of the materials used in various fields of engineering. Such materials have appeared in aviation, the automotive industry, the naval industry, in the construction of robots and manipulators, in medicine, etc. However, these materials may contain voids accidentally created during manufacturing, in certain situations, may significantly alter the properties of the materials. Sometimes, these voids can be introduced intentionally or appear controlled, if a certain technological process is used. This has caused the interest of researchers in this type of materials [1].

Because to the difficulty in developing appropriate models, most studies are experimental. Tests have been carried out to ascertain how the mechanical properties change when voids appear [2,3,4]. Theoretical approaches are still at an unsatisfactory level, due to the numerous parameters that must be taken into account [5,6,7]. For this reason, it has been emphasized that there is a need to develop numerical procedures to solve these problems [8]. A first attempt to develop models using appropriate numerical methods is made in [9,10,11].

In the manufacturing process of some materials, especially composite materials, voids may appear, which are generally undesirable phenomena. These voids are generally filled with air. However, they may be filled with other gases, depending on the manufacturing process used. For manufacturers, the appearance of voids which are unexpected is an undesirable thing; a defect. These voids will have an influence on the physical properties of the respective material. Their existence can lead to the decrease in some important mechanical properties for the user and in the worst cases can cause the destruction of the material or the manufactured part. A void that appears can cause the initiation of cracks and ultimately reduce the life span.

It is also possible that other substances with a negative final influence on the material appear in these voids. The appearance of voids can also influence the isotropy of the material. In the study of composite and polymeric materials in general, voids negatively influence behavior. Therefore, it is necessary to develop suitable models with numerical solution capabilities for the study of polymeric materials used in engineering applications. There are many reasons for the appearance of voids and they relate to the manufacturing process, the wear of the machines, and the training of the personnel involved in the manufacture of the material. One example, the case of adhesive bonding, is presented in [12]. Another example, where voids can significantly influence material properties in the case of a polymer thin film, is presented in [13]. There are also a few cases in which the existence of voids in materials can be beneficial. Such examples are presented in specific applications [14,15].

Over time, as a result of the engineering requirements regarding the study of these materials, models that can predict their mechanical behavior have been developed and researched [16,17,18,19]. The cases solved so far are very simple. The development of complex models is necessary, and in this context the use of FEM could represent a main direction of analysis. The problem with voids is that, even in small concentrations, their existence will cause large variations in mechanical and physical properties [20,21,22,23,24,25,26]. The first studies in the field by Cowin and Goodman [27] resulted in the creation of the granular theory. Applications immediately appeared in various fields, first of all in the research of soils and rocks, and then of human bones [28].

In general, voids appear in materials that also exhibit microstructure (MS) properties, especially polymeric composite materials [29,30,31,32]. Thus, a theoretical support has been provided for the analysis of these materials. Unfortunately, a theoretical support is not enough for the analysis of practical engineering cases. This difference is mainly determined by the complexity of the equations that arise. Approximation methods and numerical methods are needed to solve some real-world applications. In this context, the FEM proves to be the main tool that makes the analysis of these structures possible [33].

The wealth of developed finite elements and well-verified procedures within FEMs make this method one that can successfully address such materials [34,35,36,37,38]. For composite, bi- or multiphase materials, useful results are presented in [39,40,41,42]. In this paper, the FEM is used for the dynamic analysis of such bodies. Basic notions for the application of FEM to hollow micropolar bodies are proposed in [43,44], and some developments are made in [45,46].

This paper aims to provide a model for determining the elastic behavior of a solid body with voids and micro-structure, using a formalism provided by Analytical Mechanics. For this, the FEM is used together with an analytical method to establish the equations of motion. The main difficulty is the complexity of the analytical formulas. To facilitate the analysis, the authors propose the use of another method (as an alternative to Lagrange’s Equations (LE), the Gibbs–Appell (GA) equations, introducing the rarely used notion of acceleration energy.

In most engineering applications, the analysis method for writing the equations of motion in the FEM is represented by Lagrange’s formalism. After performing this operation, using classic assembly methods, applying the boundary conditions, and introducing forces and moments corresponding to the studied problem, the system of equations that describe the movement of the system is finally obtained. The choice of a suitable finite element for modeling the studied mechanical system will determine the shape functions used; based on these, after the integrations, the coefficients of the motion equations will be obtained. The hypothesis of small deformations is considered to be valid. A number of papers have addressed the subject of the motion of systems with elastic elements using one-dimensional finite elements [47,48,49,50,51,52,53], and also more sophisticated finite elements [54,55,56].

LEs are usually preferred to determine the evolution equations for an element [57,58,59,60]. Recent works have begun to use and reconsider other methods of Analytical Mechanics in the analysis of such systems. GA equations will be the authors’ choice as an alternative solution to LE to study the mechanics of a material with micro-structure and voids, having a general motion. Using this method has several advantages. So, the number of operations used to determine the system of differential equations (and, consequently, the time and cost of designing using these materials) can decrease.

The GA equations represent a method for solving non-holonomic systems that presents some advantages over the other methods offered by Analytical Mechanics. The method was proposed in 1879 [61] by Gibbs and rediscovered by Appell in 1899. The notion introduced for the application of the method is that of the energy of accelerations. The method essentially is an application of Gauss’s principle of minimum constraint. The obtained results are equivalent to the results obtained by applying other equivalent formalisms provided by Analytical Mechanics. The use of these equations in engineering applications has begun to be reconsidered in recent years, offering advantages in terms of the number of calculations involved and reducing the cost of the modeling stage. Extension to systems with micro-structures and voids is a natural development of the method in the current context of the existing computational tools.

The use of robots has become a general phenomenon in all industries. It is a trend of the last few decades which requires an adaptation of the design of these systems to ensure adequate control. The latest developments have led to the design and manufacture of machines that work at high speeds and are subject to high forces. They need to be carefully designed, which is why many studies have begun to analyze various aspects of their operation. For such systems, the accelerations become important and can greatly influence operation. In this context, the use of the GA method becomes important, introducing the notion of acceleration energy. At this point, the method starts to become interesting for various engineering problems [62,63].

The method has proven to be interesting for a wide range of applications. An example of a flexible robot is presented in [64]. The theoretical results obtained by the GA method have been validated experimentally. Another application [65] concerns an N-flexible-link manipulator used in industrial applications. In [66], using the GA formalism, the evolution equations are established and, separately, an algebraic system provides the forces acting on the links [67,68].

An advantage of applying this set of equations is the reduced number of operations compared to other methods. Currently, the method is beginning to be reconsidered and several researchers have used it in their works [69,70,71,72,73,74]. The FEM has demonstrated its validity and utility in numerous problems in a very wide spectrum of domains. It is a powerful tool for studying the behavior of mechanical systems. For multibody systems, FEM has been used along with Lagrange’s equations. FEM currently has computational procedures and models incorporated into available software that ensure a very wide range of applications are handled.

Many finite elements have been developed, along with powerful pre- and post-processing tools. The basic assumption for the application of this method remains the hypothesis of small deformations. In this paper, we address the problem of using the FEM together with an alternative Analytical Mechanics method to obtain the motion equations, namely, the GA method. For systems of such complexity, the use of micro-structured materials with voids requires the use of procedures in the analysis that require as few calculations as possible. We mention that, in the last twenty years, the classical methods of Analytical Mechanics have begun to be reconsidered. This is because engineering applications have started to become more and more complex. In this context, in which the developed models are very sophisticated and involve a huge number of arithmetic calculations, some advantages offered by the particular form of some of the mechanical principles can have a major influence on the cost and the time required to achieve a high-performance design. The work is part of this trend of revaluing classical knowledge and reevaluating the advantages it can offer.

The accuracy issues of the obtained results are not presented in this paper. They are a field of study in themselves. In this research, we were mainly concerned with presenting an equivalent method for determining the set of equations. The justification for this is the opportunity to obtain some computational advantages, namely, a shorter modeling time and, consequently, lower costs for obtaining numerical results. The use of this method must also ensure relevant results in order to make sense from a technical point of view.

2. Energy of Accelerations and Basic Notions

The system is considered a mechanical system made up of N material points whose position is defined by n independent coordinates. These coordinates will be noted with $q_i$, where $i=\overline{1,n}$. If the studied system is a system with scleronoma links, the acceleration of one of the material particles, let us say i, of the system is expressed by the relation

$${\overline{a}}_i={\displaystyle \sum _k{\displaystyle \sum _j\frac{{\partial }^2{\overline{r}}_i}{\partial q_k\partial q_j}{\dot{q}}_k{\dot{q}}_j}}+{\displaystyle \sum _k\frac{\partial {\overline{r}}_i}{\partial q_k}{\ddot{q}}_k} ; k,j=\overline{1,n} ; i=\overline{1,N}.$$

(1)

In this expression, ${\overline{r}}_i$ is the position vector of the point i. In the following, this will be denoted with a dot placed over a scalar, vector, or matrix size of the derivative of the considered size with respect to time. The energy of acceleration denoted by S is a basic notion used in the further considerations [72]:

$$S=\frac{1}{2}{\displaystyle \sum _i{m}_i{a}_i^2} ; i=\overline{1,N}.$$

(2)

Formal Equation (2) is similar to the relation used to define the kinetic energy. However, the mechanical significance is obviously different. Considering a solid elastic body, Equation (2) for a system of material points can be extended for the solid under the for:

$$S=\frac{1}{2}{\displaystyle {\int }_V{a}^{2}dm}.$$

(3)

If we consider the acceleration of a point arbitrary chosen, this has the expression

$$\overline{a}={\displaystyle \sum _k{\displaystyle \sum _j\frac{{\partial }^2\overline{r}}{\partial q_k\partial q_j}{\dot{q}}_k{\dot{q}}_j}}+{\displaystyle \sum _k\frac{\partial \overline{r}}{\partial q_k}{\ddot{q}}_k} ; k,j=\overline{1,n}.$$

(4)

Using Equations (3) and (4) we get

$$\begin{gathered} S=\frac{1}{2}{\displaystyle \underset{V}{\int }\rho a^{2}dV}=\frac{1}{2}{{\displaystyle \underset{V}{\int }\rho \left({\displaystyle \sum _k{\displaystyle \sum _j\frac{{\partial }^2\overline{r}}{\partial q_k\partial q_j}{\dot{q}}_k{\dot{q}}_j}}+{\displaystyle \sum _k\frac{\partial \overline{r}}{\partial q_k}{\ddot{q}}_k}\right)}}^{2}dV= \\ {E}_{ao}(\dot{q})+E_{a1}(\dot{q},\ddot{q})+E_{a2}(\ddot{q}) , k,j=\overline{1,n}. \end{gathered}$$

(5)

where it is denoted by

$$E_{ao}(\dot{q})=\frac{1}{2}{\displaystyle \sum _k{\displaystyle \sum _j{\displaystyle \sum _l{\displaystyle \sum _m\left({\displaystyle \underset{V}{\int }\rho \frac{{\partial }^2\overline{r}}{\partial q_k\partial q_j}\frac{{\partial }^2\overline{r}}{\partial q_l\partial q_m}dV}\right)}}}}{\dot{q}}_k{\dot{q}}_j{\dot{q}}_l{\dot{q}}_m ; j,k,l,m=\overline{1,n},$$

(6)

which represents the part of the energy of accelerations that do not contains accelerations; by

$$E_{a1}(\dot{q},\ddot{q})=\frac{1}{2}{\displaystyle \sum _k{\displaystyle \sum _j{\displaystyle \sum _l\left({\displaystyle \underset{V}{\int }\rho \frac{{\partial }^2\overline{r}}{\partial q_k\partial q_j}\frac{\partial \overline{r}}{\partial q_l}dV}\right)}}}{\dot{q}}_k{\dot{q}}_j{\ddot{q}}_l ; k,j,l=\overline{1,n}.$$

(7)

which represents the part of the energy of accelerations that contains the linear terms in accelerations; and by

$$E_{a2}(\ddot{q})=\frac{1}{2}{\displaystyle \sum _k{\displaystyle \sum _j\left({\displaystyle \underset{V}{\int }\rho \frac{\partial \overline{r}}{\partial q_k}\frac{\partial \overline{r}}{\partial q_j}dV}\right)}}{\ddot{q}}_k{\ddot{q}}_j ; k,j=\overline{1,n}.$$

(8)

which represents the part that contains the quadratic terms in accelerations.

The classic forms of the GA equations are:

$$\frac{\partial S}{\partial {\ddot{q}}_i}=Q_i i=\overline{1,n}.$$

(9)

In Equation (9), $Q_i$ represents the generalized force corresponding to the generalized coordinate $q_i$.

Now, we consider an elastic body that can present, in its volume, a series of voids. The density $\rho (x,y,z)$ of the body changes in every point, and can be expressed by the relation [34]

$$\rho (x,y,z)=v(x,y,z) \gamma (x,y,z).$$

(10)

Here, $\gamma (x,y,z)$ is the mass density of the solid if the voids do not exist, and $v(x,y,z)$ represents the volume fraction of the existing material, ($0<\nu \le 1$). The volume fraction represents the fraction of the volume of voids over the total volume. We can consider the situation at an initial moment:

$${\rho }_o(x,y,z)={\nu }_o(x,y,z) {\gamma }_o(x,y,z).$$

(11)

Let us consider that the moving reference system, related to the studied finite element, has the angular velocity $\overline{\omega }$, the angular acceleration $\overline{\epsilon }$, the velocity ${\overline{v}}_o$ of the origin $O(X_o,Y_o,Z_o)$, and ${\overline{a}}_o$ the acceleration of the point O. In the theory of the elastic body with voids, the position of an elementary element is described by the coordinate of the point around which is considered the elementary element and the third rotation of the element. The approximation made in FEM links the displacements of a certain point denoted M with the independent coordinates ${\left\{\delta \right\}}_L$, representing the displacement of the nodes of the finite elements, and ${\left\{\varphi \right\}}_L$, representing the rotation of a reference frame situated in the point. The local kinematics of the finite element are therefore defined through the relation [34]

$${\left\{u\right\}}_L=\left[N_{\delta }\right]{\left\{\delta \right\}}_L ; {\left\{\varphi \right\}}_L=\left[N_{\varphi }\right]{\left\{\delta \right\}}_L.$$

(12)

and the distribution of the percentage of voids can be written as

$${\nu }_L=\left[N_{\nu }\right]{\left\{\delta \right\}}_L.$$

(13)

In these relations, $\left[N_{\delta }\right]$, $\left[N_{\varphi }\right]$, and $\left[N_{\nu }\right]$ are shape function matrices. The quantity ${\nu }_L$ is scalar, and the shape matrix $\left[N_{\nu }\right]$ is a line matrix.

If ${\left\{u\right\}}_L$ represents the displacement of the point M that moves into M′, it can be written

$${\left\{r_{M^{\prime }}\right\}}_G={\left\{r_O\right\}}_G+\left[R\right]\left({\left\{r\right\}}_L+{\left\{u\right\}}_L\right) .$$

(14)

In Equation (14), $\left[R\right]$ represents the rotation matrix which realizes the transformation of a vector from the local coordinate system L to the global coordinate system G. Additionally, ${\left\{r_{M^{\prime }}\right\}}_G$ is the position vector of M′ having the components expressed in a global reference frame. The velocity of M′ is

$${\left\{v_{M^{\prime }}\right\}}_G={\left\{{\dot{r}}_{M^{\prime }}\right\}}_G= {\left\{{\dot{r}}_O\right\}}_G+\left[\dot{R}\right]\left({\left\{r\right\}}_L+\left[N\right]{\left\{\delta \right\}}_L\right)+ \left[R\right]\left[N\right]{\left\{\dot{\delta }\right\}}_L.$$

(15)

Using Equations (12) and (13), we obtain

$${\left\{\dot{\varphi }\right\}}_L=\left[N_{\varphi }\right]{\left\{\dot{\delta }\right\}}_L,$$

(16)

$${\dot{\nu }}_L=\left[N_{\nu }\right]{\left\{\dot{\delta }\right\}}_L.$$

(17)

With Equations (15)–(17), the kinetic energy of a finite element, considered as a part of a micro-structured material with voids, is:

$$E_c=\frac{1}{2}{\displaystyle \underset{V}{\int }{\rho }_0{\left\{v_{M^{\prime }}\right\}}_G^T{\left\{v_{M^{\prime }}\right\}}_G}dV+\frac{1}{2}{\displaystyle \underset{V}{\int }{\rho }_0{\left\{{\dot{\varphi }}_{M^{\prime }}\right\}}_G^T\left[Y\right]{\left\{{\dot{\varphi }}_{M^{\prime }}\right\}}_G}dV+\frac{1}{2}{\displaystyle \underset{V}{\int }{\rho }_0{\kappa }_{\nu }{\dot{\nu }}_{M^{\prime }}^2{}_G}dV.$$

(18)

It is denoted with ${\kappa }_{\nu }$ the so called “equilibrated inertia” [34], and with $\left[Y\right]$, the matrix of the coefficients of inertia [1,2,3,5,10,12].

The classical form of the elastic potential energy is [58]

$$E_p=\frac{1}{2}{\displaystyle {\int }_V{\left\{\sigma \right\}}^T\left\{\epsilon \right\}}dV.$$

(19)

The classical Hooke’s law is [54]

$$\left\{\sigma \right\}=\left[H\right]\left\{\epsilon \right\}.$$

(20)

and it exists the following relations:

$$\left\{\epsilon \right\}=\left[b_1\right]\left\{u\right\}+\left[b_2\right]\left\{\varphi \right\}=\left(\left[b_1\right]\left[N_{\delta }\right]+\left[b_2\right]\left[N_{\varphi }\right]\right){\left\{\delta \right\}}_L=\left[N_s\right]{\left\{\delta \right\}}_L ,$$

(21)

where $\left[b_1\right]$ and $\left[b_2\right]$ are two differential operators, and

$$\left[N_s\right]=\left[b_1\right]\left[N_{\delta }\right]+\left[b_2\right]\left[N_{\varphi }\right]$$

(22)

(see [54]).

So, the potential energy can be written

$$E_p=\frac{1}{2}{\left\{\delta \right\}}_L^T\left({\displaystyle {\int }_V {\left[N_s\right]}^T{\left[H\right]}^T\left[N_s\right]dV}\right){\left\{\delta \right\}}_L.$$

(23)

Equation (11) can be expressed as

$$E_p=\frac{1}{2}{\left\{\delta \right\}}_L^T\left[k\right]{\left\{\delta \right\}}_L.$$

(24)

In Equation (24), $\left[k\right]$ represents the stiffness matrix.

The term of energy due to the voids is

$$E_{pv}=\frac{1}{2} {\left\{\delta \right\}}_L^T\left[k_v\right]{\left\{\delta \right\}}_L.$$

(25)

In reference [53], $E_{pv}$ is represented as

$$E_{pv}=\frac{1}{2} \left[\begin{array}{ccc}\frac{\partial \nu }{\partial x}& \frac{\partial \nu }{\partial y}& \frac{\partial \nu }{\partial z}\end{array}\right]\left[\begin{array}{ccc}{d}_{11}& d_{12}& d_{13}\\ d_{21}& d_{22}& d_{23}\\ d_{31}& d_{32}& d_{33}\end{array}\right]\left\{\begin{array}{c}\frac{\partial \nu }{\partial x}\\ \frac{\partial \nu }{\partial y}\\ \frac{\partial \nu }{\partial z}\end{array}\right\}.$$

(26)

Using Equations (12) and (13), it gives

$$\begin{array}{l}\frac{\partial \nu }{\partial x}=\frac{\partial }{\partial x}\left(\left[N_{\nu }\right]{\left\{\delta \right\}}_L\right)=\left[N_{\nu ,x}\right]{\left\{\delta \right\}}_L , \\ \frac{\partial \nu }{\partial y}=\frac{\partial }{\partial y}\left(\left[N_{\nu }\right]{\left\{\delta \right\}}_L\right)=\left[N_{\nu ,y}\right]{\left\{\delta \right\}}_L , \\ \frac{\partial \nu }{\partial z}=\frac{\partial }{\partial z}\left(\left[N_{\nu }\right]{\left\{\delta \right\}}_L\right)={\left[N_{\nu ,z}\right]}_L{\left\{\delta \right\}}_L ,\end{array}$$

(27)

or

$$\left\{\begin{array}{c}\frac{\partial \nu }{\partial x}\\ \frac{\partial \nu }{\partial y}\\ \frac{\partial \nu }{\partial z}\end{array}\right\}=\left[\begin{array}{c}\left[N_{\nu ,x}\right]\\ \left[N_{\nu ,y}\right]\\ \left[N_{\nu ,z}\right]\end{array}\right]{\left\{\delta \right\}}_L ,$$

(28)

or:

$$\left[\begin{array}{ccc}\frac{\partial \nu }{\partial x}& \frac{\partial \nu }{\partial y}& \frac{\partial \nu }{\partial z}\end{array}\right]={\left\{\delta \right\}}_L^T\left[{\begin{array}{ccc}{\left[N_{\nu ,x}\right]}^T& {\left[N_{\nu ,y}\right]}^T& \left[N_{\nu ,z}\right]\end{array}}^T\right].$$

(29)

and

$$E_{pv}=\frac{1}{2}{\left\{\delta \right\}}_L^T\left[{\begin{array}{ccc}{\left[N_{\nu ,x}\right]}^T& {\left[N_{\nu ,y}\right]}^T& \left[N_{\nu ,z}\right]\end{array}}^T\right]\left[d\right]\left[\begin{array}{c}\left[N_{\nu ,x}\right]\\ \left[N_{\nu ,y}\right]\\ \left[N_{\nu ,z}\right]\end{array}\right]{\left\{\delta \right\}}_L.$$

(30)

Considering a homogeneous and isotropic body, the coefficients $d_{ij}$ are equal, and an estimate for them is presented in [34]. The matrix [d] is the unit matrix premultiplied with a coefficient for the most practical cases [34].

The work produced by concentrated load ${\left\{q_t\right\}}_L$ and moments ${\left\{q_r\right\}}_L$ is

$$W^c={\left\{q\right\}}_L^T{\left\{\delta \right\}}_L.$$

(31)

where it is denoted

$${\left\{q\right\}}_L^T={\left\{q_t\right\}}_L^T+{\left\{q_r\right\}}_L^T.$$

(32)

and that produced by the volume forces and moments is

$$W^d={\left\{q^{*}\right\}}_L^T{\left\{\delta \right\}}_L.$$

(33)

It is noted as

$${\left\{q^{*}\right\}}_L^T=\left[{\displaystyle \underset{V}{\int }{\left\{p\right\}}_L^T\left[N_{\delta }\right]dV}\right]{\left\{\delta \right\}}_L+\left[{\displaystyle \underset{V}{\int }{\left\{m\right\}}_L^T\left[N_{\varphi }\right]dV}\right]{\left\{\delta \right\}}_L.$$

(34)

The potential energy due to the existence of voids is

$$E_{pv}=\frac{1}{2} {\left\{\delta \right\}}_L^T\left[k_v\right]{\left\{\delta \right\}}_L.$$

(35)

In what follows, the method of LE is presented to determine the evolution equations [75].

The Lagrangean [21] has the classical form

$$L=E_c-E_p+W^d+W^c+W^v.$$

(36)

In this formula, $W^d$ is the work of distributed forces, $W^c$ is the work of concentrated forces, and $W^v$ is the work determined by the existence of voids.

We use the notations

$$\left[m\right]={\displaystyle \underset{V}{\int }{\left[N_{\delta }\right]}^T\left[N_{\delta }\right]{\rho }_{o}dV}; \left[J\right]={\displaystyle \underset{V}{\int }{\left[N_{\varphi }\right]}^T\left[Y\right]\left[N_{\varphi }\right]{\rho }_{o}dV}; \left[m_{\nu }\right]={\displaystyle \underset{V}{\int }{\left[N_{\nu }\right]}^T\left[N_{\nu }\right]\kappa {\rho }_{o}dV};$$

(37)

$$\left[c\right]=2{\displaystyle \underset{V}{\int }{\left[N_{\delta }\right]}^T{\left[\dot{R}\right]}^T\left[R\right]\left[N_{\delta }\right]{\rho }_{o}dV}; \left[k(\epsilon )+k({\omega }^2)\right]={\displaystyle \underset{V}{\int }{\left[N_{\delta }\right]}^T\left({\left[\ddot{R}\right]}^T\left[R\right]\right)\left[N_{\delta }\right]{\rho }_{o}dV}.$$

(38)

$$\left[m_O^i\right]={\displaystyle \underset{V}{\int }{\left[N_{\delta }\right]}^T{\rho }_{o}dV}; \left\{q^i(\epsilon )\right\}+\left\{q^i({\omega }^2)\right\}={\displaystyle \underset{V}{\int }{\left[N_{\delta }\right]}^T{\left[R\right]}^T\left[\ddot{R}\right]{\left\{r\right\}}_L{\rho }_{o}dV}.$$

(39)

If we use the LE, it results in

$$\begin{gathered} \left(\left[m\right]+\left[J\right]+\left[m_{\nu }\right]\right){\left\{\ddot{\delta }\right\}}_L+\left[c\right]{\left\{\dot{\delta }\right\}}_L+\left(\left[k\right]+\left[k_{\nu }\right]+\left[k(\epsilon )\right]+\left[k({\omega }^2)\right]\right){\left\{\delta \right\}}_L \\ ={\left\{q\right\}}_L^T+{\left\{q^{*}\right\}}_L^T+{\left\{q_v\right\}}^T-\left[m_O^i\right]{\left\{{\ddot{r}}_O\right\}}_L-\left\{q^i(\epsilon )\right\}-\left\{q^i({\omega }^2)\right\}. \end{gathered}$$

(40)

Considering the static case, it obtains

$$\left(\left[k\right]+\left[k_{\nu }\right]+\left[k(\epsilon )\right]+\left[k({\omega }^2)\right]\right){\left\{\delta \right\}}_L={\left\{q\right\}}_L^T+{\left\{q^{*}\right\}}_L^T+{\left\{q_v\right\}}^T.$$

(41)

3. Results

In the following, the GA method is presented to obtain the evolution equations for a solid composed of such a material. Using the GA equations, it is necessary to write the energy of the accelerations. The acceleration of a point of the solid is

$${\left\{a_{M^{\prime }}\right\}}_G={\left\{{\ddot{r}}_O\right\}}_G+\left[\ddot{R}\right]{\left\{r\right\}}_L+\left[\ddot{R}\right]\left[N\right]{\left\{\delta \right\}}_L+2\left[\dot{R}\right]\left[N\right]{\left\{\dot{\delta }\right\}}_L+ \left[R\right]\left[N\right]{\left\{\ddot{\delta }\right\}}_L.$$

(42)

and the energy of the accelerations is given by the relation

$$E_c=\frac{1}{2}{\displaystyle \underset{V}{\int }{\rho }_0{\left\{a_{M^{\prime }}\right\}}_G^T{\left\{a_{M^{\prime }}\right\}}_G}dV+\frac{1}{2}{\displaystyle \underset{V}{\int }{\rho }_0{\left\{{\ddot{\varphi }}_{M^{\prime }}\right\}}_G^T\left[Y\right]{\left\{{\ddot{\varphi }}_{M^{\prime }}\right\}}_G}dV+\frac{1}{2}{\displaystyle \underset{V}{\int }{\rho }_0{\kappa }_{\nu }{\ddot{\nu }}_{M^{\prime }}^2{}_G}dV.$$

(43)

The second order derivative of Equation (12) are used in Equation (41):

$${\left\{\ddot{\varphi }\right\}}_L=\left[N_{\varphi }\right]{\left\{\ddot{\delta }\right\}}_L.$$

(44)

Performing all the necessary calculations within the expression (43) using the notations (37)–(39) and applying Equation (9), we finally obtain Relations (40) and (41). The number of operations used to obtain them is less than in the case of using LE.

An example will illustrate our theoretical consideration. We consider a circular bar in which gaps of material will be inserted at percentages of 1, 2, 3, and 4% (Figure 1). The bar has 1 m length and a diameter equal to 1 cm. To achieve this, consider the discretized bar in three-dimensional finite elements. In a cross section of the bar, with a circular shape, 32 finite elements are considered (Figure 2). Statistically, one element is extracted from this structure to assure a percentage of 1, 2, 3, or 4%. The extraction of these elements to create the voids is done statistically. For example, if the structure is divided into 3200 elements, 32 elements will be extracted to achieve 1% voids. Thus, with this procedure, a bar with 1% voids is first made, then the number of voids is successively increased until the last variant containing 4% voids is obtained. For each of these variants, the strain and stress calculation is performed to obtain the maximum values. The studied material is considered to have the Young’s modulus of 125 GPa, with $\nu =0.3$. The force value F is 100 N. The bar is clamped at both ends.

Considering the distribution of voids in Figure 2, a calculation of stresses and deformations was performed in the four considered cases and the results were compared with the situation of a homogeneous bar without voids. Figure 3 contains the maximum displacements of the bar in the studied cases. Figure 4 contains the maximum von Mises stresses of the bar in the studied cases.

In

Table 1. The displacements and von Mises stresses.

	Homogeneous Bar	1% Voids	2% Voids	3% Voids	4% Voids
Displacements [mm]	0.5126	0.5322	0.5524	0.5586	0.5723
Variation of displacements [%]	----	3.82	7.76	8.97	11.64
von Mises stress [GPa]	22.995	23.541	24.213	24.805	25.323
Variation of von Mises stress [%]	----	2.37	5.39	7.86	10.12

, these results are synthetized for a global image of the displacements and stresses.

A calculation of the eigenfrequencies of the system was made in two cases: for 0% voids and for 4% voids. The results are presented in

Table 2. The eigenmodes of the bar with voids.

Mode	Representation	Frequency for 0% Voids [Hz]	Frequency for 4% Voids [Hz]
1		139.0319	139.9093
2		139.0319	139.9093
3		382.3166	382.3166
4		382.3166	382.3166
5		747.0574	747.0574
6		747.0574	747.0574
7		1229.853	1229.853
8		1229.853	1229.853
9		1828.136	1828.136
10		1828.136	1828.136

. It can be concluded that the existence of a percentage of voids in the material will have an insignificant influence on the vibrations of the studied body.

4. Discussion

If LE are applied to determine the set of the equations of motion, there will be fewer mathematical differentiation operations if the GA equations are applied. The use of these equations makes the calculations easier, which for large systems can be significant. It is difficult to estimate the number of operations, this depends on the size of the system and the type of finite elements that was chosen, which determines the number of degrees of freedom that will be taken into account.

The LE is

$$\frac{d}{dt}\left\{\frac{\partial L}{\partial {\dot{\delta }}_e}\right\}-\left\{\frac{\partial L}{\partial {\delta }_e}\right\}=0$$

(45)

From a qualitative point of view, it can be observed that it is necessary to perform four differentiation operations for the expressions $\left\{\partial L/\partial {\dot{\delta }}_e\right\}$, another four for the expressions $\left\{\partial L/\partial {\delta }_e\right\}$, and then the derivatives in relation to time of the terms $\left\{\partial L/\partial {\dot{\delta }}_e\right\}$.

However, if we use the Gibbs–Appell equations, it is necessary to differentiate the energy of the accelerations. In that case, it is necessary to make only five operations. These operations in both cases are related to the components of a vector; for this reason, they will also depend on the number of degrees of the system. In the case of the Gibbs–Appell formalism, only five differentiations will be made.

It can therefore be concluded that if the Gibbs–Appell equations and the acceleration energy are applied, the number of operations that must be performed is less than if the Lagrange equations and the kinetic energy are used. In most applications, Lagrange’s method is used. This is due to the fact that most researchers are very familiar with the notion of kinetic energy and less so with the notion of acceleration energy. This method can therefore also be used when studying a mechanical system that has in its composition elements made of materials with micro-structure and voids. The use of the Gibbs-Appell formalism has the advantage of reducing the calculation time and, consequently, the cost of the analysis stage of such a mechanical system. In general, it can be considered that the number of operations can be halved. In the case of complex structures, of large dimensions and operating in difficult conditions, this reduction can be important.

Previous studies have shown that the existence of voids, even in small percentages, can cause a significant variation in the mechanical properties [75,76] (see

Table 3. The mechanical effect of the existence of voids in structure.

Increase in Percentage of Voids	Reduction in Flexural Rigidity/[Reference]
1%	10%/[21]
0% to 5%	25%/[21]
2.5%	20%/[20]
1%	10%/[20]
3%	17%/[77]
1%	35%/[78]
6%	40%/[78]

).

The results obtained in our analysis are presented in

Table 4. The results obtained for the studied model.

Increase in Percentage of Voids	Variation of Displacement [%]	Variation of Von Mises Stress [%]
0%	0	0
1%	3.82	2.37
2%	7.76	5.39
3%	8.97	7.86
4%	11.64	10.12

.

From this table, it can be seen that voids can lead to a significant reduction in the mechanical properties. It follows that the study of their influence in different situations has become a very important thing. So, the study of these materials becomes important for practical applications.

5. Conclusions

In this paper, the authors sought to extend the application of the Gibbs–Appell formalism to analyze systems consisting of materials with micro-s and voids. In a previous paper, the authors extended the area of FEM to apply it to MBS bodies with microstructure and voids. The method used to do this was the Lagrange Equations. As for other mechanical systems, usually made by rigid bodies, there are a lot of papers proving that using other methods such as the Lagrange formalism can be useful in some cases. We extended this observation to MBS with microstructure and voids. We found that this is possible with little effort. As a result, it is also possible for this type of material to apply the Gibbs–Appell equations and consequently benefit from the advantages of using these equations. The developed model takes into account both the microstructural properties of the material and the embedded voids.

Advantages of using the Gibbs–Appell method are mentioned in a series of papers [61,63,65,67]. The present work aims to expand this area of application of the method. If most of the works that approach this formalism refer to bodies with rigid elements, the authors tried to extend the method to elastic mechanical systems. The models used for these systems lead to a large number of operations and, as a result, any method to shorten the calculation time (and costs) is extremely beneficial. Obviously, at this stage, the work can form the basis of further developments for the study of mechanical engineering systems used in practice.