Keywords: Hybrid variational formulation, visco-plasticity, Raviart-Thomas vector space, energy conserving and decaying schemes, upscaling of visco-plastic material via Bayesian update Hybride Variationsformulierung, Visko-Plastizität, Raviart-Thomas-Vektorraum, Energieerhaltungs-und Zerfallsschemata, Hochskalierung von viskoplastischem Material über Bayes-Update . . . . . . . . . Updating PCE coefficients of input parameters q(q K, q G ) . . . .

In this thesis, a novel approach is developed for visco-plasticity and nonlinear dynamics problems. In particular, variational equations are elaborated following the Helligner-Reissner principle, so that both stress and displacement fields appear as unknown fields in the weak form. Three novel finite elements are developed. The first finite element is formulated for the axisymmetric problem, in which the stress field is approximated by low-order polynomials such as linear functions. This approach yields accurate solutions specifically in incompressible and stiff problems. In addition, a membrane and plate bending finite element are newly designed by discretizing the stress field using the lowest order Raviart-Thomas vector space RT 0 . This approach guarantees the continuity of the stress field over an entire discrete domain, which is a significant advantage in the numerical method, especially for the wave propagation problems. The developments are carried out for the viscoplastic constitutive behavior of materials, where the corresponding evolution equations are obtained by appealing to the principle of maximum dissipation. To solve the dynamic equilibrium equations, energy conserving and decaying schemes are formulated correspondingly. The energy conserving scheme is unconditional stable, since it can preserve the total energy of a given system under a free vibration, while the decaying scheme can dissipate higher frequency vibration modes. The last part of this thesis presents procedures for upscaling of the visco-plastic material behavior. Specifically, the upscaling is performed by stochastic identification method via Baysian updating using the Gauss-Markov-Kalman filter for assimilation of important material properties in the elastic and inelastic regimes.

This thesis is jointly supported between the research project SELT-TUM by Agence Nationale de la Recherche, France and the research fund from the Institut für Wissenschaftliches Rechnen-Technische Universität Braunschweig, Germany. Without these funding programs, the works in this thesis would be scarcely materialized and completed.

I am deeply grateful for having had the opportunity to pursue my doctoral program under the supervision of Prof. Ibrahimbegovic and Prof. Matthies, at Université de technologie de Compiègne & Technische Universität Braunschweig, where I benefited myself by colloquium and its academic environment. The indispensable guidance, lectures, initiatives of Prof. Ibrahimbegovic and Prof. Matthies not only enlarged my knowledge horizon but also have been motivating me throughout the research. I would like to thank the members of the thesis committee Prof. Ouahsine, Prof. Jänicke, Prof. Brank, Prof. Aparicio, Prof. Kowalsky, Prof. Hadzalic and Prof. Hajdo for their kind acceptances as committee members of my thesis defense.

The period of two years studying at Compiègne in France, is one of the most memorable time, thanks to the intellectual and academic sensibilities of Dang Phong, Pablo and Ivan. I was the beneficiary of an excellent team and vigorous academic discussions. I also appreciate the friendships with Tuan Anh, Duc Phong, Nhu, Linh, Yunbo, Tea, Adela, Samir, Dali, Khang, Hoi ... and many others, those together we shared joys and other moments of life.

For the period of one year in Braunschweig of Germany, I would like to express my gratitude to Dr. Vinh and Simona for their kindness of sharing knowledge. Their generosity encouraged me, provided prompt suggestions such that I could familiarize myself with the academic settings, gain the access to facilities of research within a short period of time. In addition, my daily life also had plenty of enjoyable moments thanks to Tobias, Vien, Anh and Sharana.

Finally, my deepest acknowledgment is attributed to my family and my love. They have remotely accompanied me with their unconditional love and care despite the burden of life they are encountering because of global viral contagion. Especially, to my parents, without their precious lifelong dedication and sacrifice, I would hardly achieve this academic milestone.

Résumé

Dans cette thèse, une nouvelle approche est développée pour les problèmes de viscoplasticité et de dynamique non linéaire. En particulier, les équations variationnelles sont élaborées selon le principe de Helligner-Reissner, de sorte que les champs de contrainte et de déplacement apparaissent comme des champs inconnus sous la forme faible. Trois nouveaux éléments finis sont développés. Le premier élément fini est formulé pour le problème axisymétrique, dans lequel le champ de contraintes est approximé par des polynômes d'ordre inférieur tels que des fonctions linéaires. Cette approche donne des solutions précises spécifiquement dans les problèmes incompressibles et rigides. De plus, un élément fini de flexion de membrane et de plaque est nouvellement conçu en discrétisant le champ de contraintes en utilisant l'espace vectoriel de Raviart-Thomas d'ordre le plus bas RT 0 . Cette approche garantit la continuité du champ de contraintes sur tout un domaine discret, ce qui est un avantage significatif dans la méthode numérique, notamment pour les problèmes de propagation des ondes. Les développements sont effectués pour le comportement constitutif visco-plastique des matériaux, où les équations d'évolution correspondantes sont obtenues en faisant appel au principe de dissipation maximale. Pour résoudre les équations d'équilibre dynamique, des schémas de conservation et de décroissance de l'énergie sont formulés en conséquence. Le schéma de conservation de l'énergie est inconditionnellement stable, car il peut préserver l'énergie totale d'un système donné sous une vibration libre, tandis que le schéma décroissant peut dissiper des modes de vibration à plus haute fréquence. La dernière partie de cette thèse présente les procédures d'upscaling du comportement des matériaux visco-plastiques. Plus précisément, la mise à l'échelle est effectuée par une méthode d'identification stochastique via une mise à jour baysienne en utilisant le filtre de Gauss-Markov-Kalman pour l'assimilation des propriétés importantes des matériaux dans les régimes élastique et inélastique.

Mots-clés:

Formulation variationnelle hybride, viscoplasticité, espace vectoriel de Raviart-Thomas, schémas de conservation et de décroissance de l'énergie, mise à l'échelle du matériau visco-plastique via la mise à jour bayésienne.

Zusammenfassung

In dieser Arbeit wird ein neuartiger Ansatz für Viskoplastizitäts-und nichtlineare Dynamikprobleme entwickelt. Insbesondere werden Variationsgleichungen nach dem Helligner-Reissner-Prinzip aufgestellt, so dass sowohl Spannungs-als auch Verschiebungsfelder als unbekannte Felder in schwacher Form erscheinen. Drei neuartige finite Elemente werden entwickelt. Das erste finite Element wird für das axialsymmetrische Problem formuliert, bei dem das Spannungsfeld durch Polynome niedriger Ordnung wie lineare Funktionen angenähert wird. Dieser Ansatz liefert genaue Lösungen speziell bei inkompressiblen und steifen Problemen. Darüber hinaus wird ein finites Element mit Membran-und Plattenbiegung neu entworfen, indem das Spannungsfeld unter Verwendung des Raviart-Thomas-Vektorraums RT 0 niedrigster Ordnung diskretisiert wird. Dieser Ansatz garantiert die Kontinuität des Spannungsfeldes über eine gesamte diskrete Domäne, was insbesondere bei Wellenausbreitungsproblemen ein wesentlicher Vorteil der numerischen Methode ist. Die Entwicklungen erfolgen für das viskos-plastische Materialverhalten, wobei die entsprechenden Evolutionsgleichungen unter Berufung auf das Prinzip der maximalen Dissipation erhalten werden. Um die dynamischen Gleichgewichtsgleichungen zu lösen, werden entsprechende Energieerhaltungs-und -zerfallsschemata formuliert. Das energieerhaltende Schema ist bedingungslos stabil, da es die Gesamtenergie eines gegebenen Systems unter einer freien Schwingung erhalten kann, während das abklingende Schema höherfrequente Schwingungsmoden zerstreuen kann. Im letzten Teil dieser Arbeit werden Verfahren zum Upscaling des viskoplastischen Materialverhaltens vorgestellt. Konkret erfolgt die Hochskalierung durch stochastische Identifikationsverfahren mittels des Baysian Update unter Verwendung des Gauss-Markov-Kalman-Filters zur Assimilation wichtiger Materialeigenschaften im elastischen und inelastischen Bereich.

List of Notations

Introduction

In this first chapter, the motivation of the current research is presented along with the research objectives. Afterwards, the most notable research works relevant to this thesis topic is briefly reviewed. Then, the methodology section is followed by the outline of the thesis, which is at the end of this chapter.

Motivation

During the recent decades, the computational mechanics has been developed enormously so that many tough mathematical, physic and mechanic problems can be now simulated and solved numerically. Some conventional research topics such as wave propagation [START_REF] Kadapa | A novel semi-implicit scheme for elastodynamics and wave propagation in nearly and truly incompressible solids[END_REF], [START_REF] Kwon | An analysis of implicit time integration schemes for wave propagations[END_REF], [START_REF] Quintal | Sensitivity of s-wave attenuation to the connectivity of fractures in fluid-saturated rocks[END_REF], fluid-structure interaction (see [START_REF] Bazilevs | Computational fluid-structure interaction: methods and applications[END_REF], [START_REF] Ohayon | Nonlinear model reduction for computational vibration analysis of structures with weak geometrical nonlinearity coupled with linear acoustic liquids in the presence of linear sloshing and capillarity[END_REF], [START_REF] Huang | Modeling, simulation and validation of supersonic parachute inflation dynamics during mars landing[END_REF], [START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF], [START_REF] Souli | Ale formulation for fluid-structure interaction problems[END_REF], [START_REF] Cai | Moving immersed boundary method[END_REF], [START_REF] Cai | Implicit immersed boundary method for fluid-structure interaction[END_REF], [1], [START_REF] Ding | Cfd/dem coupled approach for the stability of caisson-type breakwater subjected to violent wave impact[END_REF], [START_REF] Walhorn | Fluid-structure coupling within a monolithic model involving free surface flows[END_REF], [START_REF] Zorrilla | An embedded finite element framework for the resolution of strongly coupled fluid-structure interaction problems. application to volumetric and membrane-like structures[END_REF], [START_REF] Boujleben | An efficient computational model for fluid-structure interaction in application to large overall motion of wind turbine with flexible blades[END_REF]) have been still getting special interests by many researchers with an aim to solve the problem more efficient and with less computational cost. Additionally, there have been some newly-emerged research branches such as isogeometric analysis [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF][START_REF] Nguyen | Isogeometric analysis: an overview and computer implementation aspects[END_REF][START_REF] Nguyen-Thanh | Rotation free isogeometric thin shell analysis using pht-splines[END_REF], virtual element method [START_REF] Wriggers | A virtual element method for contact[END_REF][START_REF] Hudobivnik | A low order 3d virtual element formulation for finite elasto-plastic deformations[END_REF]2], topological optimization (see [START_REF] Bendsoe | Topology optimization: theory, methods, and applications[END_REF], [START_REF] Haertel | Design of thermal systems using topology optimization[END_REF], [3], [START_REF] Paulino | Polydyna: a matlab implementation for topology optimization of structures subjected to dynamic loads[END_REF], [START_REF] Wall | Isogeometric structural shape optimization[END_REF]), multi-physics (see [START_REF] Kadapa | A robust and computationally efficient finite element framework for coupled electromechanics[END_REF], [START_REF] Deü | Piezoelectric shunt vibration damping of structural-acoustic systems: finite element formulation and reduced-order model[END_REF], [START_REF] González | Partitioned vibration analysis of internal fluid-structure interaction problems[END_REF], [124], [START_REF] Pérez-Aparicio | Multiphysics and thermodynamic formulations for equilibrium and non-equilibrium interactions: nonlinear finite elements applied to multi-coupled active materials[END_REF], [START_REF] Moreno-Navarro | Plasticity coupled with thermo-electric fields: Thermodynamics framework and finite element method computations[END_REF], [START_REF] Vu | Numerical modelling of non-linear electroelasticity[END_REF], [START_REF] Mehnert | A complete thermo-electroviscoelastic characterization of dielectric elastomers, part i: Experimental investigations[END_REF], [START_REF] Cramer | Coupled chemical and mechanical processes in concrete structures with respect to aging[END_REF], [START_REF] Kowalsky | Modeling of multifield processes regarding the life-cycle of concrete structures[END_REF], [START_REF] Tu | Computational homogenization of the electro-chemically coupled multi-scale transport of li-ions in structural battery electrolytes[END_REF], [START_REF] Kaessmair | Variationally consistent computational homogenization of chemomechanical problems with stabilized weakly periodic boundary conditions[END_REF], [START_REF] Palma | Non-linear and hysteretical finite element formulation applied to magnetostrictive materials[END_REF], [START_REF] Pérez-Aparicio | Finite element analysis and material sensitivity of peltier thermoelectric cells coolers[END_REF], [START_REF] Dinkler | Finite-macro-element for modelling thin structures with integrated piezoelectric sensor-actor-elements[END_REF]), localized failure, and crack propagation (see [START_REF] Brancherie | Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures: Part i: theoretical formulation and numerical implementation[END_REF], [START_REF] Ibrahimbegovic | Combined hardening and softening constitutive model of plasticity: precursor to shear slip line failure[END_REF], [START_REF] Saksala | Combined continuum damage-embedded discontinuity model for explicit dynamic fracture analyses of quasi-brittle materials[END_REF], [START_REF] Saksala | Thermal shock weakening of granite rock under dynamic loading: 3d numerical modeling based on embedded discontinuity finite elements[END_REF], [START_REF] Saksala | Cracking of granitic rock by high frequency-high voltagealternating current actuation of piezoelectric properties of quartz mineral: 3d numerical study[END_REF], [START_REF] Nguyen | Modelling hydraulic fractures in porous media using flow cohesive interface elements[END_REF], [START_REF] Stanić | Crack propagation simulation without crack tracking algorithm: embedded discontinuity formulation with incompatible modes[END_REF]), phase-field method (see [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF], [START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. part ii. coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids[END_REF], [START_REF] Mandal | Fracture of thermo-elastic solids: Phase-field modeling and new results with an efficient monolithic solver[END_REF], [START_REF] Nguyen | Phase field simulation of earlyage fracture in cement-based materials[END_REF], [START_REF] Do | Dynamics framework for 2d anisotropic continuum-discrete damage model for progressive localized failure of massive structures[END_REF], [START_REF] Bharali | Computational homogenisation of phase-field fracture[END_REF]), nano-composite simulation (see [START_REF] Bach | Size effect in nanocomposites: Xfem/level set approach and interface element approach[END_REF], [START_REF] Le | Stochastic continuum modeling of random interphases from atomistic simulations. application to a polymer nanocomposite[END_REF], [START_REF] Le | Stochastic continuum modeling of random interphases from atomistic simulations. application to a polymer nanocomposite[END_REF]), automatic code generation [START_REF] Korelc | Automatic generation of finite-element code by simultaneous optimization of expressions[END_REF], and bio-cells modeling [START_REF] Soleimani | Red blood cell simulation using a coupled shell-fluid analysis purely based on the sph method[END_REF]. The up-rising trend with the aim to apply the uncertainty quantification and analysis into the engineering practices has also been evolving recently (see [START_REF] Hoang | Machine learning-based conditional mean filter: a generalization of the ensemble kalman filter for nonlinear data assimilation[END_REF], [START_REF] Wu | Parameter identification for phase-field modeling of fracture: a bayesian approach with samplingfree update[END_REF], [START_REF] Noii | Bayesian inversion for anisotropic hydraulic phase-field fracture[END_REF], [START_REF] Matthies | Inverse problems in a bayesian setting[END_REF], [START_REF] Stabile | Bayesian identification of a projection-based reduced order model for computational fluid dynamics[END_REF], [START_REF] Dunić | A thermo-mechanically coupled finite strain model considering inelastic heat generation[END_REF], [START_REF] Hoang | Stochastic multiscale model of mems stiction accounting for high-order statistical moments of non-gaussian contacting surfaces[END_REF], [START_REF] Ibrahimbegovic | Reduced model of macro-scale stochastic plasticity identification by bayesian inference: Application to quasi-brittle failure of concrete[END_REF], [START_REF] Park | Uncertainty reduction for model error detection in multiphase shock tube simulation[END_REF], [START_REF] Dodwell | Multilevel monte carlo simulations of composite structures with uncertain manufacturing defects[END_REF]). Beside that, several works are dedicated on other challenging topics such as parallel/higher-performance computing and digital twins. Several advanced research and industrial applications are shown in Fig. 1.1, e.g. aircraft futuristic development of DLR1 , Orion parachute of NASA2 , and wind turbine tower. Starkly, the research society is continuously looking for better and more efficient solutions to both conventional and newly-rising problems in the field of computational mechanics. In the early days, the conventional finite element method is dominant thanks to its feasibility and accuracy in many mechanics problems. However there still exists several problems, e.g. notably incompressible problems, in which this method is not capable to yield precise results. And the former method can not guarantee the continuity of stress field under dynamics analysis. To overcome those difficulties, the mixed finite element approach, e.g. Helligner-Reissner and Hu-Washizu methods, has been applied widely. The mixed approach is proved to be very effective for both compressible and incompressible mechanics problems. Relevant researches using mixed variational formulation for small deformation are applied in electrostatics and magnetostatics [START_REF] Moreno-Navarro | Multi-field variational formulations and mixed finite element approximations for electrostatics and magnetostatics[END_REF], laminated beams [START_REF] Trinh | A mixed inverse differential quadrature method for static analysis of constant-and variable-stiffness laminated beams based on hellinger-reissner mixed variational formulation[END_REF], incompressible elasticity [START_REF] Viebahn | A simple and efficient hellingerreissner type mixed finite element for nearly incompressible elasticity[END_REF], and stabilized quadrilateral Reissner-Mindlin plate element [START_REF] Gruttmann | A stabilized one-point integrated quadrilateral reissner-mindlin plate element[END_REF]. Likewise, other works for finite deformation using mixed formulation can be listed as finite rotation using Biot's stress [START_REF] Kožar | Finite element formulation of the finite rotation solid element[END_REF], higher-order elements in morphoelasticity [START_REF] Kadapa | On the advantages of mixed formulation and higher-order elements for computational morphoelasticity[END_REF], improved shell element [START_REF] Wagner | An improved quadrilateral shell element based on the hu-washizu functional[END_REF], discontinuous Galerkin method [START_REF] Terrana | A hybridizable discontinuous galerkin method for both thin and 3d nonlinear elastic structures[END_REF], structure of polyconvex stored energy functions [START_REF] Betsch | A mixed variational framework for the design of energy-momentum schemes inspired by the structure of polyconvex stored energy functions[END_REF], cable elements with a continuous or discontinuous axial force field [START_REF] Crusells-Girona | A mixed formulation for nonlinear analysis of cable structures[END_REF], quasi-incompressible finite elasticity [START_REF] Schröder | A new mixed finite element based on different approximations of the minors of deformation tensors[END_REF], and first order system of hyperbolic equations [START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part i: Total lagrangian isothermal elasticity[END_REF][START_REF] Gil | A first order hyperbolic framework for large strain computational solid dynamics. part ii: Total lagrangian compressible, nearly incompressible and truly incompressible elasticity[END_REF][START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part iii: Thermo-elasticity[END_REF]. Examples for the above method in finite deformation are shown in Fig. 1.2. In the same main stream, this thesis concentrates in the development of two-field mixed finite element method for visco-plasticity and elasto-dynamics problems. Additionally, one last chapter is devoted to parameter identification method for inelastic problem and multi-scale with visco-plasticity problems. A more comprehensive literature review relevant to thesis contents are given in the following Section 1.2.

.2 -Mixed finite element method in finite deformation: (a) Displacement magnitude's contour in rolling-up beam using algorithm of [START_REF] Kožar | Finite element formulation of the finite rotation solid element[END_REF], (b) Pressure in twisted beam [START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part i: Total lagrangian isothermal elasticity[END_REF], (c) Displacement magnitude in upward-bending of Lily flower [START_REF] Kadapa | On the advantages of mixed formulation and higher-order elements for computational morphoelasticity[END_REF] 1

.2 Literature review

There are several typical material's inelastic behaviors, e.g. plasticity, damage, and visco-plasticity. The main difference lays on the underlying mathematical models and also the yield functions. Many research works relevant to material's inelastic response are solved thoroughly via the conventional displacement-based finite element method. Nevertheless these type of problem can be also formulated via mixed finite element method with an aim to improve the accuracy of stress field. Namely, the Helligner-Reissner approach can be a potential choice to derive the weak variational form. An advantage of Helligner-Reissner variational approach is that both displacement and stress fields can be considered as independent fields. This approach guarantees the accurate computation for stress, see [START_REF] Auricchio | The dimensional reduction modelling approach for 3d beams: Differential equations and finite-element solutions based on hellinger-reissner principle[END_REF], [START_REF] Auricchio | The dimensional reduction approach for 2d non-prismatic beam modelling: a solution based on hellingerreissner principle[END_REF], and [START_REF] Viebahn | A simple and efficient hellingerreissner type mixed finite element for nearly incompressible elasticity[END_REF]. Expanding further into dynamics problem, one can use a cluster of stable time-stepping schemes which are successfully developed for long-term computation. The unconditional sta-bility property of some straightforward time-stepping schemes, which can support the computation for a large number of time steps, are proved with a sound theoretical basis in [START_REF] Bathe | Conserving energy and momentum in nonlinear dynamics: a simple implicit time integration scheme[END_REF], [START_REF] Bathe | On a composite implicit time integration procedure for nonlinear dynamics[END_REF], [START_REF] Hulbert | Automatic time step control algorithms for structural dynamics[END_REF], [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF], [START_REF] Erlicher | The analysis of the generalized-α method for non-linear dynamic problems[END_REF], and [START_REF] Leyendecker | Objective energy-momentum conserving integration for the constrained dynamics of geometrically exact beams[END_REF]. Other examples are the second-order accurate scheme [START_REF] Artioli | Second-order accurate integration algorithms for von-mises plasticity with a nonlinear kinematic hardening mechanism[END_REF] and β 1 /β 2 -Bathe method with L-stable state for wave propagation problem [START_REF] Malakiyeh | New insights into the β1/β2-bathe time integration scheme when l-stable[END_REF]. Another research effort is given at formulating a so-called energy decaying scheme with an aim to decay the contribution of higher-frequency modes in a coarse finite element mesh where they are not full resolved, see [START_REF] Ibrahimbegovic | Energy conserving/decaying implicit time-stepping scheme for nonlinear dynamics of three-dimensional beams undergoing finite rotations[END_REF]5]. Specifically, the variation of velocities and stress computation within a time interval is modified so that it reduces the participation of high-frequency modes in controllable manner via dissipation coefficients α and β. A relevant work, see [START_REF] Malakiyeh | The bathe time integration method revisited for prescribing desired numerical dissipation[END_REF], is given in developing a time integration schemes with capability of prescribing the numerical dissipation.

Rather than using the Helligner-Reissner approach, a condensed two-field variational formulation can be derived directly from the regularized three-field variational formulation [START_REF] Hughes | On drilling degrees of freedom[END_REF] by eliminating the independent rotation field, see relevant works [START_REF] Hughes | Numerical assessment of some membrane elements with drilling degrees of freedom[END_REF][START_REF] Ibrahimbegovic | Mixed finite element with drilling rotations for plane problems in finite elasticity[END_REF][START_REF] Ibrahimbegovic | A robust quadrilateral membrane finite element with drilling degrees of freedom[END_REF]. The new-formed two-field variational formulation is capable of offering an improvement of stability properties. One of the critical point in this variation formulation is that the non-symmetric part of stress tensor is assumed. Consequently the stress is not restricted to be symmetric, it enables us to interpolate stress field by Whitney's approximation [START_REF] Moreno-Navarro | Multi-field variational formulations and mixed finite element approximations for electrostatics and magnetostatics[END_REF][START_REF] Bossavit | Whitney forms: A class of finite elements for three-dimensional computations in electromagnetism[END_REF], which is successful deployed for scalar field problems, or by the lowest order of Raviart-Thomas (RT 0 ) vector space, see [START_REF] Dubois | Dual raviart-thomas mixed finite elements[END_REF][START_REF] Ervin | Computational bases for rtk and bdmk on triangles[END_REF]. The latter approximation can improve the computation accuracy and maintain the continuity of stress traction over element's edges, which is an advantageous point for dynamics. Particularly, an enhanced finite element can be developed from the traditional constant stress/strain element, where three more nodes are added at mid-edges to place the degrees of freedom for stress field. Here, the stress degrees of freedom are actually components of normal traction at a corresponding element edge with the use of RT 0 discretization. In additional to the Raviart-Thomas interpolation, there exist other choices of stress interpolation, such as Brezzi-Douglas-Marini (BDM) [START_REF] Ervin | Computational bases for rtk and bdmk on triangles[END_REF][START_REF] Braess | Error indicators for mixed finite elements in 2-dimensional linear elasticity[END_REF], though this type of interpolation is complex to implement. To improve solution's accuracy in dynamics problems, several research works opt for the discontinuous Galerkin method [START_REF] Johnson | Numerical solution of partial differential equations by the finite element method[END_REF] or by isogeometric analysis [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF] to improve solution's accuracy in dynamics problems.

Due to the important role in engineering practice, the plate bending element has been still receiving many attentions of research and development works. In general, there exists thick and thin plate theories. The main difference between the two models concerns the treatment of shear strain and its work-conjugate shear forces. In the case of the Kirchhoff plate model, the plate fiber is assumed to remain perpendicular to mid-surface upon deformation, which eliminates shear strains, contrary to the Reissner-Mindlin plate theory that accounts for shear force by considering fiber to remain straight but not necessarily perpendicular to plate mid-surface. Several early researches are pioneered by [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF], [START_REF] Mindlin | Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates[END_REF], and [START_REF] Reissner | On the theory of transverse bending of elastic plates[END_REF]. Numerous research works on Reissner-Mindlin plate aim to alleviating the locking phenomenon, such as by using enhanced assumed strain [START_REF] Hughes | Finite elements based upon mindlin plate theory with particular reference to the four-node bilinear isoparametric element[END_REF], [START_REF] Bathe | A four-node plate bending element based on mindlin/reissner plate theory and a mixed interpolation[END_REF], [4], [START_REF] Ibrahimbegovic | Quadrilateral finite elements for analysis of thick and thin plates[END_REF], [START_REF] Felippa | Construction of optimal 3-node plate bending triangles by templates[END_REF], [START_REF] Mostafa | A solid-shell corotational element based on andes, ans and eas for geometrically nonlinear structural analysis[END_REF]. Other research works employ mixed finite element, see [START_REF] Batoz | A study of three-node triangular plate bending elements[END_REF], [START_REF] Bathe | A four-node plate bending element based on mindlin/reissner plate theory and a mixed interpolation[END_REF], [START_REF] Brezzi | Mixed-interpolated elements for reissner-mindlin plates[END_REF], [START_REF] Zienkiewicz | Plate bending elements with discrete constraints: new triangular elements[END_REF], [START_REF] Auricchio | A shear deformable plate element with an exact thin limit[END_REF], [START_REF] Auricchio | A triangular thick plate finite element with an exact thin limit[END_REF], [START_REF] Gruttmann | A stabilized one-point integrated quadrilateral reissner-mindlin plate element[END_REF], [START_REF] Wagner | An improved quadrilateral shell element based on the hu-washizu functional[END_REF]. Additional approaches are the sub-element method [START_REF] Clough | A refined quadrilateral element for analysis of plate bending[END_REF], the reduced integration method [START_REF] Hughes | A simple and efficient finite element for plate bending[END_REF], the series solution [START_REF] Brank | On boundary layer in the mindlin plate model: Levy plates[END_REF] and the alternative alpha finite element method [START_REF] Nguyen | An alternative alpha finite element method with discrete shear gap technique for analysis of isotropic mindlin-reissner plates[END_REF]. A more detailed review of many important developments on plate bending problem can be found in [START_REF] Bathe | Finite element method[END_REF] and [START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF]. Some more recent review and comparison works be-tween common plate and shell elements are reported in [START_REF] Katili | A comparative formulation of t3γs, dst, dkmt and mitc3+ triangular plate elements with new numerical results based on s-norm tests[END_REF], [START_REF] Batoz | On the formulation and evaluation of old and new efficient low order triangular plate bending elements with shear effects[END_REF], [START_REF] Lavrenčič | Hybrid-mixed low-order finite elements for geometrically exact shell models: Overview and comparison[END_REF], and [START_REF] Nguyen | Notable highlights on locking-free techniques of reissner-mindlin plate finite elements in elastostatics[END_REF]. For nonlinear analysis of plate bending problem, some significant works are given on fracture computation [START_REF] Dolbow | Modeling fracture in mindlin-reissner plates with the extended finite element method[END_REF], plasticity [START_REF] Papadopoulos | An analysis of inelastic reissner-mindlin plates[END_REF][START_REF] Wagner | An adaptive strategy for the multi-scale analysis of plate and shell structures with elasto-plastic material behaviour[END_REF], visco-plasticity [START_REF] Dujc | On stress resultant plasticity and viscoplasticity for metal plates[END_REF], a von Kármán plate model [START_REF] Dongi | Active panel flutter suppression using self-sensing piezoactuators[END_REF], and mixed convection boundary-layer flow over a flat plate in a porous medium [START_REF] Guedda | Analytical and chpdm analysis of mhd mixed convection over a vertical flat plate embedded in a porous medium filled with water at 4 o c[END_REF]. Several works relevant to further refinement of inelastic behavior in the bending plate problem can be found in [START_REF] Owen | Elasto-viscoplastic analysis of anisotropic laminated plates and shells[END_REF], [START_REF] Dinis | Elastic-viscoplastic analysis of plates by the finite element method[END_REF], [START_REF] Klosowski | Dynamics of elasto-viscoplastic plates and shells[END_REF], and [START_REF] Dujc | On stress resultant plasticity and viscoplasticity for metal plates[END_REF].

The field of uncertainty quantification has been gaining many attentions from the engineering society of computational and numerical methods, e.g. dynamical systems in the frequency domain [START_REF] Römer | An adaptive sparse grid rational arnoldi method for uncertainty quantification of dynamical systems in the frequency domain[END_REF], low-frequency electromagnetic devices [START_REF] Galetzka | A multilevel monte carlo method for high-dimensional uncertainty quantification of low-frequency electromagnetic devices[END_REF], sampling-free linear Bayesian update for fracture phase-field modeling [START_REF] Wu | Parameter identification for phase-field modeling of fracture: a bayesian approach with samplingfree update[END_REF], Bayesian inference of heterogeneous visco-plastic material parameters [START_REF] Janouchová | Bayesian inference of heterogeneous viscoplastic material parameters[END_REF], stochastic multiple crack patterns [START_REF] Gerasimov | Stochastic phase-field modeling of brittle fracture: computing multiple crack patterns and their probabilities[END_REF], polynomial chaos expansions for Maxwell's source problem [START_REF] Georg | Conformally mapped polynomial chaos expansions for maxwell's source problem with random input data[END_REF], polynomial chaos in evaluating failure probability [START_REF] Janouchová | Polynomial chaos in evaluating failure probability: A comparative study[END_REF], polynomial chaos for stochastic differential equations [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF], adaptive Gaussian process for optimization problem [START_REF] Kim | Quantile surrogates and sensitivity by adaptive gaussian process for efficient reliability-based design optimization[END_REF], inverse problem in piezoelectric material to identify damage [START_REF] Rus | Damage identification inverse problem for a piezoelectric material[END_REF], structural fragility analysis underground motion [START_REF] Yi | Gaussian mixture-based equivalent linearization method (gm-elm) for fragility analysis of structures under nonstationary excitations[END_REF], and random fluctuations in material behavior of metals with visco-plastic and damage [START_REF] Kowalsky | Random fluctuations of material behaviour in fe-damage-analysis[END_REF]. There are several other researches adopting neural networks and machine learning into the engineering practices, e.g solving phase-field modeling of fracture [START_REF] Goswami | Transfer learning enhanced physics informed neural network for phase-field modeling of fracture[END_REF], predicting energy consumption [START_REF] Bui | An artificial neural network (ann) expert system enhanced with the electromagnetism-based firefly algorithm (efa) for predicting the energy consumption in buildings[END_REF], detecting damage in wind turbine tower [START_REF] Nguyen | Damage identification of windturbine tower using modal properties-based artificial neural networks[END_REF][START_REF] Nguyen | Vibration-based damage detection in wind turbine towers using artificial neural networks[END_REF], tracking high-dimensional non-Gaussian state models [START_REF] Hoang | Machine learning-based conditional mean filter: a generalization of the ensemble kalman filter for nonlinear data assimilation[END_REF], modeling a digital twin for composite materials [START_REF] Ghanem | Probabilistic learning and updating of a digital twin for composite material systems[END_REF], predicting load-bearing capacity of concrete-filled steel tubular [START_REF] Le | Development of user-friendly kernel-based gaussian process regression model for prediction of load-bearing capacity of square concrete-filled steel tubular members[END_REF], computational homogenization [START_REF] Nguyen | A surrogate model for computational homogenization of elastostatics at finite strain using high-dimensional model representation-based neural network[END_REF], and calibration of nonlinear mechanical models [START_REF] Mareš | Artificial neural networks in the calibration of nonlinear mechanical models[END_REF][START_REF] Kučerová | Soft computing-based calibration of microplane m4 model parameters: Methodology and validation[END_REF]. A cluster of research works is relevant to parameter identification, e.g. using conditional expectation [START_REF] Matthies | Parameter estimation via conditional expectation: a bayesian inversion[END_REF], applying Gauss-Markov-Kalman filter for the Bayesian estimation [START_REF] Matthies | Bayesian parameter estimation via filtering and functional approximations[END_REF], developing sampling-free non-linear Bayesian update [START_REF] Matthies | Inverse problems in a bayesian setting[END_REF], considering stored energy and dissipation in multi-scale analysis [START_REF] Sarfaraz | Bayesian stochastic multi-scale analysis via energy considerations[END_REF], applying the reduced model to quasi-brittle failure of concrete [START_REF] Ibrahimbegovic | Reduced model of macro-scale stochastic plasticity identification by bayesian inference: Application to quasi-brittle failure of concrete[END_REF], estimating effective parameters in anisotropic hydraulic phase-field fracture [START_REF] Noii | Bayesian inversion for anisotropic hydraulic phase-field fracture[END_REF], estimating model coefficients of a novel turbulent flow model over porous media [START_REF] Friedman | Bayesian calibration of model coefficients for a simulation of flow over porous material involving svm classification[END_REF], developing an efficient computational method to sample the posterior random variables [START_REF] Hoang | An efficient computational method for parameter identification in the context of random set theory via bayesian inversion[END_REF], stabilizing reduced order models in computational fluid dynamics problems [START_REF] Stabile | Bayesian identification of a projection-based reduced order model for computational fluid dynamics[END_REF]. The idea to bring parameter assimilation into the multi-scale problem is rather challenging. Some research works on multi-scale problems can be listed as probabilistic analysis of localized failure [START_REF] Ibrahimbegovic | Probabilistic multiscale analysis of inelastic localized failure in solid mechanics[END_REF], efficient code-coupling strategy [START_REF] Ibrahimbegovic | Code-coupling strategy for efficient development of computer software in multiscale and multiphysics nonlinear evolution problems in computational mechanics[END_REF], embedded discontinuity capable of interpreting full set of 3D failure modes for heterogeneous materials [START_REF] Ibrahimbegovic | Multiscale model with embedded discontinuity discrete approximation capable of representing full set of 3d failure modes for heterogeneous materials with no scale separation[END_REF], programming and computational procedure [START_REF] Rukavina | Ed-fem multiscale computation procedure for localized failure[END_REF], poro-viscoelastic substitute model [START_REF] Jänicke | A poro-viscoelastic substitute model of fine-scale poroelasticity obtained from homogenization and numerical model reduction[END_REF], modelling of micromorphic continua [START_REF] Jänicke | Two-scale modelling of micromorphic continua[END_REF], elastic/viscoelastic compounds [START_REF] Schüler | Multi-scale modelling of elastic/viscoelastic compounds[END_REF], micro-scale or nano-scale composites [START_REF] Chatzigeorgiou | Multiscale modelling for composites with energetic interfaces at the micro-or nanoscale[END_REF], curing processes [START_REF] Klinge | Modeling of curing processes based on a multi-field potential. single-and multiscale aspects[END_REF], large eddy simulation of turbulent flow [START_REF] Gravemeier | An algebraic variational multiscale-multigrid method for large eddy simulation of turbulent flow[END_REF], and Nitsche-type extended variational formulation for incompressible two-phase flow [START_REF] Schott | A face-oriented stabilized nitsche-type extended variational multiscale method for incompressible two-phase flow[END_REF].

Research objectives and methodology

The work in this thesis aim at providing the following contributions:

• The Helligner-Reissner approach is formulated for the conventional viscoplasticity problem, in which the stress field is interpolated by the low-order polynomials.

• The regularized variational formulation is derived for plane membrane and also plate bending problems. The lowest-order Raviart-Thomas vector space RT 0 is employed to approximate stress field in the former problem, and force resultants (moments and shear forces) fields in the latter problem. With this type of interpolation, it is advantageous to maintain the continuity of stress traction or force resultant's projection over each element's edges, specifically in dynamics and wave propagation problems.

• The energy conserving and decaying schemes are developed to guarantee the stability of those numerical computation over a long period of time.

• The Bayesian updating and Gauss-Markov-Kalman filter methods are applied to identify critical material's properties in the inelastic problem and the multiscale problem with visco-plasticity.

Thesis outline

The thesis outline is as follows. The mixed variational formulation using low-order Polynomial functions for stress approximation is presented in Chapter 2. This chapter is devoted to the numerical implementation of visco-plasticity problem.

In addition, the energy conserving and decaying schemes are formulated to control the overall energy of a system and dissipate higher frequency modes. Several selective simulations of the visco-plasticity problem in a thick-walled pipe under dynamic loading scenarios are simulated. In Chapter 3, the mixed variational formulation using lowest-order Raviart Thomas vector space for stress approximation is presented for the plane membrane problem. The numerical implementation of the linear elastodynamics problem including the chosen discrete approximation of displacement and stress fields are given. In addition, the energy conserving time integration scheme is developed to maintain the total energy of the system over computational period. Several illustrative simulations of the elastodynamics problem under static and dynamic loading scenarios are presented. In Chapter 4, the mixed variational formulation using lowest-order Raviart Thomas vector space for stress resultants approximation is developed. The theoretical and variational formulations of Mindlin-Reissner plate in linear elastodynamics using a regularized functional is presented. The numerical implementation of linear elastodynamics problem including the discretization of displacement (lateral deflection and rotations) and force resultants fields is shown. The energy conserving/decaying time integration scheme is reconfigured to maintain computational stability over a long period of time for the vibration of the plate bending problem. Several selective simulations of Mindlin-Reissner plate under dynamic loading scenarios are illustrated. In Chapter 5, the method of parameter identification using the Bayesian updating and the Gauss-Markov-Kalman filter methods are presented. A set of principle concepts and methods are laid out. Several principle material's properties in some highlight inelastic and multi-scale problems are selected to be assimilated via the two proposed methods. In detail, the multi-scale problem with the visco-plasticity is proceeded with an aim to upscale material's properties in the macro-scale via data from the micro-scale. The micro-scale is simulated via the Lattice model, while the macroscale is simulated via an enhanced triangular element developed in Chapter 3 with additional capabilities of linear isotropic and kinematic hardening visco-plasticity.

The conclusions and perspectives are given in Chapter 6.

2

Hellinger-Reissner formulation for the visco-plasticity problem

In this chapter, the theoretical formulation of the visco-plasticity problem is presented by using Helligner-Reissner variational formulation.

Then the interpolation for stress and displacement fields are proceeded by low-order polynomial functions. Afterwards, the energy conserving and decaying are developed for the visco-plasticity problem

Theoretical formulation

The Hellinger-Reissner variational principle is employed for the visco-plasticity problem at small strains. Let Ω be a bounded domain in R 3 with piece-wise smooth boundary ∂Ω, satisfying

∂Ω = ∂Ω u ∪ ∂Ω t and ∂Ω u ∩ ∂Ω t = ∅, (2.1)
where ∂Ω u denotes the displacement boundary and ∂Ω t the stress boundary. The linear total strain (u) tensor in small strain theory is defined by the symmetric gradient of the displacement displacement u field as in Eq. (2.2)

(u) := symm ∇u = 1 2 ∇u + (∇u) T . (2.
2)

The symbol ∇u denotes the gradient of a displacement field, e.g. ∇u = [u i,j ]. The above total strain tensor is assumed as a combination of an elastic strain e tensor and a visco-plastic strain vp tensor, hence we have the relation e = -vp . The problem of visco-plasticity at small strains is characterized by the strain energy functional Π int (u) as follows

Π int (u) = Ω ψ e ( -vp ) dV. (2.3)
where the elastic free energy density function ψ e ( -vp ) is a function of the elastic strain tensor

ψ e ( -vp ) = 1 2 ( -vp ) : C : ( -vp ) , ( 2.4)
with C as the fourth-order elasticity tensor. From such elastic free energy function, the formulation of stress σ tensor can be derived as

σ := ∂ψ e (•) ∂ = C : ( -vp ) .
The yield function φ(σ) governs the evolution of the visco-plastic strain vp along with time. The non-positive value of the yield function indicates the elastic regime without the development of the visco-plastic strain ( ˙ vp = 0). Any positive value of the yield function implies the visco-plastic regime, with the corresponding change of the visco-plastic strain ( ˙ vp = 0) according to

φ(σ) := sgn F |F | m ; F = 3J 2 /σ 2 y -1; φ(σ) ≤ 0 ⇒ ˙ vp = 0, φ(σ) > 0 ⇒ ˙ vp = 0, (2.5)
where σ y is the yield stress and J 2 is the stress invariant. The dot in ˙ vp denotes the derivative of the visco-plastic strain with respect to time, or the rate of the visco-plastic strain. In the following, the parameter value is set to be m = 1 in Eq. (2.5), thus the classical visco-plasticity model (e.g. [START_REF] Ibrahimbegovic | Classical plasticity and viscoplasticity models reformulated: theoretical basis and numerical implementation[END_REF]) is recovered, where the corresponding yield function accepts all admissible values in R as follows

φ ≡ F = 3J 2 /σ 2 y -1; J 2 = 1 2 dev(σ) : dev(σ), (2.6)
where the operator dev(σ) stands for the deviatoric part of stress tensor. The external load potential Π ext (u) and kinetic energy K( u) are written as follows

Π ext (u) = Ω f • u dV + ∂Ω 2 t • u dA ; K( u) = 1 2 Ω ρ u • u dV, (2.7)
with the body force f in Ω and the traction vector t on ∂Ω 2 , ρ as the material density, and the overhead dot in u denotes the derivative with respect to time. The dot between two given vectors denotes the dot product, e.g. a • b = a i b i . The stress tensor is considered as an additional independent variable by appealing to the Legendre transformation for the free energy function in Eq. (2.4). The complementary stored energy χ(σ) is a function of the stress field, in which C -1 is the inverse of the fourth-order elasticity tensor

ψ e ( -vp ) = σ : ( -vp ) -χ(σ); χ(σ) = 1 2 σ : C -1 : σ. (2.8)
From Eq. (2.3) and Eq. (2.8), the Hellinger-Reissner functional Π HR (σ, u) is now written as functional of the stress and displacement fields

Π HR (σ, u) = Ω [σ : ( -vp ) -χ(σ)] dV -Π ext (u) + K( u).
(2.9)

The least action principle can be used to obtain corresponding variational equations of the functional Π HR (σ, u) above, which can be written as equivalent d'Alembert equations (e.g. [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]) in terms of the virtual displacement δu and the virtual stress δσ as follows

G u (u, σ; δu) = Ω δ : σ dV + Ω δu • ρ ü dV - Ω δu • f dV - ∂Ω 2 δu • t dA = 0; G σ (u, σ; δσ) = Ω [δσ : ( -vp ) -δσ : ∂ σ χ(σ)] dV = 0,
(2.10) where the virtual strain δ tensor is defined as ∇ s δu in the same manner of Eq. (2.2) and the symbol ∂ σ χ(σ) denotes the derivative of the complementary stored energy with respect to the stress field. The corresponding Euler-Lagrange equations can be obtained directly from the variations in Eq. (2.10) as

G u (u, σ; δu) = - Ω δu • (div σ + f -ρ ü) dV = 0 ⇒ div σ + f = ρ ü; G σ (u, σ; δσ) = Ω δσ : ( -vp -∂ σ χ(σ)) dV = 0 ⇒ -vp = ∂ σ χ(σ).
(2.11) The second law of thermodynamics by means of the Clausius-Duhem inequality is applied to guarantee the thermodynamical admissibility of the constitutive relations

D vp = σ : ˙ -ψe = σ : ˙ vp ≥ 0, (2.12)
where D vp denotes the visco-plastic dissipation. In principle, all the stress states are admissible, but those outside the elastic domain are penalized by an additional term P (•) directly proportional to the penalty factor 1/η (with η as the viscosity coefficient depending on material property). The modified visco-plastic dissipation D vp η is now written as follows

-D vp η (σ) = -σ : ˙ vp + 1 η P (φ(σ)). (2.13)
The simplest choice of the penalty term is a quadratic functional P (•) placing higher penalty on the stress states which are further outside the elastic domain

P (φ(σ)) =      1 2 φ(σ) 2 ; φ(σ) > 0. 0; φ(σ) ≤ 0. ⇒ d dφ P (φ(σ)) =< φ(σ) >, (2.14)
where < • > is the Macauley parenthesis, which is defined as

< φ >:= (φ(σ) + |φ(σ)|)/2 = φ(σ); φ(σ) > 0. 0; φ(σ) ≤ 0. (2.15)
The evolution of the visco-plasticity can be derived by means of the principle of maximum visco-plastic dissipation. From all possible states σ * satisfying the yield criterion φ (σ * ) ≥ 0, the one σ that maximizes the visco-plastic dissipation D vp η is the solution. This problem can also be written as constrained minimization problem

D vp η (σ) = max φ(σ * )≥0 D vp η (σ * ) ⇔ -D vp η (σ) = min φ(σ * )≥0 -D vp η (σ * ) . (2.16)
The Kuhn-Tucker optimality conditions for the minimization problem defined in the equation above lead to the evolution equation of the visco-plastic strain

0 = ∂ -D vp η (σ) ∂σ = ∂[-σ : ˙ vp + 1 η P (φ(σ))] ∂σ = -˙ vp + 1 η < φ > ∂φ ∂σ . (2.17)
Hence, the rate of the visco-plastic strain is written as

˙ vp = 1 η < φ > ∂φ ∂σ =: β(σ); ∂φ ∂σ = 3 σ 2 y dev(σ). (2.18)

Numerical implementation

In order to provide the solution for the system of differential equations in Eqs. (2.11) and (2.18) above, the physical domain is discretized via the finite element approximation, while the nonlinear dynamic response is obtained by means of a proposed time-stepping scheme. The main computational task of such a procedure is to calculate the values of the nodal displacements and their time derivatives, along with the stress interpolation parameters over a single time interval, which is discussed further in this section.

Finite element approximation

The stress and displacement fields are interpolated independently. At time t n+1/2 , the displacement field u n+1/2 and acceleration field ün+1/2 are interpolated from nodal displacement d n+1/2 and nodal acceleration a n+1/2 vectors via shape function N matrix, respectively

u n+1/2 = N d n+1/2 ; ün+1/2 = N a n+1/2 .
The corresponding strain approximation can be derived from the displacement field

n+1/2 = Bd n+1/2 ; B = ∇ s N .
Meanwhile, the stress field in each element is interpolated directly using a complete linear polynomial S matrix, which is similar to the stress interpolation [START_REF] Pian | Rational approach for assumed stress finite elements[END_REF][START_REF] Wilson | Use of incompatible displacement modes for the calculation of element stiffnesses or stresses[END_REF] in terms of stress interpolation parameters τ n+1/2 matrix

σ n+1/2 = Sτ n+1/2 .

Discrete weak form

The index n+1/2 for the obvious terms is only shown where it is needed. From Eq. (2.10), the discretized weak forms for a typical element Ω e are written in the matrix forms as follows It is noted that the set of motion equations above can be computed with any selected value of α ∈ [0, 1]. However, the second-order accuracy is achieved only for the midpoint scheme, by selecting α = 1/2

G e u = δd T Ω e B T Sτ dV + Ω e ρN T N a dV - Ω e N T f dV - ∂Ω e 2 N T t dA = 0; G e σ = δτ T Ω e S T (Bd -vp -∂ σ χ(σ)) dV = 0. ( 2
σ n+α = (1 -α)σ n + ασ n+1 ; α = 1/2; ˙ vp n+α = β(σ n+α ) ⇒ vp n+1/2 = vp n + ∆t 2 β(σ n+1/2 ); ∂ σ χ (σ) = C -1 : σ.
(2.23)

By exploiting the auxiliary results above, the variation ∆G e σ (•) can be written explicitly as follows

∆G e σ (∆τ , ∆u) = δτ T Ω e S T B dV ∆d n+1/2 - Ω e S T D vp n+1/2 S dV ∆τ n+1/2 ; D vp n+1/2 := ∆t 2 β (σ n+1/2 ) + C -1 ; β (σ n+1/2 ) = ∂β(σ n+1/2 )/∂σ n+1/2 .
(2.24) In the same manner, for any arbitrary virtual stress parameter values δτ T , an additional local equation valid independently inside each element Ω e is obtained, which can be written as follows (2.25)

h e,(i+1) n+1/2 (τ , u) = h e,(i) n+1/2 (τ , u) + ∆h e,(i) n+1/2 (∆τ , ∆u) = 0; h e,(i) n+1/2 (•) := Ω e S T B dV d (i) n+1/2 - Ω e S T vp
Thus, one can establish the full system of governing equations with both nodal displacements and element stress parameters are unknown variables

A nel e=1 F e,T η M e -H e n+1/2 F e   ∆τ e,(i) n+1/2 ∆d e,(i) n+1/2   = -A nel e=1   r e,(i) u,n+1/2 (•) h e,(i) n+1/2 (•)   .
(2.26)

The system of equations above can be condensed, so that only nodal displacement is the unknown field. The stress field would be computed correspondingly from the displacement fields. From Eq. (2.26b), the stress state in each element can be rewritten as = 0.

∆τ e n+1/2 = H e,-
(2.28) Assembling each element contributions over the entire discrete domain yields the governing equation of a system in matrix form

A nel e=1 F e,T n+1/2 H e,-1 n+1/2 F e + η M e ∆d e,(i) n+1/2 = -A nel e=1 r e,(i) u,n+1/2 + F e,T H e,-1 n+1/2 h e,(i) n+1/2 ⇔ K ef f,(i) n+1/2 ∆d (i) n+1/2 = -R (i) n+1/2 .
(2.29) In each iteration, the step of displacement is computed via Eq. (2.29). Subsequently, the step of stress is computed at the level of each element from equation Eq. (2.27). The stress and displacement fields are correspondingly updated after each computational iteration of Newton-Raphson method

τ (i+1) n+1/2 = τ (i) n+1/2 + ∆τ (i) n+1/2 ; d (i+1) n+1/2 = d (i) n+1/2 + ∆d (i) n+1/2 .
(2.30)

Energy conserving scheme

An energy conserving (EC) scheme is developed to guarantee the stability of the computation over long period for the visco-plasticity problem. To achieve the secondorder accurate O(∆t 3 ) solution [START_REF] Artioli | Second-order accurate integration algorithms for von-mises plasticity with a nonlinear kinematic hardening mechanism[END_REF], the mid-point time integration scheme is employed. For enforcing the second-order accuracy in the visco-plasticity, we must simultaneously integrate the visco-plastic deformation with both displacement and stress, as opposed to the operator split algorithm ( [START_REF] Ibrahimbegovic | Classical plasticity and viscoplasticity models reformulated: theoretical basis and numerical implementation[END_REF], [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]) that separates the integration of these two, and thus reduces the time-stepping scheme to first order (e.g. [START_REF] Chorin | Product formulas and numerical algorithms[END_REF]). Accordingly, the nodal velocity v n+1/2 and acceleration a n+1/2 at the time step t n+1/2 can be formulated in the form of nodal displacement and velocity

v n+1/2 = (d n+1 -d n )/∆t; a n+1/2 = (v n+1 -v n )/∆t.
(2.31)

The increment of displacement field within an interval of time step is computed by using the mid-point approximation

d n+1 -d n = 1 2 ∆t(v n+1 + v n ). (2.32)
This equation can be formed in an alternative expression showing the simple update of the displacement vector

d n+1 = d n + u; u := 1 2 ∆t(v n+1 + v n ).
(2.33)

The stress state at t n+1/2 is selected in the following algorithmic form. This proposal, which is left open at the moment, will be explained shortly

τ alg n+1/2 := 1 2 (τ n+1 + τ n ). ( 2
(d n+1 -d n ) T f int n+1/2 = (d n+1 -d n ) T f ext n+1/2 .
The internal work includes two obvious components including kinetic and potential energy

(d n+1 -d n ) T f int n+1/2 = (d n+1 -d n ) T Ω B T Sτ n+1/2 dV +(d n+1 -d n ) T Ω ρN T N a n+1/2 dV.
The work done by inertia force can be simplified thanks to results obtained in Eq. (2.31b) and (2.32)

(d n+1 -d n ) T M a n+1/2 = 1 2 (v n+1 + v n ) T M (v n+1 -v n ) = 1 2 v T n+1 M v n+1 - 1 2 v T n M v n = K n+1 -K n , (2.35)
where K n+1 is the kinetic energy of a system at time t n+1 . It is noted that the mass matrix M is symmetric, so this identity holds v T n+1 M v n -v T n M v n+1 = 0. In the following, the subscript n + 1/2 for matrix F and H is dropped for shorter expression. The work done by the potential force can be also expressed as a combination of 2 groups

(d n+1 -d n ) T F T τ n+1/2 = 1 2 (d n+1 -d n ) T F T (τ n+1 + τ n ) = 1 2 d T n+1 F T τ n+1 - 1 2 d T n F T τ n + 1 2 d T n+1 F T τ n - 1 2 d T n F T τ n+1 L .
(2.36) Since the scalar L has a mix variables from both time step t n and t n+1 , further simplification of L should be elaborated. Via exploiting the equation Eq. (2.25b) with a test stress function δτ T = τ n+1 -τ n , the following identity is derived based on the fact that the matrix H is always symmetric, even when the system goes into the visco-plastic regime

(τ n+1 -τ n ) T F d n+1/2 -Hτ n+1/2 =0 ⇒ 1 2 d T n+1 F T τ n+1 - 1 2 τ T n+1 Hτ n+1 - 1 2 d T n F T τ n + 1 2 τ T n Hτ n = 1 2 d T n+1 F T τ n - 1 2 d T n F T τ n+1 L .
(2.37) The matrix H can be decomposed in two matrices, the first one H vp n+1/2 is a function of the visco-plastic strain rate ˙ vp n+1/2 and the other H e is a function of the material elastic modulus. It is clear that the matrix H vp n+1/2 vanishes under the elastic regime, where the visco-plastic strain does not evolve ( ˙ vp n+1/2 = 0). With the last identity, Eq. (2.36) can be rewritten

(d n+1 -d n ) T F T τ n+1/2 = 1 2 d T n+1 F T τ n+1 - 1 2 d T n F T τ n + 1 2 d T n+1 F T τ n+1 - 1 2 τ n+1 Hτ n+1 - 1 2 d T n F T τ n + 1 2 τ n Hτ n = d T n+1 F T τ n+1 - 1 2 τ T n+1 H e τ n+1 P n+1 -(d T n F T τ n - 1 2 τ T n H e τ n ) Pn + ( 1 2 τ T n H vp n+1/2 τ n - 1 2 τ T n+1 H vp n+1/2 τ n+1 D vp n+1/2 ) = P n+1 -P n + D vp n+1/2 ,
(2.38) where P n+1 is the potential energy of the system at time t n+1 and D vp n+1/2 is the visco-plastic dissipation within an interval of time step. From the above derivation, one can express the balance of the internal energy and the external work in following equation

(K n+1 + P n+1 ) E n+1 -(K n + P n ) En +D vp n+1/2 = (d n+1 -d n ) T f ext n+1/2 ∆W n+1/2 ; E n+1 -E n + D vp n+1/2 =∆W n+1/2 .
(2.39)

The symbols E n and E n+1 denote the total energy of the system at time step t n and t n+1 , respectively. We can conclude that the proposed algorithmic constitutive equation in Eq. (2.34) would conserve the total energy E n+1 = E n + ∆W n+1/2 for any bounded external loading as long as the system undergoing the elastic regime (D vp n+1/2 = 0). This property is viewed as "unconditionally stable" [START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF]. Otherwise, once the system undergoes the visco-plastic regime, the total energy can not be conserved anymore due to the existence of the visco-plastic dissipation (D vp n+1/2 > 0). However, this non-negative visco-plastic dissipation is always bounded via the controllable penalty factor as shown in Section 2.1.

Energy decaying scheme

An energy decaying (ED) scheme is an alternative approach to control potential instability problems of time-stepping algorithms, with an additional advantage of a providing superior accuracy in stress computations by eliminating the highfrequency noise especially for a set of stiff equations (e.g. see [START_REF] Ibrahimbegovic | Energy conserving/decaying implicit time-stepping scheme for nonlinear dynamics of three-dimensional beams undergoing finite rotations[END_REF]). The latter implies a large difference between max/min eigenvalues in tangent stiffness, which can also occur in the present case when one part of the domain has turned viscoplastic and the rest still remains elastic. In this case, we will enforce a so-called modified constitutive algorithm with the ability of decaying contribution of higher frequencies over each interval of time step. Such algorithmic updates for of displacement and stress at time step t n+1/2 are represented by modifying the conservative terms, as defined in Eqs. (2.32) and (2.34), with the corresponding dissipative terms:

d n+1 -d n := ∆t 2 (v n+1 + v n ) + α∆t(v n+1 -v n ); τ alg n+1/2 := 1 2 (τ n+1 + τ n ) + β(τ n+1 -τ n ).
(2.40)

where α and β are the dissipation coefficients, which are typically chosen as α, β ∈ [0 0.1]. The corresponding dissipative correction terms α∆t(v n+1 -v n ) and β(τ n+1 -τ n ) remain typically small, or practically vanish, for low frequency modes with a small difference of successive values within a typical time step; the latter establishes the consistency of the present modification in the spirit of the EC scheme, for the case where high frequency modes are not triggered. Meanwhile, the correction terms are less likely to vanish in high frequency modes, where the corresponding successive values within a time step can be quite significant. Therefore, the energy decaying scheme is only introducing a slight perturbation, which is not likely to fully reduce the scheme performance to the first-order accuracy. With such a choice of algorithmic update of displacement and stress updates, the velocity and stress at time step t n+1 take a modified form

v n+1 = 1 (1/2 + α)∆t (d n+1 -d n ) - 1/2 -α 1/2 + α v n ; τ n+1 = 1 1/2 + β τ alg n+1/2 - 1/2 -β 1/2 + β τ n .
(2.41)

One can apply the same implementation as presented in Section 2.2.3. The work done by inertia force takes a new form

(d n+1 -d n ) T M a n+1/2 = 1 2 v T n+1 M v n+1 - 1 2 v T n M v n + α(v n+1 -v n ) T M (v n+1 -v n ) D K,n+1/2 = K n+1 -K n + D K,n+1/2 .
(2.42) where D K,n+1/2 is the numerical dissipation of kinetic energy within an interval of time step. In the following, the subscript n + 1/2 for matrix F and H is dropped for shorter expression. The work done by potential force takes a combination form of 3 groups

(d n+1 -d n ) T F T τ n+1/2 = 1 2 (d n+1 -d n ) T F T (τ n+1 + τ n ) = 1 2 d T n+1 F T τ n+1 - 1 2 d T n F T τ n + 1 2 d T n+1 F T τ n - 1 2 d T n F T τ n+1 L + β(τ n+1 -τ n ) T H(τ n+1 -τ n ) D P,n+1/2 .
(2.43) where D P,n+1/2 is the numerical dissipation of potential energy within an interval of time step. The similar identity to the one in Eq. (2.37) can be written as

(τ n+1 -τ n ) T F d n+1/2 -Hτ n+1/2 =0; ⇒ 1 2 d T n+1 F T τ n+1 - 1 2 τ T n+1 Hτ n+1 - 1 2 d T n F T τ n + 1 2 τ T n Hτ n + β(τ n+1 -τ n ) T H(τ n+1 -τ n ) = 1 2 d T n+1 F T τ n - 1 2 d T n F T τ n+1 L .
(2.44) With this identity on hand, the work done by potential force takes a new form

(d n+1 -d n ) T F T τ n+1/2 = d T n+1 F T τ n+1 - 1 2 τ T n+1 H e τ n+1 P n+1 -(d T n F T τ n - 1 2 τ T n H e τ n ) Pn + ( 1 2 τ n H vp n+1/2 τ n - 1 2 τ n+1 H vp n+1/2 τ n+1 ) D vp n+1/2 + β(τ n+1 -τ n ) T H(τ n+1 -τ n ) D P,n+1/2 = P n+1 -P n + D vp n+1/2 + D P,n+1/2 .
(2.45) From the above derivation, the balance of the internal energy and external work can be formed in following equation

(K n+1 + P n+1 ) E n+1 -(K n + P n ) En +D vp n+1/2 + (D K + D P ) D KP,n+1/2 = (d n+1 -d n ) T f ext n+1/2 ∆W n+1/2 ; E n+1 -E n + D vp n+1/2 + D KP,n+1/2 =∆W n+1/2 .
(2.46) where D KP,n+1/2 is the numerical dissipation of total energy due to the dissipative constitutive algorithm in Eq. 2.41. It is obvious that switching between energy conserving/decaying schemes can be executed flexibly without changing the main code of user-defined element, since the update of displacement and velocity updates can be programmed in a separate subroutine for the time-stepping integration algorithm. Additionally, the amount of numerical dissipation is fully controllable via an appropriate selection of the dissipation coefficients α and β.

Numerical examples

In order to illustrate the favorable properties of the newly proposed EC and ED schemes for the visco-plastic problem, several numerical simulations are presented and compared against a conventional visco-plasticity approach. All the computations are programmed as a user-defined element in computer code FEAP v8.4, see [START_REF] Taylor | Feap-finite element analysis program[END_REF]. The output data is then post-processed via Matlab scripts for plotting. The computations are performed for both quasi-static and dynamic pressure loads applied to a thick-walled pipe by using the proposed finite element formulation based on Hellinger-Reissner variational pinciple, hybrid-stress interpolations and time-stepping schemes developed for the visco-plasticity in the previous section. The computed results are compared against either the exact (analytic) solution in the elastic regime or the available results in the visco-plasticity regime [START_REF] Lubliner | Plasticity theory[END_REF], as well as against the previous study results [START_REF] Hughes | Unconditionally stable algorithms for quasistatic elasto/visco-plastic finite element analysis[END_REF]. Additionally, the performance of the proposed element is compared against the quadrilateral rectangle element embedded the B-bar method (Q4/P1).

Quasi-static loadings

The elasticity and visco-plasticity problems of a pipe under a quasi-static internal pressure is considered. In quasi-static analysis, the inertial term, proportional to the mass matrix, is neglected in Eq. (2.26). The pipe is assumed very long compared to the outer radius, so that it is subjected to a plane strain constraint with zero deformation along the pipe axis. Thanks to the axisymmetry of this 3D problem, the domain can be discretized by finite elements using the shape functions varying only in the r direction. The inner and outer radii of the pipe are chosen as r 1 = 1 and r 2 = 2. The properties of the pipe material are listed as: elastic modulus E = 30.10 6 , Poisson's ratio ν = 0.3, yield stress σ y = 30.10 3 and 1/η = 10 -8 . Regarding the boundary conditions, on the inner end, we apply pressure, while the outer end of the pipe is free of stress. The domain is discretized into 8 such elements, each with 3 nodes and 2 Gauss points as in Fig. 2.5(a) in r direction, which provide the second order displacement accuracy. This interpolation is handled with isoparametric elements using quadratic shape functions, see [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]. Thus, the displacement u matrix and strain matrix are functions of nodal displacement d matrix, see [START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF]

u = 3 i=1 N i u i → u = N d; r θ z = ∂u ∂r u r 0 → = Bd.
(2.47)

The elasticity matrix in Eq. (2.23c) is rewritten for axisymmetric and plane strain conditions

C = E (1 + ν)(1 -2ν)    1 -ν ν ν ν 1 -ν ν ν ν 1 -ν    .
(2.48)

The interpolation of stress field is selected based on the equilibrium equations. The latter, expressed in polar coordinates with zero body forces and plane strain con-straint, reduces to

∂σ rr ∂r + 1 r (σ rr -σ θθ ) = 0; ∂σ θθ ∂θ = 0; ∂σ zz ∂z = 0. (2.49)
Since the quadratic shape functions for the displacement field in this element result with the linear strain field, we choose similarly the corresponding stress approximations. In particular, the stress σ rr should be a linear function to avoid the trivial case that it is a constant and equal to σ θθ . If the stress σ zz and σ θθ are constants in order to satisfy the Eqs. (2.49b,c), they can not satisfy the condition for plane strain (σ rr + σ θθ )ν = σ zz . As the result, all of three stress components are selected as linear functions. When the element gets smaller, Eq. (2.49) would be satisfied in a "weak" sense. The chosen approximation for stress field can be interpolated from 6 stress parameters, which become part of unknowns in this problem:

   σ rr σ θθ σ zz    =    1 r 0 0 0 0 0 0 1 r 0 0 0 0 0 0 1 r              τ 1 τ 2 τ 3 τ 4 τ 5 τ 6           ⇒ σ = Sτ . (2.50)
All the ingredients of chosen discrete approximations are readily replaced into the theoretical formulation presented in Section 2.2. The analytical solution of stress field for the pipe under internal pressure p, see [START_REF] Lubliner | Plasticity theory[END_REF], is given as follows

σ r = -p (r 2 /r 1 ) 2 -1 r 2 2 r 2 -1 ; σ θ = p (r 2 /r 1 ) 2 -1 r 2 2 r 2 + 1 ; σ z = 2νp (r 2 /r 1 ) 2 -1 .
(2.51)

(a) A pipe under internal pressure p 0 As shown in Fig. 2.2(a), the pipe under the elastic regime is analyzed by the quasi-static analysis, in which the internal pressure p 0 load varies from 0 to 5.10 3 within 100 time steps. A coarse mesh with 8 by 8 elements is required by using enhanced quadrilateral element (Q4/P1) and the result is extracted at Gauss points laying on the radial line AB, as shown in Fig. 2.4(a). Meanwhile, a mesh with only 8 elements is adequate by using proposed element, as shown in Fig. 2.5(a). The solution from the proposed element matches exactly the analytic solution, while the stress computed by Q4/P1 element shows a staggered trend, as shown in Figs. 2.2(bd). The total time of CPU (by a standard portable computer) to solve this problem is 0.38s by proposed element and 0.86s by Q4/P1 element. It is clear that the proposed element has an advanced computational cost with respect to the mesh generation and also the computational time. To alleviate the high incompressibility constraint and enhance the accuracy of two-dimensional elements in this type of problem, one can also employ other methods such as enhancing displacement gradient by the F-bar method (see [START_REF] Neto | Design of simple low order finite elements for large strain analysis of nearly incompressible solids[END_REF][START_REF] Neto | F-bar-based linear triangles and tetrahedra for finite strain analysis of nearly incompressible solids. part i: formulation and benchmarking[END_REF]) and/or generate a finer mesh. As shown in Fig. 2.3, the pipe under the visco-plastic regime is analyzed by quasistatic analysis, in which the internal pressure p 0 load stay constant at the value of 21.10 3 . After 15 seconds, the visco-plastic zone develops to the middle of the pipe. This current result is simulated with time step ∆t = 10 -4 . The zone with F +1 > 1 equivalent to φ > 0 , which has the non-zero visco-plastic deformation, corresponds to the region with r 1 ≤ r ≤ r p = 1.4r 1 , as shown in Fig. 2.3(a). Meanwhile the elastic zone is located at r p ≤ r ≤ r 2 . The current solution matches with previous numerical study by Hughes and Taylor, see [START_REF] Hughes | Unconditionally stable algorithms for quasistatic elasto/visco-plastic finite element analysis[END_REF], which is in agreement with the exact solution given by Prager and Hodge with the ratio r p /r 1 = 1.4, see page 271 of [START_REF] Lubliner | Plasticity theory[END_REF]. Compared to the elastic analysis, the similar trend of the radial stress can be observed, as shown in Fig. 2.3(b). Meanwhile, the circumferential and longitudinal stresses show clearly two different zones due to the existence of the visco-plasticity, as shown in Fig. 2.3(c,d). The stress field is also compared to that from Q4/P1 element under perfectly plastic regime. The location r p of the separation line is the same as that from the proposed element, as shown in Fig. 2.4(b). It is also noticed that the solution from Q4/P1 element shows a staggered trend, as shown in Fig. 2.3(b-d). For the proposed element, the visco-plastic problem is solved within 3 iterations under a quadratic convergence rate (energy residual up-to 10 -15 ). The proposed element shows an advanced smooth and accurate stress field compared to the Q4/P1 element in the pipe.

Dynamic loadings

The visco-plastic problem inside selected thick-walled pipes is simulated under dynamic loads. There are two models with different radius for the simulation. The first one has r 1 = 1 and r 2 = 2 referred as "soft" model, while the other has r 1 = 5 and r 2 = 10 referred as "stiff" model. The difference between two concerns the corresponding ratio between the smallest and largest eigenvalues of the tangent stiffness. The same number of elements is used to perform the computations for both models. The results presented subsequently for cases 1-3 belong to the "soft" model and for cases 4-6 belong to the "stiff" model. The chosen time step in all simulations is ∆t = 10 -5 , which is able to capture high frequency response.

Regarding the energy conserving (EC) scheme, only one triangle pulse with a peak of pressure p 0 is introduced to both models as shown in Fig. 2.5(a), which subsequently leads to a free vibration of the system. Several cases with different pulse loading duration are examined, as listed in Table 2.1. Regarding the energy decaying (ED) scheme, an additional pulse of short duration is later introduced into the system. This pulse is applied on end-nodes of each element with opposite direction, as shown in Fig. 2.5(b), in order to activate high frequency modes of the system. The results are extracted at nodes 1, 9 and 17 as shown in Fig. 2.5(a) and at elements 1, 4 and 8 as shown in Fig. 2 The results in Fig. 2.6 show that the total internal energy of the system is conserved exactly under a free vibration. The case 1-3 is from soft model and case 4-5 is from stiff model. After releasing the external load, the total energy of the system, including kinetic and potential energies, remains constant and coincides with the external energy introduced to system by the pulse. For example, the total energy of the system is around 0.09 in case 1 and 0.02 in case 4, as shown in Figs. 2.6(a,b). Longer duration of loading leads to higher total energy gained by the system. It is noted that these cases are fully under the elastic regime without any evolution of the visco-plasticity strain. The dynamic response of the system including displacement, velocity and stress are recorded at node number 1, 9 and 17, which are, respectively, at the inner radius, in the middle and at the outer radius of the mesh, as shown in Fig. 2.5(a). As shown in Fig. 2.7, the amplitude of displacement at node 1 is higher than those at other nodes, since the pulse is applied at this node. It is noted that the displacement field in soft model is dominated by one frequency, whereas this is no longer the case in stiff model. A similar observation is made regarding the velocity and stress fields, as shown in Figs. 2.8 and 2.9, where higher values of absolute displacements lead to higher stress. During the elastic regime, the norm of deviatoric part of largest stress || dev(σ)|| always remains lower than the critical value 2 3 σ y ≈ 24.5 × 10 3 (see Eq. (2.6)), which is equivalent to φ < 0. The same numerical examples are further examined under the visco-plastic regime. Here, higher values of the pressure pulse have to be introduced in order to activate the visco-plastic regime. In particular, the stiff model requires a higher pulse amplitude compared to the soft model. During the visco-plastic regime, the total energy of the system is not conserved due to the evolution of the visco-plastic dissipation, as shown in Fig. 2.10. However, since the visco-plastic dissipation remains rather small, the total energy of the system is approximately comparable to the total external work after releasing the pulse. For example, the total external work is around 38 in case 1 and 800 in case 4, as shown in Figs. 2.10(a,b). The viscoplastic dissipation increases overall; however it stays constant sometimes when there is no excursion in the visco-plastic regime. This trend can be explained clearly by corresponding variations of the stress field, as in Fig. 2.11.

The stresses are extracted from the first Gauss point of elements 1, 4 and 8 (numbered starting from the inner radius), as shown in Fig. 2.5(b). During the period when every element has the value || dev(σ)|| < 2 3 σ y (or equivalent to φ < 0), the whole system is under the elastic regime and there is no evolution of the viscoplastic strain, which keeps the corresponding visco-plastic dissipation constant. By contrast, the visco-plastic dissipation evolves when the yielding criterion is satisfied, precisely || dev(σ)|| > 2 3 σ y . This energy dissipation depends on a type of material in which the viscosity coefficient η is taken into account. In Fig. 2.12, we show the total energy of the system in the elastic regime, computed by ED scheme. The only dissipation sources are introduced by energydecaying scheme, where the higher frequency dissipation is controlled by dissipation factors α, β that regulate potential and kinetic dissipation energies. In both soft and stiff models, there are 3 cases with different dissipation factors (α, β) ∈ {(0.01, 0.01), (0.05, 0.05), (0.1, 0.1)}. It is obvious that the dissipative energy increases proportionally with respect to the chosen values of dissipation factors (α, β) as shown in Figs. 2.12(a,b). In case 3 of soft model, the total energy decreases slightly due to the dissipative energy over time as shown in 2.12(c). The same pattern can be better recognized after the injection of high frequency modes contribution produced by the second pulse. This high frequency pulse does not dissipate neither kinetic nor potential energy of the soft model as much as it dissipates in case 6 of the stiff model. The noise generated by the second pulse is dissipated out quickly for both kinetic and potential energies, as shown in Fig. 2.12(d). When dissipation factors (α, β) approach zeros, Figs. 2.12(a,b) will agree with Figs. 2.6(e,f) in the period before introducing the second pulse at ∆t 2 = 30.10 -3 .

In Fig. 2.13, the stresses are taken again in elements 1, 4 and 8 at the Gauss points as in the previous simulation. As shown in Figs. 2.13(a,c,e), the second pulse in soft model does not bring much noise to the stress, except for an obvious effect that it increases the magnitude of stresses. However, the perturbation of high frequency modes can be viewed from cases 4-6 in the stiff model. As shown in Fig. 2.13(f), high frequency stresses are quickly filtered out of the system with α = β = 0.1, which is faster than the other cases with smaller values of α and β, as shown in Fig. 2.13(b,d). As shown in Fig. 2.14, the total energy of the system under the visco-plastic regime are computed by ED scheme. There are in total 3 sources of dissipation, namely visco-plastic dissipation, numerical dissipation of energy-decaying scheme for kinetic and for potential energy. The direct correlation between the higher values of α and β versus resulting numerical dissipation is obviously captured in Figs. 2.14(a) and 2.14(b). The visco-plastic dissipation increases in time whenever the system satisfies the yield criterion, as shown in Figs. 2.14(c) and 2.14(d). It is recognized that numerical dissipation is much higher than for the visco-plasticity. From the results in Fig. 2.12(a,b) and Fig. 2.14(a,b), the trend and relative dissipation of the system due to ED scheme is comparable between elastic and visco-plastic regimes.

The stresses are taken again in elements 1, 4 and 8 at the first Gauss points, as presented in Fig. The stiff model is discretized via 16 elements as in Figs. 2.16(a,c) and via 32 elements as in Figs. 2.16(b,d), respectively. The energy-decaying scheme results obtained with dissipation factors (α, β) ∈ {(0.01, 0.01), (0.1, 0.1)} are then evaluated with the two above discrete models, where only the elastic regime is considered. In Figs. 2.16(a,c), the higher values of α and β are chosen, the faster higher frequency modes are filtered out for the stiff system with coarse mesh. In contrast, the lower values of α and β are chosen, the slower higher frequency modes are filtered out from the system with fine mesh, as in Figs. 2.16(b,d). Moreover, any finer mesh can better resolve high frequencies, and thus produces smaller numerical dissipation and less significant on amplitude decay due to higher-frequency modes. From the above observations, the most proper choice for the numerical dissipation factors α and β depends both on the type of the system (soft vs. stiff), and on the element size in the mesh. This best choice would be left to the users depending upon the kind of results that are the most acceptable for them. Practically, it is recommended that higher values of dissipation factors (α, β) should be chosen for a system with coarser mesh than those for a system with finer mesh.

The visco-plastic problem itself dissipates energy through the visco-plastic strain.

The energy decaying scheme can be employed if the user expects that the external forces can lead to a significant high frequency modes, in which the contribution from these modes should be dissipated especially on a coarse mesh where they can not be properly resolved. Meanwhile, the energy conserving scheme can be selected within any case, regardless of elastic or visco-plastic regimes. It is recorded that the dynamic problem coupled with the visco-plastic regime is solved within 3 iterations under the quadratic convergence rate (energy residual up-to 10 -15 ). 3

Mixed variational formulation for the elasto-dynamics membrane problem

In this chapter, a regularized variational formulation of linear elastodynamics is derived following the work of [START_REF] Hughes | On drilling degrees of freedom[END_REF]. Then the displacement is interpolated by low-order polynomial functions while the stress field are interpolated by the lowest-order Raviart-Thomas vector space. Afterwards, the energy conserving scheme is applied for the dynamics computation. Several numerical examples are presented to validate the performance of proposed algorithms.

Regularized variational formulation

In this section, we follow closely the work of Hughes and Brezzi [START_REF] Hughes | On drilling degrees of freedom[END_REF], see also [START_REF] Ibrahimbegovic | A novel membrane finite element with an enhanced displacement interpolation[END_REF][START_REF] Ibrahimbegovic | Mixed finite element with drilling rotations for plane problems in finite elasticity[END_REF][START_REF] Ibrahimbegovic | A robust quadrilateral membrane finite element with drilling degrees of freedom[END_REF]. The problem is limited to linear elastodynamics with small strain theory.

Let Ω ⊂ R 3 be an open set with piece-wise smooth boundary.

∂Ω = ∂Ω u ∪ ∂Ω t and ∂Ω u ∩ ∂Ω t = ∅, (3.1)
where ∂Ω u denotes the displacement boundary and ∂Ω t the stress boundary. Both stress and displacement fields are considered as independent fields. The key point here is that the stress is not assumed to be symmetric, but the symmetry of stress tensor is enforced through corresponding moment equilibrium equation of the strong form of the problem ∀x ∈ Ω div σ + f = ρü;

(1) skew σ = 0;

(2) symm σ = C : symm ∇u;

(3) u = ū| ∂Ωu and σ • n = t| ∂Ωt , (4)

(3.2)
where ( 1) to ( 4) are, respectively, the equilibrium equations, the symmetry conditions for stress, the constitutive equations and the boundary condition. The secondrank stress tensor is decomposed into symmetric and skew-symmetric parts

σ = symm σ + skew σ, ( 3.3)
where

symm σ = 1 2 σ + σ T ; skew σ = 1 2 σ -σ T . (3.4)
The displacement gradient can be additionally decomposed into symmetric part defining the infinitesimal strain tensor (u)

(u) := symm ∇u = 1 2 ∇u + ∇u T , ( 3.5)
and the skew-symmetric part defining the infinitesimal rotation

ω(u) := skew ∇u = 1 2 ∇u -∇u T . (3.6)
The following identity between arbitrary symmetric and skew-symmetric tensors is exploited for several forthcoming equations

symm A : skew B = 0; ∀A, B ∈ R 2 . (3.7)
For the isotropic material, the fourth-order elastic constitutive modulus tensor C is conventionally defined as

C =λ1 ⊗ 1 + 2µI s ; C ijkl =λδ ij δ kl + 2µ 1 2 (δ ik δ jl + δ il δ jk ); i, j, k ∈ [1, 2, 3], (3.8)
where

λ = νE (1 + ν)(1 -2ν) ; µ = E 2(1 + ν) , ( 3.9)
with E and ν as Young's modulus and Poisson's ratio, respectively. For the weak form of the boundary value problem, the mixed form of variational functional [START_REF] Hughes | On drilling degrees of freedom[END_REF] with three independent fields is given as follows

Π(u, ω, σ) = - 1 2 Ω symm σ : C -1 : symm σ dV + Ω σ : (∇u -ω) dV -Π ext (u), (3.10
) where u ∈ U, ω ∈ W and σ ∈ S are trial displacement, rotation and stress fields, respectively. The variational statement in Eq. (3.10) above requires that the stress tensor σ and the infinitesimal rotation tensor ω belong to the space of square integrable functions over the region Ω, L 2 (Ω) as follows

S = {σ|σ ∈ L 2 (Ω)} ; W = {ω|ω ∈ L 2 (Ω); symm ω = 0} ; (3.11)
The space of trial displacements, however, must belong to the functions whose generalized derivative belongs to L 2 (Ω), i.e. a subset of the Sobolev space H 1 (Ω) denoted as H 1 0 (Ω) which satisfies the boundary condition on

∂Ω u U = u|u ∈ H 1 0 (Ω) . ( 3.12)
The regularized functional (see [START_REF] Hughes | On drilling degrees of freedom[END_REF]) is derived by adding an extra term to improve the ellipticity of the standard functional

Π (u, ω, σ) = - 1 2 Ω symm σ : C -1 : symm σ dV + Ω σ : (∇u -ω) dV - 1 2 -1 Ω | skew σ| 2 dV -Π ext (u). (3.13)
Such a regularization allows us to make any convenient choice of finite element discrete approximation upon regularization. The external virtual work G ext (v) is formulated in the spirit of d'Alembert principle with an additional term from inertia force, which appears in the elasto-dynamics problem

G ext (v) = - Ω ρü • v dV + Ω f • v dV + ∂Ωt t • v dA. (3.14)
The corresponding Euler-Lagrange equations are obtained by taking the variations with respect to the displacement δu, stress δσ and rotation δω field

G u (u, ω, σ; δu) := ∂Ωt δu • σ • n -t dA - Ω δu • (div σ + f -ρü) dV = 0; G ω (u, ω, σ; δω) := Ω δω : skew σ dV = 0; G σ (u, ω, σ; δσ) := Ω symm δσ : C -1 : symm σ -symm ∇u dV - Ω δσ : skew ∇u -ω --1 skew σ dV = 0. (3.15)
It is noted that the Gauss divergence theorem is employed to obtain the result in Eq. (3.15a). Since the trial displacement δu, stress δσ and rotation δω are arbitrary, all equations corresponding to the strong form in Eqs. (3.2) can be recovered, plus an additional relationship from Eq. (3.15c) showing that ω = skew ∇u --1 skew σ. (3.16) Here, we choose to eliminate the rotation field in mixed formulation above by using the additional relationship from Eq. (3.15c). This result is plugged back into the regularized functional in Eq. (3.13), and thus we recover a two-field variational formulation with displacement and stress as independent fields by exploiting the identity Eq. (3.7), which can be written as

Π (u, σ) = - 1 2 Ω symm σ : C -1 : symm σ dV + Ω symm σ : symm ∇u dV + 1 2 -1 Ω | skew σ| 2 dV -Π ext (u).
(3.17) The corresponding variational equation of the functional Π (u, σ) above with respect to the trial displacement field v can be written as follows:

G u (u, σ; δu) = Ω symm ∇δu : σ dV - Ω δu • (f -ρü)dV - ∂Ωt δu • t dA = 0.
(3.18) Similarly, by using again the identity Eq. (3.7), the corresponding variational equation of the functional Π (u, σ) above with respect to the trial stress field δσ can be written as follows

G σ (u, σ; δσ) = Ω symm δσ : symm ∇u -C -1 : symm σ + skew δσ : 1 skew σ) dV = Ω δσ : symm ∇u -C -1 : symm σ - 1 skew σ dV = Ω δσ : symm ∇u -Ĉ-1 : σ dV = 0. (3.19)
with an optimal choice = 2µ [START_REF] Hughes | On drilling degrees of freedom[END_REF]. It is noted that the new elasticity tensor in the above equation is formed by regrouping the symmetric and skew-symmetric part of stress tensor into a single stress tensor σ, which can be written as follows

Ĉ-1 = C -1 : I s -C -1 : I a , ( 3.20)
where the two fourth-order tensors I s and I a are given as follows (for more details see Appendix A)

I s = 1 2 (δ ik δ jl + δ il δ jk ); I a = 1 2 (δ ik δ jl -δ il δ jk ). (3.21)

Numerical implementation

In the following, we will simplify discussion to a 2D plane strain case. In defining the discrete approximation, the variational formulation will first be restated in a matrix notation. The details of mapping from tensor to matrix form are given in Appendix A. The mapping operation L is introduced by taking advantages of the plane strain condition. The second order tensor of stress and symmetric part of displacement gradient are transformed into a 4 by 1 matrix

σ = [σ ij ] (i,j)∈[1,3] →L(σ) := σ h =      σ 11 σ 22 σ 12 σ 21      ; symm ∇u = [ ij ] (i,j)∈[1,3] →L(symm ∇u) := h =      11 22 12 21      . (3.22)
Finally, the fourth order tensor of elasticity is transformed into a 4 by 4 matrix:

Ĉ-1 = Ĉ-1 ijkl (i,j,k,l)∈[1,3] → L( Ĉ-1 ) := D h 4×4 .
(3.23)

Finite element approach

For the space discretization, the well-known constant strain/stress triangle element (CST) is selected as the element to enhance its performance, in which the physical meaning for the degree of freedom of stress field is additionally discussed. The conventional discretization method is employed for the displacement field [START_REF] Hughes | The finite element method: Linear static and dynamic finite element analysis[END_REF], while the Raviart-Thomas vector space with the lowest order (RT 0 ) [START_REF] Dubois | Dual raviart-thomas mixed finite elements[END_REF][START_REF] Ervin | Computational bases for rtk and bdmk on triangles[END_REF] is introduced to approximate the stress field. For time approximation, the appropriate timestepping integration is consequently formulated to preserve the total energy of a system in dynamics. The resulting discrete problem is solved in a single level of computation, with the global set of equilibrium equations at time step n+1/2. Let us consider a domain Ω ⊂ R 2 , represented by the discrete triangularized mesh as a set T of triangle Ω e such that Ω = ∪ Ω e ∈T Ω e . The common parts between two neighboring triangles are either an edge, two corner nodes, a midpoint node or empty. In the parent coordinate system, the coordinate of vertex node 1, 2 and 3 are (0,0), (0,1) and (1,0), respectively. These nodes contain the degrees of freedom for the displacement field. Meanwhile the mid-edge nodes 4, 5 and 6 are (1/2,1/2), (0,1/2) and (1/2,0). These nodal degrees of freedom correspond to the stress field. The shape functions to interpolate the displacement field at any point in the given element Ω e is given as follows

N 1 = 1 -ξ -η; N 2 = ξ; N 3 = η.
At any time-step n+1/2, the displacement and strain fields are interpolated from nodal displacement values d, by using linear shape functions N

u h n+1/2 | Ω e = N d h n+1/2 | Ω e ; h n+1/2 | Ω e = Bd h n+1/2 | Ω e ; B = ∇ s N . (3.24)
To interpolate the stress field, the lowest order Raviart-Thomas basis function (RT 0 ) is utilized, with the normal basis functions given as

Φ1 (ξ, η) = √ 2 ξ η ; Φ2 (ξ, η) = ξ -1 η ; Φ3 (ξ, η) = ξ η -1 . (3.25)
For any midpoint ĝj along edge j ∈ [1, 2, 3], we can show that the independence property of normal basis function

Φi ĝj • ni = 1; if j = i. 0; if j = i. (3.26)
where ni is the normal vector along edge i ∈ [1, 2, 3], which is opposite node i as shown in Fig. 3

-→ Φ i (x) = F ( Φi )(x) := |e k | |ê k | J T |J T | Φi (x); Φ i (g j ) • n i = 1; if j = i. 0; if j = i. (3.27)
The independence property of basis function in parent space is maintained by the Piola transformation to physical space as in Eq. (3.27b). The stress field is now written as a combination of these independent basis functions. Let us consider two faces in x and y direction of a infinitesimal element, on which the stress components could be combined to define the stress vector at the corresponding face

σ h 11 σ h 12 = σ h 1 | Ω e = τ 1 1 Φ 1 + τ 2 1 Φ 2 + τ 3 1 Φ 3 ; σ h 21 σ h 22 = σ h 2 | Ω e = τ 1 2 Φ 1 + τ 2 2 Φ 2 + τ 3 2 Φ 3 .
(3.28)

Casting the above terms in the correct order in Eq. (3.22a), the stress interpolation can thus be formulated in matrix form as follows

σ h n+1/2 | Ω e = Sτ h n+1/2 | Ω e . (3.29)
By using the results in Eq. (3.28) and Eq. (3.27b), the traction due to the stress tensor at the edge i can be expressed as follows

t h i = σ h • n h i = σ h 1 • n h i σ h 2 • n h i ≡ τ i 1 τ i 2 . (3.30)
Thus, the couple (τ i 1 , τ i 2 ) in Eq. (3.28) actually represents the components of the traction vector at the edge i. Since τ i 1 , τ i 2 are also the degrees of freedom, chosen same on a given edge which is shared by a couple of neighboring triangular elements, this type of interpolation guarantees the continuity of a normal traction vector on an edge between two neighboring elements.

Discrete weak form

With the chosen interpolation of displacement and stress fields, the discrete weak forms in Eq. (3.18) and Eq. (3.19) for a typical element Ω e can be written in a matrix notation as follows In the same manner in Chapter 2, the visco-plasticity is implemented with the yield function using von Mises stress tensor norm J 2 . The set of linearized discrete governing equations in a matrix notation is now established with unknown nodal displacement and stress fields at each iteration:

G e u (u, σ; δu) =δd T Ω e ρN T N a dV + Ω e B T Sτ dV - Ω e N T f dV - ∂Ω e t N T t dA = 0; G e σ (u, σ; δσ) =δτ T Ω e S T Bd dV - Ω e S T DSτ dV = 0. ( 3
A nel e=1 η M e F e,T F e -H e n+1/2   ∆d e,(i) n+1/2 ∆τ e,(i) n+1/2   = -A nel e=1   r e,(i) u,n+1/2 (•) h e,(i) n+1/2 (•)   . (3.36)
Having obtained the solution at each iteration by using of Newton-Raphson iterative method, the stress and displacement nodal values are updated

τ (i+1) n+1/2 = τ (i) n+1/2 + ∆τ (i) n+1/2 ; d (i+1) n+1/2 = d (i) n+1/2 + ∆d (i) n+1/2 .
(3.37)

Energy conserving scheme

In this section, an energy conserving (EC) scheme for the linear elastodynamics problem, is developed to guarantee the stability of the computation over long period. The mid-point time integration scheme is employed to maintain the second-order accurate O(∆t 3 ) solution, see [START_REF] Artioli | Second-order accurate integration algorithms for von-mises plasticity with a nonlinear kinematic hardening mechanism[END_REF]. Accordingly, the nodal velocity and acceleration at the time step t n+1/2 can be formulated in the form of nodal displacement and velocity

v n+1/2 = (d n+1 -d n )/∆t; a n+1/2 = (v n+1 -v n )/∆t. (3.38)
The increment of displacement field within an interval of time step is computed by using the mid-point approximation:

d n+1 -d n = 1 2 ∆t(v n+1 + v n ). (3.39)
This equation can be formed in an alternative expression showing the simple update of a position vector:

d n+1 = d n + u; u := 1 2 ∆t(v n+1 + v n ). (3.40)
The stress state at t n+1/2 is selected in the following algorithmic form, which is left open to be explained shortly:

τ alg n+1/2 := 1 2 (τ n+1 + τ n ). (3.41)
By choosing a test displacement vector δd T = d n+1 -d n , a work done by both the internal f int and external f ext can be elaborated from the weak form of balance equation in Eq. (3.32b)

(d n+1 -d n ) T f int = (d n+1 -d n ) T f ext .
The internal work includes two components including kinetic and potential energies

(d n+1 -d n ) T f int = (d n+1 -d n ) T Ω B T Sτ n+1/2 dV +(d n+1 -d n ) T Ω ρN T N a n+1/2 dV.
The work done by inertia force can be simplified thanks to results obtained in Eq. (3.38b) and Eq. (3.39)

(d n+1 -d n ) T M a n+1/2 = 1 2 (v n+1 + v n ) T M (v n+1 -v n ) = 1 2 v T n+1 M v n+1 - 1 2 v T n M v n = K n+1 -K n , (3.42)
where K n+1 is the kinetic energy of the system at time t n+1 . It is noted that the mass matrix M is symmetric, so this identity holds:

v T n+1 M v n -v T n M v n+1 = 0.
The work done by the potential force can also be expressed as a combination of 2 groups

(d n+1 -d n ) T F T τ n+1/2 = 1 2 (d n+1 -d n ) T F T (τ n+1 + τ n ) = 1 2 d T n+1 F T τ n+1 - 1 2 d T n F T τ n + 1 2 d T n+1 F T τ n - 1 2 d T n F T τ n+1 L .
(3.43) Since the scalar L has variables from both time step t n and t n+1 , further simplification of L should be elaborated. Via exploiting the result in Eq. (3.34b) with a test stress function chosen as δτ T = τ n+1 -τ n , the following identity is derived based on the fact that the matrix H is always a symmetric matrix:

(τ n+1 -τ n ) T F d n+1/2 -Hτ n+1/2 =0; ⇒ 1 2 d T n+1 F T τ n+1 - 1 2 τ T n+1 Hτ n+1 - 1 2 d T n F T τ n + 1 2 τ T n Hτ n = 1 2 d T n+1 F T τ n - 1 2 d T n F T τ n+1 L .
(3.44) With the last identity on hand, Eq. (3.43) can further be rewritten as

(d n+1 -d n ) T F T τ n+1/2 = 1 2 d T n+1 F T τ n+1 - 1 2 d T n F T τ n + 1 2 d T n+1 F T τ n+1 - 1 2 τ n+1 Hτ n+1 - 1 2 d T n F T τ n + 1 2 τ n Hτ n = d T n+1 F T τ n+1 - 1 2 τ T n+1 H e τ n+1 P n+1 -(d T n F T τ n - 1 2 τ T n H e τ n ) Pn = P n+1 -P n , (3.45)
where P n+1 is the potential energy of the system at time t n+1 . From the above derivation, one can express the balance of the internal energy and the external work in following equation

(K n+1 + P n+1 ) E n+1 -(K n + P n ) En = (d n+1 -d n ) T f ext n+1/2 ∆W n+1/2 ; ⇔E n+1 -E n =∆W. (3.46)
We can conclude that the algorithmic constitutive equation Eq. (3.41) would conserve the total energy E n+1 = E n + ∆W for any bounded external loading. This property is viewed as "unconditionally stable", see [START_REF] Hughes | Unconditionally stable algorithms for quasistatic elasto/visco-plastic finite element analysis[END_REF][START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF].

Numerical examples

Several numerical simulations are presented in this section to illustrate the performance of the enhanced triangular element and the corresponding time-stepping integration scheme. All the computations are programmed as so-called user-defined element and subroutines in FEAP v8.4 (see [START_REF] Taylor | Feap-finite element analysis program[END_REF]). The output data is then postprocessed via Matlab scripts for plotting.

Patch test

The patch test can reveal any spurious modes which may exist in the element. In this test, model 1 (1 by 1 in size) is denoted by a domain with a homogeneous material while model 2 (2 by 1 in size) is denoted by a domain with a non-homogeneous material including 2 different properties as shown in Fig. 3.2 (see [START_REF] Auricchio | Isogeometric collocation for elastostatics and explicit dynamics[END_REF]). In model 1, the regular and irregular mesh corresponds to case 1 and case 2, respectively. The same material is applied for either mesh where elastic modulus E = 200 and ν = 0.499 for incompressible condition. The numerical displacement at point A is (0.01, -9.95.10 -3 ), matching with the 'exact' solution. This solution is obtained by using σ 22 ≈ 0, the details of stress-strain relation are given in Appendix A

∆l x l x = 11 = (1 -ν 2 ) E σ 11 ⇒ ∆l ux(A) = 1 -ν 2 E σ 11 p =0.01; ∆l y l y = 22 = (-ν -ν 2 ) E σ 11 ⇒ ∆l uy(A) = -ν -ν 2 E σ 11 p = -9.95.10 -3 .
(3.47)

The displacement field in case 1 is shown in Fig. 3.2(a-d). In both cases, the vertical component of traction τ y vanishes at a midpoint of any given edge while the horizontal counterpart's magnitude τ x on vertical edge dominates those on the horizontal and inclined edges, see Fig. 3.2(e,f). The stress σ 11 is constant as expected over the entire domain for both cases. In model 2, the material properties of 2 regions are chosen as follows: (E 1 = 768, ν = 0.2) and (E 2 = 1000, ν = 0.25). The numerical displacement at point A computed by the eCST is (2.1875.10 -3 , -3.125.10 -4 ), matching with a solution from a higher order approximation by an isogeometric element, see [START_REF] Auricchio | Isogeometric collocation for elastostatics and explicit dynamics[END_REF]. The corresponding displacement field is shown in Fig. 3.4(a,b). The vertical component of traction τ y vanishes at a midpoint of any given edge, while the horizontal counterpart's magnitude τ x on vertical edge dominates those on the horizontal and inclined edges, see Fig. 3.4(c). The stress σ 11 is constant as expected over the entire domain for both cases. Finally, we conclude that the proposed element passes the patch test.

Higher order patch test

A simple beam 10 by 1 under a force couple load p = 1 applied at each end is subjected to a pure bending state, which represents a higher order patch test as shown in Fig. 3.5(a). The material properties are chosen as: elastic modulus E = 100 and ν = 0. The deflection at point M(5,1) computed by the proposed element (eCST) is compared against the counterpart obtained by the conventional constant stress/strain triangle element (CST) and four-node quadrilateral element (Q4) under several mesh configurations. The 'exact' solution is computed by a very fine mesh of nine-node quadrilateral element (Q9). As shown in Fig. 3.5(b), the rate of convergence and accuracy of the eCST element is better than all the others. Detailed results are given in two cases. A coarse mesh with 6 by 1 element in each direction is denoted in case 1, while a finer mesh with 12 by 6 element is denoted in case 2. The displacement field in case 1 is shown in Fig. 3.6(a,c). The finer mesh in case 2 yield the better contour of displacement field as in Fig. 3.6(b,d). The deformed shape is shown in Fig. 3.6(e,f). The illustration of traction field is depicted in a finer mesh. The horizontal and vertical vertical components of traction (τ x ,τ y ) for case 2 are shown in Fig. 3.7(a,b), respectively. The values are normalized with respect to the corresponding maximum. The traction vector on midpoint of each edge is given in Fig. 3.7(c). It is expected that there are 2 different zones including tension and compression as depicted in Fig. 3.7(d). We conclude that the proposed element passes the higher order patch test. Detailed results are given in two cases. A coarse mesh with 4 by 1 element in each direction is denoted in case 1, while a finer mesh with 12 by 6 element is denoted in case 2. The displacement field in case 1 is shown in Fig. 3.9(a,c). The finer mesh in case 2 yields the smooth contours of the displacement field, as in Fig. 3.9(b,d). The deformed shape is shown in Fig. 3.9(e,f). The illustration of traction field is depicted with finer mesh. The horizontal and vertical components of traction (τ x ,τ y ) for case 2 are shown in Fig. 3.10(a,b), respectively. The values are normalized with respect to the corresponding maximum. The traction vector on the midpoint of each edge is given in Fig. 3.10(c). It is expected that there are two different zones including tension and compression as depicted in Fig. 3.10(d).

Cook's membrane

Cook's membrane test proposed in [START_REF] Cook | Improved two-dimensional finite element[END_REF] presents the shear dominated behavior and also displays the effects of mesh distortion, as shown in Fig. 3.11. The material properties are chosen as: elastic modulus E = 1, ν = 0.33 for quasi-incompressible condition and ν = 0.499 for an incompressible condition. A total vertical load p = 1 applied at the right end. The deflection at point 3 [START_REF] Dolbow | Modeling fracture in mindlin-reissner plates with the extended finite element method[END_REF][START_REF] Ghanem | Probabilistic learning and updating of a digital twin for composite material systems[END_REF] computed by the proposed element (eCST) is compared against the counterpart obtained by the conventional constant stress/strain triangle element (CST) and four-node quadrilateral element (Q4) with several mesh configurations. The 'exact' solution is computed by a very fine mesh of the nine-node quadrilateral element (Q9). As shown in Fig. 3.12, the rate of convergence and accuracy of the eCST element is better, especially for the incompressible case with ν = 0.499. Detailed results are given in two cases with ν = 0.33 to represent the incompressible condition. A coarse mesh with 2 by 2 elements in each direction is denoted as case 1, while a finer mesh with 8 by 8 elements is denoted as case 2. The displacement field in case 1 is shown in Fig. 3. 13(a,c). The finer mesh in case 2 yields the smooth contours of displacement field, as shown in Fig. 3.13(b,d). The deformed shape is shown in Fig. 3.13(e,f). The illustration of traction field is depicted with finer mesh. The horizontal and vertical components of traction (τ x ,τ y ) for case 2 are shown in Fig. 3.14(a,b), respectively. The values are normalized with respect to the corresponding maximum. The traction vector on the midpoint of each edge is given in Fig. 3.14(c). It is expected that there are 2 different zones including tension and compression as depicted in Fig. 3.14(d).

Dynamic loadings

A simple cantilever beam with width 10 and height 1, as shown in Fig. 3.16(a), is selected to simulate a vibration under a dynamic pulse. There are two models with different load setting: axial and shear pressure. The same fine mesh (10 by 2 elements in each direction) is generated in both cases to allow a wave propagation in the beam under introduced pulse. The properties of both model materials are chosen as: elastic modulus E = 200, Poisson's ratio ν = 0.33 and ρ = 0.785. The time step chosen for the simulation is ∆t = 10 -i , i ∈ [1,2,3,4]. A triangle pulse with a peak of pressure p 0 = 1 is introduced to both models as shown in Fig. 3.16(b), which leads subsequently to free vibrations of the system. These loading scenarios are chosen as shown in Table 3.1. The energy conserving scheme (EC scheme) is employed to compute the evolution of displacement and traction fields. Regarding model 1, the amplitude of displacement in both directions at node 123 is shown in Fig. 3.17(a,c,e,g). It is noted that the displacement in the horizontal direction is larger than its counterpart in vertical direction, since the axial pulse is introduced. As shown in Fig. 3.17(b,d,f,h), the similar trend can be observed for the first several steps in the traction at node 82. As shown in Fig. 3.18, the results show that the total internal energy of the system is conserved exactly under free vibration in every case in model 1. It is noted that longer duration of loading leads to higher total energy gained by the system. These cases are fully under the elastic regime. As shown in Fig. 3.19, the results show that the total internal energy of the system is conserved exactly under the free vibration in every case for model 2. It is noted that a longer duration of loading leads to higher total energy gained by the system. These cases are fully under the elastic regime. Regarding model 2, the amplitude of displacement in both direction at node 123 is shown in Fig. 3.20(a,c,e,g). Since the shear pulse is introduced, the displacement in the vertical direction is larger than its counterpart in the horizontal direction. The traction at node 82 is shown in Fig. 3.20(b,d,f,h).

To implement the visco-plasticity, the same material is employed with additional properties including yield stress σ f = 0.5 and viscosity coefficient η = 1/10 -8 . The results from case 3 in model 1 are computed for the visco-plastic simulation, as shown in Fig. 3.21. The corresponding nodes and elements can be found in Fig. 3.16(a).

4

Mixed variational formulation for the visco-plasticity plate bending problem under dynamics

In this chapter, a regularized variational formulation is derived in the same manner as presented in Chapter 3 for the visco-plasticity plate bending problem. Then the displacement is interpolated by low-order polynomial functions, while the stress-resultants (moment and shear forces) are approximated by the lowest-order Raviart-Thomas vector space. Afterwards, the energy and decaying conserving scheme are configured for the dynamics computation. The performance of the proposed element is validated via several numerical simulations.

Regularized variational formulation

We first briefly review the conventional model of the Reissner-Mindlin plate theory.

Then the hybrid-stress variational formulation of a governing equation is given with the moment and shear force fields introduced as additional independent variables by appealing to the Hellinger-Reissner method, and its regularized format, which is the most suitable for the Raviart-Thomas discrete approximation, see [START_REF] Nguyen | Hybrid-stress triangular finite element with enhanced performance for statics and dynamics[END_REF] for the elastodynamics problem.

Let Ω in Fig. 4.1 be a bounded flat body in R 3 representing a plate of thickness t, in which the thickness t is small compared to the other two dimensions, see Eq. (4.1a). The corresponding piece-wise smooth boundary ∂Ω, including imposed deflection boundary ∂Ω u and imposed stress resultants boundary ∂Ω t , satisfying Eq. (4.1b). The xy plane of chosen coordinate system is placed in the plate mid-surface (denoted as A), and the surface loading p(x, y) is applied normal to the mid-surface, as shown in Fig. 4.1. The Reissner-Mindlin plate theory assumes that a straight fiber normal to the mid-surface before applying the loading remains straight but not necessarily normal to the plate mid-surface. We choose the corresponding rotation components of this fiber, denoted as φ x and φ y .

Conventional Reissner-Mindlin plate model

In this section, we reformulate the classical Reissner-Mindlin plate model, where the governing equations are recast in tensor notation, which is needed to provide the regularized version starting from 3D continuum mechanics model. The Reissner-Mindlin plate kinematics considers that for any point C (see Fig. 4.1) placed outside of the plate mid-surface, the displacement field components u x , u y and u z can be expressed in terms of the lateral deflection w at the plate mid-surface and the fiber rotation vector φ = (φ x , φ y ), which can be written as follows

u x = -zφ x (x, y); u y = -zφ y (x, y); u z = w(x, y). (4.2)
The sign convention of the rotation vector θ = (θ x , θ y ) about the x and y directions is shown in Fig. 4.1, (see [START_REF] Ibrahimbegovic | Quadrilateral finite elements for analysis of thick and thin plates[END_REF]) hence the relation between two rotation vectors components is given as follows

φ x = -θ y ; φ y = θ x . (4.
3)

The gradient of the displacement field can be written in a tensor notation, or rather as corresponding 3D matrix representation in the Euclidean setting (e.g. [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]), conjugate to stress tensor σ assumed to be non-symmetric. In agreement with Eq. (4.2), we can write

∇u =    u x,x u x,y u x,z u y,x u y,y u y,z u z,x u z,y u z,z    =    -zφ x,x -zφ x,y -φ x -zφ y,x -zφ y,y -φ y w ,x w ,y 0    =    zθ y,x zθ y,y θ y -zθ x,x -zθ x,y -θ x w ,x w ,y 0   
(4.4) This allows to redefine the generalized strain measure for the Reissner-Mindlin plate model in terms of the curvature tensor and the shear strain. The corresponding matrix representation of the curvature tensor κ of the plate Reissner-Mindlin model, can be written as

κ = κ xx κ xy κ yx κ yy = θ y,x 1 2 (-θ x,x + θ y,y ) 1 2 (-θ x,x + θ y,y ) -θ x,y . (4.5)
The curvature tensor κ and the shear strain γ are work-conjugates to the bending moment tensor M and the shear force Q, respectively. The latter are internal forces in the plate expressed in terms of stress resultants, which are defined by integrating the normal and shear stresses through the thickness of the plate, leaving out the component σ zz . However, we will not enforce the symmetry of the bending moment tensor, keeping the corresponding moment equilibrium equation in the strong form (more details in next section), so that we can write: In numerical computations we further switch to the corresponding vector representation of the defined tensors κ and M . Firstly, the curvature tensor κ components for the Reissner-Mindlin plate model can be written as

M = M xx M xy M yx M yy =   t/2 -t
κ = κ xx , κ yy , κ xy , κ yx T = θ y,x , -θ x,y , 1 2 (-θ x,x + θ y,y ), 1 2 (-θ x,x + θ y,y ) T . (4.7)
The work-conjugated bending moment tensor M can also be rewritten in the vector form as follows

M = [M xx , M yy , M xy' , M yx ] T = t/2 -t/2 σ xx z dz, t/2 -t/2 σ yy z dz, t/2 -t/2 σ xy z dz, t/2 -t/2 σ yx z dz T . (4.8)
The shear strain components for the Reissner-Mindlin plate are placed in the vector γ written as

γ = γ xz , γ yz T = w ,x + θ y , w ,y -θ x T . ( 4.9)
The work-conjugated shear force Q can also be written in the vector form as follows

Q = [Q xz , Q yz ] T = t/2 -t/2 σ xz dz, t/2 -t/2 σ yz dz T . (4.10)
For the isotropic linear elastic constitutive behavior of a plate material, the constitutive relation between moment tensor components and corresponding curvature tensor components can be written as

       M xx M yy M xy M yx        = D      1 ν 0 0 ν 1 0 0 0 0 (1 -ν)/2 0 0 0 0 (1 -ν)/2      C b        κ xx κ yy κ xy κ yx        → M = C b κ, ( 4.11)
where

D = Et 3 12(1 -ν 2 )
is the flexural rigidity of the plate with the Young's modulus E and the Poisson's ratio ν, whereas the bending elastic matrix C b is derived from the fourth order tensor C and integrated over thickness t of the plate. The constitutive equation for the shear force versus the shear strain relation can be written in the vector notation as follows

Q x Q y = kGt 1 0 0 1 C s γ xz γ yz → Q = C s γ, (4.12)
where

G = E 2(1 + v)
is shear modulus and k is the correction factor. The usual value k = 5/6 is used for the homogeneous plate and the parabolic distribution of the transverse shear stress in the z-direction. We note in passing that the compliance matrices for bending and shear constitutive relations are computed simply as the inverse of corresponding constitutive matrices D b = C b,-1 and D s = C s,-1 . We further assume a more general constitutive behavior of the plate material in terms of the visco-plasticity model [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF], where we regroup all stress resultants components within the generalized stress resultant vector s, and also regroup the corresponding strain components in the vector χ; both are simply defined as follows

χ = [κ xx , κ yy , κ xy , κ yx , γ xz , γ yz ] ; s = [M xx , M yy , M xy , M yx , Q xz , Q yz ] . (4.13)
For the case of visco-plasticity at small strains, the strain components can be decomposed additively into elastic and visco-plastic strain measures for the plate with χ = χ e + χ vp . This implies that both the curvature tensor and shear strain vectors are assumed to be split additively into elastic strains (κ e , γ e ) and visco-plastic strains (κ vp , γ vp ) parts. Hence, we can enforce the following relations κ = κ e + κ vp and γ = γ e + γ vp .

Stress resultant-based variational equations

The strong form of the governing equations for the plate model is rewritten in regularized form. We follow closely the work in [START_REF] Hughes | On drilling degrees of freedom[END_REF] for the 3D continuum problem, which is proposed to accommodate any most convenient discrete approximations. From this starting point, we switch to the hybrid-stress version where only stress and displacement fields are retained as independent fields. We then use the development in previous section to recast such formulation in terms of stress resultants. In particular, this allows to provide the moment vector and shear continuity across the edge between the adjacent plate elements, which is of special interest for dynamics to capture (smooth) wave propagation. The regularized formulation for plate problem is recast for plate bending problem as follows

Π (w, θ, M , Q) = - 1 2 A symm M • D b symm M + Q • D s Q dA + A (symm M • κ + Q • γ) dA + 1 2 -1 A 12 t 3 | skew M | 2 dA -Π ext (u).
(4.14) The variational equations corresponding to the functional in Eq. (4.14) above for plate bending problem can be written as follows

G u (w, θ, M , Q; δw, δθ) = A δκ • M dA + A δγ • Q dA - ∂Ωt δwp dA - Ω δwρ ẅ dV G ext (•) = 0; G M ,Q (w, θ, M , Q; δM , δQ) = A δM • κ - Db M + δQ • (γ -D s Q) dA = A δs • χ -C -1 s dA = 0, (4.15
) in which the inverse elasticity matrix C -1 = Db ⊕ D s . We also note that the external virtual work G ext (w; δw) is formulated in the spirit of d Alembert principle with an additional term from the inertia force for the dynamics problem. In the case of visco-plastic constitutive behavior of the plate material, Eq. (4.15b) is rewritten by exploiting the additive decomposition of strain fields into elastic and visco-plastic components

G M ,Q (w, θ, M , Q; δM , δQ) = A {δM • κ -κ vp + Db M +δQ • [γ -(γ vp + D s Q)]} dA = A δs • χ -χ vp + C -1 s dA = 0. (4.16)
The corresponding form of updated elastic compliance matrix Db related to the bending moment field in Eq. (4.15) is presented in Appendix B. The final point in defining the nonlinear behavior with visco-plasticity concerns the selection of the yield criterion, which triggers the evolution of the visco-plastic deformation. For the present plate bending problem, we follow the choice made in [START_REF] Ibrahimbegovic | An efficient implementation of stress resultant plasticity in analysis of reissner-mindlin plates[END_REF] by selecting the yield function φ(s) as a quadratic form in stress resultant, which can be written as follows

φ (s) = 1 2 s • As -1; A = 1 m 2 y P ⊕ 1 q 2 y I 2 , (4.17)
where I 2 is the 2 × 2 identity matrix. The chosen yield criterion is obtained by recasting the classical von Mises yield criterion in the stress resultant form, which explicitly defines the deviatoric matrix in the stress resultant form P , see [START_REF] Ibrahimbegovic | An efficient implementation of stress resultant plasticity in analysis of reissner-mindlin plates[END_REF] for details. The corresponding ultimate bending moment m y and ultimate shear force q y , each defined as function of the yielding stress σ y of a background continuum problem are

m y = σ y t 2 4 ; q y = σ y t √ 3 . (4.18)
Note that for the chosen visco-plasticity model, the yield function φ(s) accepts all values of stress in the space of real number R. This is the crucial choice for development of the second-order time-integration scheme, where both global motion equations and local visco-plastic strain evolution equations have to be integrated simultaneously [START_REF] Nguyen | Visco-plasticity stress-based solid dynamics formulation and time-stepping algorithms for stiff case[END_REF] (contrary to the operator split procedure of rate-independent plasticity; e.g. [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]). The evolution of visco-plastic generalized strain χ vp is obtained by appealing to the penalized version of maximum plastic dissipation [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF], resulting with the standard format that involves the derivatives of the corresponding yield function φ(s). The visco-plastic flow rule is triggered by positive values of the yield function φ(s) > 0, whereas the elastic regime is characterized with non-positive values of the yield function φ(s) ≤ 0, which can be written as follows

φ(s) ≤ 0 ⇒ χvp = 0; φ(s) > 0 ⇒ χvp = 0. (4.19)
Here, we provide more detailed developments of these ideas. First, we can appeal to the second law of thermodynamics in terms of the Clausius-Duhem inequality [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF], in order to provide the corresponding visco-plastic dissipation D vp in terms of plate model stress resultants

D vp = s • χ -ψe = s • χvp ≥ 0, (4.20)
where s denotes the stress resultants and χ vp denotes generalized visco-plastic strains for the Reissner-Mindlin plate model, which contains the visco-plastic curvature and shear strains. For our choice of the visco-plasticity constitutive model, all the stress resultant values are admissible, but those outside the elastic domain are penalized by an additional penalty function P (φ). For simplicity, we choose this penalty function as a quadratic form in φ(s), which will penalize more the stress states far outside the elastic domain

P (φ(s)) =      1 2 φ(s) 2 ; φ(s) > 0. 0; φ(s) ≤ 0. (4.21)
The principle of maximum visco-plastic dissipation [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF] can then be applied to obtain the corresponding evolution equation for the visco-plastic strain

∂P (φ) ∂φ =< φ >, (4.22)
where < • > is the Macauley parenthesis (e.g. [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]), which is defined as

< φ >:= (φ(s) + |φ(s)|)/2 = φ(s); φ(s) > 0. 0; φ(s) ≤ 0. (4.23)
The penalty function is co-jointed with the factor 1/η, in which the scalar η is the viscosity coefficient characterized by the corresponding material. The modified visco-plastic dissipation D vp η is now written as follows

-D vp η [s] = -s • χvp + 1 η P (φ(s)). (4.24)
The evolution of visco-plasticity can be derived by means of the principle of maximum visco-plastic dissipation. From all possible states s * satisfying the yield criterion φ (s * ) > 0, the stress resultant state s that maximizes the visco-plastic dissipation D vp η is the solution. This problem can also be written as constrained minimization problem

D vp η (s) = max φ(s * )>0 D vp η (s * ) ⇔ -D vp η (s) = min φ(s * )>0 -D vp η (s * ) . ( 4.25)
The Kuhn-Tucker optimality conditions for the minimization problem defined in the above equation lead to the evolution equation of the visco-plastic resultant strain

0 = ∂ -D vp η (s) ∂s := ∂[-s • χvp + 1 η P (φ(s))] ∂s = -χvp + 1 η < φ > ∂φ ∂s . ( 4

Numerical implementation

In this section, we present the discrete approximation for displacement (including lateral deflection and rotations), moment and shear force fields. The proposed element is developed based on the three-node triangle element. The linear discretization is employed for both lateral deflection and rotation fields [START_REF] Hughes | The finite element method: Linear static and dynamic finite element analysis[END_REF]. Moreover, the lowest order Raviart-Thomas vector space (RT 0 ) [START_REF] Ervin | Computational bases for rtk and bdmk on triangles[END_REF] is introduced to approximate moment and shear force fields. Having defined the discrete approximation of all independent fields, we can obtain the corresponding semi-discretized weak form. For dynamics framework studied herein, we then carry on to present the time-integration scheme, which is used to fully discretize the problem. Here, we propose a scheme that can conserve the total energy (in a free vibration) and its extension that decays the en-ergy of higher frequency modes. The resulting discretization of the problem results with a linear set of motion equations at time step t n+1/2 for elastic plate behavior, and accompanying evolution equations for the visco-plastic plate behavior, to be solved at each step of computations.

Discrete approximation for deflection and rotation fields

The following development is based on the previous development of quadrilateral plate elements applicable to the analysis of the thick Reissner-Mindlin plate, see [START_REF] Ibrahimbegovic | Quadrilateral finite elements for analysis of thick and thin plates[END_REF]. This development is adapted to a triangle element with low-order shape functions N employed for interpolation. More precisely, the deflection and the rotation fields are constructed by generalizing the Timoshenko beam interpolations, see [START_REF] Ibrahimbegovic | Quadrilateral finite elements for analysis of thick and thin plates[END_REF]. The discrete approximation of the rotation field θ = (θ x , θ y ) at a given point in a discrete element Ω e is written in terms of nodal rotation θ i = (θ i x , θ i y )

θ x θ y = 3 i=1 N i (ξ, η) θ i x θ i y → θ = N θ θ d . (4.28)
The corresponding curvature components are set in a matrix form which is derived by substituting Eq. (4.28) into Eq. (4.7) to obtain

κ = 3 i=1 B i (ξ, η)d i → κ = B κ d, (4.29)
where

B i (ξ, η) =      0 0 ∂N i (ξ, η)/∂x 0 -∂N i (ξ, η)/∂y 0 0 -1 2 ∂N i (ξ, η)/∂x 1 2 ∂N i (ξ, η)/∂y 0 -1 2 ∂N i (ξ, η)/∂x 1 2 ∂N i (ξ, η)/∂y      ; d i =    w i θ i x θ i y   .
(4.30)

The discrete approximation of the lateral deflection w is written in terms of nodal values w i as follows

w = 3 i=1 N i (ξ, η)w i → w = N w w d . (4.31)
In the same manner, the corresponding acceleration of the lateral deflection field ẅ is computed from nodal values of acceleration a as follows ẅ =

3 i=1 N i (ξ, η)a i → ẅ = N w a. (4.32)
The discrete approximation of the lateral deflection w and acceleration ẅ are employed to compute the virtual work of external forces. To alleviate the shear locking, a bi-linear distribution for the assumed shear strain is proposed (e.g. [START_REF] Ibrahimbegovic | Quadrilateral finite elements for analysis of thick and thin plates[END_REF]), which can be written as follows

γ xz γ yz = 1 3 γ 1 + γ 2 + γ 3 → γ = N γ γ d . (4.33)
where the nodal shear strains γ i (see Fig. 4.2) are computed in the manner that ensures the constant shear strain distribution along each element edge, and thus the continuity between adjacent plate elements. By enforcing the equality between the projection of the nodal shear strain along the edge of two neighbors, the corresponding matrix of nodal shear strain is formed as a function of nodal unknown matrix d. The form of nodal shear strain γ is written as a function of lateral deflection and rotation fields as follows

γ = N γ γ d = N γ Λd = B γ d. (4.34)
With the chosen interpolations in Eq. (4.29) and Eq. (4.33), the proposed element inherits the enhanced performance with the ability to overcome the shear locking phenomena from the previous developed quadrilateral element [START_REF] Ibrahimbegovic | Quadrilateral finite elements for analysis of thick and thin plates[END_REF].

Discrete approximation for moment and shear force fields

As shown in Fig. 4.3, the element corner nodes 1, 2 and 3 are assigned the degrees of freedom for the deflection w and rotation θ i = (θ i x , θ i y ) fields. Meanwhile the nodal degree of freedom corresponding to bending moment m i = (m i x , m i y ) and shear force q i fields are located at mid-edge nodes 4, 5 and 6. For the corresponding parent element, the normal basis functions Φ i (ξ, η) of the lowest order Raviart-Thomas vector space (RT 0 ) are defined as follows

Φ1 (ξ, η) = √ 2 ξ η ; Φ2 (ξ, η) = ξ -1 η ; Φ3 (ξ, η) = ξ η -1 . (4.35)
For any midpoint ĝj along edge j ∈ [1, 2, 3], the orthogonality property of normal basis functions can be expressed as follows

Φi ĝj • ni = 1; if j = i. 0; if j = i. , (4.36)
where ni is the unit exterior normal at edge i ∈ [1, 2, 3], placed opposite node i as shown in Fig. 4.3. The basis function Φi (ξ, η) is mapped back to physical coordinate space by the Piola transformation [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF] represented with an affine function F . Since one common edge between two elements in physical space can be mapped into two different edges (with different lengths) in parent coordinate, the terms |e k | and |ê k | (defining length of considered edges in the physical and parent space respectively) have to be added in order to maintain Eq. (4.36) in the physical coordinate, see [START_REF] Ervin | Computational bases for rtk and bdmk on triangles[END_REF]. Hence, the basis function Φi at edge i in the physical coordinate system is now written as follows

Φi

-

→ Φ i (x) = F ( Φi )(x) := |e k | |ê k | J T |J T | Φi (x); Φ i (g j ) • n i = 1; if j = i. 0; if j = i. , (4.37)
where the notation J T and |J T | are the Jacobian matrix and its determinant of the corresponding affine function. The independence property of basis function in parent space is maintained by the Piola transformation to physical space as indicated in Eq. (4.37b).

The moment field discrete approximation can now be written as a function of the independent basis functions

M h xx M h xy = M h x | Ω e = m 1 x Φ 1 + m 2 x Φ 2 + m 3 x Φ 3 ; M h yx M h yy = M h y | Ω e = m 1 y Φ 1 + m 2 y Φ 2 + m 3 y Φ 3 , (4.38)
or expressed in corresponding more compact form by rearranging terms in the right order to write

M h n+1/2 | Ω e = Sm h n+1/2 | Ω e . ( 4.39)
The matrix S denotes the interpolation matrix that contains the basis functions for moment field. From the results in Eq. (4.37b) and Eq. (4.38), the projection of the moment tensor upon a unit normal vector n i at an edge i is written as

m i,h = M h x • n h i M h y • n h i ≡ m i x m i y . (4.40)
Such a projection in Eq. (4.40) above is further referred to as the moment vector m i,h (m i x , m i y ), since it is equivalent to the traction or stress vector obtained by projecting the stress tensor (through Cauchy principle) upon the unit vector of exterior normal n i . The moment vector m i,h (m i x , m i y ) is called the projection vector of moment tensor at the edge i, and constitutes the nodal degrees of freedom of bending moment interpolations. These nodal degrees of freedom (m i x , m i y ) are the same on a given edge i, which is shared by two neighboring triangular plate elements. Thus, the chosen discrete approximation of the bending moment guarantees the continuity of moment components orthogonal and parallel to element edge separating two adjacent plate elements. In other words, we obtain the continuity of the moment vector when crossing the plate element boundaries, which is of great interest for dynamics to ensure a smooth wave propagation.

The discrete approximation of the shear force field can be written in the same manner to ensure the continuity as follows

Q h xz Q h yz = Q h | Ω e = q 1 Φ 1 + q 2 Φ 2 + q 3 Φ 3 .
(4.41)

The corresponding compact matrix form can be written as follows

Q h n+1/2 | Ω e = Xq h n+1/2 | Ω e . ( 4

.42)

The matrix X denotes the interpolation matrix that contains the basis functions for shear force field. From the results in Eq. (4.37b) and Eq. (4.41), the projection of shear force field on a normal vector n i at an edge i is written as

q i,h = Q h • n h i . (4.43)
This nodal degree of freedom q i is the same on a given edge i, which is shared by a couple of neighboring triangular elements. From Eq. (4.40) and Eq. (4.43), it follows that the lowest order Raviart-Thomas interpolation (RT 0 ) guarantees the continuity of projection m h of moment and q h shear force fields over edges of each element.

Discrete weak form

The governing equations of the weak form are set up at time step t n+1/2 , which is consistent with the choice of second-order energy conserving time integration scheme. To alleviate notation, the index n + 1/2 indicating mid-step is only shown when it is needed. With the chosen interpolation of lateral deflection, rotations, moment and shear force fields, the discrete weak forms in Eq. (4.15) for a typical element Ω e can be written in a matrix notation as follows with the plate load p representing a general lateral distributed load. The mid-point scheme is used to discrete the visco-plastic compliance as follows The set of governing equations to be solved at time step n + 1/2 is written for all unknown parameters as follows The nodal degrees of freedom and visco-plastic strain resultants are updated at each iteration leading to current iterative values

D vp n+1/2 = ∂χ n+1/2 ∂s n+1/2 = C -1 + ∆t 2 ∂β(s n+1/2 ) ∂s n+1/2 =:   D vp,b n+1/2 D vp
A nel e=1     η M e F e,
d (i+1) n+1/2 =d (i) n+1/2 + ∆d (i) n+1/2 ; m (i+1) n+1/2 =m (i) n+1/2 + ∆m (i) n+1/2 ; q (i+1) n+1/2 =q (i) n+1/2 + ∆q (i) n+1/2 ; χ (i),vp n+1/2 =χ vp n + ∆t 2 β (i) n+1/2 .
(4.52)

Energy conserving scheme

If we seek to ensure the robust performance of a time-integration scheme enforcing the computations stability, one very interesting choice is an energy conserving scheme [START_REF] Ibrahimbegovic | Nonlinear dynamics of fexible beams in planar motion: formulation and time-stepping scheme for stiff problems[END_REF]. The energy conserving (EC) time stepping scheme proposed herein is able to guarantee the stability of the computation over long period of time for the linear elastodynamics. At the time step t n+1/2 , the nodal velocity v n+1/2 and the acceleration a n+1/2 can be computed from the nodal deflection and velocity at two subsequent steps

v n+1/2 = (d n+1 -d n )/∆t; a n+1/2 = (v n+1 -v n )/∆t. (4.53)
Here the increment of deflection field within a time step is computed by using the mid-point approximation

d n+1 -d n = 1 2 ∆t(v n+1 + v n ). (4.54)
This equation can be recast in an alternative format showing the simple update of deflection vector

d n+1 = d n + u; u := 1 2 ∆t(v n+1 + v n ). (4.55)
Next, the lateral deflection, rotations, moment and shear force states at time step t n+1/2 are selected in optimal form with respect to proposed algorithm in order to enforce energy conservation

d alg n+1/2 := 1 2 (d n+1 + d n ); m alg n+1/2 := 1 2 (m n+1 + m n ); q alg n+1/2 := 1 2 (q n+1 + q n ).
(d n+1 -d n ) T f int = (d n+1 -d n ) T f ext .
The internal work includes two components including kinetic and potential energies

(d n+1 -d n ) T f int = (d n+1 -d n ) T M a n+1/2 + G T q n+1/2 + F T m n+1/2 .
Furthermore, the work done by inertia force can be simplified thanks to results in Eq. (4.53b) and Eq. (4.54) to reduce this contribution to increment in the kinetic energy

(d n+1 -d n ) T M a n+1/2 = 1 2 (v n+1 + v n ) T M (v n+1 -v n ) = 1 2 v T n+1 M v n+1 - 1 2 v T n M v n = K n+1 -K n , (4.57)
where K n+1 is the kinetic energy at time step t n+1 . We note that the mass matrix M is symmetric, the following identity holds v T n+1 M v n -v T n M v n+1 = 0. Moreover, the work done by the moment field can be decomposed into 2 parts

(d n+1 -d n ) T F T m n+1/2 = 1 2 (d n+1 -d n ) T F T (m n+1 + m n ) = 1 2 d T n+1 F T m n+1 - 1 2 d T n F T m n + 1 2 d T n+1 F T m n - 1 2 d T n F T m n+1 L 1 . (4.58)
Similarly, the work done by the shear force field can be decomposed into 2 parts, increment in the potential energy P q and extra term denoted by L 2

(d n+1 -d n ) T G T q n+1/2 = 1 2 (d n+1 -d n ) T G T (q n+1 + q n ) = 1 2 d T n+1 G T q n+1 - 1 2 d T n G T q n + 1 2 d T n+1 G T q n - 1 2 d T n G T q n+1 L 2 . (4.59)
The L 1 and L 2 scalars are further simplified to separate the terms from time step t n and t n+1 . By appealing to the equation Eq. (4.48b) when using the test functions for moments δm = m n+1 -m n , the following identity can be derived, given that the matrix H is symmetric

(m n+1 -m n ) T F d n+1/2 -Hm n+1/2 = 0; ⇒ 1 2 d T n+1 F T m n+1 - 1 2 m T n+1 Hm n+1 - 1 2 d T n F T m n + 1 2 m T n Hm n = 1 2 d T n+1 F T m n - 1 2 d T n F T m n+1 L 1 . (4.60)
The similar identity related to scalar L 2 can be derived by using Eq. (4.49b), given that the matrix Z is symmetric

(q n+1 -q n ) T Gd n+1/2 -Zq n+1/2 = 0; ⇒ 1 2 d T n+1 G T q n+1 - 1 2 q T n+1 Zq n+1 - 1 2 d T n G T q n + 1 2 q T n Zq n = 1 2 d T n+1 G T q n - 1 2 d T n G T q n+1 L 2 . (4.61)
By exploiting the identity in Eq. (4.60), Eq. (4.58) can be rewritten as an increment in the moment field contribution to the potential energy

(d n+1 -d n ) T F T m n+1/2 = 1 2 d T n+1 F T m n+1 - 1 2 d T n F T m n + 1 2 d T n+1 F T m n+1 - 1 2 m T n+1 Hm n+1 - 1 2 d T n F T m n + 1 2 m T n Hm n = d T n+1 F T m n+1 - 1 2 m T n+1 Hm n+1 P m,n+1 -(d T n F T m n - 1 2 m T n Hm n ) Pm,n
= P m,n+1 -P m,n , (4.62) where P m,n+1 is the potential energy of the system related to the moment field at time step t n+1 . Similarly by using Eq. (4.61), Eq. (4.59) can be rewritten as increment in the shear force contribution to the potential energy

(d n+1 -d n ) T G T q n+1/2 = 1 2 d T n+1 G T q n+1 - 1 2 d T n G T q n + 1 2 d T n+1 G T q n+1 - 1 2 q T n+1 Zq n+1 - 1 2 d T n G T q n + 1 2 q T n Zq n = d T n+1 G T q n+1 - 1 2 q T n+1 Zq n+1 P q,n+1 -(d T n G T q n - 1 2 q T n Zq n )
Pq,n = P q,n+1 -P q,n , (4.63) where P q,n+1 is the potential energy of the system related to the shear force field at time step t n+1 . From the above derivation, one can express the balance of the internal energy and the external work as indicated next

(K n+1 + P m,n+1 + P q,n+1 ) E n+1 -(K n + P m,n + P q,n ) En = (d n+1 -d n ) T f ext ∆W n+1/2 ; E n+1 -E n =∆W n+1/2 .
(4.64)

The symbols E n and E n+1 denote the total energy of the system at time step t n and t n+1 , respectively. We can conclude that the algorithmic constitutive equations, which are proposed in Eq. (4.56), would conserve the total energy E n+1 = E n for any zero contribution of external loading. This property can also be viewed as the proof of "unconditionally stable" scheme, see [START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF]. In the case of visco-plasticity, we can decompose the following matrices into elastic and visco-plastic components

H = H el + H vp ; Z = Z e + Z vp . (4.65)
The equation Eq. (4.60) is now updated as follows

(m n+1 -m n ) T F d n+1/2 -Hm n+1/2 -Y T q n+1/2 = 0; ⇒ 1 2 d T n+1 F T m n+1 - 1 2 m T n+1 H el m n+1 - 1 2 d T n F T m n + 1 2 m T n H el m n + 1 2 m T n H vp m n - 1 2 m T n+1 H vp m n+1 -(m n+1 -m n ) T Y T q n+1/2 D vp,1 n+1/2 = 1 2 d T n+1 F T m n - 1 2 d T n F T m n+1 L 1 . ( 4.66)
The equation Eq. (4.61) is now updated as follows

(q n+1 -q n ) T Gd n+1/2 -Y m n+1/2 -Zq n+1/2 = 0; ⇒ 1 2 d T n+1 G T q n+1 - 1 2 q T n+1 Z el q n+1 - 1 2 d T n G T q n + 1 2 q T n Z el q n + 1 2 q T n Z vp q n - 1 2 q T n+1 Z vp q n+1 -(q n+1 -q n ) T Y m n+1/2 D vp,2 n+1/2 = 1 2 d T n+1 G T q n - 1 2 d T n G T q n+1 L 2 . (4.67)
The correspondingly updated form of the balance of the internal energy and external work, see Eq. (4.64), for the visco-plastic regime coupled with the energy conserving scheme is written as

E n+1 -E n + D vp,n+1/2 = ∆W n+1/2 ; D vp,n+1/2 = D vp,1 n+1/2 + D vp,2 n+1/2 . (4.68)

Energy decaying scheme

The energy decaying (ED) scheme presented here-in is also an approach to control potential instability problems of time-stepping algorithms, but with an additional advantage of providing superior accuracy of stress computations. The latter is established by eliminating the high-frequency modes contribution that can not be resolved on a coarse mesh, especially for a set of stiff equations (e.g. see [START_REF] Ibrahimbegovic | Energy conserving/decaying implicit time-stepping scheme for nonlinear dynamics of three-dimensional beams undergoing finite rotations[END_REF]). The latter implies a large difference between max/min eigenvalues in tangent stiffness, which would require a very fine mesh. In order to handle this case, we will enforce a modified algorithm constitutive equation with the ability of decaying contribution of higher frequencies over each time step. Such algorithmic updates for of all degrees of freedom at time step t n+1/2 are represented by modifying conservative terms, as defined in Eqs. (4.54) and (4.56), with the corresponding dissipative terms

d n+1 -d n := ∆t 2 (v n+1 + v n ) + α∆t(v n+1 -v n ). (4.69)
where α is the dissipation coefficient, which is typically chosen as α ∈ [0, 0.1]. The corresponding dissipative correction term α∆t(v n+1 -v n ) remains typically vanishingly small for low frequency modes with a small difference of successive velocity values within a typical time step; the latter establishes the consistency of the present modification in the spirit of the EC scheme, for the case where high frequency modes are not triggered. Meanwhile, the correction terms are less likely to vanish for high frequency modes, where the differences in successive values within a time step can be quite significant. Therefore, the energy decaying scheme is only introducing correction for higher order modes, which is not likely to reduce the scheme performance from the second-order accuracy. With such a choice of algorithmic update, the velocity at time step t n+1 takes a modified form

v n+1 = 1 (1/2 + α)∆t (d n+1 -d n ) - 1/2 -α 1/2 + α v n . (4.70)
One can follow the same procedure as presented in Section 4.2.4 to cast the work done by the inertia force as an increment in kinetic energy and dissipation terms

(d n+1 -d n ) T M a n+1/2 = 1 2 v T n+1 M v n+1 - 1 2 v T n M v n + α(v n+1 -v n ) T M (v n+1 -v n ) D K,n+1/2 = K n+1 -K n + D K,n+1/2 , ( 4
.71) where D K,n+1/2 is considered as the numerical dissipation of kinetic energy within an interval of time step due to the introduction of dissipative correction term α∆t(v n+1 -v n ) in Eq. (4.69). Note that dissipation is always a positive scalar since the mass is a positive-definite matrix. The work done by potential force is the same as derived in the energy conserving scheme. From the above derivation, the balance of the internal energy and external work can be formed in following equation

(K n+1 + P n+1 ) E n+1 -(K n + P n ) En +D K,n+1/2 = (d n+1 -d n ) T f ext n+1/2 ∆W n+1/2 ; E n+1 -E n + D K,n+1/2 =∆W n+1/2 , (4.72)
It is obvious that switching between energy conserving/decaying schemes can be executed flexibly without changing the main code of user-defined element since the update of displacement and velocity updates can be programmed in a separate subroutine for time-stepping integration algorithm. Additionally, the amount of numerical dissipation is fully controllable via an appropriate selection of dissipation coefficient α.

The algorithmic updates can be applied to moment and shear force projections in the same manner

m alg n+1/2 := 1 2 (m n+1 + m n ) + β(m n+1 -m n ); q alg n+1/2 := 1 2 (q n+1 + q n ) + β(q n+1 -q n ), (4.73)
where β is also the dissipation coefficient. The projections of moment and shear forces at time step t n+1 take a modified form

m n+1 = 1 1/2 + β m alg n+1/2 - 1/2 -β 1/2 + β m n ; q n+1 = 1 1/2 + β q alg n+1/2 - 1/2 -β 1/2 + β q n . (4.74)
The work done by potential force related to the moment field can be written as follows

(d n+1 -d n ) T F T m n+1/2 = 1 2 (d n+1 -d n ) T F T (m n+1 + m n ) + β(d n+1 -d n ) T F T (m n+1 -m n ) = 1 2 d T n+1 F T m n+1 - 1 2 d T n F T m n + β(m n+1 -m n ) T H(m n+1 -m n ) + 1 2 d T n+1 F T m n - 1 2 d T n F T m n+1 L 1 . (4.75)
The similar approach to the derivation of L 1 in Eq. ( 4.60) can be employed. Correspondingly, the work done by the potential force related to the moment field takes a new form as follows

(d n+1 -d n ) T F T m n+1/2 = (d T n+1 F T m n+1 - 1 2 m T n+1 Hm n+1 ) P n+1 -(d T n F T m n - 1 2 m T n Hm n ) Pn + β(m n+1 -m n ) T H(m n+1 -m n ) D P m,n+1/2 = P n+1 -P n + D P m,n+1/2 . (4.76)
where D P m,n+1/2 is the numerical dissipation of the potential energy related to the moment field within an interval of time step. In the same manner as Eq. (4.75), the work done by the potential force related to the shear forces field can be written as follows

(d n+1 -d n ) T G T q n+1/2 = 1 2 (d n+1 -d n ) T G T (q n+1 + q n ) + β(d n+1 -d n ) T G T (q n+1 -q n ) = 1 2 d T n+1 G T q n+1 - 1 2 d T n G T q n + β(q n+1 -q n ) T Z(q n+1 -q n ) + 1 2 d T n+1 G T q n - 1 2 d T n G T q n+1 L 2 . (4.77)
The similar approach to the derivation of L 2 in Eq. (4.60) can be employed. Corre-spondingly, the work done by the potential force takes a new form as follows

(d n+1 -d n ) T G T q n+1/2 = (d T n+1 G T q n+1 - 1 2 q T n+1 Zq n+1 ) P n+1 -(d T n G T q n - 1 2 q T n Zq n ) Pn + β(q n+1 -q n ) T Z(q n+1 -q n ) D P q,n+1/2
= P n+1 -P n + D P q,n+1/2 , (4.78) where D P q,n+1/2 is the numerical dissipation of the potential energy related to shear forces field within an interval of time step. From the above derivation, the balance of the internal energy and the external work can be formed in following equation

(K n+1 + P n+1 ) E n+1 -(K n + P n ) En + (D K,n+1/2 + D P m,n+1/2 + D P q,n+1/2 ) D n+1/2 = (d n+1 -d n ) T f ext n+1/2 ∆W n+1/2 ; E n+1 -E n + D n+1/2 = ∆W n+1/2 . (4.79)
where D n+1/2 is the numerical dissipation of the total energy due to the dissipative constitutive algorithm in Eq. (4.74). The correspondingly updated form of the balance of the internal energy and the external work, see Eq. (4.79), for the viscoplastic regime coupled with energy decaying scheme is written as follows

E n+1 -E n + D n+1/2 + D vp,n+1/2 = ∆W n+1/2 , (4.80)
in which the visco-plastic dissipation energy D vp,n+1/2 is the same as the one in Eq. (4.68).

Numerical simulations of plates under dynamic loads

The proposed energy conserving and decaying time stepping schemes are employed to compute the lateral deflection, rotations, moment and shear force responses in several illustrative examples of plate structures vibrations. The examples are simulated with a time step ∆t = 10 -2 , for the total number of time steps n = 800. Two triangle pulses T P 1 and T P 2 are introduced to the system to let each model vibrate freely, see Fig. 4.6(a). The first pulse T P 1 acts within a time interval 0 ≤ t ≤ t 1 , in which t 1 = 80 × ∆t, peaking in half of the interval. The second pulse T P 2 acts within a time interval t 2 ≤ t ≤ t 3 , in which t 2 = 250 × ∆t and t 3 = 260 × ∆t, peaking in the middle of this interval. We choose the second pulse T P 2 in order to introduce high-frequency modes contribution to the total energy of the system. Three different plate structures are selected to present the numerical performance including a cantilever, square and circular plate, which are denoted as model 1, 2 and 3, respectively. Regarding the elastic regime coupled with energy conserving scheme, each model will run with a higher value of yield stress σ y and zero value of visco-plastic parameter γ so that the total energy of the system is always conserved. Meanwhile, the visco-plasticity energy decaying scheme is conducted via five numerical simulations, with remaining parameters including α, β and σ y , γ = 1/η are given in Table 4.1. For simplicity, we select the same values for both α, β to show a smoothening effect of higher value of dissipation parameters. All cases from 1 to 5 indicate the using of current simple algorithms. Meanwhile case 6, which is used for comparison in Fig. 4.12(a), indicates the elastic plate element proposed by [START_REF] Auricchio | A shear deformable plate element with an exact thin limit[END_REF] along with the Newmark time integration scheme.

Cantilever rectangular plate

In model 1, a free vibration of a 10 × 1 rectangular cantilever plate with thickness t = 0.1 is selected and investigated via the time history of degrees of freedom along with stress resultants at the free end, the total energy and relevant dissipations. The vibration of cantilever plate is illustrated via the corresponding plate response, time history of degrees of freedom along with stress resultants at the free tip, energy and dissipation over the entire domain. The plate is fixed at the left end. The regular mesh is constructed by 512 triangular plate elements, by splitting each of 32 × 8 square sub-domains into 2 triangle elements. The material properties are selected as follows: elastic modulus E = 30 × 10 6 , Poisson's ratio ν = 0.25 and mass density ρ = 0.785. The first pulse T P 1 has the peak magnitude p 1 = p, while the second pulse T P 2 has the peak magnitude 10 times bigger p 2 = 10 × p.

The first pulse T P 1 is created by a vertical distributed force p = 0.5 applied on the right end. The second pulse T P 2 is applied by using distributed forces at six different locations along the plate in order to incite high-frequency modes. Those locations along the x axis are 27/32 28/32 29/32 30/32 31/32 1 with the corresponding direction as +1 -1 +1 -1 +1 -1 . The couple (+1, -1) represents a positive and negative direction, respectively. Firstly, the accuracy of the proposed plate element is compared with respect to robust non-locking Reissner-Mindlin plate elements (e.g. see [START_REF] Auricchio | A shear deformable plate element with an exact thin limit[END_REF]), as shown in Fig. 4.4. Only the first pulse T P 1 is introduced to the system, which leads to the free vibration after releasing the pulse. With respect to the time integration scheme, the proposed plate element is accompanied by the energy conserving scheme, meanwhile the plate element of [START_REF] Auricchio | A shear deformable plate element with an exact thin limit[END_REF] is accompanied by the Newmark scheme with standard parameters. The different regular mesh settings using for the comparison in Fig. 4.4 include mesh 1 (128 elements), mesh 2 (512 elements) and mesh 3 (2048 elements). It is recognized that the proposed element delivers a high accuracy (for both lateral deflection w at the free end and total reaction moment Myy at the boundary) even with the coarse mesh, whereas the plate element of [START_REF] Auricchio | A shear deformable plate element with an exact thin limit[END_REF] does not fully converge even for the finest mesh with vibration magnitude increasing in a number of elements. As shown in Fig. 4.5, the moment M yy along the plate is presented at time step n = 150 and n = 250. At each time step, the moment M yy is highly agreed from the two different meshes. To test the energy conserving scheme, we select a high value of yield stress σ y and a zero value of visco-plastic parameter γ so that the cantilever plate always remains in the elastic regime. The total energy of the plate is conserved during the free vibration phase by using the energy conserving scheme, see Fig. 4.6(b). We note that the total energy gets high value at the beginning due to the external load, and then it stays constant after the end of the first pulse T P 1 until the introduction of T P 2 . The total energy still conserves after the second pulse T P 2 . Correspondingly, the cantilever plate response at time step n = 250 is presented in Fig. 4.6(c,d). The lateral deflection w and rotation θ y are dominated by plate bending along the y-axis. The maximum absolute values of the lateral displacement w and the rotation θ y are recorded at 0.14 and 0.02, respectively. Next, we introduce the visco-plasticity for testing energy decaying scheme via five numerical simulations. The remaining parameters including α, β, σ y and γ = 1/η are given in Table 4.1. The results are extracted at point A(10, 0.5) for displacement degrees of freedom and B(10, 0.4375) moment, shear degrees of freedom. In the first several steps, the amplitude response from case 1 is comparable with that from case 6 (see [START_REF] Auricchio | A shear deformable plate element with an exact thin limit[END_REF]) under the Newmark scheme. Soon after the introducing of higher frequency modes via the second pulse T P 2 , we can see the insignificant differences between cases as shown in Fig. 4.7(a). The higher value of dissipation parameters (α, β) get, the corresponding response get smoother. The ED acts as a high frequency filter which is capable of smoothening the response curve. The clear illustration of smoothening effect can be seen in in Fig. 4.7(b) by comparing case 1 with (α, β) = (0, 0) and case 5 with a much higher value of (α, β) = (0.08, 0.08). It is expected that the computed response gets smoother with the higher value of dissipation parameters. This choice can be adjusted by the users for particular problems or even in each particular simulation. Within time period t ≤ t 2 , all cases yield the similar stress resultant responses from an element located at the middle right edge, see Fig. 4.8. The second pulse T P 2 generates more noises to the stress resultant response, triggered by participation of higher frequency modes. It is clear that the smoothening effect of the ED scheme can be observed in these cases. The bending moments M xx and M yy get smooth evolution rather fast with a higher value of dissipation parameters (α, β) = (0.08, 0.08), see Fig. 4.8(a,b). It is noted from Fig. 4.8(c,d) that the bending moments M xy and M yx are comparable although the approximation by RT 0 can not guarantee the symmetric of stress tensor. The shear force evolution also gets smooth rather faster for higher values of dissipation parameters, see Fig. 4.8(e,f). With the time step n ≤ 250, the total energy, algorithmic and visco-plastic dissipations due to the ED scheme and visco-plasticity are presented in Fig. 4.9. After introducing the first load pulse T P 1 , the system oscillates freely, the total energy reduces negligibly due to the existence of both numerical and visco-plastic dissipation, see Fig. 4.9(a). With the zero dissipation parameters (α, β), there is no algorithmic dissipation in case 1, which is equivalent to the EC. Meanwhile, it is clear that the higher value of dissipation parameters (α, β) are chosen, the higher value of algorithmic dissipation gain, see Fig. 4.9(b). The algorithmic dissipation due to the ED scheme rises linearly in those five simulations. Regarding the visco-plastic dissipation D vp , the dissipation develops moderately after t 1 /2, see Fig. 4.9(c). It develops slowly at some periods, which can be due to the state of stress resultants yielding a low value of φ(s). Within time step n ≤ 250, the dissipation is dominated by algorithmic dissipation. The total dissipation including both algorithmic and visco-plasticity is shown in Fig. 4.9(d). At n = 250, the maximum dissipation approaches a value of 8 × 10 -3 in case 5 and a value of 1 × 10 -3 in case 1. To illustrate the whole computational period n ≤ 800, full history of both energy and dissipation are presented in Fig. 4.10. The second pulse T P 2 injects more higher frequencies energy to the system and then dissipates gradually along with time, see Fig. 4.10(a). Again, with the zero dissipation parameters (α, β), there is no algorithmic dissipation in case 1. It is clear depicted in Fig. 4.10(b) that much more higher frequencies energy are dissipated during t ≥ t 3 compared to the first period t ≤ t 2 in Fig. 4.9(b). The same correlation between the dissipation parameters (α, β) and algorithmic dissipation can be observed as in Fig. 4.9. The visco-plastic dissipation still increases after t 3 , see Fig. 4.10(c). The total dissipation is given in Fig. 4.10(d). At the end of computational period n = 800, the maximum dissipation approaches 0.085 in case 5 and the one of case 1 gains a value of 0.04.

Square plate

In model 2, the free vibration of a 1 × 1 square plate with thickness t = 0.1 is chosen and illustrated via the time history of degrees of freedom along with stress resultants at the center, the total energy and relevant dissipation. The plate is clamped along the outer boundary. The regular mesh is constructed by 512 triangular plate elements, by splitting each of 16 × 16 square sub-domains into 2 triangle elements. The material properties are selected with elastic modulus E = 10.92 × 10 6 , Poisson's ratio ν = 0.3 and mass density ρ = 0.785. The first triangular pulse T P 1 is assigned by only a point load P = 100, which is applied at the center. Meanwhile, the second pulse T P 2 is applied in terms of point loads at five locations in order to trigger high-frequency modes. These locations are (1/4, 1/4) (3/4, 1/4) (3/4, 1/4) (1/4, 3/4) (1/2, 1/2) with the corresponding direction as -1 +1 -1 +1 +1 . The couple (+1, -1) represents a positive and negative direction, respectively. In both pulses, the peak magnitude is the same

p 1 = p 2 .
To test the energy conserving scheme, we select a high value of yield stress σ y and zero value of visco-plastic parameter γ so that the square plate always remains in the elastic regime. The total energy of the plate is conserved under free vibration by using energy conserving scheme, see Fig. 4.11(a). We note that the higher frequency modes are introduced into the system after the second pulse T P 2 . Both kinetic and potential energies lost the smoothness after time step t 2 due to the second pulse T P 2 . Correspondingly, the square plate response at time step n = 250 is presented. The lateral deflection w, rotation θ x and rotation θ y are shown in Figs. 4.11(b-d), respectively. The maximum absolute value of lateral deflection is 5 × 10 -4 located at the center. Subsequently, we introduce the visco-plasticity and energy decaying scheme via five numerical simulations. The remaining parameters including α, β and σ y , γ = 1/η are given in Table 4.1. After the introduction of higher frequency modes via the second pulse T P 2 , we can see the differences between case 6 (see [START_REF] Auricchio | A shear deformable plate element with an exact thin limit[END_REF]) by using the Newmark scheme and others as shown in Fig. 4.12(a) though all magnitudes remain nearly the same. Meanwhile the response by case 6 seems to be lagging compared to the current algorithm. The response is recorded at the plate center. The higher value of dissipation parameters (α, β) get, the corresponding response get smoother. The ED takes a role as a high frequency filter which is capable of smoothening the response curve. The clear illustration of smoothening effect can be seen at moment degree of freedom m x , shown in Fig. 4.12(b), by comparing case 1 and case 5. It is expected that the curve gets smoother with the higher value of dissipation parameters. The remaining degrees of freedom are shown at Figs. 4. 12(c,d). We note that the magnitude of m y and q d.o.f get higher by t ≥ t 3 while it is not the case for m x d.o.f. Within time period t ≤ t 2 , all cases yield the similar stress resultant responses at the plate center, see Fig. 4.13. The second pulse T P 2 generates more noises to the stress resultant responses compared to the first pulse T P 1 . The noises are corresponding to higher frequency modes. The bending moments M xx and M yy get smoothening fast with a high value of dissipation parameters (α, β) = (0.08, 0.08), see Fig. 4.13(a,b). It is noted from Fig. 4. 13(c,d) that the bending moments M xy and M yx are not exactly the same due to the approximation by RT 0 . The shear forces also decay fast with higher dissipation parameters, see Fig. 4.13(e,f). All data is recorded at the plate center. With the time step n ≤ 250, the total energy, algorithmic and visco-plastic dissipations due to the ED scheme and visco-plasticity are presented in Fig. 4.14. After introducing the first pulse T P 1 , the system oscillates freely, the total energy reduces negligibly due to numerical dissipation, see Fig. 4.14(a). With the zero dissipation parameters (α, β), there is no numerical dissipation in case 1, which is equivalent to the EC. Meanwhile, it is clear that the higher value of dissipation parameters (α, β) are chosen, the higher value of numerical dissipation gain, see Fig. 4.14(b). The algorithmic dissipation due to the ED scheme rises linearly in those five simulations. Regarding visco-plastic dissipation D vp , the dissipation develops substantially after t 1 /2, see Fig. 4.14(c). It develops slowly at some periods, which can be due to less number of yielding elements. Within time step n ≤ 250, the total dissipation is dominated by algorithmic dissipation as shown in Fig. 4.14(d). At n = 250, the maximum dissipation approaches a value of 1.7 × 10 -3 in case 5 and the one of case 1 gains a value of 0.25 × 10 -3 . To illustrate the whole computational period n ≤ 800, all energy and dissipations are presented in Fig. 4.15. As expected, the second pulse T P 2 generates more higher frequencies energy to the system, see Fig. 4.15(a). Again, with the zero dissipation parameters (α, β), there is no numerical dissipation in case 1, which is equivalent to the EC. It is clearly depicted in Fig. 4.15(b) that much more higher frequencies energy are dissipated with time t ≥ t 3 compared to the first period t ≤ t 2 in Fig. 4.14(b). The same correlation between the dissipation parameters (α, β) and algorithmic dissipation during n ≥ 250 can be observed similarly to that at the beginning period with n ≤ 250. The visco-plastic dissipation also increases after t 3 due to the higher stress values in the plate, see Fig. 4.15(c). During the computational period, the total dissipation is dominated by algorithmic dissipation as shown in Fig. 4.14(d). At the end of computational period n = 800, the maximum dissipation approaches 0.1 in case 5 and the one of case 1 gains a value of 0.046.

Circular plate

In model 3, the free vibration of a circular plate with radius r = 5 and thickness t = 1 is selected to compute the time history of the total energy and relevant dissipations, along with stress resultants computations. The plate is simply supported around the circumferential boundary, and loaded in the center. The symmetry condition is applied to analyze only a quarter of the plate by imposing constraints upon the tangential rotations along radial boundaries, see [START_REF] Ibrahimbegovic | Quadrilateral finite elements for analysis of thick and thin plates[END_REF]. The mesh is divided into 3 sub-regions, with each one counting 128 triangular plate elements. The material properties are selected with elastic modulus E = 10.92 × 10 6 , Poisson's ratio ν = 0.3 and mass density ρ = 0.785. The first pulse T P 1 is assigned by only a point load P = 200 is applied at the center. The second pulse T P 2 is assigned at mesh intersection point I(2r/5, 2r/5) with the same amplitude but in the reversed direction. In both pulses, the peak magnitude is the same p 1 = p 2 .

To test the energy conserving scheme, we select a high value of yield stress σ y and zero value of visco-plastic parameter γ so that the circular plate is always under the elastic regime. The total energy of the plate is conserved under free vibration by using energy conserving scheme, see Fig. 4.16(a). We note that the total energy get high value at the beginning of each pulse and then it stays constant after the end of each pulse. It is clear that the higher frequency modes are introduced into the system after the second pulse T P 2 . Both kinetic and potential energies time histories get more noisy after time step t 2 corresponding to application of the second pulse T P 2 . The cantilever plate responses at time step n = 250 is presented. The lateral deflection w, rotation θ x and rotation θ y are shown in Figs. 4.16(b-d), respectively. The maximum absolute value of lateral deflection is 7 × 10 -4 located at the center. We now introduce the visco-plasticity and energy decaying scheme via five numerical simulations. The remaining parameters including α, β and σ y , γ = 1/η are given in Table 4.1. In the first several steps, the deflection response from case 1 highly agrees with that from case 6 (see [START_REF] Auricchio | A shear deformable plate element with an exact thin limit[END_REF]) by using the Newmark scheme. Soon after the introducing of higher frequency modes via the second pulse T P 2 , all deflection cases are still comparable as shown in Fig. 4.17(a). The response is extracted from the point at the plate center. The higher are the values of dissipation parameters (α, β), the smoother is the corresponding computed response. The ED plays a role of a high frequency filter which is capable of smoothening the response curve. The clear illustration of smoothening effect can be seen in Fig. 4.17(b) by comparing case 1 and case 8 with a much higher value of (α, β). It can be seen as expected that the curve gets smoother with the higher value of dissipation parameters. Within t ≤ t 2 , all cases yield the similar stress resultant responses at the plate center, see Fig. 4. 18(a-d). The second pulse T P 2 again creates a lot of high frequency noises to the stress resultant responses compared to those in the plate deflection. The bending moments M xx and M yy get smoothened fast with dissipation parameters (α, β) = (0.08, 0.08), see Fig. 4. 18(a,b). The bending moments M xy and M yx have the same decaying trend but large differences in magnitudes, see Fig. 4. 18(c,d). The shear forces are plotted in Fig. 4.18(e,f). They decay fast within the computational period. All data is recorded at the point at the plate center. With the time step n ≤ 250, the total energy, algorithmic and visco-plastic dissipations due to the ED scheme and visco-plasticity are presented in Fig. 4.19. After introducing the first load pulse, the system vibrates freely, the total energy nearly remains the same due to small values of both numerical and visco-plastic dissipation, see Fig. 4. 19(a). With the zero dissipation parameters (α, β), there is no numerical dissipation in case 1, which is equivalent to the EC. Meanwhile, it is clear that the higher value of dissipation parameters (α, β) are chosen, the higher value of numerical dissipation gain, see Fig. 4. 19(b). Regarding the visco-plastic dissipation, the dissipation starts during the pulse introduction and then later get stable around a constant, see Fig. 4.19(c), which reflects the fact that all elements in the system are under the elastic regime. During this period of time, the total dissipation is dominated by algorithmic dissipation as shown in Fig. 4. 19(d). At n = 250, the maximum dissipation approaches a value of 1.2 × 10 -3 in case 5 and the one of case 1 gains a value of 0.2 × 10 -3 . To illustrate the whole computational period n ≤ 800, all energy and dissipations are presented in Fig. 4.20. The second pulse T P 2 injects more higher frequencies energy to the system and then dissipates gradually along with time, see Fig. 4.20(a). Again, with the zero dissipation parameters (α, β), there is no numerical dissipation in case 1. It is clearly depicted in Fig. 4.20(b) that much higher frequencies energy are dissipated during t ≥ t 3 compared to the first period t ≤ t 2 in Fig. 4. 19(b). The same correlation between the dissipation parameters (α, β) and algorithmic dissipation can be observed as in the previous computational period t ≤ t 2 . The visco-plastic dissipation also increases greatly after remaining constant with time t ≥ t 3 due to the higher stress values in the plate, see Fig. 4.20(c). The total dissipation is dominated by algorithmic dissipation as shown in Fig. 4. 20(d). At the end of computation period n = 800, the maximum dissipation approaches 0.035 in case 5 and the one of case 1 gains a value of 0.02.

We also notice that there is a correlation between dissipation parameters and the evolution of visco-plastic dissipation. In contrast to algorithmic dissipation, the visco-plastic dissipation gets lower with higher value of (α, β). The dissipation coefficients (α, β) have independent smoothening effects on computed response for both displacements and stress-resultants. The higher value of this coefficient, the more significant is the effect of smoothening in computed structural response. The selection of different values between (α, β) are still left open. In these simulations, the total dissipation is dominated by algorithmic dissipation, see Fig. 4.10(d), Fig. 4.15(d) and Fig. 4.20(d). This phenomenon can be explained due to the fact that the higher frequency modes are dominant after introducing the second pulse T P 2 , hence the numerical dissipation is dominant part of the total dissipation. The selected value of visco-plastic parameter γ is typically small (see [START_REF] Hughes | Unconditionally stable algorithms for quasistatic elasto/visco-plastic finite element analysis[END_REF]) hence the visco-plastic dissipation is comparable low to its counterpart. Here, we also examine the same simulation with higher value of visco-plastic parameter γ at 1×10 -4 and 1 as shown in Fig. 4.21(a) and Fig. 4.21(b), respectively. The visco-plastic dissipation gains much higher value, but still remains low compared to algorithmic dissipation. However, we can increase the part of visco-plastic dissipation by choosing lower value of yield stress σ y (which means plate crosses into visco-plastic regime early) and fixing visco-plastic parameter γ = 1 × 10 -8 as shown in Fig. 4. 22(a,b). For such lower value of yield stress σ y , there are more plate elements satisfying yield condition and thus the visco-plastic dissipation increases; see Fig. 4. 22(a). The increase of visco-plastic dissipation is not significant compared to those in Fig. 4.21 due to the fact that the visco-plastic dissipation depends much on the rate of viscoplasticity β(s) as shown in Eq. (4.47), and this rate of visco-plasticity β(s) increases much with higher value of visco-plastic parameter γ via relation in Eq. (4.27). It is also noted that the algorithmic dissipation remains practically the same with the change of yield stress σ y . We further examine the simulation with higher value of visco-plastic parameter γ and lower value of dissipation coefficients (α, β) while keeping other parameters the same as the beginning of this section. As shown in Fig. 4.23(a), when the visco-plastic parameter γ takes a value of 1 × 10 -4 , the algorithmic dissipation with the input α = β = 1 × 10 -8 is lower than the visco-plastic counterpart. While in Fig. 4.23(b), when the visco-plastic parameter γ takes a value of 1, the algorithmic dissipation with the input α = β = 1 × 10 -4 is lower than visco-plastic counterpart. With the lower values of dissipation coefficients (α, β), the time history response gets smoothening slower as observed in the previous Figs. (4.17,4.18). From those observations, the visco-plastic dissipation can dominate with a proper selection of visco-plastic parameters γ and dissipation parameters (α, β). Final results are seeking to illustrate the smoothness of computed moment field. With the same simulation for cantilever and circular plates, we present the plots of degree of freedom of moment over entire plates as shown in Fig. 4.24(a) and Fig. 4.24(b) under the EC scheme and the elastic regime at time step n = 250, respectively. The max magnitude of vector m(m x , m y ) is recorded at 13.6 and 128.1 in case of cantilever and circular plates, respectively. For cantilever plate simulation results in Fig. 4.24(a), it is observed that there is a main vectorial stream along the horizontal direction which dominates the rest. In circular plate simulations in Fig. 4.24(b), it is observed that there are two main vectorial streams including centripetal and tangential directions. The current approximation of stress-resultant fields by the lowest Raviart-Thomas vector space allows the continuity of those degree of freedom over each element's edges. Thus, we consider that the stress-resultant waves can be smoothly propagated throughout the entire domain of interest, which is in contrast with the conventional finite element method where the stress-resultant fields are discontinuous from element to element.

5

Parameter identification for inelastic and multi-scale problems

In this chapter, the parameter identification for inelastic and multiscale problems are presented. Firstly, the theoretical background of uncertainty quantification is reviewed. In details, several key definitions including random variables, random fields, Karhunen-Loève expansion (KLE), Bayesian theorem, and Gauss-Markov-Kalman filter are briefly summarized. An example to show a capability of the parameter identification using the Bayesian updating method is conducted in a one dimensional bar problem to assimilate fracture energy G f . Finally, the parameter identification using the Gauss-Markov-Kalman filter is employed for a multi-scale problem to identify bulk and shear moduli and other material properties in a macro-scale with the data from a microscale as quantities of interest (QoI). Equivalently, the problem can be viewed as upscaling homogenization.

Theoretical background

The uncertainty quantification has been widely utilized in the field of computational mechanics. Many research works use this method to quantify the source of uncertainty in the numerical models. There are several sources of uncertainties in a numerical simulation, such as material and geometry uncertainties. In general, there are two common definitions, namely aleatory and epistemic uncertainties. The aleatory uncertainty is the inherent variation in a quantity while the epistemic uncertainty is due to the lack of knowledge of an analyst. There are two popular problems in uncertainty quantification. The first one is the so-called forward problem, which is used to propagate the uncertainty through the model of interest and investigate the outputs, so that some of key features or critical effects of inputs onto the model can be discovered. Meanwhile, the backward or inverse problem is employed to approximate the parameter of interest via the observed outputs of a given numerical model. In the following multi-scale problem, the data from the micro-scale can be used to assimilate data in the macro-scale, such process is named as an upscaling process. In contrast, the inversed one is named as a downscaling process.

Principle concepts and methods

The sample space Θ comprises all possible elementary outcomes θ of a random experiment. By tossing a fair coin two times, we can have the sample space for this experiment is Θ = [HH, HT, T H, T T ], in which the notation H stands for the coin head and the notation T stands for the coin tail.

A random variable (RV) is a variable taking different value whenever it is observed, e.g. a length of a ruler. In the rest of this chapter, we deal with a real-valued random variable which is a mapping X : Θ → R, e.g. the material elastic modulus, hardening and softening moduli. An observed value of a random variable is called realization x = X(θ). Some of random variables can be described by a probability density function. For a continuous random variable, the probability density function f x is a derivative of the cumulative distribution function

F X F X (b) -F X (a) = b a f X (x) dx.
(5.1)

A mean or expected value of a discrete random variable can be computed by using the sum of all possible values divided by the number of values, while the mean value of a continuous random variable is defined as

µ x = E[x] = ∞ -∞ xf X (x) dx. (5.2)
The variance of a random variable is defined as

V[X] := E (X -E[X]) 2 = E[X 2 ] -E[X] 2 . (5.3)
Correspondingly, the standard deviation is defined as a root square of the variance, e.g σ = V[X].

There are some common distributions, e.g Uniform, Normal or Gaussian, Log-Normal and Beta distributions. One of the most widely used distributions is the Normal or Gaussian distribution which is denoted as X ∼ N (µ, σ). This distribution is fully defined by the expected value µ and the standard deviation σ. The probability density function is given as

f X (x) = 1 σ √ 2π e -(x-µ) 2 2σ 2 .
(5.4)

A Log-Normal random variable X ∼ LN (µ, σ) is defined as X = e Z with Z ∼ N (µ, σ). The corresponding probability density function is given as A random field (RF) is a collection of random variables and it explicitly depends on space and takes different value whenever it is observed. For example, twodimensional sources of water or oil underground can be viewed as random fields. It can be stated that a random field is a mapping α : (R n , Θ) → R with n ≥ 2. Similar as the random variable, the mean value function is the most basic characteristic of a random field

f X (x) = 1 xσ √ 2π e -(ln(x)-µ) 2 /(2σ 2 ) . ( 5
µ α (x) = E[α(x, ω)] = Ω α(x, ω)P(dω). (5.6)
In a practical approach, the corresponding oscillating part α(x, ω) -µ α (x) is commonly represented via a linear combination of deterministic functions and coefficients, in which the latter can be random variables. There are several methods to separate the representation as above, such as the Karhunen-Loève expansion (KLE) and the Proper Orthogonal Decomposition (POD). In the following, only the Karhunen-Loève expansion (KLE) is discussed since it is one of the most common methods in the field of uncertainty quantification. Let us consider a random field in a two-dimensional domain Ω. It is clear that there is an infinite number of random variables at each point x ∈ Ω. The target is to reduce the number of these random variables by smoothening a given random field. In other words, it can be represented as a finite series of products of spatial functions multiplied by scalar random variables, in which the error with respect to the original random field can be neglected. Given a Gaussian stochastic field, the Karhunen-Loève Expansion, see [START_REF] Smith | Uncertainty quantification: theory, implementation, and applications[END_REF] and [START_REF] Chiachio-Ruano | Bayesian Inverse Problems: Fundamentals and Engineering Applications[END_REF], is defined as

Θ(x, ω) = µ Θ (x) + ∞ i=1 σ i ψ i (x)θ i (ω).
(5.7)

It can be seen clearly that Eq. (5.7) is composed of the first term as a mean value function µ Θ (x), which depends only on space, and the second term as an oscillating part. The spatial functions ψ i (x) are square integrable and orthogonal. θ i (ω) is the independent standard Gaussian N (0, 1) random variable, and the σ i are multipliers.

The spatial functions ψ i (x) can be discretized in the spatial domain in a form as

ψ i (x) ≈ ψ h,i (x) = n j=1 Ψ j (x)v ji = Ψ V . (5.8)
in which Ψ is a matrix of spatial basis functions as shape functions in the finite element method. The optimal spatial functions are the solution of the generalized eigenvalue problem

GCGv i = λ i Gv i = σ 2 i Gv i , ( 5.9)
in which λ i and v are the generalized eigenvalues and eigen functions. The G(Ψ i , Ψ j ) is the Gramian matrix of the basis functions, which can be assembled in the same manner as the mass matrix ∀x ∈ Ω as

G(Ψ ) = Ω Ψ i (x)Ψ j (x) dΩ. (5.10)
Meanwhile, the covariance matrix C is computed via Matérn function between each couple of points, e.g. C i,j = C νc (d ij /l c ). The Matérn function is defined as

C νc (d ij /l c ) = σ 2 2 1-νc Γ (ν c ) √ 2ν c d c νc K νc √ 2ν c d c , ( 5.11)
where ν c is a non-negative parameter, σ is also an input parameter, d ij is the distance between the two input points (i, j), and l c is the correlation length. The Gaussian Θ(x) field can be approximated (5.12) where L ≤ n is the truncated number of the eigen functions. The condition to get L is defined as

Θ(x, ω) ≈ µ Θ + L i=1 σ i ψ i (x)θ i (ω) = µ Θ + Ψ (x)V Sθ(ω),
ρ L = L i=1 λ i n j=1 λ j , ( 5.13)
in which n is the total possible eigenvalues in a given system. The relative cumulative factor ρ L can be predefined, normally ρ L ≥ 0.9.

Algorithm 1 Template to generate realization of Gaussian distribution variables To generate the realization of a random variable, one can use several available methods, such as the Monte-Carlo method, the Metropolis algorithm, the Metropolis-Hastings algorithm, and the Gibbs sampling algorithm. Some of those methods are available in the open-source SGLIB library coded in MATLAB, see [START_REF] Zander | Sglib[END_REF]. An example to generate realization using the Monte-Carlo method for two random variables, namely var1 and var2, is given in Alg. 1 by using the SGLIB library. The SimP aramSet is a general-purpose object which is capable of executing several tasks by its programmed methods. After initialization of the object SimP aramSet in Step 1, one can define several prior distribution by a method socalled add. Next, the type of distribution with corresponding default parameters are defined as the method inputs. The realization then can be generated in Step 3 via a method so-called sample. The realizations of random variable x ∼ N (10 3 , 10) using SimP aramSet is shown in Fig. 5.2 with different number of samples. The larger the number of samples, the closer the histogram gets to the Gaussian distribution. In the final step, it depends on the working procedure of an analyst to set up all the necessary files in which the new values of realization are written.

Bayesian updating procedure

The non-intrusive and intrusive methods are common approaches for the uncertainty quantification analysis. The former is proceeded by treating the available finite element model, such as in-house codes or commercial computational programs, as black boxes. The main work is to generate the realization of input parameters and then apply these parameters into the legacy code to run and extract the interest output data. The output data is then jointed with the experiment data to access the property of QoI, e.g. input parameters. Thanks to the feasibility and simplicity, this method is widely accepted in many research works. By contrast, the latter method heavily reformulates the weak form or variational formulation of the problem and introduce the uncertainty into the system explicitly. Hence, this approach is more complicated and sometimes it is impossible with a number of uncertainty parameters due to the high complexity.

Only non-intrusive methods are investigated to identify material properties in this chapter. Specifically, the Bayesian updating method is adopted for the parameter identification. The inverse problem of identifying or calibrating the material properties in a given numerical model is addressed in the framework of the Bayesian estimation, which lead to a computation of the conditional expectation. The Bayesian theorem is considered as the consistent way to update a probabilistic description when new data in the form of observation y is available. In such case, it is possible to state the conditional density π (Q|Y ) (θ | y) of θ given y, see [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Estimation[END_REF] and [START_REF] Matthies | Parameter estimation via conditional expectation: a bayesian inversion[END_REF]

π (Q|Y ) (θ | y) = π (Q,Y ) (θ, y) π Y (y) = π (Y |Q) (y | θ) π Y (y) π Q (θ), (5.14)
in which The flowchart of the Bayesian updating method is shown in Fig. 5.3. Firstly, the realization of each interest input variable is generated via its corresponding prior distribution. The prior probability density function can be selected from several available distributions so that it is mathematically feasible to generate the realization of a given random variable. Secondly, the target physic or mechanic problem can be cast in a general form of mathematical representation M(x), which can be built via commercial computational programs or in-house codes. It can be also truncated and approximately represented via surrogate model, such as general Polynomial chaos (gPC) with a fast computational capability. In the following works, the model M(x) is solved by finite element method, which is programmed as a FORTRAN user-defined element in the FEAP computational program. Then the quantity of interest (QoI) depends on the nature of the problem and the selection of an analyst. However, it is normally selected based on the criterion that it can be measured explicitly with the aid of experimental tools. The common examples of QoI are displacement, velocity, acceleration, strain, reaction force, stored energy, and even dissipation. Finally, with the data from the experiment, all essential ingredients are completed in the preparation for the Bayesian updating. In the following, the Bayesian updating method is conducted to update the fracture energy G f under the softening regime in a one-dimensional coupled elastoplasticity. A cantilever bar is fixed on the left and an imposed displacement ū = 0.01 is applied within 500 time steps on the right. The reaction force and elastic energy at the last time step are selected as the quantities of interest. The assumed true values of the material properties are given in Table 5.1, e.g. cross-section area A, Young elastic modulus E, yield stress σ y , linear hardening modulus K h , ultimate stress σ u , linear softening modulus K s and fracture energy G f in the case of exponential softening regime. Table 5.1 -One-dimensional inelastic bar material properties It is assumed that the prior distribution given to the fracture energy G f is following the Normal distribution N (µ Gf , σ Gf ), in which the standard deviation σ Gf = 52.5 is set to be 15% of a mean value µ Gf = 3.5 × 10 2 . The realization of G f is generated by Monte-Carlo method with n s = 10 3 samples. With each value of the fracture energy G f , the problem is executed once. The responses from reaction and elastic energy are given in Fig. 5.4(a) and Fig. 5.5(a). As shown in Fig. 5.4(b) and Fig. 5.5(b), it can be seen that those updated standard deviations are much smaller than the prior counterpart. This means that after using the Bayesian updating method the updated value of the fracture energy is getting more accurate and reliable than the prior mean value of µ Gf . The collection of prior and posterior values of the fracture energy G f is summarized in Table 5.2.

π Y (y) = Q π (Q,Y ) (θ, y) dθ = Q π (Y |Q) (y | θ)π Q (θ) dθ, ( 5
Property A E σ y K h σ u K s G f

Gauss-Markov-Kalman filter

A mechanics problem itself can be viewed as a mapping M, which can be a linear or nonlinear mapping depending on the nature of the given problem, from the input data x to the output data y via relation y = M(x). To obtain more insights into the given problem, the input data x can be formulated as random variables or also as random fields. After going through the mapping M, it is clear that the output data y is no more deterministic but rather a stochastic or random one. The uncertainty in the input x can propagate thorough out the mapping M. More precisely the uncertainty propagation refers to mean value, variance, probability distribution and other quantities. Not only the forward propagation problem but also backward or parameter identification are of research interest. In practices, there are several methods to formulate the surrogate model of the mapping M. One of the most popular is so-called Polynomial Chaos Expansion (PCE), which is discussed briefly and employed for the process of updating in this chapter. More detailed discussion on PCE and generalized Polynomial Chaos can be found in [START_REF] Xiu | Numerical methods for stochastic computations[END_REF], [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF] and [START_REF] Römer | Methods of Uncertainty Analysis and Quantification[END_REF]. The Polynomial Chaos Expansion (PCE) for a random variable y, which has finite mean value and variance (E[y], V[y]) < ∞, and corresponding probability density function f y is defined as follows (5.16) in which the involving polynomial must be orthogonal with respect to the probability density function If the random variable y is represented via a Gaussian distribution basis function or germ ξ(ω) ∼ N (0, 1) then the corresponding polynomials are the Hermite functions since it is orthogonal with respect to the Gaussian distribution. Several computations can be executed via numerical integration instead of direct sampling. The coefficients of the PCE of the model can be obtained using output from the numerical model y (ξ j ) = M (ξ j ) as follows

y = ∞ i=0 q i Φ i (ξ),
f x E ξ [Φ i Φ j ] = R Φ i (α)Φ j (α)f ξ (α) dα = δ ij , (5.17
q i ≈ n int j=1 y (ξ j ) Φ i (ξ j ) w j E [Φ i (•) 2 ] , ( 5.19)
in which (ξ j , w j ) are numerical integration points j and corresponding weights. The extension to the PCE of multi-variable y = M(x 1 , x 2 , ..., x n ) is in the same manner with more germs and integration points. The essential mapping from germs ξ to input parameter is omitted sometimes if it is an identity map, nevertheless, the map is an injective function in general.

The second ingredient for the Gauss-Markov-Kalman filter is Kalman filter, which is presented in the following. Let us consider again the mathematical or numerical model, such as finite element model, with y = M(x) with x ∈ X . The task is now to calibrate or identify input parameter x a from a forecast x f , representing initial knowledge, and an observation ŷ for the random variable z = y + .

The filtered or assimilated random variable x a after the observation ŷ is given as (5.20) in which the term so-called innovation x i is actually an orthogonal error. It is observed that x a is unbiased and linear in x f , while E X (x f | ŷ) is linear and optimal.

x a = x f + (E X (x f | ŷ) -E X (x f | z)) = x f + x i ,
The problem turns into finding the solution of the minimization equation, see [START_REF] Matthies | Bayesian parameter estimation via filtering and functional approximations[END_REF], as

x f -(K(z) + a) 2 X = min L,b x f -(L(z) + b) 2 X , (5.21)
in which a = x f -K(z). The Kalman gain K satisfies the above minimization problem, hence the final form of the Kalman filter is given as

x a = x f + K(ŷ -z(ω)) = x f + K (ŷ -(M (x f (ω)) + (ω))) , ( 5.22)
in which the Kalman gain is computed as K = C x f z C -1 z and C x f z is the covariance of x f and z, while C z is the auto-covariance of z.

With the brief discussion on the Polynomial Chaos Expansion (PCE) and the Kalman filter, the Gauss-Markov-Kalman filter is presented briefly, see [START_REF] Chiachio-Ruano | Bayesian Inverse Problems: Fundamentals and Engineering Applications[END_REF][START_REF] Rosić | Sampling-free linear bayesian update of polynomial chaos representations[END_REF]169] for more detailed discussion and [START_REF] Marsili | Parameter identification via gpce-based stochastic inverse methods for reliability assessment of existing structures[END_REF][START_REF] Kumar | Bayesian calibration of volume averaged rans model parameters for turbulent flow simulations over porous materials[END_REF] for relevant approaches. The key idea of this filter is using the PCE method to approximate all relevant random variables and observation, and then the Kalman filter can be applied with these PCE coefficients in order to update the PCE coefficients of the targeted random variables. Hence, this method is proceeded without any sampling work. The input parameter, the error and the measurement are written in the PCE format

q = l i=0 c q i Φ q i (ξ); = n i=0 c e i Φ e i (η); y = m i=0 c y i Φ υ i (ξ).
(5.23)

The input parameter q is introduced to generalize the problem with an additional map from q and x via x = f (q), which can be essential for the latter updating non-negative material parameters. The input parameter q and measurement y are estimated via the same basis functions ξ, while the error is estimated via the different one η. Hence, it is required to fuse them all together so all considered variables are approximated by the same set of basis functions (ξ, η). This transformation is given as follows

q = h i=0 ĉq i Φi (ξ, η) = Q Φ; = h i=0 ĉe i Φi (ξ, η) = Ê Φ; y = h i=0 ĉy i Φi (ξ, η) = Υ Φ.
(5.24) With the current PCE form of all random variables, the Kalman filter can be applied directly as

Q Φ = Q Φ + K( Ẑ Φ -( Υ Φ + Ê Φ)). ( 5

.25)

By removing the common basis functions Φ, the new filtered or assimilated PCE coefficients of input variables are written as

Q = Q + K( Ẑ -( Υ + Ê)), (5.26)
in which the Kalman gain is computed via following equation

K = C QY (C Y + C E ) -1 , ( 5.27)
and corresponding covariance matrices are obtained directly from the PCE coefficients.

Parameter identification in multi-scale problem using Gauss-Markov-Kalman filter

In this section, the parameter identification is configured for a multi-scale problem, e.g a micro-scale is simulated via the Lattice model while a macro-scale is simulated via a mixed triangular element with Raviart-Thomas interpolation. The geometry selected for the multi-scale problem is 2 × 10 in length and height as shown in There are 738 elements in the micro-scale and only 8 elements in the macroscale. As shown in Fig. 5.7, the micro and macro-scales are denoted correspondingly as domain Ω m and Ω M . One-way connection from Ω m and Ω M is the QoI, which is computed from the micro-scale. In general, both bulk and shear moduli (K m (x, y, ω), G m (x, y, ω)) can be considered as random fields in the micro-scale Ω m . In the case of vanishing Poisson ration, the shear modulus is a function of the Young modulus, hence only the Young modulus E m (x, y, ω) is considered as a random field. This random field is then truncated via the Karhunen-Loève method to remove the high-order modes with its correspondingly low eigenvalues. Then one can use any sampling method to generate a set of N m realizations for E m (x, y, ω). The number of simulations is exactly the number of realizations N m . All of the QoI data, which function is used to compute the correlation between a couple of two given points with respect to correlation length l c . Specifically, to compute the Matérn covariance, the values of the non-negative parameter µ c and correlation length l c are selected, respectively, as 1 and 0.4×min(l, h) in which (l, h) are length and height of geometry. Following steps 4 and 5, all computed eigenvalues are shown in Fig. 5.9(a), and the corresponding relative cumulative sum of eigenvalues in Fig. 5.9(b). The first L = 35 eigenvalues contribute to more than 90% (ρ L = 0.9) of the total sum of all eigenvalues. These eigenvalues are selected for coming steps in the algorithm of the Karhunen-Loève Expansion. The truncated form of exponential parameters in random fields are now written as

Θ E (x, ω) = µ Θ E (x) + L i=1 σ i ψ i (x)θ E,i (ω)/ρ L .
(5.29)

The relative cumulative ration ρ L is added in Eq. 5.29 to counter the removed 10% contribution. The first nine eigen functions are shown in Fig. 5.8. Each eigen function covers a full domain of the given geometry. In the first three responses, they look similar to low-order bi-harmonic functions. From the computed eigen functions, the random fields including bulk and shear moduli are generated after truncation, see Eq. 5.12. One realization of the Young modulus E m (x, y, ω) is in Fig. 5.10. It is noted that the realization of Young modulus is computed with the low variance var(Θ). The realizations of truncated Young modulus is inserted into corresponding N m input files for simulation. The geometry, mesh, and boundary conditions remain the same throughout all input files.

Upscaling in the elastic regime

The updating processes are divided into two main phases. The first phase is to identify bulk and shear moduli (K M (ω), G M (ω)) in the elastic regime, while the second one is to assimilate yield stress σ y (ω), linear hardening modulus H h (ω), kinematic hardening modulus H k (ω), and visco-plasticity parameter η(ω) and under the inelastic regime. In the macro-scale Ω M , both bulk and shear moduli are considered as random variables (K M (ω), G M (ω)), which means that they are constant over the entire domain in each simulation for each realization. Firstly, the surrogate model of QoI (stored energy) is prepared as a function of simultaneously both bulk and shear moduli (K M (ω), G M (ω)). To ensure the positive value of both moduli, they are configured as exponential functions of input parameters (q K (ω), q G (ω)), which the latter can be modeled as the Gaussian distribution. With this approach, the distributions of (K M (ω), G M (ω)), in fact, are of Lognormal distribution. The mapping is written as follows

K M (ω) = e Θ K (ω) ; Θ K (ω) = µ Θ K + σ Θ K q(ω); G M (ω) = e Θ G (ω) ; Θ G (ω) = µ Θ G + σ Θ G q(ω).
(5.30)

Firstly, the surrogate model of QoI, which is the stored energy Ê at each time step, is formed via the PCE approach as shown in Alg. 3. The data of QoI is stored in the form of y = Ê1 Ê2 . . . Ên , with n as the number of time steps in each simulation.

Algorithm 3 Generating surrogate model of QoI via PCE 1:

Step 1: Define prior distribution of (q K , q G ) as Gaussian distributions, e.g. K = e q K and G = e q G 2:

Step 2: Specify orthogonal basis polynomials, e.g. multi-variable Hermite functions Φ q α (ξ) Step 5: Update or assimilate PCE coefficients Q of input parameter (q K , qG ) after filtering via K

The numerical experiments are set up to identify the bulk and shear moduli in the macro-scale. Firstly, the shear test is proceeded to identify both moduli. This test would yield a good assimilation value of the shear modulus, since it brings a lot of information of the shear modulus. Afterwards both bulk and shear moduli are considered as random variables in the second experiment where the posterior values of both bulk and shear moduli from previous test are input as prior values for the current expansion test. In the shear test, the load is applied via the imposed displacement within 20 time steps at the top of the domain. Then, the expansion test is conducted via applying imposed displacement perpendicular to both right and top edges. The max imposed displacement in both cases is ū = 0.1. Basically, there are two updating approaches. The first type is to use simultaneously all stored energies and update only once, this method can be called as the simultaneous update. Meanwhile, the second type is named as the sequential update. In detail, only stored energy at the first time step from the micro-scale model is inserted into the Gauss-Markov-Kalman filter as QoI. In the next update, the newly assimilated or the posterior PCE coefficients of the input parameter Q are introduced as the prior values to the next update with corresponding stored energy at the second time step. This procedure is repeated until all stored energies are used. Both approaches are illustrated in Fig. 5.12. In the following, the second approach is employed. The prior values of bulk and shear moduli in the macro-scale are given in Table 5.3. After updating, the assimilated of shear modulus get more confident with a low value of variance, see Fig. 5.14(b). However, it is not the case for bulk modulus, see Fig. 5.14(a). This phenomenon can be explained since the shear test brings less information on the bulk modulus K M . Upscaling an ensemble of E m (x, y, ω) on the micro-scale Ω m To propagate more uncertainty in the micro-scale Ω m , the Young modulus E m is generated as a random field. In micro-scale Ω m , there are N m numerical simulations in total. In each simulation, the stored energy Ê is saved into the data base. At the beginning, the shear test is conducted. The upscaling procedure using the Gauss-Markov-Kalman filter is applied for each set of data from each simulation in mirco-scale. After completing all updates in the shear test, those posterior data are employed as prior data in the expansion test. The same procedure is repeated to update both moduli in the expansion test. The workflow is shown in Fig. 5.18. After repeating the procedure within N m times, the final posterior pdf of both bulk and shear moduli are computed via averaging all assimilated pdf as follows

Property

fK M = 1 N m N m i=1 f i,K M ; fG M = 1 N m N m i=1 f i,G M . ( 5.31)
From these average pdf, the corresponding mean and variance ( μ, σ2 ) are computed.

For the discrete pdf, the formulation of mean and variance are given as

μ = N m i=1 fi x i ; σ2 = N m i=1 fi (x i -μ) 2 .
(5.32)

It is noted that the average pdf f is normalized so that N m i=1 fi = 1.

Figure 5.18 -Flow of data in both tests

Since the shear test gives little information on the bulk modulus K M , the pdf remains nearly the same in every update. Meanwhile the variance of shear modulus G M shrinks in each update. The posterior from each update in the previous shear test is used as corresponding prior pdf for the expansion test. The average pdf of all updates is given in Table 5.6, in which both mean and variance are computed from the average pdf. To zoom in the average pdf of both moduli, their values are scaled up 5 times in Fig. 5.19. Compared to the updates of (K M , G M ) from using one realization of E m , the new variances are considered bigger but the mean values remain nearly comparable. Table 5.6 -Bulk and shear moduli (K M , G M ) after both tests

Upscaling in the inelastic regime

To begin with, the inelastic models used in both scales are briefly described. Additionally, random fields on the micro-scale and targeted random variables on the macro-scale are selected. Only important equations are listed in the following section.

Inelastic model in the micro-scale Ω m On the micro-scale, the Lattice model is simulated via the Timoshenko beam capable of linear isotropic hardening plasticity and nonlinear kinematic hardening viscoplasticity, see [START_REF] Hadzalic | Theoretical formulation and seamless discrete approximation for localized failure of saturated poro-plastic structure interacting with reservoir[END_REF][START_REF] Hadžalić | Analysis of pore pressure influence on failure mechanisms in structural systems[END_REF]. It is noted that the additive decomposition of the total strain is applied only to axial and shear strains ( , γ), not for the curvature κ. The strain energies due to moment, axial and shear forces are given as

ψ m (κ) = 1 2 κIA e κ; ψ a ( , vp , ξ a ) = 1 2 ( -vp ) EA e ( -vp ) + 1 2 ξ a H a h A e ξ a ; ψ s (γ, γ vp , ξ s ) = 1 2 (γ -γ vp ) k c GA e (γ -γ vp ) + 1 2 ξ s H s h A e ξ s , (5.33)
in which (ξ a , ξ s ) are internal hardening variables of axial and shear forces, and (H a h , H s h ) are hardening moduli for axial and shear forces. It is assumed that there are independent plasticity mechanisms activated by axial and shear forces. The independent yield functions (φ a , φ s ) for axial and shear forces (N, Q) are

φ a (N, χ a , q a ) = |N -χ a A e | -(N y -q a A e ) ≤ 0; φ s (Q, χ s , q s ) = |Q -χ s A e | -(Q y -q s A e ) ≤ 0, (5.34)
in which (χ a , χ s ) are back-stress variables. The yield axial and shear forces (N y , Q y ) are computed from the Lattice yield stress σ y . For linear hardening, the stress-like isotropic hardening variables (q a , q s ) are given as where the relevant internal variables are listed as ξ vp and ζ vp . We also consider additional hardening effects as a general case by including the corresponding potentials Ξ vp 1 and Ξ vp 2 . For those potential energies, we select a quadratic form as follows

q a = -H a h ξ a ; q s = -H s h ξ s . ( 5
Ξ vp 1 (ξ vp ) = 1 2 ξ vp K vp ξ vp ; q vp = - ∂Ξ vp ∂ξ vp ; Ξ vp 2 (ζ vp ) = 1 3 H vp ζ vp : ζ vp ; α vp = - ∂Ξ vp ∂ζ vp , ( 5.39)
in which the terms K vp and H vp are the linear hardening and kinematic hardening moduli. Correspondingly, the terms q vp and α vp are the stress-like variable and back-stress tensor in the visco-plastic model. The yield criterion is now written as φ vp (σ, q vp , α vp ) := dev (σ) + α vp -(σ y -q vp ) ≤ 0.

(5.40)

The total dissipation can be written in the following form 0 < D = σ : ˙ vp + q vp ξvp + α vp : ζvp . The list of material properties considered as random fields is given in Table 5.9, e.g. yield stress, linear hardening modulus, nonlinear kinematic hardening modulus, and visco-plastic parameter. These material properties are chosen to be updated via the Gauss-Markov-Kalman filter method. The workflow in the elastic regime is adapted here.

Property

Random variables Yield stress σ M y (ω) Linear hardening modulus H M h (ω) Nonlinear kinematic hardening modulus H M k (ω) Visco-plastic parameter γ M (ω) Table 5.9 -Macro-scale material properties as random variables

Upscaling one realization of (σ m y,t (ω), σ m y,c (ω), σ m y,s (ω), H m h (ω), H m k (ω), η m (ω)) on the micro-scale Ω m The boundary condition and loads are illustrated in Fig. 5.21(a). The bottom edge is fixed, while an imposed displacement ū is applied on the top edge. The first cycle for an imposed displacement u 1 takes place within 0 ≤ t ≤ t 4 , see Fig. 5.21(b). The second and third ones are defined via t 4 ≤ t ≤ t 8 and t 8 ≤ t ≤ t 12 , respectively. The time step in each cycle is proportional to its loading magnitude so that the rate of loading is maintained, e.g. ∆t 3 = 1.4∆t 1 and ∆t 2 = 1.2∆t 1 . In both cases, the time step in the first cycle is selected with a value of ∆t 1 = 0.1. The load magnitudes are given as ū1 = 0.04, ū2 = 0.08 and ū3 = 0.12, respectively. The response in the micro-scale is shown in Fig. 5.21(c). The max stored energies and dissipations in three cycles are selected as QoI for the identification process, see Table The prior values of all targeted material variables in the macro-scale are given in Table 5.11. After updating, it is observed that the visco-plastic parameter γ = 1/η is successfully updated with a smaller variance. The pdf of other parameters remains nearly the same, which means that the current experimental test does not contain much information about these parameters.

Property

Conclusions

Chapter 2 illustrates a detailed picture on the new approach to solve the viscoplasticity problem, in which the Hellinger-Reissner variational principle and the hybrid-stress finite element approximation are utilized to deliver a high accuracy of the stress field. Both displacement and stress fields are formulated as independent variables, which increases the problem size with respect to the conventional displacement-based finite element approach. However, the approach is capable of providing high stress accuracy even with coarse mesh settings. Other remarks of the provided approach lay at the desirable smoothness of the stress field, which is of special interest for capturing the wave propagation phenomena. The energy conserving and decaying schemes are developed with a feature to control the numerical computation in the dynamics problem. The proposed algorithm presents satisfying performances for a cluster of numerical simulations. It is clear that the computation is unconditionally stable by using the energy conserving scheme under the elastic phase. Meanwhile, the energy decaying scheme can dissipate high frequencies modes which can not be resolved properly on a coarse finite element mesh. The dissipation parameters α and β in decaying scheme can be selected with respect to a mesh size. Practically, the higher values of α and β dissipate faster inappropriate highfrequency modes in a coarse-mesh domain. On the other hand, these parameters can be chosen relatively lower for a fine-mesh domain which has fewer unresolved high frequencies and thus naturally reduce the numerical dissipation. Chapter 3 presents the mixed variational formulation suitable for constructing hybrid-stress discrete approximation with independent displacement and stress fields. Particularly, the reduced variational formulation is derived by eliminating the rotation field from the regularized functional [START_REF] Hughes | On drilling degrees of freedom[END_REF] for drilling rotations. Hence, the new variational format is able to include both stress and displacement as independent fields. The lowest Raviart-Thomas vector space is introduced to discretize the stress field with an aim to maintain the continuity of traction vector across element boundaries. Meanwhile the displacement field is interpolated via linear shape functions. The proposed triangular finite element for two-dimensional planar plate has better performance than some classic low-order finite elements for either elastic or visco-plastic behaviors. In the dynamics framework, the energy-conserving scheme is applied to control the overall stability of numerical computation over a long period of time. The proposed formulation and algorithm are examined under many numerical simulations and found with the satisfying performance in both statics and dynamics problems.

Chapter 4 performs an attempt to adapt the weak variational formulation in Chapter 3 into the plate bending problem. The same approach is applied to derive the reduced variational formulation. In this effort, all displacement (including lateral deflection and rotations), moments and shear forces fields are considered as independent variables at the expense of increasing the problem size with respect to conventional displacement-based finite element approach. Nevertheless, it also provides the possibility to maintain the desirable projection continuity of moments and shear forces fields over element's edges, which is of special interest for capturing the wave propagation phenomena. The energy conserving and decaying schemes are also reformulated to adapt with the stated plate bending problem. The same positive observation on those time integration schemes is gained via several numerical simulations. The proposed formulation has a potential on other elasto-plastic and dynamic problems.

Chapter 5 depicts the feasibility study on the application of parameter identification methods, including the Bayesian update and the Gauss-Markov-Kalman filter, into inelastic problems. In particular, the Bayesian update method is implemented to identify the fracture energy G f in the simple one dimensional problem. Then, the Gauss-Markov-Kalman filter is applied in the multi-scale problem. The Lattice model is used to simulate the micro-scale with a random field such as Young modulus E m (x, y, ω), while the mixed triangular membrane element is employed for the macro-scale. The finite elements using in both scales are embedded with the visco-plasticity. The main idea here is to assimilate the material's properties in the macro-scale Ω M by using the stored energy and dissipation as quantity of interests from the micro-scale Ω m . This procedure can be also viewed as an upscaling homogenization procedure.

Perspectives

A new shell finite element can be formed by adding membrane and plate bending finite elements together, which are developed in Chapter 3 and 4 respectively. A thorough investigation on the performance and accuracy of that shell finite element can be an worth-trying topic. In addition, the application of Raviart-Thomas vector space for stress discretization can be applied on quadrilateral, tetrahedral and brick finite elements for nonlinear and inelastic problems, specially for wave propagation problems. The weak form can be derived using the regularized mixed variational functional as presented in this thesis. Beside that, the application of parameter identification process is still open for many other problems in the field of computational mechanics. The Gauss-Markov-Kalman filter can be implemented along with neural networks to compare the efficiency between the two methods.

A

Mapping constitutive equation into matrix notation

For non-symmetric stress tensor, the constitutive equation is

   ij=32    1 0 0 1 0 0    ij=33                           (A.3)
Regarding fourth-order tensor I s = 1 2 (δ ik δ jl + δ il δ jk )

I s =                              1
   ij=33                           (A.5)
After executing an operation of double contraction, the constitutive equation symmσ = C : is now written as follows.

B

Explicit form of matrix Db relevant to moment field

The explicit form of matrix D b symm M in Eq. (4.15) can be written as follows

D b symm M = 12(1 -ν 2 ) Et 3             -1 ν 2 -1 ν ν 2 -1 0 0 ν ν 2 -1 -1 ν 2 -1 0 0 0 0 1 1 -ν 0 0 0 0 1 1 -ν                       M xx M yy 1 2 (M xy + M yx ) 1 2 (M yx + M xy )           (B.1)
With the selection = 2µ hence 1 = 1 + ν E , the explicit form of matrix 1 12 t 3 skew M in Eq. (4.15) can be written as follows

1 12 t 3 skew M = 12(1 -ν 2 ) Et 3             -1 ν 2 -1 ν ν 2 -1 0 0 ν ν 2 -1 -1 ν 2 -1 0 0 0 0 1 1 -ν 0 0 0 0 1 1 -ν                       0 0 1 2 (M xy -M yx ) 1 2 (M yx -M xy )           (B.2)
Finally, the explicit form of considered term Db can be revealed in the following equation

D b symm M - 1 12 t 3 skew M = 12 t 3             1 E -ν E 0 0 -ν E 1 E 0 0 0 0 0 1 + ν E 0 0 1 + ν E 0             Db        M xx M yy M xy M yx        := Db M (B.3)

C

Implementation of matrices in Gauss-Markov-Kalman filter

The random variables (q K (ω), q G (ω)) are written as

q (ξ 1 , ξ 2 ) = q K q G = log K log G = α∈J (2) l c q α Φ q α (ξ 1 , ξ 2 ) ; c q α ∈ R 2×|J l | , (C.1)
for an expansion up to degree l in two random variables (ξ 1 , ξ 2 ) ∈ N (0, 1). Take an expansion in J 3 , i.e. l = 3 with α = (α 1 , α 2 ), then J 3 = {(0, 0), (1, 0), (0, 1), (2, 0), (1, 1), (0, 2), (3, 0), (2, 1), (1, 2), (0, 3)}, (C.2) and the corresponding matrix Φ q α ∈ R 10×1 containing basis functions or germs is written as

Φ q α = 1 ξ 1 ξ 2 ξ 2 1 -1 ξ 1 ξ 2 ξ 2 2 -1 ξ 3 1 -3ξ 1 (ξ 2 1 -1)ξ 2 ξ 1 (ξ 2 2 -1) ξ 3 2 -3ξ 2 T . (C.3)
The full coefficient matrix is c q α ∈ R 2×10 . However, at the beginning the random variables (q K (ω), q G (ω)) are of Gaussian distributions as q K = c q (0,0),1 + c q (1,0),1 ξ 1 ; q G = c q (0,0),2 + c q (0,1),2 ξ 2 , (C.4)

or precisely, the first prior only the coefficients c q α are non-zero at these indices α = (α 1 , α 2 ) = (0, 0), (1, 0), and (0, 1). (C.5)

After the updating process, the coefficients c q α which are zero at the beginning may be updated to nonzero values. The number of QoI (n) corresponds to the number of time steps, at which the stored energy is computed. The synthetic errors (ω) of numerical measurement of QoI are written in PCE format as white noise for an expansion up to degree 1 in n random variables η ∈ N (0, 1). Take an expansion in J 1 , i.e. with β = (β 1 , β 2 , ..., β n ), then J 1 = {(0, 0, ..., 0), (0, 1, ..., 0), ..., (0, 0, ..., 1)}, (C.7)

and the corresponding matrix Φ β ∈ R (n+1)×1 containing basis functions or germs is written as

Φ β = 1 η 1 η 2 . . . η n T .
(C.8)

The corresponding explicit matrix (η) is written in PCE form as follows

=       1 2 . . . n       =       0 σ 1 0 • • • 0 0 0 σ 2 • • • 0 . . . . . . . . . 0 0 0 • • • σ n                1 η 1 η 2 . . . η n         
.

(C.9)

The full coefficient matrix is c β ∈ R n×(n+1) . The stored energy y ∈ R n×1 (QoI) is written in PCE format, with the same germs as of input parameters, as follows in which υ ≡ c y α , e.g. υ i,1 = c y (0,0),i , υ i,2 = c y (1,0),i and similarly υ i,10 = c y (0,3),i . Before applying the PCE method, all random variables should be written in a same set of basis functions, which is established via fusing Φ q α (ξ) and Φ β (η). The new form of combined basis Φ(ξ, η) ∈ R (n+10)×1 containing all PCE polynomials is written as

y (ξ 1 , ξ 2 ) =
Φ = Φ q α ∪ Φ β = 1 ξ 1 ξ 2 . . . ξ 3 2 -3ξ 2 η 1 η 2 . . . η n T .
(C.12)

The corresponding coefficient matrices Ê, Q and Υ with the same germs are written as Q = 0 σ q K 0 0 . . . 0 0 . . . 0 0 0 σ q G 0 . . . Consequently, the Kalman gain is obtained via Eq. (5.26). Then, the new filtered or assimilated PCE coefficients of input variables can be achieved via Eq. (5.27).

Contents 1 . 1 2 1. 2 3 1. 3 5 1. 4

 11223354 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . Literature review . . . . . . . . . . . . . . . . . . . . . . . Research objectives and methodology . . . . . . . . . . . Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . 6

Figure 1 . 1 -

 11 Figure 1.1 -Advanced numerical simulation and applications: (a) Research concept on future aircraft development (courtesy of dlr.de), (b) NASA's Orion parachute validation (courtesy of nasa.gov), (c) Wind turbine tower at Sainte-Colombe, France.

Figure 2 . 1 -

 21 Figure 2.1 -Domain Ω and its boundary ∂Ω

Figure 2 . 2 -

 22 Figure 2.2 -Computed stress under elastic regime at p 0 = 5.10 3

Figure 2 . 3 -

 23 Figure 2.3 -Computed stress under visco-plastic regime by constant p 0 = 21.10 3

Figure 2 . 4 -

 24 Figure 2.4 -A quarter of a pipe by Q4/P1 element

Figure 2 . 6 -

 26 Figure 2.6 -Energy conserving under elastic regime (EC scheme)

Figure 2 . 7 -

 27 Figure 2.7 -Displacement under elastic regime (EC scheme)

6 Figure 2 . 8 -

 628 Figure 2.8 -Velocity under elastic regime (EC scheme)

6 Figure 2 . 9 -Figure 2 . 10 -

 629210 Figure 2.9 -Stress under elastic regime (EC scheme)

6 Figure 2 . 11 - 1 Figure 2 . 12 -

 62111212 Figure 2.11 -Stress under visco-plastic regime (EC scheme)

6 Figure 2 . 13 -Figure 2 . 14 -

 6213214 Figure 2.13 -Stress under elastic regime (ED scheme)

  2.15. Similar to the previous simulation, the second pulse bring only little noise to the stress in soft model, as shown in Figs. 2.15(a,c,e). The same observation on the energy-decaying on high frequency mode's contribution no longer applies for the stiff model, as shown in Figs. 2.15(b,d,f).

6 Figure 2 . 15 -

 6215 Figure 2.15 -Stress under visco-plastic regime (ED scheme)

Figure 2 . 16 -

 216 Figure 2.16 -Energy under elastic regime (stiff model-ED scheme)

Contents 3 . 1 36 3. 2 39 3. 2 . 1 39 3. 2 . 2 41 3. 2 . 3 42 3. 3 44 3. 3 . 1 45 3. 3 . 2 47 3. 3 . 3 50 3. 3 . 5

 313623921392241234234431453247335035 Regularized variational formulation . . . . . . . . . . . . Numerical implementation . . . . . . . . . . . . . . . . . . Finite element approach . . . . . . . . . . . . . . . . . . . Discrete weak form . . . . . . . . . . . . . . . . . . . . . . Energy conserving scheme . . . . . . . . . . . . . . . . . . Numerical examples . . . . . . . . . . . . . . . . . . . . . . Patch test . . . . . . . . . . . . . . . . . . . . . . . . . . . Higher order patch test . . . . . . . . . . . . . . . . . . . Cantilever beam . . . . . . . . . . . . . . . . . . . . . . . 48 3.3.4 Cook's membrane . . . . . . . . . . . . . . . . . . . . . . Dynamic loadings . . . . . . . . . . . . . . . . . . . . . . 53

Figure 3 . 1 -

 31 Figure 3.1 -Mapping of isoparametric triangle Ω to physical triangle Ω e

  (a) Model 1: Homogeneous material (b) Model 2: Non-homogeneous material

Figure 3 . 2 -

 32 Figure 3.2 -The patch test

(a) Case 1 :Figure 3 . 3 -

 133 Figure 3.3 -Model 1-Homogeneous material

  Figure 3.4 -Model 2-Non-homogeneous material

  Figure 3.5 -A simple beam under bending state

  Figure 3.6 -Displacement field and deformed configuration

  (a) Traction τ x , max=5.72 (b) Traction τ y , max=5.58 (c) Traction τ , max=7.02 (d) Hydrostatic pressure

Figure 3 . 7 -3. 3 . 3

 3733 Figure 3.7 -Traction field and hydrostatic pressure-Case 2

  Figure 3.8 -A cantilever beam

  (a) Case 1: Displacement u x (b) Case 2: Displacement u x (c) Case 1: Displacement u y (d) Case 2: Displacement u y (e) Case 1: Deformed shape (f) Case 2: Deformed shape

Figure 3 . 9 -

 39 Figure 3.9 -Displacement field and deformed configuration

  Figure 3.10 -Traction field and hydrostatic pressure-Case 2

Figure 3 .

 3 Figure 3.11 -Cook's membrane

(a) Case 1 :Figure 3 . 13 -

 1313 Figure 3.13 -Displacement field and deformed configuration

  (a) Distribution of yield function φ(σ) (b) Hydrostatic pressure

Figure 3 .

 3 Figure 3.15 -Cook's membrane under visco-plastic regime

Figure 3

 3 Figure 3.16 -Dynamic test

(a) Case 1 :Figure 3 . 17 -

 1317 Figure 3.17 -Model 1-Displacement and traction

(a) Case 1 : 4 Figure 3 . 18 -

 14318 Figure 3.18 -Model 1-Total energy

(a) Case 1 : 4 Figure 3 . 19 -

 14319 Figure 3.19 -Model 2-Total energy

(a) Case 1 :

 1 Displacement at node 123 (b) Case 1: Traction at node 82 (c) Case 2: Displacement at node 123 (d) Case 2: Traction at node 82 (e) Case 3: Displacement at node 123 (f) Case 3: Traction at node 82 (g) Case 4: Displacement at node 123 (h) Case 4: Traction at node 82

Figure 3 . 20 -

 320 Figure 3.20 -Model 2-Displacement and traction

Figure 3 . 21 -

 321 Figure 3.21 -Model 1-Case 3: Dynamics under visco-plastic regime

Figure 4 . 1 -, + t 2 ,

 412 Figure 4.1 -Plate and sign convention

Figure 4 . 2 -

 42 Figure 4.2 -Assumed shear strain γ

Figure 4 . 3 -

 43 Figure 4.3 -Mapping of isoparametric triangle Ω to physical triangle Ω e

N

  T p dA = 0; G M,Q (•; δM, δQ) =δm T Ω e S T B κ ddA -Ω e S T κ vp dA + Ω e S T Db SmdA +δq T Ω e X T B γ ddA -Ω e X T γ vp dA + Ω e X T D s XqdA = 0(4.44) Since the variations δd T , δm T and δq T are arbitrary, the term inside each bracket in the above equation should be zero. The linearization of the equation Eq. (4.44a) yields r e,(i+1) d,n+1/2 (d, m, q) = r e,(i) d,n+1/2 (d, m, q) + ∆r e,(i) d,n+1/2 (∆d, ∆m, q) T N dV ; η := ∂a n+1/2 /∂d n+1/2 ; G e := Ω e X T B γ dA; F e := Ω e S T B κ dA; f ext =

(4. 56 )

 56 By choosing a test deflection vector δd = d n+1 -d n ≡ u, the work done by both internal forces f int and external forces f ext can be obtained from the weak form of balance equation in Eq. (4.45b)

  (a) Deflection w (b) Reaction moment Myy

Figure 4 . 4 -

 44 Figure 4.4 -Model 1-Plate response under only pulse T P 1 (proposed element vs [11])

Figure 4 . 5 -

 45 Figure 4.5 -Model 1-Moment M yy along plate under only pulse T P 1

Figure 4 . 6 -

 46 Figure 4.6 -Model 1-Plate response under EC and elastic regime

Figure 4 . 7 -

 47 Figure 4.7 -Model 1-Time history of plate degrees of freedom

Figure 4 . 8 -

 48 Figure 4.8 -Model 1-Time history of stress resultants

Figure 4 . 9 -

 49 Figure 4.9 -Model 1-Total energy and dissipation due to T P 1 (n ≤ 250)

Figure 4 . 10 -

 410 Figure 4.10 -Model 1-Total energy and dissipation due to T P 1 , T P 2 (n ≤ 800)

  Rotation θ x (d) Rotation θ y

Figure 4 . 11 -

 411 Figure 4.11 -Model 2-Plate response under EC and elastic regime

  Shear d.o.f. q

Figure 4 . 12 -

 412 Figure 4.12 -Model 2-Time history of plate degrees of freedom

Figure 4 . 13 -

 413 Figure 4.13 -Model 2-Time history of stress resultants

Figure 4 . 14 -

 414 Figure 4.14 -Model 2-Total energy and dissipation due to T P 1 (n ≤ 250)

Figure 4 . 15 -

 415 Figure 4.15 -Model 2-Total energy and dissipation due to T P 1 , T P 2 (n ≤ 800)

  Rotation θ x (d) Rotation θ y

Figure 4 . 16 -

 416 Figure 4.16 -Model 3-Plate response under EC and elastic regime

  Shear d.o.f. q

Figure 4 . 17 -

 417 Figure 4.17 -Model 3-Time history of plate degrees of freedom

Figure 4 . 18 -

 418 Figure 4.18 -Model 3-Time history of stress resultants

Figure 4 . 19 -

 419 Figure 4.19 -Model 3-Total energy and dissipation due to T P 1 (n ≤ 250)

Figure 4 . 20 -

 420 Figure 4.20 -Model 3-Total energy and dissipation due to T P 1 , T P 2 (n ≤ 800)

Figure 4 . 21 -

 421 Figure 4.21 -Model 3-Visco-plastic dissipation with higher γ

Figure 4 . 22 -

 422 Figure 4.22 -Model 3-Visco-plastic and algorithmic dissipation with different σ y

8 0

 8 γ = 1 × 10 -4 and α = β = 1 × 10 -γ = 1 and α = β = 1 × 10 -4

Figure 4 . 23 -

 423 Figure 4.23 -Model 3-Dissipation with higher γ and lower α, β

Figure 4 . 24 -

 424 Figure 4.24 -Continuity of d.o.f. m x and m y over entire domain (EC and elastic regime)

. 5 )

 5 Several plots for those common distribution are shown in Fig. 5.1. The Normal distributions with different values of mean value µ and standard deviation σ is shown Fig. 5.1(a), while the Log-Normal distribution is shown in Fig. 5.1(b).

( a )Figure 5 . 1 -

 a51 Figure 5.1 -Normal and Log-Normal distributions

  (a) n s = 10 2 (b) Histogram with n s = 10 2 (c) n s = 10 3 (d) Histogram with n s = 10 3

Figure 5 . 2 -

 52 Figure 5.2 -Generate realization of random variable x ∼ N (10 3 , 10) via SimP aramSet object

. 15 )

 15 where π Y (y) is the probability density function of the random variable y (the evidence) and π Q (θ) is the prior probability density function θ, and π (Y |Q) (y | θ) is the likelihood of y = Y (θ, ) given θ. The prior probability density function π Q (θ) describes a belief of an analyst about the possible population characteristics of random variable θ. The likelihood π (Y |Q) (y | θ) describes our belief that the observations y if we know θ is true. And the posterior π (Q|Y ) (θ | y) describes our belief on the possible population characteristics of random variable θ after observing y.

Figure 5 . 3 -

 53 Figure 5.3 -Flow of Bayesian updating method

Value 1 30 ×

 30 10 6 30 × 10 3 15 × 10 6 60 × 10 3 7.5 × 10 6 4 × 10 2

Figure 5 . 4 -Figure 5 . 5 - 2 -

 54552 Figure 5.4 -Bayesian updating G f using reaction R x

) and the polynomials are normalized via R Φ 2 i

 2 (α)f ξ (α) dα = 1.(5.18)

Fig. 5 .

 5 6(a). The corresponding mesh in the micro-scale and the macro-scale are shown in Fig.5.6(b,c).

Figure 5 . 6 -

 56 Figure 5.6 -Multi-scale setting: (a) Geometry, (b) Micro-scale: Lattice model, (c) Macro-scale: Enhanced triangular model

Figure 5 . 8 -

 58 Figure 5.8 -First nine eigen functions

Figure 5 . 9 - 9 Algorithm 2 1 :

 59921 Figure 5.9 -Eigenvalues and desired relative contribution ρ L = 0.9

Figure 5 . 10 -

 510 Figure 5.10 -One realization of random field Young modulus E m (x, y, ω)

Figure 5 . 11 -

 511 Figure 5.11 -Experiment setup in elastic regime to identify (K M , G M )

Figure 5 . 12 -

 512 Figure 5.12 -Simultaneous and sequential update approaches

Figure 5 . 13 -

 513 Figure 5.13 -One ensemble of E m : Stored energy responses

Figure 5 . 19 -

 519 Figure 5.19 -Posterior and average pdf after both tests

8 -

 8 Micro-scale: One realization values of (σ m y,t , σ m y,c , σ m y,s , H m h , H m k , η m ) (a) Imposed displacement ū = 0.1 (b) Hysteresis loop

Figure 5 . 20 -

 520 Figure 5.20 -One realization of (σ m y,t , σ m y,c , σ m y,s , H m h , H m k , η m ): Micro-scale responses in inelastic regime by imposed displacement

(5. 41 )

 41 Following the penalty method as in Chapter 3, the evolution equations of this viscoplastic model are written as follows0 = ∂D vp η (σ, •) ∂σ = -˙ vp + 1 η vp < φ vp > ∂φ ∂σ ⇒ ˙ vp = 1 η vp < φ vp > n = f ; 0 = ∂D vp η (σ, •) ∂α vp = -ζvp + 1 η vp < φ vp > ∂φ ∂α vp ⇒ ζvp = 1 η vp < φ vp > n = f ζ ; 0 = ∂D vp η (σ, •) ∂q vp = -ξvp + 1 η vp < φ vp > ∂φ ∂q vp ⇒ ξvp = 1 η vp < φ vp >= f ξ ,(5.42)in which the notation n is the unit normal tensor n = dev (σ) + α vp dev (σ) + α vp . The numerical implementation of this model is straightforward as the one in Chapter 3.

  .

  5.10.

Figure 5 . 21 -

 521 Figure 5.21 -Experiment setup with imposed displacement ū in inleastic regime to identify (σ M y , H M h , H M k , γ M )

Figure 5 . 3 Table 5 . 10 -

 53510 Figure 5.22 -Micro-scale Ω m responses

  (a) Yield stress σ M y (b) Modulus H M h (c) Modulus H M k (d) Visco-plastic parameter γ M

Figure 5 . 23 -

 523 Figure 5.23 -Prior and posterior pdf using both stored energy and dissipation as QoI

  Upscaling an ensemble of (σ m y,t (x, y, ω), σ m y,c (x, y, ω), σ m y,s (x, y, ω), H m h (x, y, ω), H m k (x, y, ω), η m (x, y, ω)) on the micro-scale Ω m On the micro-scale, the parameters (σ m y,t , σ m y,c , σ m y,s , H m h , H m k , η m ) are generated as lowvariance random fields in the micro-scale Ω m . The mean values of these random variables are given in Table 5.8. The procedure takes place as shown in Fig. 5.18, where the stored energy and dissipation are used as QoI. The result yields agreement with updated values using one realization of (σ m y,t , σ m y,c , σ m y,s , H m h , H m k , η m ).

  Visco-plastic parameter γ M

Figure 5 . 24 -

 524 Figure 5.24 -Prior and posterior pdf using an ensemble of random fields

symmσ = C : ⇔ 1 2 ((A. 2 )

 22 σ ij + σ ji ) = C ijkl kl (,j + u j,i ); i, j ∈ [1 : 3] C = λ1 ⊗ 1 + 2µI s ⇔C ijkl = λδ ij δ kl + 2µ 1 2 (δ ik δ jl + δ il δ jk ); i, j, k ∈ [1:3]Regarding fourth-order tensor 1 ⊗ 1 = δ ij δ kl

  -order elasticity tensor C = λ1 ⊗ 1 + 2µI s

1 2 (σ 21 + σ 12 ) 7 ) 8 ) 1

 12112781 The above equation leads to the form of inverse of elasticity matrix C -1 .Applying the same mapping L for the second-order tensor 1 skew σ with = 2µ = E second-order tensor D in Eq. (3.19) can be written as followsL(C -1 symmσ -1 skewσ) = C -1

1 c

 1 β Φ β (η) ; c β ∈ R n×(n+1) , (C.6)

  α (ξ 1 , ξ 2 ) ; c y α ∈ R n×|J l | , (C.10)or in explicit matrix form, with l = 3

  .1. The basis functions Φi are mapped back to physical coordinate space by the Piola transformation represented with an affine function F . Since one common edge between two elements in physical space can be mapped into two different edges (with different lengths) in the parent coordinate, the terms |e k | and |ê k | (length of considered edges in the physical and parent spaces, respectively) are added into Eq. (3.27) to maintain Eq.(3.26) in the physical coordinate, see[START_REF] Ervin | Computational bases for rtk and bdmk on triangles[END_REF]. The notation J T and |J T | are the Jacobian matrix and its determinant of the corresponding affine function. Hence, the basis function at edge i in the physical coordinate system is now written as follows Φi

Table 3 .

 3

1 -Load groups for EC scheme

  /2 σ xx z dz

	t/2 -t/2 σ yx z dz	t/2 -t/2 σ xy z dz t/2 -t/2 σ yy z dz	  .	(4.6)

Table 4 .

 4 , 10 -8 2.10 3 , 10 -8 2.10 3 , 10 -8 2.10 3 , 10 -8 2.10 3 , 10 -8 1 -Parameters for EC and ED schemes in elastic or visco-plastic regimes

	Case	1	2	3	4	5
	α, β	0,0	0.01,0.01	0.02,0.02	0.03,0.03	0.08,0.08
	σ y , γ 2.10 3 Scheme EC	ED	ED	ED	ED

3 :

 3 Step 3: Compute squared norms and generate integration points with weights (ξ, w) 4: Step 4: Compute measurable response QoI by FEM solver at integration points 5: Step 5: Obtain PCE coefficients Υ of QoI y After getting the PCE model of QoI, the Gauss-Markov-Kalman filter can be proceeded with step by step shown in Alg. 4. The explicit form of relevant matrices are illustrated in Appendix. C.

	Algorithm 4 Updating PCE coefficients of input parameters q(q K , q G )
	1: Step 1: Define error model (η) in PCE form
	2: Step 2: Fuse different germs into a unified form Φ(ξ, η)
	3: Step 3: Update forms of random variables Ê, Q and Υ compatible with corre-
	sponding new germs Φ(ξ, η)
	4: Step 4: Compute Kalman gain K via Eq. (5.27)
	5:

Table 5 .

 5 5 -Updated bulk and shear moduli (K M , G M ) after both tests

	Prior µ	Prior σ Posterior μ Posterior	σ
	Bulk modulus K M 1.66 × 10 4	182	1.63 × 10 4	1
	Shear modulus G M 7.69 × 10 3	124	4.54 × 10 3	1.33

  In the following, all linear hardening moduli are assumed to be one unified value of H h to reduce the total number of input variables. Similarly, the values of (H a nk , H s nk ) take one unified value H k . The list of material properties considered as random fields are given in Table5.7, e.g. yield stress in tension, yield stress in compression, yield stress in shear, linear hardening modulus, nonlinear kinematic hardening, and viscosity coefficient. The procedure to generate the random field (RV) is similar in Section 5.2.1. ) are given in Table.5.8. It can be seen that the stored energy and dissipation keep increasing in magnitude under loading scenarios. However, they remain nearly the same in unloading scenarios.

	Property		Random field
	Yield stress in tension	σ m y,t (x, ω)
	Yield stress in compression	σ m y,c y(x, ω)
	Yield stress in shear	σ m y,s (x, ω)
	Linear hardening modulus	H m h (x, ω)
	Nonlinear hardening modulus H m k (x, ω)
	Viscosity coefficient	η m (x, ω)
	Table 5.7 -Micro-scale material properties as random fields
	An example of this Lattice model undergoing a simple shear test ū = 0.1 with
	one realization of (σ m y,t , σ m y,c , σ m y,s , H m h , H m k , η m ) is shown in Fig. 5.20. The material
	properties of (σ m y,t , σ m y,c , σ m y,s , H m h , H m k , η m
						.35)
	The Fredrick-Armstrong [6] nonlinear kinematic hardening law is employed as
	χa = H a k ˙ νp -H a nk	ξa χ a ; χs = H s k	γvp -H s nk	ξs χ s .	(5.36)
	The total dissipation can be written in the general form
	0 < D = σ : ˙ -ψ.			(5.37)

  Visco-plastic parameter γ M 1 × 10 -5 1.4 × 10 -7 2.7 × 10 -5 5.4 × 10 -8

		Prior µ	Prior σ	Posterior μ Posterior	σ
	Yield stress σ M y	1	3.2 × 10 -3	0.8	6 × 10 -4
	Modulus H M h	4 × 10 4	282.8	6.6 × 10 4	227.3
	Modulus H M k	600	34.6	548.9	31.9

Table 5 .

 5 11 -Prior and posterior values of (σ M y , H M h , H M k , γ M ) using both stored energy and dissipation as QoI

  (σ 12 + σ 21 ) 1 2 (σ 13 + σ 31 ) 11 + λ 22 + λ 33 µ 12 + µ 21 2µ 12 µ 13 + µ 31 µ 21 + µ 12 λ 11 + (λ + 2µ) 22 + λ 33 µ 23 + µ 32 µ 31 + µ 13 µ 32 + µ 23 λ 11 + λ 22 + (λ + 2µ) 33Applying the plane strain condition 13 , 23 , 33 = 0, the mapping L from tensor to matrix notation yields

				σ 11
			   	σ 22 1 2 (σ 12 + σ 21 )
		  	1 2 1 σ 11 2 (σ 21 + σ 12 ) 1 2 (σ 31 + σ 13 ) 1 2 (σ 32 + σ 23 ) σ 22	1 2 (σ 23 + σ 32 ) σ 33	  
			(λ + 2µ) 
	=	   			   
					(A.6)

  • • • υ 1,10 0 0 . . . 0 υ 2,1 υ 2,2 • • • υ 2,10 0 0 . . . 0 . . . n,1 υ n,2 • • • υ n,10 0 0 . . . 0

								
	. . .	. . .	. . .	. . .	. . .	. . .	. . .	    	; Υ ∈ R n×(n+10) .	(C.15)

υ

Deutsch Luft und Raumfahrt

National Aeronautics and Space Administration

Acknowledgements

Remerciements Danke

can be stored energy Ê and dissipation D, at each time step are saved for further updating procedure in the macro-scale Ω M . Figure 5.7 -Flow of data and update for upscaling in elastic regime Once the data QoI computed in the micro-scale is successfully processed, the updating procedure by the Gauss-Markov-Kalman filter in the macro-scale Ω M can be started. This process can be viewed as a homogenization work, in which both bulk and shear moduli (K M (ω), G M (ω)) are parameters to be identified in the macro-scale Ω M . It is noted that those parameters are implemented as random variables, which are constant over the entire domain Ω M . The surrogate model using the polynomial chaos expansion (PCE) is applied for QoI in the macro-scale Ω M , which is necessary for the updating process by the Gauss-Markov-Kalman filter. The QoI from both scales are delivered into the procedure of the Gauss-Markov-Kalman filter with an aim to update the bulk and shear moduli of current macro-scale (K M , G M ). All key tasks for this multi-scale problem are drawn in Fig. 5.7. In order to update other parameters in the macro-scale, the same procedure can be applied.

Setting the random fields on the micro-scale Ω m

The algorithm in Alg. 2 shows step by step to truncate a random field via the KLE method. It is applied to truncate the Young modulus E m (x, y, ω). The mean value of this random field is selected by a value of µ E m = 1 × 10 4 . To guaranteed the positive value of this modulus, the transformation is set up as follows

(5.28)

The full series of Young modulus E m (x, y, ω) is truncated via the KLE method.

With the given micro-scale mesh Fig. 5.6(b), there exists a Delaunay triangular mesh, in which all nodes in the Lattice mode mesh stand on. This hidden mesh with corresponding bilinear shape function Ψ (x) is key component to compute the Gramian matrix at each hidden triangular element and then assemble over the entire domain to gain the full Gramian matrix G as in step 2. In step 3, the Matérn The prior data in Table . 5.3 and posterior data in Table . 5.4 are assembled in Table . 5.5 to present the successful update of both moduli via the Gauss-Markov-Kalman filter method after using both tests.

Conclusions and Perspectives

In this chapter, the concluding remarks and perspectives are presented. To begin with, several highlights from the previous chapters are summarized. Then the possible future works and perspectives are drawn and discussed.