# Methodology of combined application of directional derivatives and the extended finite element method (X-FEM) for solving vibrations eigenvalue problems

Daria Serbichenko, Pascal Cosson, Patrick Rozycki

## ► To cite this version:

Daria Serbichenko, Pascal Cosson, Patrick Rozycki. Methodology of combined application of directional derivatives and the extended finite element method (X-FEM) for solving vibrations eigenvalue problems. European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2016), Jun 2016, Chersónissos, Greece. pp.6483. hal-01349713

HAL Id: hal-01349713

https://hal.science/hal-01349713

Submitted on 28 Jul 2016

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.
METHODOLOGY OF COMBINED APPLICATIONS OF
DIRECTIONAL DERIVATIVES AND EXTENDED FINITE ELEMENT
METHOD (X-FEM) FOR SOLVING VIBRATION EIGENVALUE
PROBLEMS

D. Serbichenko¹, P. Cosson¹ and P. Rozycki¹

¹ Institut de recherche en Génie Civil et Mécanique (GeM) - UMR CNRS 6183
Ecole Centrale de Nantes, 1 rue de la Noë, BP 92101, 44321 Nantes Cedex 3 (France)
e-mail: daria.serbichenko@ec-nantes.fr, {pascal.cosson,patrick.rozycki}@ec-nantes.fr

**Keywords:** Directional Derivatives; X-FEM; Vibrations; Eigenvalue problem.

**Abstract.** This paper presents a new methodology for solving the eigenvalue problem for time dependent structures. The time dependent structures of interest are structures with a moving discontinuity such as crack or structures with moving free (external/internal) surfaces. For the last case, they can result from a removal of material during a machining process or from a deterioration of the structure's geometry. The methodology that we developed, is based on a combination of the eXtended Finite Element Method (X-FEM) and the Directional Derivatives method. X-FEM enables to overcome the drawbacks of conformity and remeshing: indeed, using standard FEM, a moving discontinuity in time within a structure requires not only that the mesh must conform to the discontinuity geometry but also to fully remesh the structure as much as necessary to follow the discontinuity in time. In order to alleviate this last point, the directional derivatives are a powerful tool because they allow to estimate the evolution of quantities from on reference domain to another one. In our case, they will allow to estimate the solutions of the eigenvalue problem. We suggest on the first sections to remind the main keys of both methods and we present then the combined methods in order to solve an eigenvalue problem. The application will be done on a one-dimensional eigenvalue problem and the numerical results will be presented to demonstrate the accuracy and the advantages of selected approaches. We conclude on the future prospects of the current work that mainly consist of to develop the methodology at the second order in order to increase the accuracy and to find a criteria in order to automatize the combined methods.
# 1 INTRODUCTION

The eigenvalue problems are of a great importance and play a crucial role in many engineering fields. In the aeronautical field for instance, the identification of natural frequencies and natural shapes during the design of a plane is compulsory before its commercialization to ensure the security of passengers. A modal analysis is usually done for the whole plane as well as for some structural parts such as engines or blades. In the field of machining, we can also cite the problem of quality insurance for the dimensions of machined pieces: the reliability of the processing machines must be ensured in order to avoid some critical interactions of their possible vibrations modes on the constraints of the dimensions of the piece. At last, in some other cases, the design of optimized structures to prevent failure due to instabilities and vibrations brings the problem of the determination of optimal physical parameters so that the load carrying capacity or the fundamental natural frequency is maximized.

Nowadays, there is a need for solving eigenvalue problems for time dependent structures. Changes of a structure can occur due to the propagation of cracks. They can also appear due to the fact that the shapes of a structure are modified because of the removal of material during its machining process or because of the deterioration of its geometry. They can also happen during design process when the structure has to be optimized with respect to frequency criteria. To investigate such issues, standard ways involve remeshing the structure in order to follow its changes and to conform to its current geometry. They also involve solving the eigenvalue problem for each configuration of the structure which has been meshed. Besides, we can also notice a possible loss of accuracy as the data is mapped from the old mesh to the new one. In other words, the FEM applications can lead to lot of limitations and complexities for engineers and can really be time consuming.

To offer an alternative solution to the above issues, we propose in this paper a new methodology combining two approaches. The first one is the Extended Finite Element Method (X-FEM) [1-2] which is based on the Partition of Unity [3]. The aim of this method is to alleviate the mesh constraints that come from discontinuities, using the level sets [4-7]. That is the main reason that X-FEM is effectively applied to solve many problems in material modeling, crack propagation problems [8] and structures with holes (or free surface) [9,10]. For our problem, the X-FEM method solves the mesh dependency of the structure due to the changing of boundaries or/and due to the discontinuities.

Since the simulation in order to follow the discontinuity in time and its consequence on the natural frequencies imposes to continuously compute the eigenvalue problem, it seems that this method is not quite optimal in terms of CPU time cost. In order to avoid computing the eigenvalue problem on each configuration where information is needed, we focus on a second approach: the Directional Derivatives [11-13]. This mathematical tool allows us to get the evolution of a quantity from one domain to another one. Applied to our problem, it becomes possible to estimate natural frequencies and natural shapes for a given configuration without solving the eigenvalue problem for this configuration.

In the next sections of this paper, we generalize joint application of the above techniques in one unique approach; we first introduce the state of eigenvalue problem, expose some key points of the X-FEM method and the basic definitions for the Directional Derivatives of first order in the general case. We then suggest applying the joint techniques (X-FEM and Directional Derivatives) for solving the eigenvalue problem in the one-dimensional case. In a last section, we discuss about the obtained results for the 1D case and we present the works and the main key-points for the application of the methodology in 2D case.
## 2 STATE OF THE PROBLEM

<figure><img src="image_1.png" /><figcaption>Figure 1: The eigenvalue problem associated with free oscillations.</figcaption></figure>

The free-vibration natural frequencies and mode shapes of a linear structural system can be computed by solving the following eigenvalue problem. With the assumption of small perturbations in the framework of the linear elasticity theory, the free oscillations of a structure (Figure 1) are governed by the following equations system:

* local equations

$$ div \bar{\bar{\sigma}}(u) = \rho \ddot{u} \quad (1) $$

* behavior law

$$ \bar{\bar{\sigma}}(u) = \bar{\bar{E}}:\bar{\bar{\epsilon}}(u) \quad \text{with} \quad \bar{\bar{\epsilon}}(u) = \frac{\nabla u + \nabla^T u}{2} \quad (2) $$

* boundary conditions

$$ u(X, t) = 0 \quad \forall X \in \partial \Omega_u \quad (3) $$

$$ \bar{\bar{\sigma}}(u)\vec{n} = 0 \quad \forall X \in \partial \Omega_{\sigma} \quad (4) $$

Where:

* $u(X, t)$ is the displacement field;

* $\epsilon(u)$ is the deformation field induced by $u(X, t)$ in the framework of the assumption of small perturbations;

* $\bar{\bar{\sigma}}(u)$ is the stress field induced by $u(X, t)$;

* $\partial \Omega_u$ is the part of boundary where the displacement field is known; on $\partial \Omega_u$, boundary conditions are essential or geometrical ones;

* $\partial \Omega_{\sigma}$ is the part of the boundary where the stress field is known; on $\partial \Omega_{\sigma}$, boundary conditions are natural ones;

* $\vec{n}$ is the external normal on the boundary $\partial \Omega_{\sigma}$;

* $\rho$ is the density;

* $\bar{\bar{E}}$ is elasticity tensor.

Usually, the solution of (1) is set in a following way:

$$ u(X, t) = U(X) \cdot q(t) $$

where $U(X)$ is a displacement field (function of the variable $X$) and $q(t)$ is a scalar function of time. We obtain:

$$ div \bar{\bar{\sigma}}(U)q(t) = \rho U \ddot{q} \quad (5) $$
For longitudinal vibrations of an uniform bar, the equation (5) becomes $EU_{,x^2}q(t) = \rho U\ddot{q}(t)$. Due to stability of solution in time, it is necessary that $\ddot{q} = -\omega^2q$ with $\omega$ constant. That leads for longitudinal vibrations: $EU_{,x^2} = -\rho\omega^2U$.

For the general case, we can write:

$$ div \bar{\sigma}(U) = - \rho \lambda U \quad \text{with} \quad \lambda \ge 0 \qquad (6) $$

If the set of kinematically admissible displacement fields is defined by:

$$ U_{ad}(\Omega) = \{ \hat{W} \in H^1(\Omega) / \hat{W}(X) = 0 \forall X \in \partial\Omega_u \} $$

The solutions of the eigenvalues problem $(\lambda, U^\lambda)$ satisfies the three following relationships:

$$ (\lambda, U^\lambda) \in \mathbb{R}^+ \times U_{ad}(\Omega) \quad ; \quad div \bar{\sigma}(U^\lambda) = - \rho \lambda U^\lambda \quad ; \quad \bar{\sigma}(U^\lambda(X))\vec{n} = 0 \forall X \in \partial\Omega_\sigma $$

and therefore $\forall \hat{W} \in U_{ad}(\Omega)$, we obtain:

$$ \int_{\Omega} \bar{\sigma}(U^{\lambda}): \bar{\epsilon}(\hat{W}) d\Omega = \lambda \int_{\Omega} \rho U^{\lambda} \hat{W} d\Omega \qquad (7) $$

Owing to the fact that the stress tensor is symmetric, we can write the relation (7) as:

$$ \int_{\Omega} \bar{\sigma}(U^{\lambda}): \nabla \hat{W} d\Omega = \lambda \int_{\Omega} \rho U^{\lambda} \hat{W} d\Omega $$

The above relation corresponds to the standard variational formulation of the eigenvalue problem that occurs in Dynamics. We obtain that the solutions $(\lambda, U^\lambda)$ satisfy:

$$ (\lambda, U^\lambda) \in \mathbb{R}^+ \times U_{ad}(\Omega) \qquad (8) $$

$$ \int_{\Omega} \bar{\sigma}(U^{\lambda}): \nabla \hat{W} d\Omega = \lambda \int_{\Omega} \rho U^{\lambda} \hat{W} d\Omega \quad \forall \hat{W} \in U_{ad}(\Omega) $$

## 2.1 The X-FEM method

To solve the eigenvalue problem for time dependent structures, we use the X-FEM method which introduces a local enrichment of the approximation space in order to treat any kind of discontinuities within a structure (material discontinuity, cracks, free surfaces). These enrichment functions can include an analytical solution or an *a priori* knowledge of the experimental solution of test results. They are a combination of nodal functions that are associated with a mesh and the product of shape functions, which describe the discontinuity. This approach allows to independently model a discontinuity from the mesh. In the case of discontinuity's propagation, the method also avoids remeshing at each step. Besides, the enrichment functions are only added locally, i.e. when the domain requires to be enriched. As a result, the algebraic system of equations consists of two types of unknowns: standard degrees of freedom and enriched degrees of freedom. All the above features provide some important advantages to the X-FEM method compared to the standard FEM method for modeling arbitrary discontinuities.

The feature of X-FEM is to add special functions to some nodes of the global mesh domain. The purpose of these additional functions is to enrich the approximation of the existing displacement fields. The description of discontinuities in the framework of the X-FEM is
usually realized by the level-set method: it allows to locate the cutting elements by the level set and to determine the nodes of interest for which the additional degrees of freedom will be added. According to the method of partition of unity [3], the X-FEM displacement approximation $u^{XFEM}(x)$ is assumed to be as the following general form [6]:

$$u_{XFEM} = u_{FEM} + u_{enr} = \sum_{\substack{i \in I \\ FEM}} N_i(x)u_i + \sum_{\substack{j \in J \\ enrichment}} N_j(x)F(x)a_j \quad (9)$$

Where:

* $N_i$ are shape functions of node i;

* $u_i$ are unknown of finite element part at node i;

* $j \in J$ is nodal subset of the enrichment;

* $N_j$ are partition of unity functions of node i;

* $F(x)$ is enrichment function;

* $a_j$ are the additional degrees of freedom.

## 2.2 Definition of the changing of configuration in two-dimensional case

As it has been said, the changing of a structure can concern, for instance, the moving of its boundaries or the growth of a crack. So, for a given structure, each position of its boundaries, or each position of a crack, allows to define one configuration of this structure. For each configuration, we can associate one eigenvalue problem. Consequently, eigenvalues and eigenshapes will depend of the configuration of structure and will change as the configuration of the structure will evolves. The main benefit of X-FEM method lies in the fact that it becomes possible to follow the discontinuities within a structure without remeshing: the mesh dependency vanishes. But using merely X-FEM, the eigenvalue problem must still be solved on every configuration where an estimation of natural frequencies is needed.

To overcome the above shortcoming that can be time consuming, we propose to use the directional derivatives. They are a universal tool to avoid multiple calculations. It can be used to estimate the evolution of a quantity with respect to the change of configuration. In our case, the directional derivatives will provide a way to estimate the derivative of eigensolutions with respect to a change of configuration of the structure. Therefore, the eigenvalue problem will only be solved for some configurations until the accuracy has decreased drastically.

In order to take into account the time dependency of the structure, we introduce a function $q(X)$ which governs the changing of the shape from a reference configuration (noted $\Omega_0$) to a current one (noted $\Omega(\tau)$). $\tau$ is a scalar parameter which allows to follow the evolution of the structure and $X$ represents the position of a point on the reference configuration. If we note $x_\tau$ its position on the current configuration, we have (10):

$$x_{\tau} = X + \vec{q}(X) \cdot \tau \quad (10)$$

Concerning the function $\vec{q}(X)$ which governs the changing of configuration, the following remarks can be done:

* **Remark 1:**

The transformation between both configurations must be a bijective one.

* **Remark 2:**

The transformation doesn't change the part of boundary where displacement conditions are imposed.
### Remark 3:

It should also be noted that the changing of configuration characterizes a changing of geometry in the initial considered domain, but the material remains the same. The density and the elasticity tensor for both configurations are the same.

<figure><img src="image_2.png" /><figcaption>Figure 2: The transformation from $\Omega_0$ to $\Omega(\tau)$.</figcaption></figure>

The transformation's gradient from the reference configuration to the current one is defined by (11):

$$ \bar{\bar{F}} = \nabla \vec{x}_{\tau} = \mathbb{I} + \tau \nabla \vec{q}(\vec{X}) \quad (11) $$

The determinant of the gradient transformation can be obtained in the following way:

$$ \det \bar{\bar{F}} = \begin{vmatrix} 1 + \tau q_{1,1} & \tau q_{1,2} \\ \tau q_{2,1} & 1 + \tau q_{2,2} \end{vmatrix} = 1 + \tau \operatorname{div} \vec{q} + \frac{\tau^2}{2} ((\operatorname{div} \vec{q})^2 - \nabla^T \vec{q} : \nabla \vec{q}) \quad (12) $$

A displacement field $v$ on the current configuration can be regarded as a displacement field on the reference configuration. We have:

$$ v(x_{\tau}) = v(X + \vec{q}(X) \cdot \tau) = V(\tau, X) \quad (13) $$

It follows that:

$$ \operatorname{grad} v(x_{\tau}) = \nabla V(\tau, X) \cdot \bar{\bar{F}}^{-1} \quad (14) $$

From relation (11), it can be readily deduced that for small values of parameter $\tau$, we have:

$$ \bar{\bar{F}}^{-1} = \mathbb{I} - \tau \nabla \vec{q}(\vec{X}) + \tau^2 \nabla^2 \vec{q}(\vec{X}) \nabla \vec{q}(\vec{X}) + \tau^2 \eta(\tau) \quad \text{with} \quad \lim_{\tau \to 0} \|\eta(\tau)\| = 0 \quad (15) $$

## 2.3 Expressions of the directional derivatives for eigenvalues and eigenshapes

The solution of eigenvalue problem, illustrated on Figure 2 for initial and current configurations, are noted (16):

$$ \begin{gathered} (\lambda^i, U^i(X)) \in \mathbb{R}^+ \times U_{ad}(\Omega_0) \\ (\lambda^i(\tau), u^i_\tau(x_\tau)) \in \mathbb{R}^+ \times U_{ad}(\Omega(\tau)) \quad \text{with} \quad u^i_\tau(x_\tau) = u^i(\tau, X) \end{gathered} \quad (16) $$

Where $U_{ad}(\Omega_0)$ and $U_{ad}(\Omega(\tau))$ are the sets of admissible functions for both reference and current configurations. With these definitions, the directional derivatives of eigenvalue $D\lambda^i[q]$ and eigenshape $Du^i[q]$ for the first order are given by:
$$D\lambda^i[\vec{q}] = \lim_{\tau \to 0} \frac{\lambda^i(\tau) - \lambda^i}{\tau} \qquad (17)$$

$$Du^i[\vec{q}] = \lim_{\tau \to 0} \frac{u^i(\tau, X) - U^i(X)}{\tau} \qquad (18)$$

Starting from expressions (17) and (18), it becomes possible to obtain an estimation of eigenvalues and eigenshapes for the current configuration as soon as the directional derivatives are known. We have:

$$\lambda^i(\tau) = \lambda^i + \tau D\lambda^i[\vec{q}] \qquad (19)$$

$$u^i(\tau, X) = U^i(X) + \tau Du^i[\vec{q}] \qquad (20)$$

The both above expressions represent the connection between previous domain and current one and express the eigenvalues at the current time from the reference time as well as the eigenshape $u^i(\tau, X)$, which depends on the parameter $\tau$ and the space. More precisely, as it can be noted, the eigenvalues and the eigenshapes on the current configuration are written in terms of the solution of the eigenvalue problem on the reference configuration and its directional derivatives, both quantities being computing on the reference configuration. Therefore, knowing the solution of the eigenvalue problem on the reference configuration and its directional derivatives, the solution of the eigenvalue problem on the current configuration can be approximated, according to (19) and (20), by varying the value of parameter $\tau$.

## 2.4 Variational formulation of eigenvalue problem in two-dimensional case

To obtain the expression of the directional derivatives of eigenvalues and eigenshapes, one can start from the variational formulation of the eigenvalue problem that has to be solved in Dynamics. On the initial configuration, this variational formulation can be set as follows:

$$\int_{\Omega_0} \bar{\bar{\sigma}}(U^i(X)) : \nabla \hat{W} d\Omega_0 = \lambda^i \int_{\Omega_0} \rho U^i(X) \hat{W} d\Omega_0 \qquad (21)$$

$$\forall \hat{W} \in U_{ad}^*(\Omega_0) \text{ and } U_{ad}^*(\Omega_0) = \{\hat{W} \in H^1(\Omega_0)/\hat{W}(X) = 0 \forall X \in \partial\Omega_u(0)\}$$

The variational formulation of the eigenvalue problem on the current configuration takes the form:

$$\int_{\Omega(\tau)} \bar{\bar{\sigma}}(u_{\tau}^{i}(x_{\tau})) : \nabla \hat{W} d\Omega = \lambda^{i}(\tau) \int_{\Omega(\tau)} \rho u_{\tau}^{i}(x_{\tau}) \hat{W} d\Omega \qquad (22)$$

$$\forall \hat{W} \in U_{ad}^*(\Omega(\tau)) \text{ and } U_{ad}^*(\Omega(\tau)) = \{\hat{W} \in H^1(\Omega(\tau))/\hat{W}(x_{\tau}) = 0 \forall x_{\tau} \in \partial\Omega_u(\tau)\}$$

As the function doesn't change the part of the boundary where displacements are imposed, we have:

$$\forall \hat{W} \in U_{ad}^*(\Omega(\tau)) \Leftrightarrow \forall \hat{W} \in U_{ad}^*(\Omega_0)$$
Starting from the relation (22), using the relation (10) between $x_τ$ and $X$ and taking into consideration the relation (14), one can note that solution $(\lambda^i(\tau), u(\tau, X))$ of the eigenvalue problem on the current configuration must satisfy:

$$ \int_{\Omega_0} \bar{\bar{\sigma}}(u^i(\tau, X)) : (\nabla \hat{W} \bar{\bar{F}}^{-1}) \det \bar{\bar{F}} d\Omega_0 = \lambda^i(\tau) \int_{\Omega_0} \rho u^i(\tau, X) \hat{W} \det \bar{\bar{F}} d\Omega_0 \quad \forall \hat{W} \in U_{ad}^*(\Omega_0) \quad (23) $$

$$ \bar{\bar{\sigma}}(u^i(\tau, X)) = \bar{\bar{E}} : (\nabla u^i(\tau, X) \bar{\bar{F}}^{-1})^S $$

It has to be noticed that all integrals of the relation (23) are established on the reference configuration. Therefore, the solution of eigenvalue problem $(\lambda^i(\tau), u^i(\tau, X))$ for configuration $\Omega(\tau)$ is an element of the set $\mathbb{R}^+ \times U_{ad}^*(\Omega_0)$ satisfying for any element $\hat{W}$ of $U_{ad}^*(\Omega_0)$:

$$ \frac{\int_{\Omega_0} [\bar{\bar{E}} : (\nabla u^i(\tau, X) \bar{\bar{F}}^{-1})^S] : (\nabla \hat{W} \bar{\bar{F}}^{-1}) \det \bar{\bar{F}} d\Omega_0}{A(\tau)} = \lambda^i(\tau) \frac{\int_{\Omega_0} \rho u^i(\tau, X) \hat{W} \det \bar{\bar{F}} d\Omega_0}{B(\tau)} \quad (24) $$

In the neighborhood of $\tau = 0$, we can write $A(\tau)$ and $B(\tau)$ as an expansion in series of the variable $\tau$. By taking into account only the zero and first orders, we obtain:

$$ A(\tau) = A(0) + \left. \frac{dA}{d\tau} \right|_{\tau=0} \cdot \tau + O(\tau^2) \quad (25) $$

$$ B(\tau) = B(0) + \left. \frac{dB}{d\tau} \right|_{\tau=0} \cdot \tau + O(\tau^2) \quad (26) $$

The identification of terms of order zero, namely $A(0) = B(0)$, leads to the definition of the eigenvalue problem on the reference configuration. The equality of terms of first order leads to:

$$ \left. \frac{dA}{d\tau} \right|_{\tau=0} = \left. \frac{dB}{d\tau} \right|_{\tau=0} $$

The derivative of $A(\tau)$, with respect to the parameter $\tau$, is a sum of four terms, each one depending of $\tau$:

$$ \frac{dA(\tau)}{d\tau} = \underbrace{\int_{\Omega_0} [\bar{\bar{E}} : (\nabla \frac{\partial u^i(\tau, X)}{\partial \tau} \bar{\bar{F}}^{-1})]^S] : (\nabla \hat{W} \bar{\bar{F}}^{-1}) \det \bar{\bar{F}} d\Omega_0}_{a_1(\tau)} + \\ \underbrace{\int_{\Omega_0} [\bar{\bar{E}} : (\nabla u^i(\tau, X) \frac{\partial \bar{\bar{F}}^{-1}}{\partial \tau})]^S] : (\nabla \hat{W} \bar{\bar{F}}^{-1}) \det \bar{\bar{F}} d\Omega_0}_{a_2(\tau)} + \\ \underbrace{\int_{\Omega_0} [\bar{\bar{E}} : (\nabla u^i(\tau, X) \bar{\bar{F}}^{-1})]^S] : (\nabla \hat{W} \frac{\partial \bar{\bar{F}}^{-1}}{\partial \tau}) \det \bar{\bar{F}} d\Omega_0}_{a_3(\tau)} + \quad (27) $$
$$ \frac{\int_{\Omega_0} [\bar{\bar{E}}:(\nabla u^i(\tau, X)\bar{\bar{F}}^{-1})^S]:(\nabla \hat{W}\bar{\bar{F}}^{-1}) \frac{\partial \det \bar{\bar{F}}}{\partial \tau} d\Omega_0}{a_4(\tau)} $$

The derivative of $B(\tau)$, with respect to the parameter $\tau$, is a sum of three terms, each one also depending of variable $\tau$:

$$ \frac{dB(\tau)}{d\tau} = \underbrace{\frac{\partial \lambda^i(\tau)}{\partial \tau} \int_{\Omega_0} \rho u^i(\tau, X) \hat{W} \det \bar{\bar{F}} d\Omega_0}_{b_1(\tau)} + $$

$$ \lambda^i(\tau) \underbrace{\int_{\Omega_0} \rho \frac{\partial u^i(\tau, X)}{\partial \tau} \hat{W} \det \bar{\bar{F}} d\Omega_0}_{b_2(\tau)} + \lambda^i(\tau) \underbrace{\int_{\Omega_0} \rho u^i(\tau, X) \hat{W} \frac{\partial \det \bar{\bar{F}}}{\partial \tau} d\Omega_0}_{b_3(\tau)} \quad (28) $$

Starting from expression (27), using relations (11), (12), (14), (15) and (20), we finally obtain, by taking the limits when $\tau$ vanishes:

$$ \left. \frac{dA}{d\tau} \right|_{\tau=0} = \int_{\Omega_0} [\bar{\bar{E}}:(\nabla Du^i[\vec{q}])^S]:(\nabla \hat{W}) d\Omega_0 + \quad (29) $$

$$ \int_{\Omega_0} [-\bar{\bar{E}}:(\nabla U^i \nabla \vec{q})^S - (\bar{\bar{E}}:(\nabla U^i)^S) \nabla^T \vec{q} + (\bar{\bar{E}}:(\nabla U^i)^S) \mathrm{div} \vec{q}]: (\nabla \hat{W}) d\Omega_0 $$

In a similar way, starting from expression (28), using relations (12), (19) and (20), it can been shown, by taking the limits when $\tau$ vanishes, that:

$$ \begin{aligned} \left. \frac{dB}{d\tau} \right|_{\tau=0} ={}& D\lambda^i[\vec{q}] \int_{\Omega_0} \rho U^i \hat{W} d\Omega_0 + \lambda^i \int_{\Omega_0} \rho Du^i[\vec{q}] \hat{W} d\Omega_0 \\ & + \lambda^i \int_{\Omega_0} \rho U^i \hat{W} \mathrm{div} \vec{q} d\Omega_0 \end{aligned} \quad (30) $$

Consequently, the expressions (29) and (30) of $\left. \frac{dA}{d\tau} \right|_{\tau=0}$ and $\left. \frac{dB}{d\tau} \right|_{\tau=0}$ lead to the following relation between the directional derivatives of eigenvalue and eigenshape. We have:

$$ \begin{aligned} & \int_{\Omega_0} [\bar{\bar{E}}:(\nabla Du^i[\vec{q}])^S]:(\nabla \hat{W}) d\Omega_0 - D\lambda^i[\vec{q}] \int_{\Omega_0} \rho U^i \hat{W} d\Omega_0 \\ & \qquad - \lambda^i \int_{\Omega_0} \rho Du^i[\vec{q}] \hat{W} d\Omega_0 \\ & = \int_{\Omega_0} [\bar{\bar{E}}:(\nabla U^i \nabla \vec{q})^S + (\bar{\bar{E}}:(\nabla U^i)^S) \nabla^T \vec{q}]: (\nabla \hat{W}) d\Omega_0 \\ & \qquad + \int_{\Omega_0} [-\bar{\bar{E}}:(\nabla U^i)^S]:(\nabla \hat{W}) + \lambda^i \rho U^i \hat{W} \mathrm{div} \vec{q} d\Omega_0 \end{aligned} \quad (31) $$

Clearly, the relation (8) defines an eigenvalue problem for a self-adjoint system. The self-adjointness comes from the symmetry of integrals and the nullity of boundary conditions. It can be ascertained through integrations by parts. It follows from this self-adjointness that the
eigenfunctions $U^i(X)$ constitute a complete orthogonal set of infinite dimension [14] which we use as basis for $U_{ad}(\Omega_0)$. Consequently, using the expansion theorem for self-adjoint systems, the directional derivative of eigenshape on the reference configuration can be set as:

$$Du^i[\vec{q}] = \sum a_{ki} U^k(X) \qquad (32)$$

In relation (32), it should be pointed out that the series $\sum a_{ki} U^k(X)$ converge in energy to $Du^i[\vec{q}]$. The expressions of the directional derivatives for eigenvalues and eigenshapes can be deduced from relation (31) using the expansion (32) of $Du^i[\vec{q}]$ and the orthogonality relations between eigenshapes on the reference configuration.

### 3 RESULTS

This section deals with the numerical calculation of the eigenvalue problem for a clamped-free bar which is illustrated on Figure 3. The length of the bar which contains two different sections is named $L$. It is discretized using 200 finite elements. In this one-dimensional case, the discontinuity corresponds to the change of the cross section of the bar: in this particular case, the level set is reduced to one unique point. It follows that the enrichment function has to be built for only one element (the one containing the cross section discontinuity).

The used material is steel for which Young's modulus $E$ is equal to $2.06 \times 10^{11}$ Pa and the density is defined by $\rho = 7860 \text{ kg/m}^3$. The length of the bar $L$ is $2 \text{ m}$. The initial position of the discontinuity (change of cross section) is noted $X_0$. The sections $S_1$ and $S_2$ are equal respectively to $12.566 \times 10^{-4} \text{ m}^2$ and $8.0425 \times 10^{-4} \text{ m}^2$. Our purpose is to obtain an estimation of natural frequencies and natural shapes when the position of the discontinuity moves, that is to say when the size of each part is varied, the whole length of the bar staying unchanged.

<figure><img src="image_3.png" /><figcaption>Figure 3: Bar with a discontinuity of the cross section.</figcaption></figure>

For our problem, the bar is a free-clamped one. As the length of the whole remains constant, we have to consider a function $q(X)$ satisfying $q(0) = q(L) = 0$.

<figure><img src="image_4.png" /><figcaption>Figure 4: The transformation for one-dimensional case.</figcaption></figure>
The calculation of the directional derivatives of the eigenvalues and the eigenshapes was
carried out for one-dimensional case. Starting from the relations (19) and (22) we have:

$$
\begin{align*}
\int_0^L ES \frac{\partial}{\partial X} (Du^i[q]) \frac{\partial \hat{W}}{\partial X} dX - \lambda^i \int_0^L \rho S Du^i[q] \hat{W}(X) dX - D \lambda^i[q] \int_0^L \rho S u^i(X) \hat{W}(X) dX &= \\
&= \int_0^L \left[ ES \frac{\partial U^i(X)}{\partial X} \frac{\partial \hat{W}}{\partial X} + \lambda^i \rho S U^i(X) \hat{W}(X) \right] q'(X) dX
\end{align*}
$$

As the above expression is verified, $\forall \hat{W} \in U_{\text{ad}}^*(\Omega_0)$, we obtain, after discretization:

$$
\begin{align}
& ([K] - \lambda^i(\tau)[M])\{Du^i[q]\}_{EF} - D\lambda^i[q][M]\{U^i\}_{EF} \nonumber \\
& \quad = ([KDq] + \lambda^i(\tau)[MDq])\{U^i\}_{EF} \tag{33}
\end{align}
$$

Where:

$$
\begin{align*}
[K] &= \sum_k \int_{\text{element } k} ES \left[ \frac{\partial N}{\partial \xi} \right]^T \left[ \frac{\partial N}{\partial \xi} \right] dX \\
[K] &= \sum_k \int_{\text{element } k} ES \left[ \frac{\partial N}{\partial \xi} \right]^T \left[ \frac{\partial N}{\partial \xi} \right] dX
\end{align*}
$$

$$
\begin{align*}
[M] &= \sum_k \int_{\text{element } k} \rho S [N(\xi)]^T [N(\xi)] dX \\
[M] &= \sum_k \int_{\text{element } k} \rho S [N(\xi)]^T [N(\xi)] dX
\end{align*}
$$

$$
\begin{align*}
[KDq] &= \sum_k \int_{\text{element } k} ES \left[ \frac{\partial N}{\partial \xi} \right]^T \left[ \frac{\partial N}{\partial \xi} \right] q'(X) dX \\
[KDq] &= \sum_k \int_{\text{element } k} ES \left[ \frac{\partial N}{\partial \xi} \right]^T \left[ \frac{\partial N}{\partial \xi} \right] q'(X) dX
\end{align*}
$$

$$
\begin{align*}
[MDq] &= \sum_k \int_{\text{element } k} \rho S [N(\xi)]^T [N(\xi)] q'(X) dX \\
[MDq] &= \sum_k \int_{\text{element } k} \rho S [N(\xi)]^T [N(\xi)] q'(X) dX
\end{align*}
$$

Due to the orthogonality of eigenshapes and after pre-multiplying left and right sides of
expression (35) by {$U^i$}$_{EF}^T$, we obtain:

$$
D\lambda^i[q] = - \frac{\{U^i\}_{EF}^T ([KDq] + \lambda^i(\tau)[MDq])\{U^i\}_{EF}}{\{U^i\}_{EF}^T [M] \{U^i\}_{EF}} \quad (34)
$$

Owing to relation (32), we assume a directional derivative {$Du^i[q]$}$_{EF}$ in the form (35).

$$
\{Du^i[q]\}_{EF} = \sum_{k \neq i} \alpha_k^i \{U^k(X)\}_{EF} \quad (35)
$$

Starting from the relation (33), introducing the expression (35) of the directional derivative
{Du^i[q]}_{EF}, pre-multiplying by {U^k}^T_{EF} with k ≠ i and using the orthogonality relation
between eigenshapes on the reference configuration leads to:

$$
\alpha_k^i = \frac{\{U^k\}_{EF}^T ([KDq] + \lambda^i(\tau)[MDq])\{U^i\}_{EF}}{\{U^k\}_{EF}^T [M] \{U^k\}_{EF} \cdot (\lambda^k - \lambda^i)} \quad (36)
$$

Thus, the approximate expression of {$u^i(\tau, X)$}$_{EF}$ can be written as:

$$
\{u^i(\tau, X)\}_{EF} = \{U^i\}_{EF} + \tau \cdot \sum_{k \neq i} \alpha_k^i \{U^k\}_{EF} \quad (37)
$$

For a bar of two different sections, the exact solution of the eigenvalue problem is given by:

$$
\begin{align*}
u_1^i(X) &= \sin W_i X \quad \text{for } 0 \le X \le l_1 \\
u_2^i(X) &= A_2^i \cos W_i X + B_2^i \sin W_i X \quad \text{for } 0 \le X \le l_2
\end{align*}
$$
where $A_2^i = \sin W_i l_1$ and $B_2^i = \frac{S_1}{S_2} \cos W_i l_1$. $W_i$ satisfy:

$$ \cos W_i (l_1 + l_2) + \frac{S_1 - S_2}{S_1 + S_2} \cos W_i (l_1 - l_2) = 0 $$

It leads the opportunity to compare exact solution of this problem and the solution obtained with the proposed method.

For studying the dependence of the eigenvalues and eigenshapes with respect to the position of the discontinuity, we have to take a function $q(X)$ which allows to move the position of the discontinuity, the total length of the bar is remained to be constant. Consequently, the function $q(X)$ has to satisfy:

* $q(0) = 0$ for the fixed extremity;

* $q(L) = 0$ in order to impose the total length of the bar unchanged;

* $q(X_0) = 1$ if we have noted $X_0$ the initial position of the discontinuity.

We here present the results that we have obtained by choosing a polynomial piecewise function for $q(X)$. Two initial positions of the discontinuity have been studied:

$$ X_0 = 0.9854m \quad \text{and} \quad X_0 = 1.1854m $$

For these two cases, the functions $q(X)$ that we have considered are presented on the Figure 5.

<figure><img src="image_5.png" /><figcaption>Figure 5: Functions $q(X)$ for both initial position cases.</figcaption></figure>

The estimated values of natural frequencies have been obtained owing to relation (19). As our formulation only uses directional derivatives of first order, the approximation of the eigenvalue is linear. Directional derivatives have been computed according to relation (34). The estimated values have been compared to exact ones. To compare both solutions, the frequency deviation can be used. It corresponds to the relative error between both methods. For the i-th natural frequency, this frequency deviation is defined by

$$ \left| \frac{f_{i, \text{estimated}} - f_{i, \text{exact}}}{f_{i, \text{exact}}} \right| $$

The graphs of Figure 6, Figure 7 and Figure 8 show the evolution of the first three frequencies with respect to the position of the discontinuity. For each natural frequency, the exact values and estimated ones have been represented for two positions of the initial discontinuity.

$$ X_0 = 0.9854m \quad \text{and} \quad X_0 = 1.1854m $$
<figure><img src="image_6.png" /><figcaption>Figure 6: Dependence of the first natural frequency on position of discontinuity (directional derivatives are calculated for $X_0 = 0.9854m$ in the first case and $X_0 = 1.1854m$ in the second case).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 7: Dependence of the second natural frequency on position of discontinuity (directional derivatives are calculated for $X_0 = 0.9854m$ in the first case and $X_0 = 1.1854m$ in the second case).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 8: Dependence of the third natural frequency on position of discontinuity (directional derivatives are calculated for $X_0 = 0.9854m$ in the first case and $X_0 = 1.1854m$ in the second case).</figcaption></figure>

As seen in these graphs, the curves of exact and estimated solutions are very close. But the exact solution has a parabolic shape. The approximate solution, obtained with relation (19), is linear. On these curves, it can be observed that when directional derivatives of eigenvalue are close to zero (for instance when $X_0 = 0.9854m$), the gap between exact natural frequencies
and estimated ones increases rapidly. It provides from the fact that according to relation (19), our approximation is linear. Therefore, when directional derivatives are close to zero, the accuracy is weak. In order to overcome this drawback, it would be useful to calculate directional derivatives of second order which would allow to define a parabolic approximation of eigenvalues and eigenshapes. In the second case, when $X_0 = 1.1854m$, directional derivatives of eigenvalues are not close to zero and the gap between both solutions decreases significantly. In this second case, the linear approximation provides more accurate results.

Graphs of the Figure 9 show the frequency deviation for the three first natural frequencies when directional derivatives are calculated in $X_0 = 1.1854m$. The relative error is less than 3% for the three natural frequencies.

<figure><img src="image_9.png" /><figcaption>Figure 9: Frequency deviation for the three first natural frequencies when directional derivatives are calculated for $X_0 = 1.1854m$.</figcaption></figure>

The evolution of the three first eigenshapes is given Figure 10 for several values of parameter $\tau$. Exact eigenshapes and estimated ones are represented. The estimated eigenshapes have been obtained according to relation (20). Just like for the eigenvalues, our approximation for eigenshapes is linear. The directional derivatives $Du^i[q]$ are calculated with 20 modes.
<figure><img src="image_10.png" /><figcaption>Figure 10: Comparison of exact and approximate solutions.</figcaption></figure>

For the eigenshapes, to make an accurate comparison between estimated values and theoretical ones, we use the Modal Assistance Criterium. When we compare the estimated and theoretical eigenshapes of the same order, the MAC is close to one and decreases as $\tau$ grows up. On top of that, the MAC between estimated and theoretical eigenshapes of different order is nearly to zero. From these two results, we can deduce that
* the correlation between exact eigenshapes and estimated ones is good;

* orthogonality properties are verified between estimated eigenshapes.

<figure><img src="image_11.png" /><figcaption>Figure 11: Modal Assistance Criterium for the first natural shape (directional derivatives are calculated for $X_0 = 0.9854m$ in the first case and $X_0 = 1.1854m$ in the second case).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 12: Modal Assistance Criterium for the second natural shape (directional derivatives are calculated for $X_0 = 0.9854m$ in the first case and $X_0 = 1.1854m$ in the second case).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 13: Modal Assistance Criterium for the third natural shape (directional derivatives are calculated for $X_0 = 0.9854m$ in the first case and $X_0 = 1.1854m$ in the second case).</figcaption></figure>
# 4 CONCLUSION

A methodology for solving the eigenvalue problem for a time dependent structures has been proposed. This methodology combines both the directional derivatives and X-FEM methods. The convenience of the use of X-FEM is the possibility to avoid remeshing. A discontinuity can be depicted with a Level Set function. This Level Set function is used to define a local enrichment for taking into account this discontinuity. Consequently, due to enrichment, the mesh is not mandatory to be conform to the geometry. The application of directional derivative allows to avoid many calculations at every time step and gives an approximation of the solutions of the eigenvalue problem. This study focuses on the fact that directional derivatives is a convenient tool for following the evolution of natural frequencies and shapes.

The main idea of our approach consisted in obtaining the expression of derivatives for eigenvalues and eigenshapes in a given direction which is described by a function $q(X)$. This function determines the transformation from a reference configuration to a current one. Starting from the variational formulation of the eigenvalue problem on the current configuration, derivatives of eigenvalues and eigenshapes can be expressed in terms of the solutions of the eigenvalue problem for the reference domain by the mean of this function $q(X)$. Once we have computed the directional derivatives of eigenvalues and eigenshapes, it becomes possible to evaluate the evolution of these values when the domain is modified. To obtain solutions at the required time step (that is to say on the current configuration), it is sufficient to know information from the initial one. The main advantage of this procedure lies in the fact that numerical calculations are done without time consuming.

The numerical results for a one-dimensional structure varying in time and the accuracy of the new technique were presented. Geometric transformation of the eigenvalue problem for one-dimensional element with discontinuity was defined. To compute the natural shapes for a current configuration, the directional derivatives of natural shapes have been performed on the reference configuration, these directional derivatives have been set in terms of natural shapes of the reference configuration.

Using the general methodology which combines of the offered methods gives possibility to illustrate the application of this technique on a two-dimensional eigenvalue problem for a plate with shapes which are moving. One of the aspects for future work is to look upon the expressions of directional derivatives of eigenvalue problem for second order to improve the existing linear approximation. Another issue which has to be solved is to find a criteria to determine when accuracy has decreased significantly which means that directional derivatives have to be computed. The numerical implementation will be tested on on real simple and more complex cases.

# REFERENCES

[1] T. Belytschko and T. Black, Elastic crack growth in finite elements with minimal remeshing. *International journal for numerical methods in engineering, 45*, 601-620, 1999.

[2] N.Moës, John Dolbow, and T. Belytschko, A finite element method for crack growth without remeshing. *International journal for numerical methods in engineering, 46*, 131-150, 1999.
[3] J.M. Melenk and I. Babuska, The partition of unity finite element method: basic theory and applications. *Computer Methods in Applied Mechanics and Engineering*, **139**, 289-314, 1996.

[4] S. Osher and J. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. *Journal of Computational Physics*, **79**, 712-749, 1988.

[5] M. Duflot, A study of representation of crack with level sets. *International journal for numerical methods in engineering*, **70**, 1261-1302, 2007.

[6] S. Bordas and B. Moran, Enriched finite elements and level sets for damage tolerance assessment of complex structures. *Engineering Fracture Mechanics*, **73**, 1176-1201, 2006.

[7] M. Stolarska, D. Chopp, Modeling thermal fatigue cracking in integrated circuits by level sets and the extended finite element method. *International Journal of Engineering Science*, **41**, 2381-2410, 2003.

[8] D.L. Chopp, N. Sukumar, Fatigue Crack Propagation of Multiple Coplanar Cracks with the Coupled Extended Finite Element/Fast Marching Method. *International Journal of Engineering Science*, **41** (8), 845-869, 2003.

[9] N. Sukumar, D. Chopp, N. Moës and T. Belytschko, Modelling holes and inclusions by level sets in the extended finite-element method. *Computer Methods in Applied Mechanics and Engineering* 190, 6183-6200, 2001.

[10] N. Moës, M. Cloirec, P. Cartraud, J.F. Remacle, A computational approach to handle complex microstructure geometries. *Computer Methods in Applied Mechanics and Engineering*, **192**, 3163-3177, 2003.

[11] E. Taroco, Shape sensitivity analysis in linear elastic fracture mechanics. *Computer Methods in Applied Mechanics and Engineering*, 697-712, 2000.

[12] K. Choi, N. Kim, *Design Sensitivity Analysis of Linear Structural Systems*. In K. Choi. N. Kim, Structural Sensitivity Analysis and Optimization 1: Linear Systems and Applications. Berlin: Springer Science & Business Media, 2006.

[13] J.P. Zolésio, *The material derivative (or speed) method for shape optimization*. *Optimization of Distributed Parameter Structures*, Sijthoff and Noordhoff, Alphen aan den Rijn, Netherlands, 1089-1151, 1981.

[14] L. Meirovitch, *Principles and Techniques of Vibrations*, Prentice-Hall International Editions, 1997.