# Low intrusive coupling of implicit and explicit integration schemes for structural dynamics: application to low energy impacts on composite structures

Teddy Chantrait, Johann Rannou, Anthony Gravouil

<figure><img src="image_1.png" /></figure>

## To cite this version:

Teddy Chantrait, Johann Rannou, Anthony Gravouil. Low intrusive coupling of implicit and explicit integration schemes for structural dynamics: application to low energy impacts on composite structures. V International Conference on Computational Methods for Coupled Problems in Science and Engineering, Jun 2013, Barcelone, Spain. pp.1373-1382. hal-04710799

HAL Id: hal-04710799

https://hal.science/hal-04710799v1

Submitted on 9 Jul 2025

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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
LOW INTRUSIVE COUPLING OF IMPLICIT AND
EXPLICIT INTEGRATION SCHEMES FOR STRUCTURAL
DYNAMICS: APPLICATION TO LOW ENERGY IMPACTS
ON COMPOSITE STRUCTURES

T. CHANTRAIT*, J. RANNOU* and A. GRAVOUIL†

*Onera, CC/DMSM/MNU
29 avenue de la division Leclerc
92322 Châtillon, France
e-mail: teddy.chantrait@onera.fr - Web page: http://www.onera.fr

†Laboratoire de Mécaniques des Contacts et des Structures (LaMCoS, CNRS UMR 5259)
Institut National des Sciences Appliquées de Lyon
20, avenue Albert Einstein
F69621, Villeurbanne cedex, France
e-mail: anthony.gravouil@insa-lyon.fr - Web page: http://lamcos.insa-lyon.fr

**Key words:** Codes coupling, Explicit, Implicit, Low intrusive, Composite, Low energy impact

**Abstract.** Simulation of low energy impacts on composite structures is a key feature in aeronautics. Unfortunately they are very expensive: on the one side, the structures of interest have large dimensions and need fine volumic meshes (at least locally) in order to capture damages. On the other side small time steps are required to ensure the explicit algorithms stability which are commonly used in these kind of simulations [4]. Implicit algorithms are in fact rarely used in this situation because of the roughness of the solutions that leads to prohibitive expensive time steps or even to non convergence of Newton-like iterative processes. It is also observed that rough phenomenons are localized in space and time (near the impacted zone). It may therefore be advantageous to adopt a multiscale space/time approach by splitting the structure into several substructures owning there own space/time discretization and their own integration schemes. The purpose of this decomposition is to take advantage of the specificities of both algorithms families: explicit scheme focuses on rough areas while smoother (actually linear) parts of the solutions are computed with larger time steps with an implicit scheme. We propose here an implementation of the Gravouil-Combescure method (GC) [1] by the mean of low intrusive coupling between the implicit finite element analysis (FEA) code *Z-set* and the explicit FEA code *Europlexus*. Simulations of low energy impacts on composite stiffened panels are presented. It is shown on this application that time step ratios up to 5000
can be reached. However, computations related to the explicit domain still remain a
bottleneck in terms of cpu time.

# 1 INTRODUCTION

Low energy impacts can be very harmful in particular for composite structures used
in the aerospace industry. In fact, they can cause significant damage (matrix cracking,
fiber failure, delamination...) inside the composite or on the side opposite to the impact.
However, the residual print it leaves on the impacted side can be almost undetectable to
the naked eye. The damage caused can therefore lead to early failure of the structure while
they can be unnoticed during a visual inspection, this is related to the concept of BVID
(Barely Visible Impact Damage). We thus understand the importance for manufacturers
to control such situations. Numerical simulations of this phenomenon can be a great help
especially to orient and to rationalize tests campaigns by the use of virtual testing. Various
researches are led in scientific and industrial communities to simulate these impacts but
these one remain currently difficult to implement on an industrial scale because they are
very difficult to manage. Indeed, sources of non-regularities introduced in the models
(contact, softening damage laws, cohesive zone models...) make convergence difficult
to achieve for implicit algorithms. Simulations performed in this context are therefore
mostly conducted though explicit time integration [3, 4]. This make it possible to better
take into account these sources of non-regularities but they require the use of time steps
depending especially on the smallest mesh element to ensure stability. Moreover, very
fine meshes are usually required (at least locally) to capture the non-linear phenomena
occuring during impact. This thus leads to a very large number of increments which can
be prohibitive. Note however that these non-linear phenomena occur on a very localized
area around the impact point. Adopting a space/time multiscale strategy thus appears
to be advantageous to solve this kind of multiscale problems. This can be performed
through domain decomposition where each subdomain owns its own time discretization.
The purpose of this decomposition is to focus on numerical computation where non-linear
phenomena appear [5]. Explicit resolution in the area close to the impact is required
because of the roughness of the solution. However, on the complementary area where the
solution is smoother, implicit integration is appropriate. Larger time steps can then be
used thus saving cpu time. The work done here is based on the Gravouil-Combescure (GC)
method [1] and a low-intrusive coupling [6] between implicit finite element analysis codes
$Z-set^1$ and explicit FEA code *Europlexus^2* has been realized. The strategy is schematically
shown in Figure 1. It shows one section of an impacted plate and a typical mesh. This
mesh is divided into two domains: an impacted domain (center) which is processed by
*Europlexus* with fine time-stepping and a complementary domain which is processed by

<footnote n="1">¹$Z-set$ is developed by Mines ParisTech, Onera and NW Numerics & Modeling</footnote>

<footnote n="2">²*Europlexus* is developed by CEA and the Joint Research Centre in Ispra, Italy</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1: Overview of a multiscale space/time approach applied in the context of low-energy impact on composite structures.</figcaption></figure>

Z-set with larger time steps.

## 2 OVERVIEW OF THE USED MULTISCALE SPACE/TIME METHOD

Different multiscale space/time approaches can be found in literature [1, 7, 8, 9]. These can be viewed as extensions of dual domain decomposition methods without overlapping, conventionally used for parallel computing [10]. Indeed, the domain decomposition methods consist in splitting spatially a structure into several subdomains and search solution on each subdomain as independently as possible (an example of structure split into two domains is shown in dotted box in Figure 2). A key point in domain decomposition methods is to determine the special boundary conditions that must be applied on the subdomains interface. This boundary condition should indeed ensure both kinematic continuity and interface equilibrium. An additional interface problem is then solved at each time step to determine this boundary condition and its size is determined by the spatial discretization of the interface. At this stage we may note that from the discrete point of view, imposing the continuity of the displacement, velocity or acceleration is not equivalent. The interface problem can then be solved directly after building interface operator (based on domains Schur complements) or be solved through iterative methods (e.g. conjugate gradient method) which do not require explicit construction of the interface operator.

The additional feature of multiscale space/time methods is that each subdomain have its own time discretization. A diagram is shown in Figure 2 where a fine time stepping is associated with the subdomain $Ω_1$ and a coarse time stepping associated with the
subdomain $\Omega_2$. Let $\Delta t$ be the fine scale time step and let $\Delta T$ be coarse scale one. $m$ is the ratio between these two time steps: $m = \frac{\Delta T}{\Delta t}$. The interface problem has to be rewritten in order to satisfy the continuity and the balance conditions in time. The GC method proposes to ensure these conditions at every fine scale time step. This allows to naturally handle non-linearities in the explicit domain which is not the case for the other methods. Interfacial velocity field of the coarse time scale is not known at each fine scale step of the fine time scale, so it is linearly interpolated [2]. The interface problem to solve is only slightly modified in comparison with a dual domain decomposition methods with a single time scale. If the problem to solve on the implicit domain is linear, the interface operator remains constant during the computation. Therefore, because of the large number of steps to achieve due to the CFL condition, the interface operator worth being explicitely build.

<figure><img src="image_4.png" /><figcaption>Figure 2: Overview of multiscale space/time methods.</figcaption></figure>

## 3 IMPLEMENTATION OF THE GRAVOUIL-COMBESCURE MULTISCALE SPACE/TIME METHOD

The particularity of the approach proposed here is to adopt a code coupling point of view rather than algebraic one. We thus achieve a structure/structure code coupling where each subdomain is associated with different computer codes which have their own characteritics. The advantage of this type of coupling is to extend the GC method to large problems while taking advantages of features existing in each FEA codes (contact algorithms, cohesive zone elements, material models...). The potential of the GC method is enhanced through application to a large low-energy impact problem (see section 4). The implementation of the method has led to the identification of the minimal entry points that codes must provide. With this minimal interface, the method can be easily implemented without being intrusive. Only the development of this interface could be
<figure><img src="image_5.png" /><figcaption>**Figure 3:** Overview of implicit/explicit code coupling realised based on the GC methode.</figcaption></figure>

few intrusive. These entry points include:

- A method to apply a nodal force boundary condition at the beginning of each increment on area of mesh.

- A method to get nodal kinematic values on area of the mesh at the end of each step (at least the velocities).

- A method to validate or invalidate an increment.

In addition, each code must provide a trace operator.

This interface was developed in the two mentioned FEA codes (*Z-set* and *Europlexus*) through a Python layer as can be schematically seen in Figure 3. Python scripts containing the GC method and tools to solve interface problem are also written. These scripts contain MPI instructions needed to exchange data between the codes. The interface problems are solved in the script attached to the code handling the finest time discretization (i.e the explicit code). This allows to minimize MPI communications. Different test cases were then carried out on 2D and 3D structures. Some proposed in [7] were reproduced in order to validate the implementation of the method.

# 4 APPLICATIONS OF THE MULTISCALE COUPLING METHOD ON STIFFENED COMPOSITE PANEL IMPACTED

Simulations with different time steps ratios were performed on a stiffened composite panel relatively representative of an aircraft's subassembly. Figure 4 illustrates the geometry, the mesh and the domain decomposition adopted. The central impacted area of the panel (in blue) as well as the impactor are processed with *Europlexus* and the complementary part which has relatively larger time step is processed with *Z-set*. The stacking sequence is [90/45/0/-45]s for the stiffeners and [90/45/0/0/-45]s for the skin. The size of elements (hexaedrons and wedges) ranges from 5mm to 0.25mm in the impact area. Note that each ply is modeled with one element in the thickness direction.

The size of the problem is about one million degrees of freedom with 90% of nodes located in the implicit domain. The size of the interface is about 3000 degrees of freedom.
<figure><img src="image_6.png" /><figcaption>Figure 4: Example of domain decomposition on stiffened composite panel split in two areas (one implicit subdomain and nine explicit subdomains).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: Elastic-damageable behavior illustrated through traction/compression in fiber and matrix direction.</figcaption></figure>
A 10J impact at 10m/s is applied. Time step of the impacted area is set up to 2.10⁻⁸s which is in accordance with the CFL condition. We then perform computations with different time step ratios with *m* ranging from 500 to 5000.

<figure><img src="image_8.png" /><figcaption>Figure 6: Matrix damage obtained on stiffened composite panel for a 10J impact and with a time step ratio equal to 500.</figcaption></figure>

A elastic-domageable material model based on Onera Progressive Failure Model (OPFM) [11] is used on the explicit domain. Its parameters was identified on T700/M21 whose behavior is shown in Figure 5 for deformation controlled loading in traction/compression. One of the results obtained with this model for a time step ratio equal to 500 is shown in Figure 6. We observe qualitatively cone-shaped matrix damage typically observed on this kind of impact.

We present in figure 7 the time evolution of the composite panel deflection for two time step ratios (m=500 and m=5000). Errors are also represented, they are defined by the absolute value of the difference between the reference solution considered (here with a time step ratio m=500) and the current solution normalized by the maximum displacement. The solutions appear to be very similar over the whole time range. However we note some differences on the error at 0.0036s and 0.0052s. These differences may be explained by the large time step ratio used for m=5000: only 60 increments are performed which is very coarse. Moreover, impactors rebounds are observed there, they are thus only poorly captured. We can see in figure 7 that the error still remains relatively low.
<figure><img src="image_9.png" /><figcaption>Figure 7: Time evolution of the stiffened composite panel deflection impacted at 10J for different time step ratios.</figcaption></figure>

# 5 CONCLUSION AND PERSPECTIVES

Structure/structure code coupling has been implemented using the Gravouil-Combescure algorithm. Both implicit (Z-set) and explicit (Europlexus) codes were involved. The minimal interface which sould be provided by the FEA codes to implement this method were identified. We have also highlighted the potential of this method on large calculations. The cpu time distributions presented in Figure 8 were measured on stiffened panel simulations for different time step ratios (50, 500 and 5000) with a linear elastic model and a domain decomposition slightly different from the one presented.

We observe that the choice to solve the interface problem on the finest time step is not too penalizing in term of cpu time with this interface size (about 3000 degrees of freedom). We also note that the overall computation cpu time decreases when the time step ratio increases. This shows that in this case the computation work is focused on the explicit area. Indeed, for time step ratios m=500 and m=5000, times spent in the implicit process are very small (less than three hours for m=500 and one hour for m=5000). However, this time saving tends to reach a limit which corresponds to the time spent in the explicit process. Work is currently underway to improve load balancing in order to reduce the computation time spent in the explicit process. Extensions that allow to use several explicit codes has been made. Calculations performed on the composite panel was made with these developments. The explicit part have been split in several subdomains as shown in Figure 4 (the explicit contains nine ring-shape subdomains centered arround the point of impact). However, they all have the same time step (m=1) which is equivalent in this case to perform a classical parallelization of the explicit code Europlexus. This
<figure><img src="image_10.png" /><figcaption>Figure 8: Computational time distribution for different time step ratios used with linear elastic stiffened composite panel impacted 10J.</figcaption></figure>

parallelization is however carried out by means of the algorithm developed (outside of the code). We also continue to enrich the model of the impacted area with feature such as cohesive zone elements or full OPFM model in order to perform experimental and numerical comparisons.

## Acknowledgements

The authors would like to thank Vincent Faucher from CEA and Steve Belon and Roland Ortiz from Onera/DADS for their help on Europlexus. This work have been performed in the framework of the Onera project "PRF transition statique/dynamique dans les matériaux composites"

## REFERENCES

[1] A. Combescure, A. Gravouil. *A time-space multi-scale algorithm for transient structural nonlinear problems*, Mécanique & Industries, Elsevier, 2(1), p. 43-55, 2001.

[2] A. Gravouil, A. Combescure. *Multi-time-step explicit-implicit method for non-linear structural dynamics*, International Journal for Numerical Methods in Engineering, Wiley Online Library, 38(1), p. 199-225, 2001.
[3] C. Lopes. Damage and Failure of non-conventional composite laminates, PhD thesis, TU Delft, Faculteit Luchtvaart-en Ruimtevaarttechniek, 2009.

[4] C. Bouvet. Etude de l'endommagement dans les structures composites, habilitation à diriger des recherches, Université de Toulouse III, Paul Sabatier, 2009.

[5] C. Dupleix-couderc. Sur une stratégie multi-échelle danalyse des grands délaminages en dynamique transitoire, Thèse, LMT-Cachan, 2011.

[6] L. Gendre, O. Allix, P. Gosselet, F. Comte. Non-intrusive and exact global/local techniques for structural problems with local plasticity, Computational Mechanics, Springer, 44(2), p. 233-245, 2009.

[7] A. Prakash, K.D. Hjelmstad, A FETI-based multi-time-step coupling method for Newmark schemes in structural dynamics, International journal for numerical methods in engineering, Wiley Online Library, 61(13), p. 2183-2204, 2004.

[8] N. Mahjoubi, A. Gravouil, A. Combescure, N. Greffet. A monolithic energy conserving method to couple heterogeneous time integrators with incompatible time steps in structural dynamics, Computer Methods in Applied Mechanics and Engineering, Elsevier, 2010.

[9] A. Ghanem, M. Torkhani, N. Mahjoubi, T. N. Baranger, A. Combescure. Arlequin framework for multi-model, multi-time scale and heterogeneous time integrators for structural transient dynamics, Computer Methods in Applied Mechanics and Engineering, Elsevier, 2012.

[10] C. Farhat, L. Crivelli, F.X. Roux. A transient FETI methodology for large-scale parallèle implicit computations in structural mechanics., International Journal for Numerical Methods in Engineering, Wiley Online Library, 37(11), p. 1945-1975, 1994.

[11] F. Laurin, N. Carrere, J.-F. Maire. A multiscale progressive failure approach for composite laminates based on thermodynamical viscoelastic and damage models, Composites Part A, Applied Science and Manufacturing, Elsevier, 38(1), p. 198-209, 2007.