Space/time global/local noninvasive coupling strategy: Application to viscoplastic structures

TL;DR

This paper extends a non-invasive global/local coupling method to large, nonlinear viscoplastic structures, proposing strategies to handle different time grids for improved accuracy and efficiency, demonstrated on 2D and 3D industrial examples.

Contribution

It introduces and compares two non-invasive coupling strategies for viscoplastic models with different time grids, enhancing simulation precision and performance.

Findings

Strategies improve accuracy over monolithic approaches

Methods are effective on 2D and 3D complex structures

Different time grids are essential for legacy code compatibility

Abstract

The purpose of this paper is to extend the non-invasive global/local iterative coupling technique [15] to the case of large structures undergoing nonlinear time-dependent evolutions at all scales. It appears that, due to the use of legacy codes, the use of different time grids at the global and local levels is mandatory in order to reach a satisfying level of precision. In this paper two strategies are proposed and compared for elastoviscoplastic models. The questions of the precision and performance of those schemes with respect to a monolithic approach is addressed. The methods are first exposed on a 2D example and then applied on a 3D part of industrial complexity.

Tables

E [MPa]	$\nu$	C [MPa]	$D$	R [MPa]	$n_f$	$K_f$	$n_s$	$K_s$
154,000	0.28	615,000	1,870	80	14	630	17.2	1,300

Strain rates $\left[s^{-1}\right]$	Abaqus default cutbacks	$\Delta p_{\max }={10}^{-3}$	$\Delta p_{\max }={10}^{-4}$
${10}^{-3}$	8.81$\%$ (2)	1.84$\%$ (9)	0.20 $\%$ (65)
${10}^{-5}$	30.40$\%$ (2)	12.36$\%$ (6)	1.80$\%$ (36)
${10}^{-8}$	30.46$\%$ (2)	30.46 $\%$ (2)	3.28$\%$ (17)

Table 3. Table 1: Influence of time integration on the computation (2D monolithic model)

Integration technique	Prediscr. #steps	Final #steps	Total time
Abaqus default cutbacks	8	12	1m 11s
$\Delta p_{\max }={10}^{-3}$	12	27	1m 33s
$\Delta p_{\max }={10}^{-4}$	44	173	7m 11s
$\Delta p_{\max }={10}^{-5}$	44	1501	1h 1m 29s

Table 4. Table 2: Error levels on the most loaded element at the end of the cycle (2D monolithic model)

	$\Delta p_{\max }={10}^{-5}$	Abaqus default cutbacks	$\Delta p_{\max }={10}^{-3}$	$\Delta p_{\max }={10}^{-4}$
	Ref. Values	Error $[\%]$	Error $[\%]$	Error $[\%]$
Mises	650.9 MPa	10.9	1.64	0.159
$p_f$	${\mathrm{1.05\hspace{0.17em}10}}^{-2}$	18.4	5.05	0.871
$p_s$	${\mathrm{1.27\hspace{0.17em}10}}^{-2}$	3.97	1.90	0.277

Table 5. Table 3: Time integration in the whole cycle according to the coupling approach (ATS: additional time steps)

Approaches	$\Delta p_{max}$	Global	Local	G/L	Total
		ATS	ATS	iterations	time
Submodeling	${10}^{-3}$	0	11	–	4min
	${10}^{-4}$	0	76	–	15.5min
Weak coupling	${10}^{-3}$	0	12	163	1h20min
	${10}^{-4}$	30	122	994	7h33min
Full coupling	${10}^{-3}$	0	13	344	1h53min
	${10}^{-4}$	0	381	5659	27h27min

Table 6. Table 4: True error on the most loaded element at the end of the cycle and CPU time for various coupling strategies

	Submodeling	Weak coupl.	Full coupl.	Reference
	error $[\%]$	error $[\%]$	error $[\%]$	$\left(\Delta p_{\max }={10}^{-5}\right)$
Mises	6.3	0.16	0.12	650.9 MPa
$p_f$	38	0.8	0.08	${\mathrm{1.05\hspace{0.33em}10}}^{-2}$
Total time	15min	7h33min	27h27min	1h11min

Equations

$$f=J_2\left({\sigma }_D-X_f\right)-R$$

$$\dot{\epsilon}^{p}=\dot{\epsilon}^{p}_{f}+\dot{\epsilon}^{p}_{s}\quad and\quad p_{i}=\int_{0}^{t}\sqrt{\frac{3}{2}\dot{\epsilon}^{p}_{i}:\dot{\epsilon}^{p}_{i}}{\color[rgb]{0,0,0}d\tau}\quad\text{ with for }i\in\{f,s\}$$

$${\dot{X}}_f=\frac{2}{3}C{\dot{ϵ}}_f^p-D{X}_f{\dot{p}}_f$$

$${\dot{p}}_f={⟨\frac{J_2({\sigma }_D-X_f)-R}{K_f}⟩}_{+}^{n_f}$$

$${\dot{p}}_s={⟨\frac{J_2({\sigma }_D)}{K_s}⟩}_{+}^{n_s}$$

$$\left(u^{GL},{\sigma }^{GL}\right)=\{\begin{array}{rl}\left(u^L,{\sigma }^L\right)& \text{in}{\Omega }_L\\ \left(u^G,{\sigma }^G\right)& \text{in}{\Omega }_C\end{array}$$

$$\text{Given }P,\text{ find }{u}^G=(u_C^G,u_{\Gamma }^G,u_A^G)\text{ such that }\left[\begin{array}{c}{f}_{int,C}^G(u_C^G,u_{\Gamma }^G)\\ f_{int,\Gamma }^G(u_C^G,u_{\Gamma }^G,u_A^G)\\ f_{int,A}^G(u_{\Gamma }^G,u_A^G)\end{array}\right]+\left[\begin{array}{c}{f}_{ext,C}^G\\ f_{ext,\Gamma }^G+P\\ f_{ext,A}^G\end{array}\right]=0$$

$${\lambda }_n^X$$

$$\text{Given }{u}_{\Gamma }^G,\text{ find }({\lambda }^A,u^A)\text{ such that }\left[\begin{array}{c}{f}_{int,\Gamma }^A(u^A;u_{\Gamma }^G)\\ f_{int,A}^A(u^A;u_{\Gamma }^G)\end{array}\right]+\left[\begin{array}{c}{f}_{ext,\Gamma }^A+{\lambda }^A\\ f_{ext,A}^A\end{array}\right]=0$$

$$\text{Given }{u}_{\Gamma }^G,\text{ find }({\lambda }^L,u^L)\text{ such that }\left[\begin{array}{c}{f}_{int,\Gamma }^L(u^L;u_{\Gamma }^G)\\ f_{int,L}^L(u^L;u_{\Gamma }^G)\end{array}\right]+\left[\begin{array}{c}{f}_{ext,\Gamma }^L+{\lambda }^L\\ f_{ext,L}^L\end{array}\right]=0$$

$$\{\begin{array}{rl}T{u}_{\Gamma }^G-u_{\Gamma }^L& =0\\ {\lambda }^G+T^T{\lambda }^L& =0\end{array}$$

$$r_i=-(T^T{\lambda }_i^L+{\lambda }_i^C)=-(T^T{\lambda }_i^L-{\lambda }_i^A+P_i)$$

$$P_{i+1}=P_i+r_i={\lambda }_{\Gamma ,i}^A-T^T{\lambda }_i^L$$

Full text

The questions of the precision and performance of those schemes with respect to a monolithic approach is addressed. The methods are first exposed on a 2D example and then applied on a 3D part of industrial complexity.

Keywords: Noninvasive coupling strategies ; viscoplasticity ; multiple time stepping

This work was supported by Safran Aircraft Engines and is part of the MAIA-MM1 research program.

**This work was published in Finite Elements in Analysis and Design, volume 156 pages 1-12 2019

doi:10.1016/j.finel.2019.01.003**

Introduction

The simulation of large structures where complex nonlinear local phenomena may occur is still a scientific and an industrial challenge. One of the main difficulties comes from the difference of length scales between the global response of the structure and the localized phenomena.

To deal with it in an effective manner, concurrent multiscale methods have been developed. They are often based on domain decomposition techniques like, among others, the primal BDD method [27], the dual FETI method [12] or the mixed Latin method [22].

These methods are not currently available in legacy codes, even though different attempts have been made. This is due both to the heavy development associated with their implementation and to their lack of robustness with respect to the variety of situations legacy codes have to deal with.

When the small wavelength phenomena are localized in space, i.e. concerning a relatively reduced part of the whole body, multiscale computations may also be based on coupling techniques as for example the Arlequin method [9]. The implementation of such methods in a legacy code is also not straightforward, mainly because the creation of the coupling operators between the two models in the transition zone requires complex integration operations.

A method often used by engineers to tackle such problems is the submodeling approach [8]: after a global computation, structural zooms are applied on the local critical zones, with details represented exactly. An advantage of this approach is that it allows to easily connect research software and commercial code, as was done for example to deal with the prediction of delamination under low velocity impacts [1]. Unfortunately, it implies neglecting the influence of the local zones on the whole structure. This may lead in turn to quite important local errors. The problem becomes more crucial when plasticity initiated locally spreads over the whole structure.

To correct the drawback of submodeling while keeping its simplicity and flexibility, a non-invasive method was proposed in order to allow exact local/global analysis embedding the same basic tools as those used in the submodeling inside an iterative procedure.

The method relies on interface coupling; its formulation and its numerical optimization have first been derived in the case of global linear models and local plasticity [15, 14]. A number of other applications and extensions have been proposed: use of XFEM at the local scale [29], the treatment of non-matching interfaces [25], coupling between a global plate model and 3D parts for bolted assemblies [18], geometrically non conforming coupling [19], multiscale time and space computation in explicit dynamic [3] with implementation in Abaqus Explicit for the analysis of delamination under impact [4], non-invasive domain decomposition approach [11]. Mesh refinement based on error estimation may also be cast in the proposed non-intrusive framework [10].

Alternative proposals exist, based on volume coupling. The volume coupling is performed either by means of a noninvasive version of the Partition of Unity method [30, 13], by using projection techniques between the local and global models [33, 20] or by means of homogenization like techniques [21].

So far the proposed method was applied for global linear models with local nonlinearities and topological changes. Safran Aircraft Engines was interested in exploring its potential for the detailed analysis of complex structures undergoing viscoplastic strains that can spread all over the structure. This situation is nowadays encountered in the case of aircraft engines. Indeed engines have undergone important improvements of their performance associated with the large increase of the working temperature and the use of optimized designs involving very thin parts. For example, micro-perforations were added in order to improve the cooling of hot parts (submitted to an air flux at $1500^{\circ}$C), like the combustion chambers and the high pressure turbine blades.

Our reference problem, see Section 5 for a full description, is an elastoviscoplastic turbine blade where two zones of interest need particular mesh refinement in order to correctly represent a complex geometry. Even if in that relatively simple case a monolithic model could be assembled and computed, our aim is to adapt to (or maybe extend) the industrial practice and only make use of a global representation of the blade (with simplified geometry in the zones of interest) and local refined representations of the zones of interest. In that case, even if they occupy less than 10% of the volume of the blade, the local meshes for the zones of interest have significantly more nodes than the global mesh (precisely three times more). We thus wish to derive an iterative coupling algorithm which converges to the reference (monolithic) solution by alternating elastoviscoplastic global and local computations.

The quality of the integration of viscoplastic models is very sensitive to the size of the time steps. When applying the non-intrusive framework with global spreading of the viscoplastic strains, several difficulties occurred using legacy codes. They appeared to be related to the management of the precision of the integration of internal variables on the different models. Indeed each model needed its own adapted time discretization, which is controlled by some automatic procedures, called cutbacks [5].

Cutbacks correspond to the subdiscretization of a given time increment when convergence difficulties or precision issues are encountered. As the Global/Local coupling relies on different models and independent computations, cutbacks are likely to generate different increment histories, i.e. different time discretizations, for the global and local models. It is thus of importance to control the precision of all computations and to provide sufficient synchronization of the adapted time grids in order to preserve the coherence of the coupling.

Based on this experience, the motivation of this paper is to propose and study different coupling strategies allowing to preserve the accuracy of the global/local coupling, making use of Abaqus’ most advanced precision controls.

The paper is organized as follow. In Section 1 the viscoplastic model is presented. A criteria for the definition of the time steps is proposed, based on the control of the increment of a sensitive variable. Then the reference model is presented in Section 2 with focus on the chosen procedure for the determination of the time discretization. In Section 3 the bases of the noninvasive global/local method are presented without addressing the problem of the coupling in time. In Section 4 two strategies for the global/local coupling in time are proposed. Those strategies are compared on the basis of a 2D example. A verification of the whole procedure on a realistic 3D industrial model involving non-matching interfaces is performed in Section 5.

1 About the elastoviscoplastic model and its sensitivity to time integration

1.1 Material model

The material model is the one proposed in [26], adapted from the Marquis-Chaboche’s behavior [6]. It is devoted to the description of a variety of phenomena which are characteristic of the elastoviscoplastic response under cycling loading. The elasticity itself is linear and isotropic, characterized by Young’s modulus $E$ and Poisson’s ratio $\nu$.

The nonlinear part of the model is ruled by the yield function based on von Mises criterion:

[TABLE]

where $J_{2}$ is the equivalent von Mises stress with $\sigma_{D}$ the deviatoric stress. The isotropic hardening is considered to be saturated at a value $R_{0}$, so that $R=R_{0}+\sigma_{y}$ where $R_{0}$ is an offset added to the elastic yield stress $\sigma_{y}$. Nevertheless, it is possible to use an exponential law in order to accurately represent the monotonic traction behavior or the cyclic stabilization.

$\epsilon^{p}$ is the plastic strain tensor and $p$ is the accumulated plasticity. These latter are split in a fast part $(f)$ and a slow part $(s)$:

[TABLE]

The kinematic hardening $X=X_{f}$ follows the classic Armstrong-Frederick’s formulation, only related to the fast potential:

[TABLE]

where $C$ and $D$ are two material parameters. The fast elastoviscoplastic and slow viscoplastic potentials follow Norton-Hoff’s laws:

[TABLE]

[TABLE]

Remark 1*.*

Even if the slow potential has no elastic threshold, it is a viscoplastic potential in its own right: it leads to plastic strains which remain after a loading/unloading cycle contrary to the case of a viscoelastic potential.

These constitutive relations allow for phenomena related to high temperature such as creep, strain rate effect, mean stress relaxation, and for all asymptotic behaviors to be represented.

The fast potential (4) dominates for strain rates in the range $\left[10^{-5},\;10^{-2}\right]s^{-1}$ whereas the slow one (5) governs in the range $\left[10^{-9},\;10^{-5}\right]s^{-1}$. The material of interest being confidential, we illustrate the type of constitutive response using the parameters given in [24] for a nickel based superalloy IN100 evolving at $800^{\circ}$C which is the mean temperature value during a flight. The corresponding values are reported in table 1.1.

In this study, the variation of temperature is taken into account only by making use of the adapted values of the material parameters. Nevertheless, some additional works were carried out with more complex laws involving thermal strains, without any trouble. For confidentiality reasons, we use classical values for the fast potential parameters, while the values for the parameters of the slow potential are chosen in order to balance the two plastic potentials for a strain rate of $10^{-5}s^{-1}$. The slow potential having no threshold, nonlinearity occurs as soon as the structure is loaded.

This law is incorporated in Abaqus Standard by the Zmat tool from Zset software [2].

1.2 Sensitivity to the time step: choice of a control parameter

The reference solver for this kind of nonlinear problem is the incremental Newton-Raphson method. Given a curve of load amplitude as Figure 2b, it is thus important to determine the load increment which ensures a good compromise between precision and CPU time. Engineering rules often permit to determine an efficient prediscretization. Anyhow it is convenient to adapt the next increment to the current properties of the computation, in particular if the solver has trouble to converge; this procedure is called cutback.

By default, Abaqus adapts the increments based on the monitoring of the speed of convergence, estimated by the evolution of the norm of the residual, compared to theoretical results and heuristics. The current increment being given, if divergence occurs then the computation is restarted with a 4 time smaller increment; if after a certain number of iterations the residual decreases too slowly, the computation is restarted with a 2 time smaller increment. If convergence is fast for two consecutive time steps, the next increment is enlarged by 50%. Abaqus also monitors the precision of the computation. In implicit dynamics analysis, the relevance of the integration is estimated by evaluating the residual in an interpolated mid-increment configuration. In quasi-static analysis, which is our case of interest, it proposes to adapt the time step in order not to allow too large evolution of internal variables. Note that the cutback procedure is highly customizable, but this process is highly dependent on the model especially for complex structures and it is difficult to set it up a priori with efficiency.

For our problem of interest, it appears that with default settings, cutbacks are only triggered out of convergence issues, and in fact accuracy is not sufficiently monitored. This is particularly problematic since it is crucial to have comparable levels of precision in the estimation of the inelastic evolution of several models.

It is thus critical to carefully tune the precision control of Abaqus which is also available through the Zmat tool of the Zset software [2]: any variable of the constitutive law can be used to master the time refinement based on a given increment of the chosen variable. After some trials, the fast accumulated plasticity was chosen as control variable since it appeared that small increments were correlated with high precision. In what follows, ${\Delta p_{\max}}$ denotes the associated plastic increment threshold: if a given increment leads to an increase of the accumulated plastic strain $\Delta p$ larger than ${\Delta p_{\max}}$ at any Gauss point then a cutback is applied: the computation is restarted with an increment of size reduced in proportion with ${\Delta p_{\max}}/\Delta p$. Please note that this process is totally automatic for the user, making it very simple to use and to calibrate on an elementary test as explained below. Such time steps added because of accuracy consideration will be called “additional time steps” (we reserve the name cutbacks to time steps added based on convergence consideration).

Sensible ranges for ${\Delta p_{\max}}$ can be determined on a tension test on a single cubic element for different strain rates, see Figure 1. In Table 1.2, the results for different thresholds are compared to an overkill solution obtained with the increment ${\Delta p_{\max}}=10^{-5}$. The results obtained using the default cutback procedure of Abaqus, starting from a single initial time increment are also displayed. Note that in practice for $\Delta p_{\max}\leqslant 10^{-4}$ no cutbacks were triggered out of convergence issues.

We observe that the solution obtained with ${\Delta p_{\max}}=10^{-4}$ is very close to the reference obtained with ${\Delta p_{\max}}=10^{-5}$, at a much lower computational cost (65 time steps instead of 500, in the case of fast loading). Note that, all errors are measured for a strain of $1\%$, which is unlikely to be reached in practice at low strain rates, so that the lower precision observed with low strain rates (3.3% for $\dot{\varepsilon}=10^{-8}s^{-1}$ vs 0.2% for $\dot{\varepsilon}=10^{-3}s^{-1}$) does not lead to loss of accuracy on structural applications.

2 Application to the reference 2D model

Before considering local/global coupling, we analyze how the increment adaptation procedure behaves on a reference 2D problem which is sufficiently small to be solved in a reasonable time using both monolithic and global/local approaches.

2.1 Reference monolithic model

The reference monolithic model is the target of the analysis; it is computed in order to evaluate the accuracy of the global/local method. It represents the structure with all the geometrical details, see Figure 2. The normal displacement of the foot of the part is prescribed to be zero, the air flux on the leading edge is modeled by a constant pressure and the rotation of the blade is taken into account by a centrifugal body force. The intensity of these external loads evolves in time according to the cycle of Figure 2 which mimics the main phases experienced by an engine during one flight. Note that all geometric nonlinearities are being neglected; in particular, the loads are evaluated in the initial configuration.

2.2 Time prediscretization over the cycle

Clearly, the definition of the load cycle by 8 time increments of Figure 2 is far from what is needed to ensure convergence and precision, or in other terms, it is far from the final time discretization resulting from cutbacks and additional time steps. In order to build a fair comparison with the global/local coupling, we propose to incorporate in the initial time discretization all the additional time steps that are deemed necessary for the resolution of the simplified model that will serve as the global model of the coupling with given $\Delta p_{\max}$ (see Subsection 3.1). In our case, the simplified model has the same behavior as the reference, but it bears no geometrical details, it is presented on Figure 5.

The time discretization resulting from this prior analysis is called prediscretization. It is then used for the reference monolithic solution (with the same $\Delta p_{\max}$).

Table 1 presents the number of time steps of the prediscretization and in the reference computation, for the Abaqus default cutback strategy and for different values of $\Delta p_{\max}$. It also gives the time spent for the reference computation. We observe that Abaqus default cutbacks introduce no new time step in the prediscretization and only 4 more for the reference computation. The control on plasticity increment leads to introducing more time steps in the prediscretization (up to 36 added time steps for $\Delta p_{\max}=10^{-4}$ or $10^{-5}$), the value of $\Delta p_{\max}$ is much more influent during the reference computation where, roughly, one order of magnitude gained in precision corresponds to the increase by one order of magnitude of the number of time steps (and of the CPU time).

Figure 3 presents the grids resulting from the prediscretization and from the reference computation. Not surprisingly, we observe that added time steps are mostly located in the phase of strongly increasing load, and to a lesser extent in the strongly unloading phase.

2.3 Accuracy of the computations

The results obtained for several thresholds are compared both in term of accuracy and cost. Due to the results of Section 1, ${\Delta p_{\max}}=10^{-5}$ serves to compute the reference solution and to measure errors. Two indicators are used, the von Mises stress and the cumulated plasticity.

As observed in Section 1, we see in Table 2 that the results obtained with ${\Delta p_{\max}}=10^{-4}$ are very close to the ones obtained with ${\Delta p_{\max}}=10^{-5}$. Table 1 shows that using ${\Delta p_{\max}}=10^{-4}$ leads to a much shorter computational time. A duration of one hour for such a simple case would correspond to very long computation time on relevant 3D industrial examples, as the one treated at the end of the paper. That is why we consider ${\Delta p_{\max}}=10^{-4}$ as a sensible compromise between accuracy and computational time.

Let us note that the error indicators are lower than the one obtained in the homogeneous case of Section 1, see Table 1.2. As said earlier, the strain rate heterogeneity (see Figure 4) leads to strain heterogeneity, so that error prone situations of large strain reached at low strain rate are never met in practice.

3 Non-invasive global/local method for local and global nonlinear models

In this section, the basic aspects of the global/local coupling are presented as if the local and global time grids were identical. The treatment of different time discretizations is discussed in Section 4.

The strategy makes use of three computational models based on the decomposition of the reference problem as explained in [17] in details for such nonlinear framework. The global model (index $G$) represents the whole structure and may bear a rough description of the zone of interest. The restriction of the global model on the zone of interest is what we call an auxiliary model (index $A$), the complement to the auxiliary model in the global model is written with the index $C$. The local model (index $L$) bears all the complexity of the zone of interest. Figure 5 illustrates the models used to set up the algorithm.

The global/local ($GL$) solution is defined on the reference model ($\Omega_{R}=\Omega_{L}\cup\Omega_{C}$) as the replacement of the global solution by the local one in the zone of interest:

[TABLE]

3.1 Global model

The global model matches the reference model on the complementary zone, that is to say everywhere except in the zone of interest where a coarser representation (named auxiliary model) can be used. The principle of the algorithm is to find an extra load $P$ to be applied on the interface $\Gamma$ of the global model between the complement and the auxiliary zones, such that the local and complement models are in balance at the interface. This is done iteratively by correcting $P$ knowing the difference between the nodal reactions of the complement and local models. For a current value of $P$ the finite element global problem can be written as:

[TABLE]

For the determination of the different nodal reactions on the interface some codes allow the shape function on one subdomain to be integrated, leading to the following formula for the nodal reaction $\lambda^{X}_{n}$ at one node $n\in\Gamma$ seen from domain $\Omega_{X}$ ($X\in\{C,A,L\}$):

[TABLE]

where $\sigma_{h}$ denotes the finite element approximation of $\sigma$, $f$ is the given body force and $F$ the given traction on $\partial_{F}\Omega_{X}$; $\phi_{n}$ is the shape function associated with node $n$.

When such a computation is not possible or too difficult (as it is the case using Abaqus), we make use of the auxiliary model. $\lambda^{A}$ is then post-processed from a computation on $\Omega^{A}$ with imposed Dirichlet conditions corresponding to the value of the current global displacement on $\Gamma$:

[TABLE]

3.2 Local model

The local model is the restriction of the reference on the zone of interest. The model is computed with the current global displacement prescribed on $\Gamma$, $u^{L}_{\Gamma}:=u^{G}_{\Gamma}$. The local system thus can be written as:

[TABLE]

The operation performed here consists in computing the local solution submitted to global displacements (Dirichlet condition) and then extracting corresponding reaction forces $\lambda^{L}$.

3.3 Non-matching interfaces

Handling non-matching interfaces provides flexibility for non-invasive approaches. It facilitates the meshing of the zone of interest (cf Figure 5). This procedure appears important especially for complex 3D meshes. Indeed two configurations can be encountered:

• •

The zone of interest is detected after a first global computation and then it is defined by a selection of global elements.

• •

The zone of interest can be defined at the level of the CAD. In that case, the local, the complement and the auxiliary models can be meshed independently so that the reference is defined with a TIE coupling [32] between the master complement and slave local models. The slave nodes are bound to the follow the kinematics of the interface of the complement domain.

In both cases the interface is well defined in all meshes and the coupling can be done on a surface, so that complex integration techniques [16, 31] are avoided.

In case of non conforming meshes the coupling equations take the following form:

[TABLE]

where $T$ is the (global to local) transfer matrix. In previous studies in the global/local framework, the transfer matrix $T$ was computed with techniques derived from the mortar method [11, 25]. Here a simpler solution is used: $T$ is the interpolation matrix which is processed using a Python script, before the computation.

3.4 Iterations

The formulation of the algorithm presented in [15] with a linear global problem is here extended to the case where both the global and local models are nonlinear. For one increment, we have:

Initialization. Set $i=0$, $P_{i}=0$. 2. 1.

Global computation. Computation with the interface load $P_{i}$, the interface displacements $u^{G}_{\Gamma,i}$ are deduced. 3. 2.

Local analysis. (Descent step) Computation with given external loads under prescribed interface displacements $Tu^{G}_{\Gamma,i}$. The reaction forces $\lambda^{L}_{i}$ are deduced. 4. 3.

Auxiliary analysis. Computation with given external loads and prescribed displacements $u^{G}_{\Gamma,i}$ applied on $\Gamma$. The reaction forces $\lambda^{A}_{i}$ are deduced. 5. 4.

Computation of the interface equilibrium residual.

[TABLE] 6. 5.

Update of $P$:

[TABLE]

Then $i\leftarrow(i+1)$ and goto 1.

The local and auxiliary analysis are computed in parallel thus the use of the auxiliary model does not add computational time to the general algorithm. The algorithm is a stationary iteration; its convergence can be proved under very general hypothesis (see [28] for a proof with weak hypothesis and [17] for the registration of the method amongst Schwarz alternating methods for which many convergence results exist). In the general case, relaxation may have to be used to ensure convergence by modifying (13) as follows: $P_{i+1}=P_{i}+\omega r_{i}$, with $\omega$ small enough. Relaxation is needed only if the auxiliary model is more compliant than the local model, which is not the case in the studied examples.

The relative norm of the residual $\|r_{i}\|_{2}/\|r_{0}\|_{2}$ is used to monitor the convergence of the method over the iterations. In practice, the tolerance to stop the iterations is set to $10^{-5}$. Figure 6 presents the evolution of the residual and of the true error compared to the monolithic computation, for different acceleration techniques. As usual Aitken’s $\delta^{2}$ leads to significant acceleration at almost null extra cost.

Remark 2*.*

The algorithm was presented for one increment of time. In order to compute a loading cycle involving many time steps, one has to reuse the current converged state to initialize the next computation. This is made using the *Restart function of Abaqus.

Remark 3*.*

Stagnation occurs when observing the true error (with respect to a monolithic computation), as can be seen on Figure 6b. This is problem very commonly encountered when using Abaqus/Standard. A probable cause is that our implementation makes use of output database (odb) files in order to extract the data necessary for the computation of the corrective load to be applied on the global model; and these database truncate nodal reaction to single precision reals in order to limit the size of the files.

4 Time coupling for the Global/Local method

The main issue, for the considered application, is that the global and local models require, in order for them to be properly integrated, different time steps that are not known in advance. The retained principle here is to make use of the automatic adaptive time stepping control method presented in Subsection 1.2). As explained in Subsection 2.2, we start from the prediscretization which correspond to the grid adapted to the global computation with given threshold $\Delta p_{\max}$. The remaining question is when to apply global/local coupling iterations, in other words at which time steps the models have to be coupled.

Let $(t^{G}_{0},t^{G}_{1},\ldots,t^{G}_{N})$ be the initial set of time steps resulting from the prediscretization. They define the initial set of time increments $\Delta t^{G}_{i}=t^{G}_{i+1}-t^{G}_{i}$.

4.1 Couplings strategies

In what follows two possible coupling strategies are presented. They are compared to the monolithic solution and to the sub-modeling procedure in terms of accuracy and cost.

• •

For the first strategy, called “weak time coupling”, the coupling is made only at the global additional time steps inserted during the local/global iterations.

• •

For the second strategy, called “full time coupling”, the coupling is ensured at each global and local additional time steps inserted during the iterations.

4.1.1 Weak time coupling

This approach corresponds to the coupling of the global and local models at all time steps of the global discretization, expecting that the additional time steps, required by the local model, are not of interest for the complementary area. A schematic view of the weak coupling is proposed on Figure 7.

• •

Initial cycle definition (blue points)

• •

Global prediscretization (green triangles $(t^{G}_{n})$)

• •

Additional global time steps (orange triangles) which are necessary because the global/local coupling introduces the extra load $P_{i}$ which increases the level of plasticity.

• •

Additional local time steps (red diamonds) are purely internal to the local computation.

The weak coupling is summed up in Algorithm 1. The principle is to try to enforce the coupling at a given time step $t_{a}=t_{c}+\Delta t_{c}$ where $t_{a}$ is the aimed time step, $t_{c}$ is current converged time step and $\Delta t_{c}$ is current value of the increment. If at any iteration of the global/local coupling a precision issue is triggered at the global level ($\Delta p^{G}>\Delta p_{\max}$ for some global computation) then the increment is reduced and $t_{a}$ redefined; the coupling is restarted. On the contrary, if a difficulty is encountered at the local level, then the local time steps are adapted until $t_{a}$ is reached; the coupling is not restarted.

We observed a difference, between loading and unloading phases, regarding the need to insert additional time steps (ATS) at the local scale, depending on the convergence of the coupling, this is illustrated on Figure 8. Indeed, during the loading phase, the global model tends to converge “from below” and the load transmitted to the local is underestimated in the first coupling iterations, this results in no local ATS in the first two iterations and up to 2 ATS in the following. On the contrary during the unloading phases, plasticity level is overestimated in the first iterations and requires temporarily up to four local ATS whereas once the coupling has almost converged, and $P_{i}$ is close to the solution, no local ATS are required.

4.1.2 Full time coupling

The full time coupling consists in performing global/local iterations for all the additional time steps required from both global and local models. It is expected to lead to high level of accuracy at a large computational cost. In Figure 7 the full time coupling would somehow correspond to backporting the local added time steps (red diamonds) to the global grid. Algorithm 2 recaps the full coupling.

For the full coupling, the unloading phases are particularly critical. Indeed, for the first iterations, the global model provides particularly inadequate boundary condition, resulting in the need of additional time steps for the local model which are no more needed when convergence improves, as evoked for the weak coupling in Figure 8b. In the full coupling, each additional local time step is transfered on the global time grid and the coupling iteration is restarted. This result in a much increased computational time without practical interest from a precision point of view.

4.2 Comparison of the different methods in 2D

In this section, submodeling and global/local couplings are compared.

Table 3 quantifies the additional time steps (ATS) and their impact on CPU time. Global and Local ATS specify which computation did not respect the $\Delta p_{\max}$ criterion and required the insertion of a time step. As explained earlier, the prediscretization was designed so that the global computation of the submodeling approach needs no additional time steps. In the full coupling, the common time discretization is always driven by the local model, which experiences the most severe plastic evolution. In the weak coupling, both global and local models trigger ATS on their own grid, but much less than the full coupling. Regarding the computational time, the full coupling is much too expensive for the objective precision of $\Delta p_{\max}=10^{-4}$.

Regarding the precision, according to the elementary tests of Subsection 1.2 and those presented in Subsection 2.2, the threshold ${\Delta p_{\max}}=10^{-4}$ is used. The couplings are assessed with respect to the reference monolithic solution integrated with ${\Delta p_{\max}}=10^{-5}$. Table 4 summarizes the results obtained for the different coupling methods. Of course the submodeling approach is attractive in term of CPU but it underestimates the level of plastic strains by $38\%$, whereas the coupling strategies maintain an error level below $1\%$. The weak coupling leads to what is considered in practice as an acceptable level of accuracy compared to the overkill solution. Figure 9 shows that, as expected, largest differences are located near the corner and near the most loaded holes.

4.3 Improvement of performance

In order to reduce the computational cost of the global/local analysis, the number of iterations performed on each time step may be reduced. As seen on Figure 6b, the accuracy is not improved after 5 iterations even if the equilibrium residual could still be decreased. The stagnation grossly corresponds to a relative residual of $10^{-3}$ which will be used (instead of $10^{-5}$) in the following computations. Furthermore, Aitken’s $\delta^{2}$ is applied. In this configuration, the total computational time is reduced by a factor 3 (2 hours and 44 minutes for the weak coupling strategy) without significantly altering the accuracy.

Remark 4* (Cost issue).*

Presently the performance of the method is strongly penalized by the fact that, for each resolution (global or local) a new independent computation of Abaqus needs to be launched. This implies three additional stages compared to the monolithic solution (apart from the first time step): verification of input data, building of the geometry, assembly of the problem. Those stages are very costly. Abaqus development team is considering evolutions to mitigate these costs.

5 Application on a 3D example of industrial complexity

5.1 Model definition

The method is now applied on a 3D model representing a high pressure turbine blade close to actual industrial problems. Figure 10 presents the meshes and gives the associated number of degrees of freedom.

Two local areas are introduced: one on the leading edge (LE) and another on the trailing edge (EE). In these two zooms, micro-perforations are added; they require meshes about ten times finer than the global one. Because the models remain of reasonable size ($5.63\,10^{6}$ degrees of freedom), a monolithic approach can be computed, with 15 processors and 70Gb of memory, and be used as a reference.

In 3D, the relative size of the interfaces, which characterizes the cost of exchanges, is quite small. On this example they correspond only to 2580 and 2520 nodes (cf. Figure 11) and the cost of the exchanges at the interfaces becomes very cheap compared to the cost of the resolution of the different models. The local and complementary areas present non matching meshes, so the method presented in Section 3.3 is applied. The reference model itself makes use the TIE coupling of Abaqus.

The material model, external loads and boundary conditions are similar to the 2D example: the blade is subjected to a pressure on the leading edge, centrifugal forces on all models and the normal component of the blade foot displacement is set to zero.

5.2 Numerical results

In order to evaluate the performance of the algorithm, the optimized settings are used and the number of coupling iterations is limited to 5. Figure 13 displays the time steps added at the local and global levels for the threshold of ${\Delta p_{\max}}=5\,10^{-4}$. The local model of the leading edge requires 62 additional times steps whereas the local model on the trailing edge does not demand any time refinement. A large plastic area, which connects all local models is obtained, see Figure 12.

The convergence rate plotted on Figure 13b is comparable to the one obtained on 2D test case. Indeed the problem might appear more complex from a domain decomposition point of view because of the complex shape of the interface, but this difficulty is largely compensated by the lesser redistribution of nonlinearity compared to the 2D case.

5.2.1 Accuracy aspects

A satisfying level of accuracy is obtained everywhere. Figure 14 reports error levels for the most loaded element. Once again the error of the submodeling approach is high, although smaller than in the 2D case. This is because the global to local redistributions are less important in this case.

5.2.2 Performance

The reference (monolithic) model is computed in 7 hours and 27 minutes and the submodeling provides a solution on the cycle in only 4 hours and 37 minutes. The weak coupling applied in its optimal configuration solves the model in 31 hours and 15 minutes, representing a ratio of 4 with the reference. As a reminder the ratio was about 24 on the 2D test case. But, as already discussed in Subsection 4.3, a fair comparison will be possible only when the issue of the model reconstruction at each time step will be solved.

Conclusion

The global/local non invasive method has been extended to allow a proper integration of the constitutive relation on each domain to be coupled. The use of different time grids over the iterations and depending on the domains seems mandatory. Two strategies have been compared and in practice the lighter one, called weak coupling, seems to offer a good compromise between cost and precision. Still, in a 3D case where a monolithic reference could be computed, the coupling strategy was 4 times slower. Let us recall here that one of the motivation is to allow the engineer to easily modify some local details without having to reconstruct the whole mesh. In fact the later operation is really cumbersome for complex industrial parts and is more and more often externalized. Moreover it will be possible to properly assess the performance of the method on Abaqus only when it will allow to construct the different models only once, a possibility offered by codes like Code_Aster for example. Departing from the context of this paper, it is also clear that using some guaranteed error estimation techniques would be a real asset to the method. It would allow to adapt the choice of the integration threshold ${\Delta p_{\max}}$ and the choice of the convergence criterion for the global/local iterations, for a given target accuracy. A huge literature devoted to error estimation is of interest here [7, 23], in particular regarding quantities of interest, but it should be ported to the context of non-invasive tools.

Acknowledgement

The authors are extremely grateful for the helpful and benevolent remarks given by the reviewers. Moreover the expertise of Reviewer 3 on Abaqus appeared to be a real asset.