Computation of the elementary stiffness matrices

Introduction

La prise en compte des mécanismes de plasticité, d’endommagement, . . . nécessite d’écrire des lois de comportement complexes fortement non linéaires. Coupler à la méthode des éléments finis, il doit être possible de décrire le comportement de structures complexes avec un degré de confiance optimal. Malheureusement, ces lois sont à l’origine de défaillances remettant en cause les résultats obtenus par simulation et faisant donc perdre tous les bénéfices issus d’une recherche fine en mécanique des matériaux. La remise en cause des résultats peut provenir de trois phénomènes :
– le premier aboutit à une absence de résultats du fait d’une divergence trop précoce du solveur ou à une remise en cause de la solution due à l’apparition de sauts de solution. L’algorithme de NewtonRaphson ou même les algorithmes de pilotage usuel (Riks 77), généralement utilisé pour résoudre les équations non linéaires issues de la discrétisation par la méthode des éléments finis, ne sont alors plus adaptés et ne pemettent pas de gérer les instabilités locales issues de l’apparition de l’endommagement ou de la plasticité.
– le second problème pouvant apparaître est une dépendance de la solution vis-à-vis des choix de modélisation (type, orientation et taille des éléments, choix des conditions aux limites, . . .). Cette défaillance peut être traduit d’un point de vue mathématique par une perte d’ellipticité des équations.
– enfin l’accroissement des temps de calcul rend les calculs trop longs voir irréalisables. En effet, l’obtention de calcul convergé nécessite généralement d’utiliser des maillages extrémement fin dans leszones ayant un comportement non élastique.
Ces trois défaillances ont fait l’objet de recherches durant les vingt dernières années. Malheureusement celles-ci apportent des solutions dans des cadres restreints nécessitant de faire des hypothèses difficilement généralisables dans le cas de strutures complexes (sur la direction de fissuration par exemple). Dans le cas des algorithmes de pilotage, Alfano (Alfano 03) a par exemple proposé de prendre en compte l’état local de la structure (dans le cas de l’utilisation d’éléments de zone cohésive nécessitant un placement a priori de ceux-ci et donc la définition du chemin de fissuration) pour piloter le chargement enpondérant l’influence des différents points de la structure (Alfano 03). Dans le cas de la dépendance aux maillages, une approche locale de la rupture a été développée. Celle-ci propose de caler les coefficients de la loi de comportement pour une taille de maille donnée et d’utiliser des maillages réglés et orientés dans la zone de fissuration ou d’endommagement (Besson 00; Besson 01; Rivalin 01). Le dernier point est plus complexe.
Ces lois de comportement peuvent être tout simplement ignorées au prix de forts surdimensionnements dus à l’utilisation de critères en post traitement de calculs quasi-élastiques. Il est aussi possible de simplifier les modélisations (les petits trous peuvent être négligés, par exemple) ou d’utiliser des résultats avec des maillages non convergés. Cela induit des risques d’erreurs souvent corrigés par le biais decampagnes expérimentales très onéreuses.
Dans cette première partie de ce manuscrit, nous allons nous attacher à proposer un schéma de résolution adapté à ces lois de comportement et ne nécessitant par d’hypothèse à priori trop forte. Cet algorithme est basé sur l’association de trois méthodes sophistiquées ayant un but précis. Ces méthodes sont :
– une modélisation non locale de l’endommagement (paragraphe 1.2). En effet, l’introduction d’un nouveau paramètre, proportionel à une longueur, permet de régulariser la dépendance aux maillages.
– un solveur à longueurs d’arc pilotant le chargement en fonction de l’état local de l’ensemble de la structure (paragraphe 1.3). Celui-ci permet de prévenir les sauts de solution et les divergences précoces tout en évitant les retours élastiques.
– une méthode de calcul par sous domaine adaptée aux deux méthodes précédentes (paragraphe 1.4) afin de réduire les temps de calcul.
Au final, l’utilisation de cette algorithme devra permettre d’effectuer des calculs « robustes » comme nous le verrons dans la dernière partie de ce chapitre (paragraphe 1.5). Il servira également de base aux autres travaux portant plus particulièrement sur les matériaux stratifiés.

Nonlocal models

The material mechanical behavior is characterized by so called “state variables”. The time evolution of these variables is usually (i.e. in the case of standard local models) expressed as a function of the variables for each material point, ~x. In the case of nonlocal models, the evolution of some of the variables at point ~x not only depends on the local state but also on the variables in the neighbourhood of ~x . In the following, one nonlocal variable, f nl , will be used which depends on a local variable f l . Based on the literature (Bazant 88; Bazant 94; Jirasek 94; Bazant 96; Jirasek 98), many local variables (equivalent strain, damage, porosity,. . .) can be chosen to define the nonlocal one (see section 1.5.5). In the pioneering work by Pijauder-Cabot (Pijaudier-Cabot 87), f nl at ~x was defined as a weighted average of f l around ~x . However this technique is difficult to implement in a finite element software. Indeed, it is necessary to have a specific algorithm (i) to built the connectivity table of the gauss points where f l is defined (They are usually defined locally in the finite element and are totally independent of their neighbors), (ii) to compute the nonlocal integral and it gradients (Baaser 03). That is why the integral relation is rewritten using gradient of f l (explicit formulation) or f nl (implicit formulation). Following the work by Peerlings (Peerlings 99) the implicit formulation should be preferred. This section describes the continuous implicit formulation and its finite element discretization.

Remarks about the local/global resolution algorithm

As said in the first and the present section, the global stiffness matrix of the problem is not symmetric. Moreover, the unknowns have not the same nature and order of amplitude. In order to overpass these difficulties and insure a good convergence of the computation, it is necessary to use specific resolution algorithms and preconditioners (see, among others, (Gosselet 03; Kincaid 03; Ayachour 03)). Because of the poor condition number of the system, an adimensionnalization of the local matrices (multiplication of Ks by the inverse of its diagonal) is mandatory.

Applications

Meshes and constitutive equations

The studied structure consists in a plate containing three holes as shown on fig. 1.5. A uniform vertical displacement is imposed and standard symmetry conditions are used. Four different meshes are used with 11, 13, 20 and 40 elements in the minimum cross section. Triangular elements are used. This leads for a standard local computations to a number of degrees of freedom varying between 30 000 and 120 000 and between 35 000 and 140 000 for the nonlocal computations.
The material is supposed to be isotropic, elastic damageable. Although simple and widely used, the simple formulation for the behavior (Saanouni 88; Ju 89) leads to the numerous numerical difficulties mentioned in the introduction. Damage.

Effects of element type and boundary conditions

Using a local damage model, results are mesh-size dependent. In this section the dependence on element type and on mesh design is also evidenced. Quadrangles (quadratic interpolation of the displacements and linear interpolation of the nonlocal damage) are used (Q8). Two meshes are used (fig. 1.13) : (i) regular radial mesh (Q8), (ii) irregular mesh with a finely meshed oriented zone (Q8-or). The global response is shown on fig. 1.12 ; results are compared with the local calculations carried out with triangles (T6) and the nonlocal calculations which depend neither on the element type nor on mesh design. It is important to outline that in the case of the oriented mesh, the crack path differs from the other cases (see fig. 1.13). In all local cases, the damage localization band width is equal to the element height. Rupture is quickest in the case of the oriented Q8 mesh because it has the smallest element width along the fracture path.
The fact that nonlocal calculations are both independent on mesh size and element type allows to use complex automatically generated meshes whereas local models require to use similar elements of fixed size along a predefined crack path (see e.g. (Rivalin 01)). Calculations have been carried out enforcing symmetry so that only 1/4 of the specimen was meshed. Calculations of the whole structure without symmetry have also been carried out. Both sets of calculations coincides in the nonlocal case ; different results are obtained for the local model as shown on fig. 1.14. The local full mesh calculation presents four distinct snap-backs which each corresponds to the failure of one ligament. Using 1/4 of the specimen, only two snap-backs are observed as only two ligaments are actually meshed. Symmetry with respect to the Oy axis is lost. In actual experiments, ligaments are expected to fail individually (i.e. four events). This could be modelled in the case of the nonlocal model using random material properties as the actual material is inhomogeneous (Besson 00). As in the local case, the localization band width corresponds to one element in the full mesh case and symmetry with respect to the Ox axis is clearly lost. Due to symmetry, the band width is actually twice as large in the other case (see fig. 1.15). 3D calculations using previous meshes extruded along the 0z axis with 1 or 3 elements in the thickness were also performed imposing null displacements along the z -direction on both sides of the mesh so that the simulation is therefore equivalent to a 2D plane strain simulation. Computation time is indeed strongly increased ; results in terms of global response and damage development are similar to those obtained in the 2D case.

Benefits of parallel computation

In previous sections, more than thirty computational results (in two or three dimensions, with complete or partial meshes), have been presented. The duration of one computation may vary, with a sequential algorithm, from one hour (for a two dimensional nonlocal model with the coarse partial mesh) to one month (for a three dimensional local model with the thin partial mesh) with a number of degrees of freedom between 30 000 and 350 000. Thus, this study could not be achieved in a non-parallel framework : the decomposition of the problem into N subproblems (fig. 1.17), enables to decrease significantly the computational cost.
In this problem, the local and nonlocal behaviors are explicit. So with or without damage the computational cost for (σ) or (σ, f l ) is always the same. Thus, it is sufficient to insure that each subdomain has the same number of elements to insure a good balance of the computational cost between subdomains4 .
Fig. 1.18 and 1.19 show the speed up (computational time for parallel simulations divided by the one for the sequential simulation) for the four methods in the case of the 5 decompositions 2D or 3D computations. The speed up of the computations are the same for local and nonlocal models. But it is important to notice that the nonlocal computations, even if they have more degrees of freedom, goes faster than the local ones : the local model induces more snap-back, which induces more increments and iterations to converge, and which increases the computational cost.
The evolution of the speed-up according to the number of domains is “always” the same. The evolution is divided in three regions. First the speed-up increases with the number of subdomains (region II), then it reaches its maximum and starts to decay (region III). Sometimes, when the size of the problem is very large, the first decompositions with a small number of subdomains increase the cost of computation due to the swap time (region I).
For the 2D-problems (fig. 1.18), the results are in the second and third regions of the previous description. The best results are obtained with the dual u -dualfnl and dualu -primalfnl methods, the little differences between the two are more likely to be induced by the computational environment (use of the network for example). The primal u -primalfnl and primal u -dualfnl method are twice more expensive : the convergence study shows that the global equilibrium requires more iterations to converge, the needed swap and computational time is then more important. In three dimensions (fig. 1.19), the results are in the first and second regions. The maximal efficiency of the method is never reached : the size of the problems is too important for the given hardwarearchitecture.
Finally, the 2D-results show that the best decomposition consists in having ten thousand degrees of freedom by subdomain for the dual u -dualfnl and dual u -primalfnl methods. The 3D-results also tends to thisvalue. For the other methods, the value is about twenty thousand degrees of freedom.

Conclusion

La prise en compte de l’endommagement dans les simulations, via le développement de lois non linéaires, doit en théorie permettre de décrire le comportement de structures complexes et de prédire leurs défaillances. Malheureusement, le couplage de ces lois à la méthode des éléments finis peut aboutir à une absence ou à une remise en cause des résultats obtenus par simulation. Il est possible d’identifier trois origines à ces défaillances : (i) les solveurs non linéaires (algorithme de Newton-Raphson ou à pilotage à longueurs d’arc globaux) ne sont plus adaptés et engendrent des sauts de solution ou une divergence précoce, (ii) les hypothèses de modélisation ne sont plus adaptées et engendrent une dépendance de la solutionobtenue vis à vis des choix de modélisation (taille et orientation des éléments par exemple), (iii) les coûts de modélisation ne permettent pas la création de maillages suffisament raffinés (taille de maille ou prise en compte précise de la géométrie de la pièce par exemple) pour obtenir des résultats satisfaisants avec des temps de calcul réalistes.

Discussion vis-à-vis des résultats expérimentaux

A l’aide de la campagne expérimentale et des modélisations associées, l’existence ou l’absence de  lien entre la longueur interne et le matériau, le modèle de comportement ou la structure sont étudiés en fonction des remarques précédentes. Dans toute cette analyse, seule la longueur interne associée aux ruptures de faisceaux de fibre est prise en compte. En effet, c’est elle qui influence le plus le comportement de la structure.
Lors des identifications, deux points importants sont apparus. Le premier est que la longueur interne est définie de manière unique pour l’ensemble des essais. Le second point concerne son ordre de grandeur à savoir 4. mm. Ces deux indications nous apportent certaines informations sur la nature de la longueur interne.
Par contre, ces indications doivent être nuancées du fait de la sensibilité relativement faible des résultats vis-à-vis de la longueur interne (en comparaison de celle de la déformation à rupture par exemple).
Lorsque la dimension de la longueur interne est comparée aux longueurs caractéristiques du T700/M21 aux différentes échelles, il est possible de remarquer que les 4. mm obtenus sont très grands devant les longueurs mesurables au niveau microscopique, mésoscopique. Par conséquent, la longueur interne ne semble pas être une grandeur matériau définissable à ces échelles.
L’unicité de la solution laisse à penser que la longueur interne n’est ni un artéfact purement numérique, ni reliable à la structure. En effet, une seule valeur est obtenue pour trois géométries différentes. De même, les types de stratifications n’ont pas donné lieu à des valeurs différentes donc une dépendance vis à-vis de l’empilement est difficilement imaginable ce qui rejette l’hypothèse d’une dépendance matériau à l’échelle macroscopique.
Au final, il ne reste plus qu’une seule solution envisageable : “la longueur interne est associable au modèle de comportement non local”. Cette définition est bien entendue très discutable. En effet, il n’est pas possible d’affirmer avec certitude une telle dépendance avec seulement six essais finalement assez semblables (même type de singularité géométrique, même matériau) et une sensibilité si “faible” des résultats de simulation par éléments finis vis-à-vis de ce nouveau paramètre. Pour aller plus loin dans cette analyse, il serait nécessaire d’envisager d’autres cas en choisissant un autre type de pli unidirectionnel, en enrichissant ou en changeant la forme de la loi de comportement, en modifiant le choix de la variable délocalisée, en changeant le type de singularité ou d’essais (essais de propagation), . . .

Conclusion

Dans ce dernier chapitre, une première confrontation expérimentale a été réalisée afin d’étudier la  validité du modèle d’endommagement non local spécifiquement développé afin de décrire le comportement de structure en composite à matrice organique.
Ce travail a consisté dans un premier temps à définir et à réaliser une campagne expérimentale afin d’acquérir des informations suffisantes sur le matériau pour identifier une loi de comportement prenant en compte les ruptures de fibres, étudier la robustesse du schéma de résolution proposé dans les trois premiers chapitres du document et étudier l’influence de la structure sur la longueur interne.
Par la suite, une loi de comportement a été développée et identifiée à partir des informations pré existantes sur le matériau T700M21 et des résultats expérimentaux.
La confrontation expérimentale a montré les difficultés liées à l’étude de matériaux complexes tels que les composites à matrice organique. Les contraintes en terme de temps et de disponibilité (du matériau et des machines d’essai) n’ont pas permis de réaliser suffisament d’essais de même nature permettant de prendreen compte la variabilité intrinsèque du T700M21.
Il en résulte que le modèle non local a fourni des résultats pertinents sur certains essais mais il a également été partiellement mis en défaut sur deux des six essais. Une analyse de sensibilité a tout de même permis de relativiser cette défaillance du modèle en montrant que les erreurs de prévision du comportement pouvaient être réduites en remettant en cause le caractère déterministe des données du modèle.
Il n’est pas apparu de dépendance entre la forme ou la stratification des éprouvettes et les longueurs internes du modèle. Ce résultat semble indiquer que cette longueur est bien intrinsèque au modèle et indépendante de la structure. Il n’est cependant pas possible de statuer définitivement sur la nature de cescoefficients.

Table des matières
Introduction
1 De la modélisation locale à la modélisation non locale 
1.1 Introduction
1.2 Nonlocal models
1.2.1 Continuous implicit gradient formulation
1.2.2 Finite element formulation
1.3 Resolution algorithms
1.3.1 Newton-Raphson algorithm
1.3.2 Arc-length algorithms
1.4 Parallel computation
1.4.1 Continuous problem
1.4.2 Finite elements application
1.4.3 Parallelization of the arc-length algorithm
1.4.4 Remarks about the local/global resolution algorithm
1.5 Applications
1.5.1 Meshes and constitutive equations
1.5.2 Local vs. nonlocal simulations
1.5.3 Size effects
1.5.4 Effects of element type and boundary conditions
1.5.5 A comment about the choice of the nonlocal variable
1.5.6 Benefits of parallel computation
1.6 Conclusion
2 Du matériau isotrope au composite stratifié 
2.1 Introduction
2.2 Composite layered material: scale and behavior
2.3 Mesoscopic nonlocal model for anisotropic materials
2.3.1 Continuous problem
2.3.2 Finite element formulation
2.3.3 Application to laminates
2.4 Nonlocal model with anisotropic internal lengths
2.4.1 Continuous problem
2.4.2 Finite element formulation
2.4.3 Simplified mesoscopic nonlocal model
2.5 Applications
2.5.1 Meshes and constitutive equations
2.5.2 Local versus mesoscopic nonlocal model
2.5.3 Simplified mesoscopic nonlocal model with anisotropic internal lengths
2.6 Conclusion
3 De la modélisation du pli à la modélisation du stratifié 
3.1 Introduction
3.2 The mechanical problem
3.2.1 Modelling scale
3.2.2 Mechanical problem
3.3 Development of a non-local layered finite element
3.3.1 The local layered element
3.3.2 Extension to the nonlocal framework
3.4 Finite element formulation of the nonlocal layered element
3.4.1 Discretization of the equilibrium equation
3.4.2 Discontinuous discretization of the ith nonlocal equation
3.4.3 Computation of the elementary stiffness matrices
3.5 Applications
3.5.1 Meshes and mechanical behaviors .
3.5.2 Layered element vs. classical element
3.6 Conclusion
4 De la modélisation à l’expérimentation 
4.1 Introduction
4.2 Campagne expérimentale
4.2.1 Matériau et drapages
4.2.2 Essais et formes des éprouvettes
4.2.3 Outils de mesure
4.2.4 Exemple de résultats: cas de l’éprouvette un trou avec un empilement [02 /90 0.5]s
4.3 Loi de comportement du pli unidirectionnel
4.3.1 Formulation
4.3.2 Cinétique d’endommagement pour les ruptures prématurées des fibres faibles
4.3.3 Cinétique d’endommagement pour les ruptures de faisceaux de fibres
4.3.4 Effet de l’endommagement
4.3.5 Loi de comportement
4.3.6 Comportement d’un VER
4.4 Pré-analyse du comportement
4.4.1 Maillages et équivalence modélisation-expérience
4.4.2 Influence des conditions aux limites dans le cas de l’éprouvette un trou[02 / ± 60 2]s92
4.4.3 Modélisation avec une loi non adoucissante dans le cas de l’éprouvette un trou
4.4.4 Modélisation avec la loi de rupture de fibre dans un cadre local dans le cas de l’éprouvette un trou
4.5 Identitication des paramètres du modèle et prédictibilité des simulations
4.5.1 Données matériaux
4.5.2 Remarques sur le choix de la valeur du module d’Young dans le sens fibre
4.5.3 Identification de la cinétique de rupture des fibres faibles
4.5.4 Identification de la cinétique de rupture des faisceaux de fibre
4.5.5 Robustesse de l’identification
4.5.6 Analyse d’une simulation complète: cas de l’éprouvette un trou [02 /90 1/2]
4.6 Longueur interne non locale: artéfact numérique, grandeur matériau ou
4.6.1 Longueur(s) caractéristique(s) dans les CMO
4.6.2 Différentes natures possibles de la longueur interne
4.6.3 Discussion vis-à-vis des résultats expérimentaux
4.7 Conclusion
Conclusions et perspectives 
Bibliographie 
A Notations
B Communications
C Implantation des outils dans un code par éléments finis
C.1 Présentation du code de calcul par éléments finis ZéBuLoN
C.2 Architecture du code de calcul de structure ZéBuLoN
C.2.1 Approches objet en calcul de structure
C.2.2 Les éléments
C.2.3 Les lois de comportement
C.2.4 Les solveurs
C.3 Implantation des éléments finis non locaux
C.3.1 Stratégie d’intégration .
C.3.2 Implantation des éléments
C.3.3 Implantation de la géométrie
C.3.4 Exemples de fichier inp
C.4 Implantation des lois de comportements non locales
C.4.1 Principe
C.4.2 Exemples de fichier matériau
C.5 Implantation des différents solveurs
C.5.1 Stratégie d’intégration
C.5.2 Classe MECHANICAL_RIKS
C.5.3 Classe RIKS_CRITERION
C.5.4 Classe RIKS_CHOICE
C.5.5 Classe RIKS_TEST
C.5.6 Classe RIKS_DECIDE
C.5.7 Exemple de fichier inp
C.6 Parallélisation du solveur à longueur d’arc
C.6.1 Remarques préliminaires
C.6.2 Parallélisation du solveur à longueur d’arc
C.6.3 Exemple de fichier inp
D Résultats expérimentaux
D.1 Empilement [02 /90 0.5]s
D.1.1 Eprouvette deux trous
D.1.2 Eprouvette quatre trous
D.2 Empilement [02 / ± 60 2]s
D.2.1 Eprouvette un trou
D.2.2 Eprouvette deux trous
D.2.3 Eprouvette quatre trous
E Transparents présentés lors de la soutenance

Lire le rapport complet