# On the computation of plate assemblies using realistic 3D joint model: a non-intrusive approach

Guillaume Guguin, Olivier Allix, Pierre Gosselet, Stéphane Guinard

► To cite this version:

Guillaume Guguin, Olivier Allix, Pierre Gosselet, Stéphane Guinard. On the computation of plate assemblies using realistic 3D joint model: a non-intrusive approach. Advanced Modeling and Simulation in Engineering Sciences, 2016, 3 (16), 10.1186/s40323-016-0069-5. http://hal-01350536

HAL Id: hal-01350536

https://hal.science/hal-01350536

Submitted on 30 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>
# On the computation of plate assemblies using realistic 3D joint model: a non-intrusive approach

Guillaume Guguin¹, Olivier Allix,¹ Pierre Gosselet¹, Stéphane Guinard²

1 LMT-Cachan, ENS-Cachan/CNRS/Université Paris-Saclay,
61 avenue du Président Wilson, 94235 Cachan, France

2 Airbus Group Innovation,
18 rue Marius Tercé Zac St Martin du Touch, 31000 Toulouse, France

July 30, 2016

## Abstract

Most large engineering structures are described as assemblies of plates and shells and they are computed as such using *adhoc* Finite Element packages. In fact their computation in 3D would be much too costly. In this framework, the connections between the parts are often modeled by means of simplified tying models. In order to improve the reliability of such simulations, we propose to apply a non-intrusive technique so as to virtually substitute the simplified connectors by a precise 3D nonlinear model, without modifying the global plate model. Moreover each computation can be conducted on independent optimized software. After a description of the method, examples are used to analyze its performance, and to draw some conclusions on the validity and limitation of both the modeling of junction by rigid connectors and the use of submodeling techniques for the estimation of the carrying capacity of bolted plates.

**keywords: assembly, non-intrusive coupling, bolt**

## 1 Introduction

The simulation of large structures undergoing complex local nonlinear phenomena is still a major scientific and industrial challenge. One of the main difficulties originates from the difference of length scale between the global response of the structure and the localized phenomena. To address those problems a first type of computational approach is based on homogenization, as $FE^2$ [1] but it works well as long as the scales are sufficiently separated. To overcome this limitation concurrent multiscale methods have been developed. They are often based on domain decomposition techniques like FETI [2], FETI-DP[3] or the LATIN multiscale method [4, 5] and its optimization for the multiscale treatment of nonlinear problems [6, 7].
Moreover most of large industrial structures are described as an assembly of plates and shells, whereas local phenomena often require 3D models to be properly analyzed. To deal with such problems, several methods have been applied or developed for the coupling of 2D and 3D models, like the Arlequin method [8, 9], transition elements [10, 11], MPCs approaches [12, 13] or Nitsche's method [14].

Most of these methods are quite demanding in terms of software development and therefore they are seldom used in industrial packages. To overcome these drawbacks, non-intrusive approaches have recently been proposed [15]. They are nowadays the subject of extensions and developments: thermoelasticity with GFEM/FEM coupling [16], crack propagation in XFEM/FEM coupling [17], stochastic simulations [18] and dynamics [19, 20].

In [21] a non-intrusive coupling between plate and 3D models was proposed in the case of linear behaviors. The present paper concerns the extension of this approach to the simulation of bolted assemblies of plates where bolts are described with full 3D nonlinear models. Such structures are good candidates for the iterative global-local non-intrusive strategy for two main reasons. First, the detailed computation of a tightened bolt with frictional contact on all surfaces is a very hard and time-consuming task to perform using commercial software. Second, the construction of the reference model, which would correspond to the assembly of a plate model and a 3D model for the bolt, would be very complex. The non-intrusive framework provides answers to both these problems. First, it allows the use of dedicated software for the local computations (in our case, COFAST a parallel software based on the LATIN domain decomposition method [22]). Second, it allows the easy coupling of a general Finite Element software (here Code\_Aster from EDF) for the plate computation, with COFAST, because the global model is unchanged during the iterative process. These properties were exploited in [23] for the simulation of damage in composite laminates at the meso and the micro scales using dedicated pieces of software.

The non-intrusive framework aims at solving the reference problem iteratively, by solving at each iteration both the global problem with prescribed residual traction at the interface and the local nonlinear problems submitted to prescribed displacement. Several techniques have been proposed to improve the convergence rate of the method by means of acceleration techniques [15, 17, 24] or improved interface conditions [25]. The method has therefore common points with so called nonlinear domain decomposition methods (or nonlinear relocalization techniques) [26, 27] which proved their efficiency and gain in robustness in the case of buckling [28], post-buckling [29] and damage analysis [30, 31]. Other proposals have been made to take into account the fact that the phenomena of interest are localized, aiming at a better representation of the target model [32, 33, 16].

The paper is organized as follows. In section 2, a summary of the reference problem corresponding to the coupling of 2D and 3D models according to [21] is presented. The non-intrusive algorithm is presented in section 3. In section 4, the results of the iterative coupling are analyzed in the case of a bolted joint.
These results are compared to those corresponding to the plate solution and to the submodeling technique. In section 5, the influence of some parameters is assessed regarding the rate of convergence of the iterative process and regarding the accuracy of the coupled model compared to a full 3D solution.

## 2 The reference problem

We consider the typical problem of two plates connected by a bolt. Starting from a simplified model of the assembly using plates elements and a simple connector (typically a beam), our aim is to perform the non-intrusive substitution of the connector by a full 3D model of the bolt. A sketch of the two models is presented in FIG. 1. In what follows lower case Greek characters are used to describe the geometry of the plates whereas capital Greek characters are used for the geometry of the 3D domains: any 3D domain $\Omega_X$ possesses a plate counterpart $\omega_X$, so does any interface $\Gamma_Y$ whose plate counterpart is written $\gamma_Y$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Left, global model: two plates assembled by one 1D connector. Right, local model: full representation of the bolt with surrounding plates.</figcaption></figure>

The target hybrid model, which will be precisely defined in this section, corresponds to plate models connected by a full 3D model of the bolt. At the convergence of the iterations, the solution of the hybrid model is obtained, as illustrated in FIG. 2. Let us note that since the problem is symmetric (in the $(e_x, e_z)$ plane) only half of the bolt is computed and shown on the figure.

### 2.1 The plate model

The global plate model is the assembly of two plates, $\omega^{\text{inf}}$ and $\omega^{\text{sup}}$ respectively the lower plate and the top plate, and a rigid connector between them $\omega^{\text{conn}}$. In FIG. 1 (left). The plates are 20 mm thick and 280 mm long. The lower plate is 160 mm wide whereas the top plate is 80 mm wide. In order to simplify the presentation, the plates are assumed to be made out of homogeneous isotropic linear elastic material with Young modulus $E = 200$ GPa and Poisson ratio $\nu = 0.3$; the handling of orthotropic composite plates is explained in [21].
<figure><img src="image_3.png" /><figcaption>Figure 2: Deformed shape of the targeted coupled problem.</figcaption></figure>

The lower plate $\omega^{inf}$ is clamped on the part $\gamma_0^{inf}$ of its boundary and a prescribed tension displacement ($e_x$ direction) is applied on the part $\gamma_L^{sup}$ of the top plate's boundary, FIG. 1:

$$ \left\{ \begin{array}{ll} \underline{U} = 0 & \text{on } \gamma_0^{inf} \\ \underline{U} = U e_x & \text{on } \gamma_L^{sup} \end{array} \right. \qquad (1) $$

Note that this tension is applied after preload (tightening) was applied to the bolt (see the description of the patch).

A classical Reissner-Mindlin plate formulation is used, and the contact is not taken into account at this stage. The kinematic is given by:

$$ \underline{U}(m, z) = \underline{V}(m) + w(m) \underline{e}_z - \underline{e}_z \wedge \underline{\theta} = \underline{V}(m) + w(m) \underline{z} - z \begin{bmatrix} \theta_y(m) \\ -\theta_x(m) \\ 0 \end{bmatrix} \qquad (2) $$

where $m \in \omega^{inf} \cup \omega^{sup}$ is a point of the mid-surface of the plates and $e_z$ the normal direction (in the thickness). The stress in the plates is described with the eight generalized forces ($\underline{N}$ tension, $\underline{M}$ bending, $\underline{Q}$ shear):

$$ \mathbb{F} = (N_{xx}, N_{yy}, N_{xy}, M_{xx}, M_{yy}, M_{xy}, Q_x, Q_y) \qquad (3) $$

Since the connector will ultimately be virtually replaced by the 3D model described in next subsection, a very crude connector is chosen: a perfectly rigid beam element surrounded by a rigid region of the size of the bolt (20 mm) in order to avoid localization around one node, see FIG. 3. Beside its simplicity, this connector has the advantage to be very rigid compared to its 3D counterpart, which is a sufficient condition to ensure the convergence of the coupling algorithm [15, 34].

## 2.2 The 3D patch

The local patch is designed to replace the connector. Thus, a full 3D representation of the bolt is used with unilateral frictional contact interfaces between each part of the assembly. The dimensions of the bolt are given in FIG. 4(a). The problem being symmetric, only half of the bolt is computed. The material
<figure><img src="image_4.png" /><figcaption>Figure 3: The definition of the connector.</figcaption></figure>

used for the screw and the nut is linear elastic with Young’s modulus $E = 300$ GPa and Poisson’s ratio $\nu = 0.3$. Coulomb’s friction coefficient is equal to 0.3 on all interfaces.

To solve this model, we use the dedicated contact solver COFAST, developed in the LMT-Cachan by L. Champaney [22]. This software uses the LATIN method [35] where interfaces between subdomains are the support of the contact conditions, see FIG. 4(b).

<figure><img src="image_5.png" /><figcaption>Figure 4: 3D representation of the bolt.</figcaption></figure>

Note that before coupling with the global problem, preload is applied by enforcing relative displacements between the nut and the screw. It is adjusted to match realistic values of tension in the screw (200 MPa).

## 2.3 Connections between the models

The plates+connector model describes the entirety of the structure and occupies the (2D+1D) domain $\omega$. The 3D model of the bolt occupies what we call the zone of interest $\Omega_I$. In the original method [15], the coarse modeling of $\omega_I = \omega \cap \Omega_I$ was simply replaced by $\Omega_I$ through iterations; in [21], it was shown that because of the edge effects which affect plate solutions, it is interesting to introduce a zone of transition between the two models.
This idea is sketched in Fig. 5: the 3D model is the only one taken into account in the inner zone of interest $\Omega_I \subset \Omega_I$, the plate model is the only one taken into account in the outer complement zone $\omega_C \subset \omega_C$, there exists an overlap $\Omega_I \cap \omega_C$ also called buffer zone where the two models are equivalent in a certain sense.

Two interfaces are defined: $\gamma_I = \partial\Omega_I \cap \omega_C$ and $\gamma_C = \partial\omega_C \cap \Omega_I$. In order to simplify the example, regular and ruled quadrilateral meshes are used along the interfaces. Thus, $\gamma_I$ and $\gamma_C$ are straight lines between two conforming meshes, see Fig. 6.

<figure><img src="image_6.png" /><figcaption>Figure 5: The domain decomposition with overlap and the interfaces.</figcaption></figure>

The substituted solution is expressed with the following form:

$$ \underline{u}^{hyb} = \begin{cases} \underline{u}^L & \text{in } \Omega_I \text{ solves the 3D equations} \\ \underline{u}^G & \text{in } \omega_C \text{ solves the plate equations} \end{cases} \qquad (4) $$

and it satisfies the following transmission conditions (in a sense which is made more precise afterwards):

$$ \left\{ \begin{array}{l} \underline{u}^L \text{ and } \underline{u}^G \text{ are continuous on } \Gamma_I \\ 3D \text{ stress and plate generalized forces are balanced on } \gamma_C \end{array} \right. \qquad (5) $$

To be more specific, let us assume that a finite element discretization is used: $\mathbf{u}^G$ is the vector of plate unknowns defined on $\omega$, $\mathbf{u}^L$ is the vector of 3D unknowns defined on $\Omega_I$, $\mathbf{f}_{ext}$ stands for the generalized forces (they include the effect of imposed displacements), $\mathbf{K}$ for the stiffness matrices, and $\mathbf{f}_{int}$ for the internal forces (work of the stress field in the finite element subspace).

The global plate model is enriched by an extra loading $\delta$ which is only non-zero on the degrees of freedom of the interface $\gamma_C$ (the values of $\delta$ are determined by the coupling algorithm). The global plate model equilibrium thus can be written as:

$$ \mathbf{K}^G \mathbf{u}^G = \mathbf{f}_{ext}^G + \delta \qquad (6) $$

Note that the assumption of a linear global model simplifies the explanation and speeds up the solving (because the factorization is done only once) but it is not a requirement of the coupling algorithm.
The local 3D equilibrium submitted to boundary conditions inherited from the plate computation can be written as:

$$ \left\{ \begin{array}{l} \mathbf{f}_{\text{int}}^{L}(\mathbf{u}^{L}) = -\mathbf{f}_{\text{ext}}^{L} \quad \text{in } \Omega_{I} \\ \mathbf{u}^{L} = \mathbf{R}\mathbf{u}^{G} + \mathbf{w}^{G}(\mathbf{u}^{G}) \quad \text{on } \Gamma_{I} \end{array} \right. \qquad (7) $$

Matrix $\mathbf{R}$ represents the classical 3D rigid section movements. It enables to express the 3D Reissner-Mindlin kinematic from the plate degrees of freedom ($\mathbf{u}^G$). Because the rigid section kinematic is a very coarse hypothesis, we supplement it with warping displacements $\mathbf{w}^G(\mathbf{u}^G)$, proportional to the plate stress state. To be more precise, the applied displacement on $\Gamma_I$ has the following form:

$$ \underline{u}^L(m, z) = \frac{\underline{u}(m) + w(m)\underline{e}_z - \underline{e}_z \wedge \theta}{\mathbf{R}\mathbf{u}^G} + \sum_{i=1}^{8} \frac{\underline{\mu}_i^I(z)\mathbb{F}_i^{\tilde{C}}}{\mathbf{w}^G(\mathbf{u}^G)} $$

with $(\mathbb{F}_i^{\tilde{C}})$ the eight generalized forces extracted along the interface $\gamma_I$ (3). $(\underline{\mu}_i^I(z))$ is the basis of warping adapted to the stacking and materials which was determined before the computation by solving preliminary problems [21]. These preliminary problems conducted on a 3D representative volume element of the composite plate enable us to compute distributions of stress and displacement (warping) associated with any solution of Saint-Venant problems. In turn, Saint-Venant problems are associated with specific components of the plate generalized stresses.

From the global and local computations, it is possible to post-process the plate and 3D nodal reactions on the $\gamma_C/\Gamma_C$ interface:

$$ \begin{aligned} \lambda^C &= \left( \mathbf{K}_{|\omega_C}^C \mathbf{u}^G - \mathbf{f}_{\text{ext}|\omega_C}^C \right)_{|\gamma_C} \\ \lambda^L &= - \left( \mathbf{f}_{\text{int}|\Omega_I}^L (\mathbf{u}^L) + \mathbf{f}_{\text{ext}|\Omega_I}^L \right)_{|\Gamma_C} \end{aligned} \qquad (8) $$

Finally the balance of the nodal reactions is enforced on the static plate quantities, thanks to the transpose operator $\mathbf{R}^T$ which enables to compute the generalized forces associated to the 3D nodal reactions $\lambda^L$:

$$ \lambda^C + \mathbf{R}^T \lambda^L = 0 \qquad (9) $$

In spite of the care taken to recover the 3D boundary conditions on $\Gamma^I$ from plate quantities, small edge effects may occur, in particular in the case of stratified plates where suppressing edges effects require high order theories which are not always available in finite element analysis software or which could difficultly be implement in a non-intrusive manner. The role of the buffer zone (the overlap between $\Omega^I$ and $\omega^C$) is then to dampen these edge effects and make the evaluation of the 3D stress reliable on $\Gamma^C$. In the example below, the size of the buffer zone is set to $b = 5$mm, which was sufficient to absorb the local
artificial effects. It corresponds to four elements between the two interfaces $\Gamma_C$ and $\Gamma_I$. The dimension $b$ is represented on FIG. 6, together with dimension $L$ which characterizes the size of the 3D domain of interest $\Omega_I$. $L$ corresponds to the size of the part of the 3D domain which is not strictly necessary to represent the bolt correctly but which was inserted as a way to keep the 3D/2D transition away from the zone dominated to by 3D effects.

<figure><img src="image_7.png" /><figcaption>Figure 6: Definition of the size of the local model and of the buffer zone.</figcaption></figure>

## 3 The non-intrusive iterative algorithm

The system (6,7,8,9) can be interpreted as finding the traction $\delta$ to be imposed to the global plate model on the inner interface $\gamma_C$ in order to generate a reaction $\lambda^C$ in balance with the reaction of the inner zone of interest submitted to the recovery of the plate displacement on its boundary $\Gamma_I$.

Starting from $\delta_0 = 0$, this can be achieved through the following iterations:

1. Run a global plate analysis with extra load $\delta$:

$$ \mathbf{u}_n^G = \mathbf{K}^{G-1} (\mathbf{f}_{ext}^G + \mathbf{\delta}) $$

2. Post-process the reaction on $\gamma_C$:

$$ \lambda_n^C = (\mathbf{K}^C \mathbf{u}_{n|\omega_C}^G - \mathbf{f}_{ext}^C)_{|\gamma_C} $$

3. Recover the 3D displacement on $\Gamma_I$:

$$ \mathbf{u}_{n|\Gamma_I}^L = \mathbf{R} \mathbf{u}_{n|\Gamma_I}^G + \mathbf{w}^G (\mathbf{u}^G)_n $$
4. Solve the local problem with imposed displacement on $\Gamma_I$:

$$ \left\{ \begin{array}{l} \mathbf{f}_{int}^{L}(\mathbf{u}_{n}^{L}) + \mathbf{f}_{ext}^{L} = 0 \\ \mathbf{u}_{n|\Gamma_I}^{L} \text{ given} \end{array} \right. $$

5. Post-process the local reaction of $\Gamma_C$:

$$ \lambda_{n}^{L} = - (\mathbf{f}_{int}^{L}(\mathbf{u}_{n}^{L}) + \mathbf{f}_{ext}^{L})_{|\Gamma_C}^{\tilde{I}} $$

6. Compute the residual on $\gamma_C$:

$$ \mathbf{r}_n = \lambda_n^C + \mathbf{R}^T \lambda_n^L $$

7. If residual is small enough then exit, else update the extra load:

$$ \delta_{n+1|\gamma_C} = \delta_{n|\gamma_C} - \mathbf{r}_n $$

and go back to 1.

One important point is that the global step 1 and the local step 4 can be processed with different software. For the test-case developed here, Code_Aster is used for the global plate problem, and COFAST3D is used for the local 3D contact problem.

In the end the local/global method is a fixed point algorithm similar to a modified Newton algorithm. The rate of convergence of this algorithm can be quite slow. However acceleration techniques can be applied:

* quasi-Newton acceleration like SR1 algorithm [25] or BFGS,

* dynamic relaxation like Aitken's algorithm ($\Delta^2$) [24],

* mixed boundary conditions [25],

* Krylov solvers (only in the linear case) [34].

Based on [34], it appears that, among those methods, SR1 quasi-Newton acceleration (5.3) technique leads to better performance while being simple to implement in a non-intrusive manner.

Note that running one global analysis (step 1) followed by a local reanalysis with given Dirichlet conditions (steps 3-4) without iterations corresponds to the industrialists' practice called submodeling (or sometimes structural zoom). Such an approach is "purely descending" in the sense that there is no feedback from the local computation towards the global scale.
# 4 Analysis of the iterative corrections at the global and local levels

In order to evaluate the convergence of the algorithm, the relative norm of the residual at iteration i is defined as:

$$r_i^{rel} = \frac{\|\lambda_i^C + \mathbf{R}^T \lambda_i^L\|_2}{\|\lambda_i^C\|_2} \qquad (10)$$

In this section, the corrections associated with the coupling along the iterations are analyzed for the fixed-point algorithm (i.e. without acceleration techniques) with convergence threshold set to $10^{-6}$ which is clearly enough for mechanical quantities to have converged.

A typical convergence curve is shown Fig. 7(a). This curve shows that the convergence is more or less independent of the level of nonlinearity in the bolt (weak for the first global time step, strong for the last one). This type of result is typical of nonlinear relocalization methods [28, 29, 30, 31]. Let us note that, as shown in Fig. 7(b), four global time steps are used for the plate computation, which appears to be sufficient, see Section 5. Anyhow, sub-stepping is used for the computation of the bolt; balance between the local and global model is only reached at the global time steps in Fig. 7(b). For more involved problems, a better control of the global time steps and of the local sub-stepping process could be obtained using error indicators on the time discretization like in [36].

<figure><img src="image_8.png" /><figcaption>Figure 7: Convergence of the fixed-point algorithm: (a) Relative norm of the residual (b) Evolution of one component of the global and local nodal reactions.</figcaption></figure>

*Remark.* With the proposed technique, for each global increment the first iteration corresponds to a classical submodeling approach: this enables us to easily measure the quality of a that approach, which is a question often raised by engineers. In this application, the level of the error of the submodeling approach is about 25%.
# 4.1 Global effect of the correction

## 4.1.1 Effect of the preload of the bolt

The global effect of the tightening of the bolt cannot be predicted using a plate theory because the associated forces are equal to zero. Nevertheless this tightening has a non-negligible global effect as can be seen in figure FIG. 8.

<figure><img src="image_9.png" /><figcaption>Figure 8: Effect of the tightening (preload): value of the $N_{xx}$ component of the plate stress tensor.</figcaption></figure>

## 4.1.2 Additional global effects of the bolt on the plate solution

After the preload of the bolt, the plate is loaded in tension, in four global time steps. The analysis of the number of global increments is presented in Section 5. The global solution is modified along the iterations to match with the 3D model of the bolt, as can be seen on FIG. 9 which shows the values of the transverse displacement in the mid-section of the plate for the initial plate solution and the corrected one.

One can notice that the correction increases with the time steps, and it becomes significant for the third step and very large for the last one. This is of course due to the fact that for the first two time steps the bolt acts more or less as a 3D rigid connector. Whereas the sliding of the bolt becomes significant for the third time step. During the last stage of the loading loosening of the bolt happens, as can be checked from the local analysis of the solution presented in the next sub-section. Let us recall that the initial plate solution also corresponds to what would be obtained by a submodeling approach since submodeling only improves the local reanalyzed area.

# 4.2 Local effect of the correction

In this section the solution within the bolt is compared for three different approaches: a reference full 3D simulation, a submodeling approach and the mixed
<figure><img src="image_10.png" /><figcaption>Figure 9: Deformation of the plates, before and after iterations.</figcaption></figure>

2D-3D model with $L = 0$ (minimal size of the local 3D model) obtained at the convergence of the iterations (convergence threshold is $10^{-6}$).

For example the figure Fig. 10 shows the comparison of the $\sigma_{xz}$ stress field between the submodeling solution, the hybrid solution and the reference one. It can be seen that the solution associated with the non-intrusive approach matches very well the 3D reference solution while the submodeling gives inaccurate results inside the screw where the shear stress is greatly underestimated.

<figure><img src="image_11.png" /><figcaption>Figure 10: Transverse shear stress distribution: reference (top), converged 2D-3D solution (middle), submodeling (bottom).</figcaption></figure>

We now compare global quantities associated to the bolt (like the amount of dissipation due to sliding or the global axial component of the force acting on the joint) as well as local quantities (like the pointwise values of the sliding) for the three models. On the figures Fig. 11(a) and Fig. 11(b), local stress quantities are extracted on a point defined on Fig. 11(d). It appears that the relative residual provides an efficient indicator for the errors on global quantities. Typically the submodeling leads to 20% of error on the dissipation. Note that
the error committed by submodeling on local quantities can be much larger. In comparison, the converged 2D-3D solution and the full 3D solution are quite close.

<figure><img src="image_12.png" /><figcaption>Figure 11: Comparison between full 3D, submodeling and 2D-3D approaches.</figcaption></figure>

Figure FIG. 12 presents the sliding obtained by the converged 2D-3D model; it matches the 3D reference (few percents of deviation on the maximal sliding). On the contrary, the submodeling approach underestimates the sliding especially for the third and fourth global time steps, as can be seen on figure FIG. 13 where the error on the maximal sliding is more than 30%.

To conclude, the coupling approach seems reliable contrarily to the submod-
eling approach for which the level of error is 20 – 25% on global quantities and
can be much more for local quantities of interest. The lack of conservatism of
the submodeling extends similar results obtained in previous studies on localized
plasticity or buckling for example.
<figure><img src="image_13.png" /><figcaption>Figure 12: Sliding at each global time step of the converged solution (mm).</figcaption></figure>

# 5 Control of some parameters of the method and acceleration technique

As already discussed several parameters can be tuned for the non-intrusive 2D-3D coupling strategy :

* the number of global time steps,

* the choice of the position of the interface in the global model (dimension L in FIG. 6),

* the width of the buffer zone (dimension *b* in FIG. 6).

The size of the buffer zone has been chosen from our experience on the 2D-3D coupling involving composite plates and orthotropic plies [21], in order to minimize the problem of artificial edge effects between the two models. Besides the parameters of the nonlinear solver used for the bolt computation have been tuned in order to ensure a high precision according to dedicated papers like [22].

## 5.1 Influence of the number of global time steps

The number of time steps has been simply chosen on the basis of the comparison of solutions, starting from one global time step and increasing this number. On the figure 14, a comparison between 1, 2, and 4 time steps is shown. The solutions with 2 and 4 time steps are considered here as quite close, whereas the one using only one time step deviates largely in the middle of the load
<figure><img src="image_14.png" /><figcaption>Figure 13: Maps of the difference of sliding between converged and submodeling solutions at each global time step (mm).</figcaption></figure>

sequence. Let us note that the use of large time steps at the global level is made possible thanks to the chosen algorithm with nonlinear localization. Indeed the possibility to use sub-stepping at the local scale greatly reduces the potential difficulties of convergence.

<figure><img src="image_15.png" /><figcaption>Figure 14: Comparison of solutions for different number of global time steps.</figcaption></figure>
## 5.2 Influence of the position of the interface between the 2D and the 3D models

The question of the position of the interface raises in fact the question of the validity of the plate theory with respect to the 3D theory and has been largely discussed in [21]. From what is known on the validity of the plate and shell theories, in the case of isotropic materials, one expects the 2D-3D model to be a good approximation of the 3D reference, for an interface situated from the bolt at a distance superior to the thickness of the plate. It is therefore interesting to analyze the influence of the position of the interface and to compare the cases of $L = 0$ mm and $L = 15$ mm (see FIG. 6).

<figure><img src="image_16.png" /><figcaption>Figure 15: Comparison of the results for $L = 0$ mm and $L = 15$ mm</figcaption></figure>

From a global point of view, the final deformed shapes of the corrected solutions in FIG. 15(c) are very close in the common plate domain. They only slightly differ in the local area of interest which are different for the two models.

From a local point of view, on figure FIG. 15(b), both cases give close results. This is confirmed when analyzing the difference of the shear stress within the
bolt $\Delta\sigma_{xz} = \sigma_{xz}^{conv} - \sigma_{xz}^{0}$, as shown in the FIG. 16. As already analyzed, if for the first two steps, most of the differences are localized around the nut, during the sliding phase the correction mostly concerns the nut itself. The FIG.17 shows the evolution of the local tangential jump with respect to the prescribed displacement. This is an interesting quantity in order to observe the initiation of the sliding. The 2D-3D models give close predictions for both values of L. They are much more closer to the reference than what is predicted by the simple submodeling.

On the figure FIG. 15(a), it appears, as expected, that a larger number of iterations has to be carried out when the interface is located on the bolt boundaries ($L = 0$). This effect of the interface location on the convergence rate is corrected when using a SR1 acceleration technique as can be seen in the next subsection.

<figure><img src="image_17.png" /><figcaption>Figure 16: Error in the local domain for two positions of the interface and at each time step.</figcaption></figure>

## 5.3 Quasi-Newton acceleration

A final component of the non-intrusive substitution is the implementation of acceleration techniques. In that context, the SR1 quasi-Newton acceleration was tested in [25], and the Aitken Delta-2 dynamic relaxation method in [17]. For this study, we only implement SR1 acceleration.

The FIG. 18 presents convergence plots. The convergence is roughly 3 times
<figure><img src="image_18.png" /><figcaption>Figure 17: Sliding between the two plates near the nut for different configurations.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 18: Effect of the SR1 acceleration on the global-local rate of convergence.</figcaption></figure>

faster with SR1. Moreover the rate of convergence is almost independent from the load step and the position of the interface. This property can be explained by the fact that the corrections induced by the SR1 acceleration technique are adapted to take into account the main differences between the 2D and 3D models.

# 6 Conclusions

In this paper, a non-intrusive coupling between plate models and 3D models [21] has been extended to deal with the precise computations of bolted plates: the plate model with simplified connector was coupled with a full 3D nonlinear model of the bolt. The flexibility of the method was exploited to easily define
the coupled model and to use of two dedicated pieces of software: Code_Aster for the plate computations, and COFAST3D for the nonlinear computation of the bolt (including many surfaces of friction).

The proposed technique enabled us to analyze the possibilities and limits of the use of plate connectors and submodeling approach in that case. It appears that, even when the bolt response is globally linear, the 3D effects induced by the bolt largely modify the plate solution itself. This shows that a rigid connector is a poor model to describe the connection between two plates. As expected, when important sliding occurs, such modeling becomes irrelevant. Moreover, the submodeling technique may lead to significant local errors and non-conservative results.

Another important feature for future applications concerning the treatment of multiple bolts in interaction, is that using an interface located at the limit of the bolt leads to acceptable results when compared to the reference solution. Works on that type of problem is in progress.

Another issue concerns the modeling of the bolt itself. In practice, some of the parameters of the bolted assembly are not precisely determined, as the preload of the nut or the friction coefficient. Such types of problems have been analyzed by dedicated techniques for the bolt computation [37, 38], including multiresolution [39]. The use of the proposed non-intrusive techniques allows to extend these studies to the case of 2D-3D structural analyses in a straightforward manner. In addition, the use of model reduction techniques for the global model itself, as proposed in [40] should lead to a very significant reduction of the computational time.

## Competing interests

The authors declare that they have no competing interests.

## Acknowledgements

This work was partially funded by the French National Research Agency as part of project ICARE (ANR-12-MONU-0002-04).

## References

[1] Feyel, F., Chaboche, J.-L.: FE² multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials. Computer methods in applied mechanics and engineering **183**(3), 309–330 (2000)

[2] Wyart, E., Coulon, D., Duflot, M., Pardoen, T., Remacle, J.-F., Lani, F.: A substructured FE-shell/XFEM-3D method for crack analysis in thin-walled
structures. International Journal for Numerical Methods in Engineering **72**(7), 757–779 (2007)

[3] Amini, A.M., Dureisseix, D., Cartraud, P.: Multi-scale domain decomposition method for large-scale structural analysis with a zooming technique: Application to plate assembly. International Journal for Numerical Methods in Engineering **79**(4), 417–443 (2009)

[4] Ladevèze, P., Dureisseix, D.: A micro / macro approach for parallel computing of heterogeneous structures. International Journal for Computational Civil and Structural Engineering **1**, 18–28 (2000)

[5] Ladevèze, P., Loiseau, O., Dureisseix, D.: A micro-macro and parallel computational strategy for highly heterogeneous structures. International Journal For Numerical Methods In Engineering **52**, 121–138 (2001)

[6] Kerfriden, P., Allix, O., Gosselet, P.: A three-scale domain decomposition method for the 3D analysis of debonding in laminates. Computational Mechanics **44**(3), 343–362 (2009)

[7] Saavedra, K., Allix, O., Gosselet, P.: On a multiscale strategy and its optimization for combined delamination and buckling simulation. International Journal for Numerical Methods in Engineering **91**, 772–798 (2012)

[8] Dhia, H.B.: Problèmes mécaniques multi-échelles: la méthode arlequin. Comptes Rendus de l'Académie des Sciences-Series IIB-Mechanics-Physics-Astronomy **326**(12), 899–904 (1998)

[9] Dhia, H.B., Rateau, G.: The arlequin method as a flexible engineering design tool. International journal for numerical methods in engineering **62**(11), 1442–1462 (2005)

[10] Gmür, T.C., Kauten, R.H.: Three-dimensional solid-to-beam transition elements for structural dynamics analysis. International journal for numerical methods in engineering **36**(9), 1429–1444 (1993)

[11] Garusi, E., Tralli, A.: A hybrid stress-assumed transition element for solid-to-beam and plate-to-beam connections. Computers & Structures **80**(2), 105–115 (2002)

[12] McCune, R.W., Armstrong, C.G., Robinson, D.J.: Mixed-dimensional coupling in finite element models. International Journal for Numerical Methods in Engineering **49**(6), 725–750 (2000)

[13] Donaghy, R.J., Cune, W.M., Bridgett, S.J., Armstrong, C.G., Robinson, D.J., McKeag, R.M., Robinson, D., M., R.: Dimensional Reduction of Analysis Models, pp. 307–320 (1996)

[14] Nguyen, V.P., Kerfriden, P., Brino, M., Bordas, S.P., Bonisoli, E.: Nitsche's method for two and three dimensional nurbs patch coupling. Computational Mechanics **53**(6), 1163–1182 (2014)
[15] Gendre, L., Allix, O., Gosselet, P., Comte, F.: Non-intrusive and exact global/local techniques for structural problems with local plasticity. Computational Mechanics **44**(2), 233–245 (2009)

[16] Plews, J., Duarte, C., Eason, T.: An improved nonintrusive global–local approach for sharp thermal gradients in a standard fea platform. International Journal for Numerical Methods in Engineering **91**(4), 426–449 (2012)

[17] Duval, M., Passieux, J.-C., Salaün, M., Guinard, S.: Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition. Archives of Computational Methods in Engineering, 1–22 (2014)

[18] Chevreuil, M., Nouy, A., Safatly, E.: A multiscale method with patch for the solution of stochastic partial differential equations with localized uncertainties. Computer Methods in Applied Mechanics and Engineering **255**(0), 255–274 (2013)

[19] Bettinotti, O., Allix, O., Malherbe, B.: A coupling strategy for adaptive local refinement in space and time with a fixed global model in explicit dynamics. Computational Mechanics, 1–14 (2013)

[20] Bettinotti, O., Allix, O., Perego, U., Oancea, V., Malherbe, B.: A fast weakly intrusive multiscale method in explicit dynamics. International Journal for Numerical Methods in Engineering **8**, 577–94 (2014)

[21] Guguin, G., Allix, O., Gosselet, P., Guinard, S.: Nonintrusive coupling of 3d and 2d laminated composite models based on finite element 3d recovery. International Journal for Numerical Methods in Engineering **98**(5), 324–343 (2014)

[22] Champaney, L., Cognard, J.Y., Ladèveze, P.: Modular analysis of assemblies of three-dimensional structures with unilateral contact conditions. Computers and Structures **73**(1–5), 249–266 (1999)

[23] Daghia, F., Ladevèze, P.: A micro–meso computational strategy for the prediction of the damage and failure of laminates. Composite structures **94**(12), 3644–3653 (2012)

[24] Liu, Y., Sun, Q., Fan, X.: A non-intrusive global/local algorithm with non-matching interface: Derivation and numerical validation. Computer Methods in Applied Mechanics and Engineering **277**, 81–103 (2014)

[25] Gendre, L., Allix, O., Gosselet, P.: A two-scale approximation of the schur complement and its use for non-intrusive coupling. International Journal for Numerical Methods in Engineering **87**(9), 889–905 (2011)

[26] Klawonn, A., Lanser, M., Rheinbach, O.: Nonlinear FETI-DP and BDDC methods. SIAM J. Sci. Comput. **36**(2), 737–765 (2014)
[27] Negrello, C., Gosselet, P., Rey, C., Pebrel, J.: Substructured formulations of nonlinear structure problems — influence of the interface condition. International Journal for Numerical Methods in Engineering online (2016)

[28] Cresta, P., Allix, O., Rey, C., Guinard, S.: Nonlinear localization strategies for domain decomposition methods: Application to post-buckling analyses. Computer Methods in Applied Mechanics and Engineering 196, 1436–1446 (2007)

[29] Hinojosa, J., Allix, O., P-A.Guidault, Cresta, P.: Domain decomposition methods with nonlinear localization for the buckling and post-buckling analyses of large structures. Advances in Engineering Software 70, 13–24 (2014)

[30] Bordeu, F., Boucard, P.A., Gosselet, P.: Balancing domain decomposition with nonlinear relocalization: Parallel implementation for laminates. Proceedings of the 1st International Conference on Parallel, Distributed and Grid Computing for Engineering 90, 46–57 (2009)

[31] Allix, O., Kerfriden, P., Gosselet, P.: A relocalization technique for the multiscale computation of delamination in composite structures. Computer Modeling in Engineering and Sciences 55(3), 271–292 (2010)

[32] Düster, A., Scholz, D., Rank, E.: pq-adaptive solid finite elements for three-dimensional plates and shells. Computer Methods in Applied Mechanics and Engineering 197(1), 243–254 (2007)

[33] Düster, A., Niggl, A., Rank, E.: Applying the hp-d version of the FEM to locally enhance dimensionally reduced models. Computer Methods in Applied Mechanics and Engineering 196(37-40), 3524–3533 (2007)

[34] Guguin, G.: Stratégie non-intrusive de couplage plaque/3D pour la simulation des assemblages de plaques composites stratifiées. PhD thesis, ENS de Cachan, 61 av. Président Wilson, 94235 Cachan, France (2014)

[35] Ladevèze, P., Simmonds, J.G.: Nonlinear Computational Structural Mechanics: New Approaches and Non-incremental Methods of Calculation. Springer, Berlin (1999)

[36] Allix, O., Kerfriden, P., Gosselet, P.: On the control of the load increments for a proper description of multiple delamination in a domain decomposition framework. International Journal for Numerical Methods in Engineering 83, 1518–1540 (2010)

[37] Gant, F., Champaney, L., Rouch, P.: Modeling of the bolted joint behavior variability with the lack of knowledge theory. In: ICCES 2010 International Conference on Computational and Experimental Engineering and Sciences, Las Vegas, USA, vol. 28, pp. 97–98 (2010)
[38] Gant, F., Rouch, P., Champaney, L.: Updating of uncertain joint models using the lack-of-knowledge theory. Computers and Structures **128**(0), 128–135 (2013)

[39] Allix, O., Vidal, P.: A new multi-solution approach suitable for structural identification problem. CMAME **191**, 2727–58 (2002)

[40] Kerfriden, P., Passieux, J.-C., Bordas, S.P.-A.: Local/global model order reduction strategy for the simulation of quasi-brittle fracture. International Journal for Numerical Methods in Engineering **89**(2), 154–179 (2012)