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

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.


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.

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 Ω X possesses a plate counterpart ω X , so does any interface Γ Y whose plate counterpart is written γ Y . 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.

The plate model
The global plate model is the assembly of two plates, ω inf and ω sup respectively the lower plate and the top plate, and a rigid connector between them ω conn Fig. 1 (left). The plates are 20 mm thick and 280 mm long. The lower plate is 160 mm wide whereas and 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 ν = 0.3; the handling of orthotropic composite plates is explained in [21].
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: where m ∈ ω inf ∪ ω 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 (N tension, M bending, Q shear): 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].

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 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  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).

Connections between the models
The plates+connector model describes the entirety of the structure and occupies the (2D+1D) domain ω. The 3D model of the bolt occupies what we call the zone of interest Ω I . In the original method [15], the coarse modeling of ω I = ω ∩ Ω I was simply replaced by Ω 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 Ω I ⊂⊂ Ω I , the plate model is the only one taken into account in the outer complement zone ω C ⊂⊂ ω C , there exists an overlap Ω I ∩ ω C also called buffer zone where the two models are equivalent in a certain sense.
Two interfaces are defined: γ I = ∂Ω I ∩ ω C and γ C = ∂ω C ∩ Ω I . In order to simplify the example, regular and ruled quadrilateral meshes are used along the interfaces. Thus, γ I and γ C are straight lines between two conforming meshes, see  The substituted solution is expressed with the following form: and it satisfies the following transmission conditions (in a sense which is made more precise afterwards): u L and u G are continous on Γ I 3D stress and plate generalized forces are balanced on γ C To be more specific, let us assume that a finite element discretization is used: u G is the vector of plate unknowns defined on ω, u L is the vector of 3D unknowns defined on Ω I , f ext stands for the generalized forces (they include the effect of imposed displacements), K for the stiffness matrices, and 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 δ which is only nonzero on the degrees of freedom of the interface γ C (the values of δ are determined by the coupling algorithm). The global plate model equilibrium thus can be written as: 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: Matrix R represents the classical 3D rigid section movements. It enables to express the 3D Reissner-Mindlin kinematic from the plate degrees of freedom (u G ).
Because the rigid section kinematic is a very coarse hypothesis, we supplement it with warping displacements w G (u G ), proportional to the plate stress state.
To be more precise, the applied displacement on Γ I has the following form: with (F C i ) the eight generalized forces extracted along the interface γ I (3). (µ 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 γ C /Γ C interface: Finally the balance of the nodal reactions is enforced on the static plate quantities, thanks to the transpose operator R T which enables to compute the generalized forces associated to the 3D nodal reactions λ L : In spite of the care taken to recover the 3D boundary conditions on Γ 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 Ω I and ω C ) is then to dampen these edge effects and make the evaluation of the 3D stress reliable on Γ C . In the example below, the size of the buffer zone is set to b = 5mm, which was sufficient to absorb the local artificial effects. It corresponds to four elements between the two interfaces Γ C and Γ I . The dimension b is represented on Fig. 6, together with dimension L which characterizes the size of the 3D domain of interest Ω 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.

The non-intrusive iterative algorithm
The system (6,7,8,9) can be interpreted as finding the traction δ to be imposed to the global plate model on the inner interface γ C in order to generate a reaction λ C in balance with the reaction of the inner zone of interest submitted to the recovery of the plate displacement on its boundary Γ I .
Starting from δ 0 = 0, this can be achieved through the following iterations: 1. Run a global plate analysis with extra load δ: Post-process the reaction on γ C : Recover the 3D displacement on Γ I : Solve the local problem with imposed displacement on Γ I : Post-process the local reaction of Γ C : 6. Compute the residual on γ C : If residual is small enough then exit, else update the extra load: 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 (∆ 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.

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: 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]. 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%.

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.

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.

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 Step 1, initial plate solution Step 1, corrected plate solution Step 2, Step 2, corrected plate solution Step 3, Step 3, corrected plate solution Step 4, Step 4, corrected plate solution initial plate solution initial plate solution initial plate solution position x (mm) for y=0, z= ± 10 mm Figure 9: Deformation of the plates, before and after iterations.
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 σ 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. 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 11: Comparison between full 3D, submodeling and 2D-3D approaches. 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 submodeling 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.

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].

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).  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 ∆σ xz = σ conv xz − σ 0 xz , 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.

L=0 L=15
Step 1 Step 2 Step 3 Step 4 ∆ σ xz ( M P a ) Figure 16: Error in the local domain for two positions of the interface and at each time step.

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   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.

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.