Mathematical analysis and numerical simulation of some nonlinear problems in solid mechanics
Full text
MATHEMATICAL ANALYSIS AND NUMERICAL SIMULATION OF SOME NONLINEAR PROBLEMS IN SOLID MECHANICS María Teresa Sánchez Rúa Departamento de Matemática Aplicada Universidade de Santiago de Compostela PhD Dissertation December 2009
Doña Peregrina Quintela Estévez, catedrática del Departamento de Matemática Aplicada de la Universidade de Santiago de Compostela, y Doña Patricia Barral Rodiño, profesora titular del Departamento de Matemática Aplicada de la Universidade de Santiago de Compostela, informan que la memoria titulada: MATHEMATICAL ANALYSIS AND NUMERICAL SIMULATION OF SOME NONLINEAR PROBLEMS IN SOLID MECHANICS fue realizada bajo su dirección por Doña María Teresa Sánchez Rúa, estimando que la interesada se encuentra en condiciones de optar al grado de Doctor en Ciencias Matemáticas, por lo que solicitan que sea admitida a trámite para su lectura y defensa pública. En Santiago de Compostela, a 29 de julio de 2009. Las directoras: Prof. Dr. Peregrina Quintela Estévez Prof. Dr. Patricia Barral Rodiño La doctoranda: María Teresa Sánchez Rúa
Contents Preface v I Aluminium casting processes: Efficient algorithms for their numerical simulation 1 1 Introduction 5 2 Mathematical model 9 2.1 Computational domain and notations . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Boundaryconditions................................... 10 2.3 Equilibriumequations .................................. 12 2.4 Behaviourlaw....................................... 12 2.5 Initialconditions ..................................... 13 2.6 Problem (P)....................................... 14 3 State of the art 15 3.1 Functionalframework .................................. 16 3.2 Assumptions on the data of Problem (P)........................ 17 3.3 Variational inequality. Problem (V P)......................... 18 3.4 Existenceofsolution................................... 18 3.4.1 Elasto-viscoplastic problem of Maxwell-Norton type with contact conditions 19 3.4.2 Thermo-elasto-viscoplastic problem of Maxwell-Norton type . . . . . . . . . 20 3.5 Formernumericalsolution................................ 21 3.5.1 Imposing the metallostatic pressure . . . . . . . . . . . . . . . . . . . . . . . 21 i
ii Contents 3.5.2 Discretizationinspace.............................. 23 3.5.3 Bermúdez-Moreno algorithm for the numerical solution of the contact problem 25 3.5.4 Bermúdez-Moreno algorithm for the numerical solution of the nonlinear constitutivelaw.................................... 27 3.6 Bermúdez-Moreno algorithm to solve the casting problem: FPM algorithm . . . . . 30 4 Improved numerical solution 33 4.1 Newton algorithm to approximate the contact multiplier . . . . . . . . . . . . . . . 34 4.2 Algorithms for the numerical solution of the nonlinear constitutive law . . . . . . . 43 4.2.1 Newton algorithm to approximate the viscoplastic multiplier . . . . . . . . . 43 4.2.2 Bermúdez-Moreno algorithm with variable parameters to approximate the viscoplasticmultiplier .............................. 51 4.3 Two new algorithms to solve the casting problem . . . . . . . . . . . . . . . . . . . 60 4.3.1 NMalgorithm................................... 60 4.3.2 NVPMalgorithm................................. 62 5 Numerical results 65 5.1 Test1:Agrowingbuttcurl ............................... 65 5.2 Test 2: A problem with large gradients . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.3 Numerical simulation of casting processes . . . . . . . . . . . . . . . . . . . . . . . 73 6 Conclusions 79 II Three-point bending tests: Mathematical study and numerical simulations 81 7 Introduction 85 8 Mathematical study of three-point bending tests 89 8.1 Mathematicalmodel ................................... 89 8.1.1 Boundaryconditions............................... 90 8.1.2 Equilibriumequations .............................. 91 8.1.3 Behaviourlaw................................... 91
Contents iii 8.1.4 Problem (P²)................................... 92 8.2 Weakformulation..................................... 92 8.2.1 Assumptions ................................... 93 8.2.2 Weak formulation of the problem . . . . . . . . . . . . . . . . . . . . . . . . 93 8.3 Existence of a solution of Problem (V P²)....................... 94 8.3.1 Properties of a solution of Problem (V P²)................... 94 8.4 Existence of a unique solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 9 Asymptotic analysis of three-point bending tests 99 9.1 Rescaling functions and changing variables . . . . . . . . . . . . . . . . . . . . . . . 100 9.2 Asymptotic expansion method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 9.3 Convergence of scaled displacements and stresses . . . . . . . . . . . . . . . . . . . 107 9.4 Returning to domain Ω²................................. 118 9.5 Some properties of the limit models . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 10 An improved formula for the MOR 125 10.1ThetheoreticalMOR................................... 125 10.2 Mathematical model for the real experiments . . . . . . . . . . . . . . . . . . . . . 126 10.3 Verifying hypotheses for two laboratory experiences . . . . . . . . . . . . . . . . . . 129 10.3.1 Description of data of the experiments . . . . . . . . . . . . . . . . . . . . . 129 10.3.2 Behaviour of applied forces . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 10.3.3 Verification of assumption (H3) ........................ 133 10.3.4 Existence and uniqueness of solution for experiments I and II . . . . . . . . 134 10.4 Obtaining a new formula for the MOR . . . . . . . . . . . . . . . . . . . . . . . . . 135 10.4.1 Analytical solution. A new formula for MOR . . . . . . . . . . . . . . . . . 137 11 Computing the MOR 141 11.1 Numerical results for Experiment I . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 11.2 Numerical results for Experiment II . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 12 Conclusions 149
iv Contents A Appendix: Modulus of rupture 151 III An elasticity problem for catalyst supports: Unfolding method 155 13 Introduction 159 14 Position of the elasticity problem 161 14.1 Computational domain and notations . . . . . . . . . . . . . . . . . . . . . . . . . . 161 14.2Elasticityproblem .................................... 164 15 Tools of the unfolding method 165 15.1 Extension of a displacement field . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 15.2 Elementary displacements in a beam . . . . . . . . . . . . . . . . . . . . . . . . . . 167 15.3 The unfolding operator Te e................................ 169 15.4 The unfolding operator Th ε............................... 172 16 The unfolded elasticity problem 175 16.1 Determination of applied forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175 16.2 Obtaining the limit unfolded problem and displacements . . . . . . . . . . . . . . . 176 16.2.1 Determinating ˆu3,U3and R3.......................... 184 16.2.2 Determinating ˆuα,Uα,Uαand U3....................... 188 16.2.3 Limit problems for the elementary displacement Ue.............. 195 16.3Strongconvergences ................................... 196 17 Conclusions 199 Acknowledgments 201 Resumen 203 Bibliography 213
Preface Nowadays, a complete understanding of phenomena in solids mechanics is fundamental in order to develop new methods or techniques which allow saving huge amounts of money and time in applied fields like aluminium or ceramics industries. Historically, engineers have considered linear models in order to simulate these processes. However, real materials often exhibit a non-linear behavior and, at this moment, more and more new materials are emerging which push the linear models to their limits. Therefore, the use of nonlinear models has become more relevant in the numerical simulations of this type of industries. In the present dissertation thesis, we distinguish three parts well defined and different but joined by a common topic: the mathematical study and numerical resolution of nonlinear problems in solid mechanics, arising from physical processes presented in the aluminium or ceramics industries. We consider real situations where different nonlinearities appear: behaviour laws of viscoplastic type, boundary conditions such as unilateral contact, etc. So, in the first part, we focus on the mathematical study of aluminium casting process, and we introduce some efficient algorithms to carry out its numerical simulation. In this problem, two nonlinearities appear: the nonlinear aluminium behaviour law and a contact condition with the block containing the metal. Then, the second part is devoted to the study of the well-known three-point bending experiment and the modulus of rupture of brittle materials. In this problem, again a contact condition introduces a nonlinearity in its formulation. Finally, the third part presents a mathematical analysis of an elasticity problem for catalyst supports. In what follows, we describe every part and chapter in more detail. Part 1. Aluminium casting processes: Efficient algorithms for their numerical simulation. In this first part we propose two algorithms in order to solve numerically the deformation suffered for a slab in the aluminium casting process. We consider that the casting process is a quasi-static problem for a thermo-elasto-viscoplastic solid with a unilateral contact condition with the bottom block. In a first stage, to deal with the plastic and contact effects, the numerical solution was based on the well-known Bermúdez-Moreno algorithm and the numerical procedure was introduced in [5]. During the stationary stage, in which the contact condition does not have a significant influence, this algorithm worked well; however, during the start stage, when two nonlinearities worked together, the convergence was worse, increasing considerably the cpu-time. The objective of this first part is to design new strategies to obtain a computationally efficient algorithm, such as: generalized Newton methods, an adimensionalization of the constitutive law, a partial Cholesky v
4TABLE OF CONTENTS 4.2.1 Newton algorithm to approximate the viscoplastic multiplier . . . . . . . . . 43 4.2.2 Bermúdez-Moreno algorithm with variable parameters to approximate the viscoplasticmultiplier .............................. 51 4.3 Two new algorithms to solve the casting problem . . . . . . . . . . . . . . . . . . . 60 4.3.1 NMalgorithm................................... 60 4.3.2 NVPMalgorithm................................. 62 5 Numerical results 65 5.1 Test1:Agrowingbuttcurl ............................... 65 5.2 Test 2: A problem with large gradients . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.3 Numerical simulation of casting processes . . . . . . . . . . . . . . . . . . . . . . . 73 6 Conclusions 79
Chapter 1 Introduction The aim of Part I of this manuscript is to obtain efficient algorithms to solve nonlinearities arising from behaviour laws of viscoplastic type, or from boundary conditions such as unilateral contact. A real situation, where those two nonlinearities appear and which is the motivation of this research, is the casting process in aluminium production. During casting processes, the liquid aluminium is poured onto a water-cooled mold which is called bottom block; when the aluminium begins to solidify the bottom block starts to descend leaving room for more liquid metal (see Figure 1.1). A detailed description of the complete process can be found in [33, 41, 75]. shell solidified 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 bottom block liquid aluminium watercooled mould water cooling jets Figure 1.1: Scheme of a casting process and a photograph of the mold. The casting process can be divided into two stages. During the start stage, the temperature field, the solidification front, and the slab shape change with time. The large thermal stresses due to the cooling jets cause the butt of the slab to bend and lose contact with the bottom block. This deformation is known as butt curl (see Figure 1.2). When the length of the slab is 1m, approximately, the temperature field has attained a steady state. In this stationary stage the solidified shell contracts inwards and the resulting cross-section has a “bone” shape. 5
6Chapter 1. Introduction x1 x2 x3 Butt curl Bottom block Figure 1.2: Butt curl deformation The deformations of butt curl and of cross-sectional could be minimized by optimizing casting parameters but first it is necessary to be able to predict the deformations with fixed casting conditions. To predict casting behaviour, the main difficulties are due to coupling effects, to the presence of nonlinear terms, to the existence of free boundaries and to the computational domain which varies with time (see [6]). Specifically, in order to simulate the butt curl deformation we must solve a problem of frictionless unilateral contact between a thermo-elasto-viscoplastic body of Maxwell-Norton type and a rigid foundation. This manuscript is the continuation of the research initiated in [5] and [9]-[13], where several mathematical aspects of this subject have been studied: •The detailed formulation as a variational inequality. •The establishment of sufficient conditions for the existence of a solution. •The numerical implementation of the metallostatic pressure on the liquidus-solidus interphase. •The numerical approximation of the variational inequality and the description of a numerical scheme to solve it. In [5], an approximate solution of this problem was obtained by using an implicit Euler scheme in time and a finite element method in space. To deal with the nonlinearities, the numerical solution was based on the Bermúdez-Moreno algorithm with the participation of two multipliers: the viscoplastic multiplier (to take into account the nonlinearity from the behaviour law) and the contact multiplier (to avoid the nonlinearity due to the contact condition). These two multipliers were approximated by using fixed point algorithms (see [9, 10, 16]). During the stationary stage, in which the contact condition does not have a significant influence, this algorithm worked well; however, during the start stage, when two nonlinearities worked together, the convergence was worse, increasing considerably the cpu-time. Moreover, this algorithm presents a strong dependence on its parameters, which could not be determined a priori. These facts made difficult the massive usage of this algorithm to help in the understanding of the different effects in casting processes.
7 The fundamental objective of this research work is to develop efficient numerical methods that incorporate adaptive procedures to accelerate the treatment of the nonlinearities and hence, be able to simulate numerically the butt curl deformation in a reasonable cpu-time. Several approaches are employed to achieve this objective: •Contact condition. Approximating the contact multiplier using generalized Newton methods together with a penalization technique to conserve the matrix’s symmetry (see [7, 14, 46, 54, 72]). The resultant algorithm for the treatment of the contact condition is fast, accurate and its convergence is independent of the algorithm parameters. Nevertheless, the stiffness matrix needs to be recalculated at each iteration. So, taking into account that usually the nonlinear boundary condition only involves a small part of the boundary, we propose to use a factorization of the stiffness matrix adapted to the problem’s geometry. •Viscoplastic law. Two options are considered: –Approximating the viscoplastic multiplier with standard Newton techniques without modification of the stiffness matrix at each iteration (see [7, 14]). –Approximating the viscoplastic multiplier with a generalization of the Bermúdez-Moreno algorithm with variable parameters automatically calculated (see [43]). When the nonlinearity coming from the viscoplastic law is solved by using Newton techniques, the strong thermal gradients of stresses force to employ some numerical techniques to obtain a good convergence: an adimensionalization technique, an Armijo rule and an optimization of the time step. We present academic problems in order to show the good behaviour of this methodology but, in the real casting process, the convergence is not always achieved. Due to this, we propose a second algorithm which is a generalization of the Bermúdez-Moreno algorithm with an automatic choice of its parameters. The outline of this Part I is as follows: In Chapter 2, we introduce the mathematical model to simulate the butt curl deformation during the casting process, giving a description of the boundary conditions, equilibrium equations and the nonlinear behaviour law associated to aluminium casting. The latter is a thermo-elasto- viscoplastic law, where the viscoplastic part is given by a nonlinear Maxwell-Norton law, and the thermal part is a generalization of Arrhenius law. In Chapter 3, we make a review of the most important results for the mathematical analysis and the numerical solution of the butt curl deformation. We propose a weak formulation for the problem in terms of suitable functional spaces, following the works of Geymotat and Suquet (see [45]) and Barral and Quintela (see [13]). Moreover, some existence results proved by Barral and Quintela in [13] and Barral et al in [8] are presented. Finally, we recall the numerical procedure proposed in [5] to approximate the solution of the problem and to simulate the butt curl deformation. This procedure is based in the Bermúdez-Moreno algorithm. With the aim of improving the poor convergence of the former algorithm in butt-curl simulation, Chapter 4 describes a new numerical procedure to approximate the solution of the problem. The contact condition is solved by means of the Bermúdez-Moreno algorithm (see [16]), in which the
8Chapter 1. Introduction contact multiplier is a fixed point of a nonlinear equation. To approximate this multiplier we propose a generalized Newton method for nonsmooth equations and a penalization technique. In order to deal with the nonlinear constitutive law we use again the Bermúdez-Moreno algorithm and we propose two methods to approximate the corresponding multiplier: a Newton algorithm together with some numerical strategies designed to improve its efficiency and a Bermúdez-Moreno algorithm with variable parameters. In Chapter 5, we present some academic examples with known analytical solutions to validate the proposed algorithms. These examples have been designed to reproduce a behaviour analogous to semicontinuous casting behaviour: large gradients of stresses with respect to time and a gap between the computational domain and the rigid foundation which grows with time. Numerical results are given to compare the efficiency of these methodologies with the former Bermúdez- Moreno algorithm used in [9, 10]. Finally, we present the numerical results for the simulation of the butt curl deformation, showing that the better results are obtained when considering the Newton algorithm for the contact, and the Bermúdez-Moreno one with variable parameters for the viscoplasticity.
Chapter 2 Mathematical model In this chapter we introduce a mathematical model to simulate the deformation suffered by an aluminium slab during the solidification process. For this purpose, we consider a quasi-static problem for a thermo-elasto-viscoplastic solid with a unilateral contact condition with the bottom block. In this manuscript, coupling between thermal and mechanical models is not considered. The temperature field is computed a priori by using the code presented in [63] and used as a data in the mechanical simulation. A detailed description of the thermal model can be found in [6, 17, 63]. Several authors have studied the mechanical deformations in casting processes. In 1982, Kristiansson and Zetterlund (see [55]) studied the deformation in steel continuous casting by means of a two-dimensional viscoelastic model, formulated on a horizontal section. Later, in 1992, Mariaux et al. (see [59]) studied the butt curl simulation in aluminium casting and, in 1995, Drezet et al. (see [34]) simulated the cross-sectional deformation; both works considered a viscoelastic model. In 2005, Sengupta et al. (see [77]) presented a coupled thermal-stress model to simulate the start-up phase of the casting process. Nevertheless, all these works present numerical results, but they do not specify the used models or algorithms. Contrarily, in 2001 Barral (see [5]) presented a mathematical analysis and a numerical study of a quasistatic evolution problem to model the thermomechanical behaviour of an aluminium casting, identifying all the governing equations as well as the suitable functional spaces and numerical algorithms. In this manuscript we are going to use this mathematical model in order to simulate the butt-curl deformation arising from casting processes. The outline of this chapter is as follows: In Section 1.1, we introduce some notations we will use hereafter in Part I, and the computational domain of the mechanical simulation. In Sections 1.2 and 1.3, we present the boundary conditions and equilibrium equations, respectively. Section 1.4 is devoted to introduce the nonlinear behaviour law; a thermo-elasto-viscoplastic law is considered, where the viscoplastic part is given by a nonlinear Maxwell-Norton law and the thermal part is a generalization of Arrhenius law. Finally, in Section 1.5, the initial conditions are given and in Section 1.6 we present the complete mathematical model. 9
10 Chapter 2. Mathematical model 2.1 Computational domain and notations Let Ox1x2x3be a fixed system of rectangular Cartesian axes associated to the slab, whose origin coincides with the center of the slab base. Hereafter in this manuscript, Latin subscripts are understood to range over the integers {1,2,3}, and Greek subscripts over the integers {1,2}. Einstein notation of summation over repeated subscripts is implied. Given two vectors u,v∈R3, with components (ui)i=1,2,3,(vi)i=1,2,3, respectively, u·vdenotes their scalar product. From now on, S3denotes the space of symmetric second order tensors endowed with the usual scalar product τ:ξ=τijξij,τ,ξ∈ S3. Given τa tensor function in S3, we denote its divergence div(τ)i=∂jτij, and its trace tr(τ) = 3 X k=1 τkk. Moreover, we introduce the following norm in S3: |τ|=ptr(τ:τ). From here on, we will use a superposed dot to denote differentiation with respect to time. Let [0, tf]be the time interval to carry out the mechanical simulation. Due to casting symmetry, Ω(t)represents a quarter of the slab at the instant t∈[0, tf]. The temperature field T(x, t)at each point x∈Ω(t)is previously computed by using the mathematical model developed in [63]. The mechanical domain at each time instant tcorresponds to the solidified part of the slab, denoted by Ωs(t), which is obtained from the solution of the thermal problem: Ωs(t) = {x∈Ω(t); T(x, t)< Tl}, where Tlis the liquidus temperature (see Figure 2.1). We assume that at the initial instant t= 0 there exists a solidified part of the slab, that is, Ωs(0) 6=∅. The mechanical problem consists of determining the displacement field u(x, t)and the stress tensor field σ(x, t)at each point x∈Ωs(t)and at each instant t∈(0, tf]. 2.2 Boundary conditions Let Γ(t)be the boundary of Ωs(t)and nthe unit outward normal vector to Ωs(t). Boundary Γ(t) is split into five disjoint parts Γ(t) = ¯ Γsl(t)∪¯ ΓC∪¯ Γsym(t)∪¯ Γ1 n(t)∪¯ Γ2 n(t), where:
2.2. Boundary conditions 11 x3 x2 GC G1 n(t) G2 n(t) Gsl(t) Gc Gsym(t) Ws(t) Bottomblock AluminiumIngot Liquidusinterface Mold Waterjet frommold Buttcurl Figure 2.1: Computational domain Ωs(t). •The upper boundary Γsl(t)is defined by the isotherm corresponding to the liquidus temperature Tl. This isotherm is the free boundary of the thermal problem and, therefore, it varies with time. Γsl(t)is subjected to the metallostatic pressure exerted by the overlying liquid metal σn=prn, pr(x, t) = ρ(T(x, t))g(x3−h(t)) on Γsl(t),0< t ≤tf, where ρis the temperature-dependent mass density of the aluminium, gthe gravitational acceleration and h(t)the length of the slab at the instant t. •ΓCis the part of the slab susceptible to be in contact with the bottom block. This is assumed to be rigid, so aluminium cannot penetrate it, that is, the normal component of displacements, un=u·n, cannot be positive: un≤0on ΓC×(0, tf]. Moreover, due to the action-reaction principle, where an effective contact exists, un= 0, the bottom block exerts a pressure in the normal direction and upwards. So, the normal stress, σn=σijnjni, verifies σn≤0. Contrarily, where a contact does not exist, un<0, the distance between the slab and the bottom block is not null and, consequently, the motion takes place freely, σn= 0. Therefore, σnun= 0 on ΓC×(0, tf]. Finally, taking into account that water is flowing down the slab sides, we assume that friction does not exist in the contact zone, that is, the tangential stress, σt=σn−σnn, is null. Summing up, to reproduce the butt curl deformation, we consider a frictionless contact condition with the bottom block, known as Signorini unilateral frictionless contact condition σt=0, σn≤0, un≤0, σnun= 0 on ΓC×(0, tf]. A detailed description of this contact condition can be found in [54].
12 Chapter 2. Mathematical model •Boundary Γ1 n(t)denotes the part of the lateral outer faces which has already solidified; so, it is free of forces σn=0on Γ1 n(t),0< t ≤tf. •The boundary Γ2 n(t)corresponds to the mushy region, which is confined by the mold; so, σt=0, un= 0 on Γ2 n(t),0< t ≤tf. •Finally, on Γsym(t)we assume usual symmetry conditions: σt=0, un= 0 on Γsym(t),0< t ≤tf. 2.3 Equilibrium equations Under the small strains assumption and in the quasistatic case, the slab behaviour is given by the equilibrium equations −div(σ) = fin Ωs(t), where frepresents the volume forces due to the gravity, that is f(x, t) = (0,0,−ρ(T)g). 2.4 Behaviour law Following the works of Drezet et al. [34], Mariaux et al. [59] and Barral et al. [9] where this type of deformations were studied, we assume that the aluminium is a nonlinear Maxwell-Norton thermoelasto-viscoplastic solid; so, the strain rate tensor is the superposition of elastic, viscoplastic and thermal components: ˙ ε(u) = ˙ εe(u) + ˙ εvp(u) + ˙ εth(u). Here the strain tensor is related to the displacement field by the usual formula εij(u) = 1 2(∂jui+∂iuj). The elastic deformations εeare related to the stress tensor σthrough Hooke’s law with material parameters dependent on temperature εe(u) = Λs(T)σ=1 + ν(T) E(T)σ−ν(T) E(T)tr(σ)I, where E,νdenote the Young’s modulus and Poisson’s ratio, respectively, and Iis the identity tensor. Reciprocally, the stress tensor σcan be written in terms of the strain tensor εe, σ= (Λs(T))−1εe(u) = λ(T)tr(εe(u))I+ 2µ(T)εe(u),
2.5. Initial conditions 13 being λand µthe Lamé coefficients related with Eand νwith the usual expressions: λ(T) = E(T)ν(T) (1 + ν(T))(1 −2ν(T)), µ(T) = E(T) 2(1 + ν(T)). A detailed description of the elastic behaviour law can be found in [28, 52]. The thermal expansion is related to the temperature by a generalized Arrhenius law ˙ εth(u) = αs(T)˙ TI, where αsis the coefficient of thermal expansion, including volume changes due to possible phase transformations. A detailed description is found in [19, 59]. The viscoplastic part of the strain rate field, ˙ εvp,is governed by Norton-Hoff’s law (see [42]). To describe this law, let us introduce the dissipation potential Φqdefined on S3by Φq(τ) = 1 qθ(T)|τ|q, q≥2is a fixed real number and θdenotes a positive material parameter depending on temperature, θ(T) = θ0e −G R(T+273) , θ0being a material parameter, Gthe activation energy and R the universal gas constant. Norton-Hoff’s law may be written as (see [42]) ˙ εvp(u) = ∇Φq(σD) = θ(T)|σD|q−2σD, where ∇Φqdenotes the gradient of Φqrelative to the inner product on S3and τDdenotes the deviatoric part of τgiven by τD=τ−1 3tr(τ)I. Summing up, the constitutive law can be formulated as ε(˙ u)(t) = ˙ z { (Λs(T)σ)(t)+(∇Φq(σD))(t) + ³αs(T)˙ T´(t)Iin Ωs(t),0< t ≤tf. 2.5 Initial conditions In order to complete the mathematical model, we consider the following initial conditions: u(x,0) = u0(x),σ(x,0) = σ0(x)in Ωs(0), where u0,σ0verify suitable compatibility conditions detailed in Chapter 4.
20 Chapter 3. State of the art 3.4.2 Thermo-elasto-viscoplastic problem of Maxwell-Norton type In [8] it was proved the existence and uniqueness of solution for a quasistatic thermo-elasto- viscoplastic problem of Maxwell-Norton type with coefficients depending on temperature when it is verified that ˙ T≥0. In this case, contact conditions were not considered. Let Ωbe a bounded, open, connected region of R3with a Lipschitz boundary. Moreover, its boundary is partitioned into two non-empty, open and disjointed parts ΓDand ΓNsatisfying Γ = ∂Ω = ¯ ΓD∪¯ ΓN, with meas(ΓD)>0. Problem (ˆ P): Given the temperature field T, find the displacement field uand the stress tensor σverifying −div(σ) = f,in (0, tf]×Ω, ε(˙ u) = ˙ z { (Λs(T)σ) + ∇Φq(σD) + αs(T)˙ TI,in (0, tf]×Ω, u=ud,on (0, tf]×ΓD, σn=h,on (0, tf]×ΓN, u(0) = u0,σ(0) = σ0,in Ω, Assumptions (d H1)The temperature field T∈W2,∞(0, tf;L∞(Ω)). Moreover it is positive and ˙ T≥0a.e. in Ω,∀t∈[0, tf]. (d H2)The applied forces satisfy f∈W2,∞(0, tf; [Lq(Ω)]3) h∈W2,∞(0, tf; [W−1 q,q(Γ)]3∩[Lq(ΓN)]3). (d H3)The elasticity tensor Λsis symmetric and the following hold true: •Λs∈[W2,∞(R)]81, therefore, there exists β1>0such that |Λs(r)τ| ≤ β1|τ|and there exists β2>0such that |Λ0 s(r)τ| ≤ β2|τ|,∀τ∈ S3,∀r∈R. •There exists β > 0such that Λs(r)τ:τ≥β|τ|2,∀τ∈ S3,∀r∈R. •Λ0 s(r)τ:τ≥0,∀τ∈ S3,∀r∈R. (d H4)The coefficient of thermal expansion αs∈W1,∞(R). (d H5)The coefficient of the viscoplastic law θ∈L∞(R)and θ(r)≥θc>0,∀r∈R.
3.5. Former numerical solution 21 (d H6)The initial stress field σ0∈Hqand the initial displacement u0∈Vpverify −div(σ0) = f(0) in Ω, σ0n=h(0) on ΓN, u0=ud(0) on ΓD, ε(u0) = Λs(T0)σ0in Ω. (d H7)The boundary Dirichlet condition satisfies ud∈W2,∞(0, tf; [W1−(1/p),p(Γ)]3). (d H8)The exponent qverifies q≥2, and its conjugate pis such that 1< p ≤2. Theorem 3.4.2. Under the above assumptions (d H1)−(d H8)there exists a unique solution (u,σ)∈ W1,2(0, tf;Vp)סW1,2(0, tf;X2)∩L∞(0, tf;Xq)¢of Problem (ˆ P). 3.5 Former numerical solution In this section we are going to review the numerical resolution of Problem (V P)defined by (3.1)- (3.3), given in [5]. The main difficulties to overcome in the numerical solution of this problem were the following: •The solidified part of the slab, which is the computational domain of the mechanical simulation, changes with time. Furthermore, on the upper boundary, which is the isotherm of the liquidus temperature and, therefore, the free boundary of the thermal problem, we must impose the metallostatic pressure due to the liquid metal. To deal with this difficulty, a fictitious domain method is considered. •To model the butt curl, we must solve a contact condition between the slab and the bottom block. To deal with this nonlinear condition, the Bermúdez-Moreno algorithm, based on maximal monotone operator techniques, was considered. •The aluminium behaviour is nonlinear and depends strongly on the temperature field. In order to avoid the nonlinearity of the constitutive law, again the Bermúdez-Moreno method was used. •The computational domain grows with time and the zone with steep gradients changes. To overcome this difficulty, the meshes were constructed in a way adapted to the problem’s geometry. 3.5.1 Imposing the metallostatic pressure In the mechanical simulation of an aluminium casting, the computational domain is the solidified part of the slab. So, its upper boundary coincides with the isotherm corresponding to the liquidus temperature, denoted by Γsl(t), which is the free boundary of the thermal problem. On this boundary, it is necessary to impose the metallostatic pressure exerted by the weight of the liquid
22 Chapter 3. State of the art metal. In order to implement numerically this condition, we are going to use the fictitious domain method introduced in [11]. It consists of solving the problem over the entire slab, assuming that the liquid metal is a very elastic material under the action of the gravity forces; so, the liquid metal does not offer resistance to solid deformations. This methodology has the advantage that remeshing is not needed at each time instant to adjust the free boundary, and so the numerical implementation is easier. Bottomblock AluminiumIngot Liquidusinterface Mold Waterjet frommold Buttcurl x3 x2 Gsl(t) Wl(t) Gl(t) Gsup(t) Gsym(t) Figure 3.1: Fictitious domain Ωl(t). Let Ωl(t)be the liquid zone of the slab at the instant t, Ωl(t) = {x∈Ω(t); T(x, t)> Tl}. The boundary of Ωl(t)is split into four disjoint parts: ∂Ωl(t) = ¯ Γsup(t)∪¯ Γsym(t)∪¯ Γl(t)∪¯ Γsl(t), where Γsup(t)denotes the upper boundary, Γsym(t)the symmetry boundary and Γl(t)the outer lateral boundary (see Figure 3.1). In the liquid domain, Ωl(t), the applied forces are the following: •To impose the weight of the liquid aluminium, we consider volume forces due to gravity, with constant density ρl=ρ(Tl). •Since the liquid must not offer resistance to solid deformations we consider that the surface forces on the upper boundary, Γsup(t), are null. •On the outer lateral boundary, Γl(t), the aluminium is confined by the mold. In order to impose that the liquid metal is a very elastic material, we consider that the behaviour law in the liquid domain is given by the Hooke’s law depending on a small parameter ², σ²(u²) = Λ−1 lε(u²) = λ² ldiv(u²)I+ 2µ² lε(u²), where the Lamé coefficients of the liquid domain, λ² l, µ² l, change with different length scale: λ² l=²βλ, µ² l=²αµ,
3.5. Former numerical solution 23 with λ, µ, α and βreal positive numbers independent of ². Assuming that α > β, this has the effect of applying the metallostatic pressure to the solidification front, as it was proved in [12] by using asymptotic expansion techniques. In order to simplify the notation, from now on, we omit the index ²in the behaviour law. With this technique, the variational inequality (3.1) can be extended to the complete slab and the integral over the thermal free boundary Γsl(t)disappears. Therefore, from now on, we will consider this new extended weak formulation over the entire slab Ω(t) = int(Ωs(t)∪Ωl(t)∪Γsl(t)): Problem (EV P) Find u∈W1,∞(0, tf;Up ad(t)) and σ∈W1,∞(0, tf;Hq(t)) verifying a.e. in (0, tf]: ZΩ(t) σ(t) : ε(v−u(t)) dx ≥ZΩ(t) f(t)·(v−u(t)) dx, ∀v∈Up ad(t),(3.4) ε(˙ u)(t) = ˙ z { (Λs(T)σ)(t)+(∇Φq(σD))(t) + ³αs(T)˙ T´(t)Iin Ωs(t), ˙ z { (Λlσ)(t)in Ωl(t), (3.5) u(0) = u0,σ(0) = σ0,in Ω(0).(3.6) Remark 3.5.1. Notice that the functional spaces are extended to the complete slab in the natural way. The initial data, u0and σ0, should be extended to all Ω(0), ensuring that the following compatibility conditions are satisfied: −div(σ0) = f(0) in Ω(0), σ0n=0on Γsup(0) ∪Γ1 n(0), σ0t=0, u0n= 0 on Γsym(0) ∪Γ2 n(0) ∪Γl(0), σ0t=0, σ0n≤0, u0n≤0, σ0nu0n= 0 on ΓC, ε(u0) = ½Λs(T0)σ0in Ωs(0), Λlσ0in Ωl(0). Remark 3.5.2. In casting processes the temperature at t= 0,T0, corresponds with the liquidus temperature Tl. 3.5.2 Discretization in space Finite element approximation to simulate the model is considered in the usual way. Displacements are discretized in space using the Lagrange tetrahedral finite element method of degree one. Stresses are assumed constant within each element. Let Th(t)be a tetrahedral mesh on the computational domain ¯ Ω(t)compatible with the boundary partition. As usual, hdenotes the maximum diameter of Th(t)and Σh(t) = {bi}1≤i≤Nh(t),
24 Chapter 3. State of the art the nodes set of the tetrahedral mesh. Moreover, Shand Zh(t)are the triangulations induced on the contact boundary ΓCand the boundaries Γ2 n(t)∪Γl(t), respectively. We construct a family of finite dimensional subspaces Vh(t)of Vp(t)by approximating the test functions vby piecewise polynomials of degree one over the tetrahedral mesh Th(t)on the computational domain ¯ Ω(t): Vh(t) = {vh∈[C0(Ω(t))]3;vh|K∈[P1(K)]3,∀K∈ Th(t)}, Uh(t) = {vh∈Vh(t); (vh)α= 0 on [xα= 0],vh|C·n|C= 0,∀C∈ Zh(t)}, where Pi(K)denotes the space of polynomials of degree idefined on Kand vh|C(resp. n|C) denotes the value of vh(resp. n) on the barycenter of C. Therefore, the discrete space of admissible displacements is Uadh(t) = {vh∈Uh(t); vh|C·n|C≤0,∀C∈ Sh}. Stresses are assumed constant within each element. We denote by Xh(t)the discretized space of stresses Xh(t) = {ξh; (ξh)ij = (ξh)ji,ξh|K∈[P0(K)]9,∀K∈ Th(t)}. In this space, we consider the norm induced by L2. Functions vh∈Vh(t)are characterized for their values in the vertices of the elements, vh(x) = Nh(t) X i=1 vh(bi)φi(x), where {φi}i=1,...,Nh(t)are the corresponding shape functions. Elements ξh∈Xh(t)are characterized for their value in the barycenter of each element, ξh(x) = X K∈Th(t) ξh|KXK(x), where XKdenotes the characteristic function of the element K. Moreover, Thdenotes the discretized temperature field, which is considered constant by element. Finally, at each time instant twe consider the operator εh:Vh(t)−→ Xh(t) vh−→ εh(vh) = X K∈Th(t) ε(vh|K)XK, and we denote by ε∗ hits adjoint operator. With this notation, the discretized variational formulation corresponding to Problem (EV P ) is the following: Problem (DV P):
3.5. Former numerical solution 25 Find uh∈W1,∞(0, tf;Uadh(t)) and σ∈W1,∞(0, tf;Xh(t)) verifying a.e. in (0, tf): ZΩ(t) σh(t) : εh(vh−uh(t)) dx ≥ZΩ(t) f(t)·(vh−uh(t)) dx, ∀vh∈Uadh(t),(3.7) εh(˙ uh)(t) = ˙ z { (Λs(Th)σh)(t)+(∇Φq(σD h))(t) + ³αs(Th)˙ Th´(t)Iin Ωs(t), ˙ z { (Λlσh)(t)in Ωl(t), (3.8) uh(0) = u0,σh(0) = σ0,in Ω(0).(3.9) To simplify the notation, from now on we omit the index hto denote the elements of the discretized spaces. 3.5.3 Bermúdez-Moreno algorithm for the numerical solution of the contact problem Both the contact condition and the viscoplastic behaviour, generate nonlinearities. A simple way to transform these nonlinearities is to use the well-known Bermúdez-Moreno lemma (see [16]). This lemma is the main idea of the algorithm developed in [9] and that we introduce in the following. In order to introduce an algorithm to deal with the nonlinearity due to the contact condition, we focus on solving a simplified problem, corresponding to Problem (DV P)suppressing the viscoplastic part in the constitutive law. Therefore, the discretized variational problem remains: ZΩ(t) σ(t) : ε(v−u(t)) dx ≥ZΩ(t) f(t)·(v−u(t)) dx, ∀v∈Uadh(t),(3.10) ε(u)(t) = (Λs(T)σ)(t) + µZT Tl αs(r)dr¶Iin Ωs(t), (Λlσ)(t)in Ωl(t), (3.11) u(0) = u0,σ(0) = σ0,in Ω(0).(3.12) Remark 3.5.3. If we remove the viscoplastic behaviour from the constitutive equation (3.8), we have ε(˙ u) = ˙ z { (Λs(T)σ) + ³αs(T)˙ T´Iin Ωs(t), ˙ z { (Λlσ)in Ωl(t). Therefore, if we integrate this equation in time, assuming that u0,σ0verify the elastic law at t= 0, we obtain the constitutive equations (3.11). In [85], it can be found an overview of different methods to solve contact problems. For numerical reasons, following [16, 27], to deal with the inequality in (3.10), we introduce a Lagrange
26 Chapter 3. State of the art multiplier. For that purpose we consider the following notations: Eh={ph∈L∞(ΓC); ph|C∈P0(C),∀C∈ Sh}, Qh={ph∈Eh;ph|C≤0,∀C∈ Sh}. In [9], it is proved that the discretized problem in space (3.10)-(3.12) can be formulated as a weak equality in the following manner: Theorem 3.5.4. Every solution of the discretized problem (3.10)-(3.12) is a solution of the problem: Find u∈W1,∞(0, tf;Uadh(t)),σ∈W1,∞(0, tf;Xh(t)) and p∈W1,∞(0, tf;Eh)such that: ZΩ(t) σ(t) : ε(v)dx +ZΓC γcunvndΓ = ZΩ(t) f(t)·vdx −ZΓC p(t)vndΓ,∀v∈Uh(t),(3.13) ε(u)(t) = (Λs(T)σ)(t) + µZT Tl αs(r)dr¶Iin Ωs(t), (Λlσ)(t)in Ωl(t), (3.14) p(t) = Gγc λc(un(t) + λcp(t)) on ΓC,(3.15) for all t∈(0, tf], and the initial conditions given in (3.12). Here λc,γcdenote real numbers such that λc>0,0≤λcγc<1and Gγc λcis the Yosida approximation of Gγc=G−γcI,G=∂IQhbeing the indicator function of the closed convex Qh. Remark 3.5.5. The Lagrange multiplier p(t)takes into account the nonlinearity due to the contact condition and it will be called contact multiplier. In order to introduce an iterative algorithm to deal with the contact condition, we discretize the problem (3.13)-(3.15) in time using the following notation: The time interval of interest is divided into Nsteps t0= 0, tj+1 =tj+ ∆t, j = 0, ..., N −1,∆t=tf N. From now on we denote by gjan approximation of a given function g(t)at time tj. It is easy to prove that the Yosida approximation of Gγcis given by the expression (see [16]) Gγc λc(ϕ) = 1 λcµϕ−ΠQhµ1 1−λcγc ϕ¶¶, ϕ ∈Eh, where ΠQhis the orthogonal projection over Qh ΠQh(ph) = ½phif ph≤0, 0if ph>0. Then, from Theorem 3.5.4, we obtain the following fixed point algorithm to solve Problem (3.10)- (3.12) introduced in [5]:
3.5. Former numerical solution 27 •Let (u0,σ0, p0)be given and γc,λc>0such that λcγc<1. •Then, for j≥0,(uj,σj, pj)known at time tj, we determine (uj+1,σj+1, pj+1)at time tj+1, an approximated weak solution of Problem (3.10)-(3.12). To do so, we propose the following iterative algorithm: i) Given (uj+1 k−1,σj+1 k−1, pj+1 k−1),we compute (uj+1 k,σj+1 k)∈Uadh(tj+1)×Xh(tj+1)by solving the variational equality ZΩj+1 σj+1 k:ε(v)dx +ZΓC γc³uj+1 k´nvndΓ = ZΩj+1 fj+1 ·vdx −ZΓC pj+1 k−1vndΓ,∀v∈Uh(tj+1), where σj+1 k= ¡Λs(Tj+1)¢−1³ε(uj+1 k)−αs(Tj+1)(Tj+1 −Tl)I´in Ωj+1 s, Λ−1 lε(uj+1 k)in Ωj+1 l. ii) The updated contact multiplier is given by relation pj+1 k=1 λc ³uj+1 k´n+λcpj+1 k−1−ΠQh ³uj+1 k´n+λcpj+1 k−1 1−λcγc . Remark 3.5.6. It can be proved that the weight of the entire slab is the total reaction force exerted by the bottom block over the slab. This property is used to initialize the contact multiplier (see [9]). 3.5.4 Bermúdez-Moreno algorithm for the numerical solution of the nonlinear constitutive law In order to avoid the nonlinearity due to the constitutive law (3.8), the numerical resolution is based on the Bermúdez-Moreno algorithm which involves a viscoplastic multiplier, which is a fixed point of a nonlinear equation. In a first stage, as in the contact case, the problem was solved in [10] by using fixed point techniques, which let to design an algorithm that we introduce in the following. As we have done in previous subsection, we are going to consider a simplified problem, corresponding to Problem (DV P ), replacing the contact boundary condition by a null Dirichlet condition. Therefore, the discretized variational problem remains: ZΩ(t) σ(t) : ε(v)dx =ZΩ(t) f(t)·vdx, ∀v∈Uadh(t),(3.16) ε(˙ u)(t) = ˙ z { (Λs(T)σ)(t)+(∇Φq(σD))(t) + ³αs(T)˙ T´(t)Iin Ωs(t), ˙ z { (Λlσ)(t)in Ωl(t), (3.17)
28 Chapter 3. State of the art u(0) = u0,σ(0) = σ0,in Ω(0),(3.18) where Uadh(t) = {vh∈Uh(t); vh=0on ΓC}. As in previous section, the thermo-elasto-viscoplastic law (3.17) is discretized in time by using an implicit Euler scheme: ε(uj+1)−ε(uj) = ¡Λs(Tj+1)σj+1 −Λs(Tj)σj¢+αs(Tj+1)(Tj+1 −Tj)I+ ∆t∇Φq³¡σj+1¢D´in Ωj+1 s, ε(uj+1)−ε(uj) =Λlσj+1 −Λlσjin Ωj+1 l. Once we have discretized the constitutive law, our objective is to obtain an explicit expression for σj+1 and to substitute it in the mixed variational formulation (3.16). However, the nonlinearity of the viscoplastic function makes difficult to compute the stress tensor in Ωs(t). In order to solve this difficulty, we are going to use again maximal monotone operator techniques. The proofs of the following results can be found in [5]. Proposition 3.5.7. At each time step tj+1, j = 0, . . . , N −1, the stress tensor in Ωj+1 sis given by the relation σj+1 =V(Tj+1)µε(uj+1)−∆tqj+1 +∆tγpsj+1 3tr(ε(uj+1))I+Fj¶,(3.19) where qj+1 is the viscoplastic multiplier at the time step tj+1 qj+1 = (∇Φq)γp³¡σj+1¢D´= (∇Φq−γpI)³¡σj+1¢D´,(3.20) and Fjis the history of the solidified metal up to time tjtogether with the expansional effects of the temperature changes at the time interval [tj, tj+1] Fj=−ε(uj)+Λs(Tj)σj−αs(Tj+1)(Tj+1 −Tj)¡1+∆tγpsj+1¢I+ ∆tγpsj+1 3µ1 sjtr(σj)−tr(ε(uj))¶I,(3.21) where γp≥0and sj=E(Tj) 1−2ν(Tj). In (3.19), the automorphism V(T)is defined by V(T)τ=e λ(T)tr(τ)I+ 2eµ(T)τ, being e λ(T) = E(T)ν(T) (L(T)−2ν(T))(L(T) + ν(T)),eµ(T) = E(T) 2(L(T) + ν(T)), and L(T) = 1 + γpE(T)∆t.
3.5. Former numerical solution 29 Lemma 3.5.8. At each time step tj+1,j= 0, . . . , N −1, the viscoplastic multiplier qj+1 is a fixed point of the equation qj+1 = (∇Φq)γp λp³¡σj+1¢D+λpqj+1´,(3.22) for λpa real positive number such that λpγp<1. In this equation, (∇Φq)γp λpdenotes the Moreau- Yosida approximation of (∇Φq)γpgiven by the expression (∇Φq)γp λp(ζ) = ζ λpµ1−1 ηj+1(1 −λpγp)¶,ζ∈ S3, where ηj+1 =η(ζ)is the unique root of the equation ηq−1−ηq−2−λpθ(Tj+1) (1 −λpγp)q−1|ζ|q−2= 0,(3.23) in the interval [1,+∞). Here we summarize the Bermúdez-Moreno algorithm to solve Problem (3.16)-(3.18) introduced in [5] and based on the results of Lemma 3.5.8: •Let (u0,σ0)be given and γp,λp>0such that λpγp<1; we compute q0= (∇Φq)γp(σD 0). •Then, for j≥0,(uj,σj,qj)known at time tj, we determine (uj+1,σj+1,qj+1)at time tj+1, an approximated weak solution of Problem (3.16)-(3.18). To do so, we propose the following iterative algorithm: i) Given (uj+1 k−1,σj+1 k−1,qj+1 k−1),we compute (uj+1 k,σj+1 k)∈Uadh(tj+1)×Xh(tj+1)by solving the variational equality, obtained by substituting expression (3.19) in equality (3.16), ZΩj+1 sV(Tj+1)·ε(uj+1 k) + ∆tγpsj+1 3tr(ε(uj+1 k))I¸:ε(v)dx +ZΩj+1 l Λ−1 lε(uj+1 k) : ε(v)dx = ZΩj+1 shV(Tj+1)∆tqj+1 k−1i:ε(v)dx −ZΩj+1 s£V(Tj+1)Fj¤:ε(v)dx+ ZΩj+1 fj+1 ·vdx ∀v∈Uadh(tj+1), where Fjis given by expression (3.21), and σj+1 k= V(Tj+1)µε(uj+1 k)−∆tqj+1 k−1+∆tγpsj+1 3tr(ε(uj+1 k))I+Fj¶in Ωj+1 s, Λ−1 lε(uj+1 k)in Ωj+1 l. ii) The updated viscoplastic multiplier is given by relation qj+1 k=³σj+1 k´D+λpqj+1 k−1 λpÃ1−1 ηj+1 k(1 −λpγp)!,
36 Chapter 4. Improved numerical solution and therefore, its Clarke generalized Jacobian is the interval ∂ˆ Gλc(0) = ·0,1 λc¸, which coincides with its subdifferential. Lemma 4.1.8. At each point s0∈R, the following approximation for ˆ Gλc(s0)holds true: ˆ Gλc(s0)∼ =ˆ Gλc(s1) + 0,if s0<0, −V s1,if s0= 0, (s0−s1)/λc,if s0>0, (4.6) for s1∈Rclose enough to s0, where Vbelongs to the subdifferential ∂ˆ Gλc(s0). Furthermore, the error of this approximation is of order O¡|s0−s1|2¢. Proof. Since ˆ Gλc, defined in (4.5), is differentiable at every point s6= 0 and strongly semismooth, we can write for all s0, s1∈Rthanks to Lemma 4.1.6 ˆ Gλc(s0) = ˆ Gλc(s1) + V(s0−s1) + O¡|s0−s1|2¢, where V∈∂ˆ Gλc(s0). Taking into account the definition of ˆ Gλc(see 4.5), we obtain (4.6). Remark 4.1.9. To solve numerically the contact problem we consider the election V= 0 in (4.6). In order to introduce an iterative algorithm to deal with the contact condition, we discretize the problem (4.1)-(4.3) in time in the same way that in Section 3.5.3. Then, from Lemma 4.1.8, at each time step tj+1, given the starting values (uj+1 0,σj+1 0, p j+1 0), successive approximations (uj+1 k,σj+1 k, p j+1 k), k ≥1, of the solution (uj+1,σj+1, p j+1)can be computed by using the recurrence formulas: ZΩj+1 σj+1 k:ε(v)dx +ZΓC pj+1 kvndγ =ZΩj+1 fj+1 ·vdx, ∀v∈Uh(tj+1),(4.7) σj+1 k= ¡Λs(Tj+1)¢−1³ε(uj+1 k)−αs(Tj+1)(Tj+1 −Tl)I´in Ωj+1 s, Λ−1 lε(uj+1 k)in Ωj+1 l. (4.8) pj+1 k=pj+1 k−1 2 + 0,if sj+1 k−1≤0, ³sj+1 k−sj+1 k−1´/λc,if sj+1 k−1>0, (4.9) where sj+1 k=³uj+1 k´n+λcpj+1 k,(4.10) pj+1 k−1 2 =ˆ Gλc(sj+1 k−1).(4.11)
4.1. Newton algorithm to approximate the contact multiplier 37 Notice that approximation (4.9) corresponds to the choice s1=sj+1 k−1,s0=sj+1 kand V= 0 for s0= 0 in (4.6). Remark 4.1.10. Note that at each iteration the variational equality corresponding to (4.7) is equivalent to the linear problem −div ³σj+1 k´=fj+1 in Ωj+1,(4.12) (uj+1 k)n= 0,(σj+1 k)t=0on Γ2 n(tj+1)∪Γl(tj+1)∪Γsym(tj+1),(4.13) σj+1 kn=0on Γ1 n(tj+1)∪Γsup(tj+1),(4.14) (σj+1 k)n=−pj+1 k,(σj+1 k)t=0on ΓC,(4.15) in the distributional sense; so the contact multiplier pj+1 krepresents the obstacle reaction on ΓC. When analyzing formulas (4.7)-(4.11), we notice that there exists a coupling between the displacements and the contact multiplier on the boundary ΓC; furthermore, from (4.15), normal stresses only depend on the contact multiplier on this boundary. In what follows, our objective is to split ΓCat each iteration into two parts: on the first part we will know the contact multiplier –and, in turn, the normal stresses–, on the second part the normal displacements will be null. This procedure will first allow us to compute displacements and stresses at each iteration, and then to compute the contact multiplier. Lemma 4.1.11. Let us consider the following sets: (Γ− C,k)j+1 ={C∈Sh;sj+1 k≤0},(4.16) (Γ+ C,k)j+1 ={C∈Sh;sj+1 k>0}.(4.17) Then, pj+1 k= 0 on (Γ− C,k−1)j+1, ³uj+1 k´n= 0 on (Γ+ C,k−1)j+1.(4.18) Proof. With the notation introduced in (4.16) and (4.17), relation (4.9) can be rewritten as pj+1 k= 0,on (Γ− C,k−1)j+1, pj+1 k−1 2 +³sj+1 k−sj+1 k−1´/λc,on (Γ+ C,k−1)j+1, since ˆ GλCis null on each C∈(Γ− C,k−1)j+1. So, we can easily update the multiplier pj+1 kby zero on (Γ− C,k−1)j+1. Nevertheless, on each face of (Γ+ C,k−1)j+1, taking into account (4.10) and (4.11), we must compute pj+1 kverifying pj+1 k=sj+1 k−1 λc +³sj+1 k−sj+1 k−1´ λc =1 λc³uj+1 k´n+pj+1 k, and this equality can only be verified if ³uj+1 k´n= 0 on (Γ+ C,k−1)j+1.
38 Chapter 4. Improved numerical solution Thus, at each iteration the boundary ΓCis split into two parts: •(Γ+ C,k−1)j+1, where the normal displacements and the tangential component of stresses are null, •(Γ− C,k−1)j+1, where the normal and tangential stresses are null. It is very important to note that this division of ΓCis the only data which is modified at each iteration. There are several ways to achieve the solution using formulas (4.7)-(4.11), but since our objective is to couple the contact effects with the viscoplastic ones, in the next subsection we propose an algorithm based on two steps. In the first one, we compute (uj+1 k,σj+1 k)using Lemma 4.1.11 and Remark 4.1.10 to remove the contact multiplier in equation (4.7). Afterwards, we compute pj+1 kand update (Γ+ C,k)j+1. Newton algorithm to solve the contact problem The proposed algorithm to solve Problem (3.10)-(3.12) is the following: •Given u0,σ0, we choose p0=−(σ0)n. •Then, for j≥0,(uj,σj, pj)known at time tj, we determine (uj+1,σj+1, pj+1)at time tj+1, an approximated weak solution of Problem (3.10)-(3.12). To do so, we propose the following iterative algorithm: Step 1: Given (uj+1 k−1,σj+1 k−1, p j+1 k−1),we compute (uj+1 k,σj+1 k)∈U∗ h(tj+1)×Xh(tj+1)by solving the variational equality ZΩj+1 σj+1 k:ε(v)dx =ZΩj+1 fj+1 ·vdx, ∀v∈U∗ h(tj+1), together with the constitutive law (4.8), where U∗ h(tj+1) = {v∈Uh(tj+1); vn= 0 on (Γ+ C,k−1)j+1}. In this first step the tangential component of σj+1 kis assumed to be null on ΓC,but normal stresses due to the contact are unknown on (Γ+ C,k−1)j+1. For the practical computation we impose the Dirichlet condition (4.18) on (Γ+ C,k−1)j+1 by using a penalty method in the variational equation ZΩj+1 σj+1 k:ε(v)dx +1 ²cZ(Γ+ C,k−1)j+1 ³uj+1 k´nvndγ =ZΩj+1 fj+1 ·vdx, (4.19) for all v∈Uh(tj+1),where ²cis a small parameter.
4.1. Newton algorithm to approximate the contact multiplier 39 Step 2: Known (uj+1 k,σj+1 k)and given that pj+1 k= 0 on (Γ− C,k−1)j+1,we update pj+1 kon (Γ+ C,k−1)j+1 thanks to the variational equality (4.7): Z(Γ+ C,k−1)j+1 pj+1 kvndγ =−ZΩj+1 σj+1 k:ε(v)dx +ZΩj+1 fj+1 ·vdx, , (4.20) for all v∈Uh(tj+1)and we update (Γ+ C,k)j+1 using (4.17). Then, taking into account expression (4.19) when updating pj+1 kon (Γ+ C,k−1)j+1, we obtain that Z(Γ+ C,k−1)j+1 pj+1 kvndγ =1 ²cZ(Γ+ C,k−1)j+1 ³uj+1 k´nvndγ, ∀v∈Uh(tj+1), that is, pj+1 k=1 ²c³uj+1 k´non (Γ+ C,k−1)j+1.(4.21) Matrix factorization To solve the matrix of equations system resulting from applying the finite element method we use a direct method of resolution. Although it might seem reasonable to employ an iterative method, such as the conjugate gradient, several factors discourage us from its use: the large gradients of temperature and stresses, the penalty term of the contact condition and the multiscale nature of the fictitious domain technique to impose the metallostatic pressure. Due to this, we decided to use a classical resolution method using a Cholesky factorization. The mean disadvantage of the Newton algorithm presented here is the need to factorize the stiffness matrix at each iteration, which implies an important increase in the computational cost with respect to the fixed point algorithm presented in Section 3.5.3. In order to avoid this difficulty, taking into account that the contact condition only involves a small part of the boundary, we propose to use a factorization of the stiffness matrix adapted to the problem’s geometry. For this purpose, we take into account that the computational domain grows with time far from the contact region. Then, there are two possibilities: •To consider a lower factorization K=LDLT–where Dis a diagonal matrix and Lis a lower triangular matrix with unitary diagonal– and number the mesh in such a way that the last nodes correspond to the contact nodes. This method is not advisable in the casting case since the computational domain grows with time far from the contact region, and the mesh is constructed by layers to model the filling process. •To consider an upper factorization K=UDUTwhere Uis an upper triangular matrix with unitary diagonal (see Figure 4.1) and number the mesh in such a way that the first nodes correspond to the contact nodes. The factorization is carried out in the following way:
40 Chapter 4. Improved numerical solution KUDUT Figure 4.1: Scheme of the upper factorization Adapted partial factorization. We consider the upper factorization K=UDUTwhere Uis an upper triangular matrix with unitary diagonal and it is easy to prove that Dand Uare characterized by the following formulae: dnn =knn, i=n−1, n −2, . . . , 1 uin =kin dnn , uij = kij − n X l=j+1 uildllujl djj , j =n−1, n −2, . . . , i + 1, dii =kii − n X l=i+1 u2 ildll. In the practice, the factorization is stored in the matrix Kitself: Dwill be the diagonal of Kand Uits upper triangular submatrix. The pseudocode of the factorization is the following: for i=n−1to 1do kin =kin/knn for j=n−1to ido for l=j+ 1 to ndo kij =kij −kilkllkjl end for if i6=jthen kij =kij/kjj end if end for end for Since with this methodology the contact rows are the first ones of the matrix, if we compute the factorization row by row from downwards to upwards and from the right to the left, during a time step we can conserve all the rows of the factorized matrix except those corresponding to the contact nodes (see Figure 4.3). So, at each iteration we must only compute and factorize the first rows of the matrix, corresponding to the contact nodes.
4.1. Newton algorithm to approximate the contact multiplier 41 Figure 4.2: Upper skyline by rows Contactrows Firstdirection Seconddirection Figure 4.3: Direction of the upper factorization Matrix storage. In order to storage the matrix K, we employ a non standard storage by means of an upper skyline by rows. So, we store the upper matrix of Kin vector form, row by row from upwards to downwards and from the left to the right. We denote by µKthe array defined by (µK)1= 0, (µK)2= 1, (µK)i+1 = (µK)i+µi−1−(i−1) + 1, where (µK)i+1 gives the position of the diagonal term Kii,µibeing the column corresponding to the last non null element of the ith row of the matrix K(see Figure 4.2). Figure 4.4 shows a flowchart of the Newton algorithm with the new matrix factorization. After initializing variables (Step 1 in Figure 4.4), a time step loop is carried out. In this loop, at the j+1 iteration, the displacement vector uj+1 is computed by constructing the fixed stiffness matrix (Step 2) and solving iteratively (Steps 5-8) equations (4.19) and (4.21), where the variable stiffness matrix is constructed at each iteration (Step 5).
42 Chapter 4. Improved numerical solution Complete factorization of fixed stiffness matrix 3. force vector Construction of 4. Start Read input data Initialization of variables fixed stiffness matrix 2. Construction of 1. 5. variable stiffness matrix Construction of 6. variable stiffness matrix Partial factorization of Stop Yes No step? Last time Yes Stop Test? No boundary & multiplier Solution of the system 7. 8.Updating of contact Time step loop Iterations loop Figure 4.4: Flowchart of Newton algorithm with adapted factorization.
4.2. Algorithms for the numerical solution of the nonlinear constitutive law 43 4.2 Algorithms for the numerical solution of the nonlinear constitutive law In order to avoid the nonlinearity due to the constitutive law (3.5), the numerical solution is based on the Bermúdez-Moreno algorithm which involves a viscoplastic multiplier, which is a fixed point of a nonlinear equation. In a first stage, as in the contact case, the problem was solved in [10] by using the Bermúdez-Moreno algorithm described in Section 3.5.4. Although this algorithm was very robust, its convergence turned out to be slow in butt curl simulation. In order to improve its efficiency, we propose two new methodologies: •To apply a Newton method to compute the multiplier. •To consider the Bermúdez-Moreno algorithm with variable parameters introduced in [43]. As we have done in previous subsection, we are going to consider a simplified problem, corresponding to Problem (DV P ), replacing the contact boundary condition by a null Dirichlet condition. Therefore, recall that the discretized variational (3.16)-(3.18) is: ZΩ(t) σ(t) : ε(v)dx =ZΩ(t) f(t)·vdx, ∀v∈Uadh(t), ε(˙ u)(t) = ˙ z { (Λs(T)σ)(t)+(∇Φq(σD))(t) + ³αs(T)˙ T´(t)Iin Ωs(t), ˙ z { (Λlσ)(t)in Ωl(t), u(0) = u0,σ(0) = σ0,in Ω(0), where Uadh(t) = {vh∈Uh(t); vh=0on ΓC}. As in Section 3.5.4, the thermo-elasto-viscoplastic law is discretized in time by using an implicit Euler scheme: ε(uj+1)−ε(uj) = ¡Λs(Tj+1)σj+1 −Λs(Tj)σj¢+αs(Tj+1)(Tj+1 −Tj)I+ ∆t∇Φq³¡σj+1¢D´in Ωj+1 s,(4.22) ε(uj+1)−ε(uj) =Λlσj+1 −Λlσjin Ωj+1 l. 4.2.1 Newton algorithm to approximate the viscoplastic multiplier In order to introduce a Newton algorithm to approximate the viscoplastic multiplier, in this section, we consider γp= 0 and λp>0a real positive number in Proposition 3.5.7 and Lemma 3.5.8. Therefore, we obtain the following expression for the stress tensor in Ωs(tj+1): σj+1 = (Λs(Tj+1))−1¡ε(uj+1)−∆tqj+1 +Fj¢,(4.23)
44 Chapter 4. Improved numerical solution where Fj=−ε(uj)+Λs(Tj)σj−αs(Tj+1)(Tj+1 −Tj)I,(4.24) is the history of the solidified metal up to time tjand the viscoplastic multiplier qj+1 verifies qj+1 = (∇Φq)λp³¡σj+1¢D+λpqj+1´,(4.25) being (∇Φq)λp(ζ) = 1 λpµ1−1 ηj+1 ¶ζ,(4.26) where ηj+1 =η(ζ)is the unique root of the equation ηq−1−ηq−2−λpθ(Tj+1)|ζ|q−2= 0,(4.27) in the interval [1,+∞). Remark 4.2.1. Since ∇Φqis a maximal monotone operator, (∇Φq)λpis contractive for λp>1 (see [16, 26, 86]). Then, the fixed point of (4.25) is unique. In order to couple the two nonlinearities due to the contact condition and the constitutive law in casting processes, a key component of our procedure must be that the stiffness matrix is preserved and we do not need to recalculate it at each iteration. At the same time we wish to achieve the efficiency of Newton techniques. To this end, assuming that the stress tensor is known, we linearize the equation (4.25) to obtain an approximation of the viscoplastic multiplier. Lemma 4.2.2. Given any σ∈ S3and λp>1, the fixed point q∈ S3of the equation q= (∇Φq)λp¡σD+λpq¢,(4.28) can be approximated by q∼ =(∇Φq)λp(σD+λpˆq) + µ1−1 η¶(q−ˆq) + ¯η¡¡σD+λpˆq¢: (q−ˆq)¢¡σD+λpˆq¢,(4.29) for ˆq ∈ S3close enough to qand provided that σD+λpˆq 6=0. Here η=η¡σD+λpˆq¢, and ¯ηis given by the expression ¯η=(q−2) (η−1) η((q−1)η−(q−2)) |σD+λpˆq|2.(4.30) The order of this approximation is O¡|q−ˆq|2¢. Proof. Since (∇Φq)λp(ξ)is twice-differentiable for ξ6=0, we can apply Taylor’s formula to approximate q. Then, for ˆ q∈ S3close enough, if σD+λpˆq 6=0, we obtain from (4.28) q= (∇Φq)λp(σD+λpˆq) + D(∇Φq)λp(σD+λpˆq) (λp(q−ˆq)) + O¡|q−ˆq|2¢.(4.31)
4.2. Algorithms for the numerical solution of the nonlinear constitutive law 45 From expression (4.26) we deduce that D(∇Φq)λp(τ)(ξ) = ξ λpµ1−1 η(τ)¶+τ λp Dη(τ)(ξ) η(τ)2,∀τ,ξ∈ S3.(4.32) Deriving implicitly with respect to ξin expression (4.27), Dη(τ)(ξ) = λpθ(T)(q−2)|τ|q−4 (q−1)η(τ)q−2−(q−2)η(τ)q−3(τ:ξ),∀τ,ξ∈ S3.(4.33) From (4.32) and (4.33), D(∇Φq)λp(τ)(ξ) = ξ λpµ1−1 η(τ)¶+θ(T)(q−2)|τ|q−4(τ:ξ) (q−1)η(τ)q−(q−2)η(τ)q−1τ. Therefore, from (4.31), q= (∇Φq)λp(σD+λpˆq) + µ1−1 η¶(q−ˆq) + ¯η¡¡σD+λpˆq¢: (q−ˆq)¢¡σD+λpˆq¢+O¡|q−ˆq|2¢, where η=η¡σD+λpˆq¢and ¯η=λpθ(T) (q−2) |σD+λpˆq|q−4 (q−1) ηq−(q−2) ηq−1. Taking into account expression (4.27), we can simplify the expression for ¯ηin the following form: ¯η=λpθ(T) (q−2) |σD+λpˆq|q−4 (q−1) ηq−(q−2) ηq−1=(q−2) |σD+λpˆq|−2(ηq−1−ηq−2) (q−1) ηq−(q−2) ηq−1= =(q−2) (η−1) η((q−1) η−(q−2)) |σD+λpˆq|2, and therefore, expression (4.30) is obtained. A detailed explanation of this procedure can be found in [73]. Remark 4.2.3. The case σD+λpˆq =0was excluded from Lemma 4.2.2 since from equation (4.25), if σD+λpq=0, then q=0is the fixed point. Summing up, if we replace expression (4.23) in equality (3.16) at time step tj+1, we obtain a variational formulation in displacements of our problem. So, we propose the following iterative algorithm: At each time step tj+1, given the starting values (uj+1 0,σj+1 0,qj+1 0), successive approximations (uj+1 k,σj+1 k,qj+1 k),k≥1, of the solution (uj+1,σj+1,qj+1)can be computed by using the recurrence formulas: ZΩj+1 s¡Λs(Tj+1)¢−1ε(uj+1 k) : ε(v)dx +ZΩj+1 l Λ−1 lε(uj+1 k) : ε(v)dx = ZΩj+1 s¡Λs(Tj+1)¢−1³∆tqj+1 k−1´:ε(v)dx −ZΩj+1 s¡Λs(Tj+1)¢−1Fj:ε(v)dx+ ZΩj+1 fj+1 ·vdx ∀v∈Uadh(tj+1),(4.34)
52 Chapter 4. Improved numerical solution where Fj=−ε(uj)+Λs(Tj)σj−αs(Tj+1)(Tj+1 −Tj)¡1 + ∆tγj psj+1¢I+ ∆tγj psj+1 3µ1 sjtr(σj)−tr(ε(uj))¶Iin Ωj+1 s,(4.45) and σj+1 k= V(Tj+1)Ãε(uj+1 k)−∆tqj+1 k−1+∆tγj psj+1 3tr(ε(uj+1 k))I+Fj!in Ωj+1 s, Λ−1 lε(uj+1 k)in Ωj+1 l, (4.46) The updated viscoplastic multiplier is given by qj+1 k=³σj+1 k´D+λj pqj+1 k−1 λj pÃ1−1 ηj+1 k(1 −λj pγj p)!in Ωj+1 s,(4.47) and ηj+1 k=ηµ³σj+1 k´D+λj pqj+1 k−1¶is the unique root of the equation ηq−1−ηq−2−λj pθ(Tj+1) (1 −λj pγj p)q−1¯¯¯¯³σj+1 k´D+λj pqj+1 k−1¯¯¯¯ q−2 = 0, in the interval [1,+∞). The main question is to determine the suitable technique to update the value of parameters γj p, λj pfrom the solution obtained at the time step tj. Taking into account that the optimal parameters to convergence verify λj pγj p= 1/2(see [16, 66]), the problem is reduced to calculate, for example, γj p. Application to plane strain case Due to the great computational demands of the three-dimensional simulation in casting processes, in the numerical simulations presented in this manuscript we consider the two-dimensional associated problem under the plane strain assumption on the symmetry plane [x1= 0]. Due to this, in this section we present the algorithm to approximate the viscoplastic multiplier under this assumption. So, at each element K∈ Th(tj), any τ∈Xh(tj)can be written in the matrix form τ|K= τ11 0 0 0τ22 τ23 0τ23 τ33 ∈ M3×3(R). Taking into account that the Norton-Hoff’s law (2.6) only involves the deviatoric part of the stress tensor, in this study we consider for each element K∈ Th(tj): VD={τ= (τij)∈ M3×3(R); τij =τji, τ12 =τ13 = 0, τ11 =−(τ22 +τ33)},
4.2. Algorithms for the numerical solution of the nonlinear constitutive law 53 the space of matrices with null trace under the plane strain assumption. It is easy to prove that a basis of the space VDis BVD=©τ1,τ2,τ3ªwith τ1= −100 0 1 0 0 0 0 ,τ2= −100 0 0 0 0 0 1 ,τ3= 000 001 010 . We recall that in Section 3.5.4 the following expression for the viscoplastic multiplier qjwas obtained from the Bermúdez-Moreno Lemma: qj= (∇Φq)γp 1/(2γp)µ¡σj¢D+1 2γp qj¶,(4.48) for γp>0and where σjis given by equation (4.46) at the time step tj. In order to simplify the notation, we denote by Hthe application of the dissipation potential H:VD−→ VD τ−→ H(τ) = (∇Φq) (τ) ; the regularized Yosida of the perturbed operator H−γpI, Hγp 1/(2γp):VD−→ VD τ−→ Hγp 1/(2γp)(τ) = ³(∇Φq)γp 1/(2γp)´(τ) ; and, at the time step tj, ˜ Hγp 1/(2γp):VD−→ VD τ−→ ˜ Hγp 1/(2γp)(τ) = ³(∇Φq)γp 1/(2γp)´µ¡σj¢D+1 2γp τ¶. So, with this notation, the viscoplastic multiplier qjcan be rewritten as: qj|K=Hγp 1/(2γp)µ¡σj¢D|K+1 2γp qj|K¶=˜ Hγp 1/(2γp)¡qj|K¢,∀K∈ Th(tj),(4.49) that is, qj|Kis obtained as the fixed point of ˜ Hγp 1/(2γp). In [43], Gallardo et al. introduced a generalized Yosida regularization and they presented a generalization of the Bermúdez-Moreno algorithm that allows the use of very general operators as parameters. Moreover, as a particular case, they analyzed the use of scalar and matrix-valued parameters, proving that the optimal choice of the parameter γj pwould be a matrix-valued function. Furthermore, in [4], Arregui et al. presented an application of this generalized algorithm with scalar parameters in one-dimensional problems. Taking into account the complexity of the real problem, in this study we propose to consider that the parameter γpis a scalar function. We will see in Chapter 5 that, although in [43] this choice did not give good results in academic cases, in the casting process the convergence improvement
54 Chapter 4. Improved numerical solution is considerable. This fact is due to that the viscoplastic behaviour varies greatly from one element to another and from one time step to another, which explains the poor convergence results when considering γpconstant. Then, the strategy proposed in [43] consists of finding a parameter γj psuch that D˜ Hγj p 1/(2γj p)(qj|K) = 0.(4.50) To simplify the notation, from now on in this section we will omit the restriction to an element K∈ Th(tj). Therefore, in order to give a scalar expression for γp, we try to minimize ρj γp, the spectral radius of D˜ Hγj p 1/(2γj p)(qj). To do that, we give the following Lemma: Lemma 4.2.5. The differential D˜ Hγp 1/(2γp)(qj)is given by the expression D˜ Hγp 1/(2γp)(qj) = "I−µ1 2I+1 2γp DH ³¡σj¢D´¶−1#, with qj= (∇Φq)γp³¡σj¢D´. Proof. In [16] it is proved the following result: Hγp λp(τ) = 1 1−λpγp Hλp 1−λpγpµ1 1−λpγp τ¶−γp 1−λpγp τ,τ∈VD, for all λp, γp>0, λpγp<1. So, for γpλp= 1/2we have that Hγp 1/(2γp)(τ) =2H1 γp (2τ)−2γpτ. By using the definition of the Yosida regularization, we obtain that Hγp 1/(2γp)(τ) =2γp"I−µI+1 γp H¶−1#(2τ)−2γpτ= 2γp"τ−µI+1 γp H¶−1 (2τ)#= 2γp"τ−µ1 2I+1 2γp H¶−1 (τ)#= 2γp"I−µ1 2I+1 2γp H¶−1#(τ). Then, applying the chain rule it results DHγp 1/(2γp)(τ) =2γp"I−Dõ1 2I+1 2γp H¶−1!(τ)#= 2γp"I−µ1 2I+1 2γp DH(η)¶−1#,(4.51)
4.2. Algorithms for the numerical solution of the nonlinear constitutive law 55 for η=µ1 2I+1 2γp H¶−1 (τ),(4.52) or equivalently, τ=µ1 2I+1 2γp H¶(η) = 1 2η+1 2γp H(η).(4.53) In particular, if we consider in (4.52) τ=¡σj¢D+1 2γp qj,(4.54) and taking into account that qjis defined as the perturbed operator of Hwith parameter γp(see (3.20)), we have that τ=¡σj¢D+1 2γp qj=¡σj¢D+1 2γphH³¡σj¢D´−γp¡σj¢Di= 1 2¡σj¢D+1 2γp H³¡σj¢D´.(4.55) From (4.53) and (4.55), we deduce that for τdefined by (4.54) the corresponding ηis η=¡σj¢D. Then, from (4.49) and (4.51), we infer that D˜ Hγp 1/(2γp)(qj) = 1 2γp DHγp 1/(2γp)µ¡σj¢D+1 2γp qj¶="I−µ1 2I+1 2γp DH ³¡σj¢D´¶−1#, and the Lemma is proved. Lemma 4.2.6. Let assume that DH ³¡σj¢D´is positive definite and consider that 0< ωj 1≤ωj 2≤ ωj 3are its eigenvalues. Then, ωj i−γp ωj i+γp are the eigenvalues of D˜ Hγp 1/(2γp)(qj)for qj= (∇Φq)γp³¡σj¢D´and its spectral radius is ρj γp= max i=1,2,3¯¯¯¯¯ ωj i−γp ωj i+γp¯¯¯¯¯.(4.56) Proof. Let µj ibe an eigenvalue of D˜ Hγp 1/(2γp)(qj); then, by using Lemma 4.2.5, there exists τ∈VD such that "I−µ1 2I+1 2γp DH ³¡σj¢D´¶−1#(τ) = µj iτ.
56 Chapter 4. Improved numerical solution By doing some simple calculations, we obtain (1 −µj i)µ1 2I+1 2γp DH ³¡σj¢D´¶(τ) = τ from where we deduce that DH ³¡σj¢D´(τ) = γp 1 + µj i 1−µj i τ. Let us denote ωj i=γp 1 + µj i 1−µj i ,(4.57) then ωj iis an eigenvalue of DH ³¡σj¢D´and we can write the eigenvalue µj iin terms of ωj iin the following manner: µj i=ωj i−γp ωj i+γp . The optimal γj pwill be that minimizing (4.56). In the following result we give an expression for γj p: Proposition 4.2.7. Under the assumptions of Lemma 4.2.6, the parameter γj pminimizing ργpis given by: γj p=qωj 1ωj 3.(4.58) Proof. The proof is similar to that given in [43] and it is based in the following property: ∀γp>0,ωj 1−γp ωj 1+γp≤ωj 2−γp ωj 2+γp≤ωj 3−γp ωj 3+γp . We have to analyze several cases depending on the position of γpwith respect to eigenvalues ωj i: •If 0< γp≤ωj 1, it is easy to prove that min 0<γp≤ωj 1 ργp=ωj 3−ωj 1 ωj 3+ωj 1 . •Analogously, min ωj 3≤γp ργp=ωj 3−ωj 1 ωj 3+ωj 1 .
4.2. Algorithms for the numerical solution of the nonlinear constitutive law 57 wj 1wj 2wj 3 r ( ) 1 p g r ( ) 2 p g r ( ) 3 p g g* p Functionsr( ) i p g (R )gp * Figure 4.5: Functions ri(γp) •Finally, let suppose that γp∈(ωj 1, ωj 3)and define ri(γp) = ¯¯¯¯¯ ωj i−γp ωj i+γp¯¯¯¯¯. Then, we have that (see Figure 4.5) ργp=(r3(γp)if ωj 1< γp≤γ∗ p, r1(γp)if γ∗ p< γp≤ωj 3, where γ∗ pis the unique element in the interval (ωj 1, ωj 3)such that r1(γ∗ p) = r3(γ∗ p). A simple calculation gives γ∗ p=qωj 1ωj 3. Notice that min ωj 1<γp<ωj 3 ργp=r1(γ∗ p) = r3(γ∗ p)≤ωj 3−ωj 1 ωj 3+ωj 1 , so the best choice for γpis precisely γ∗ p. Proposition 4.2.8. The associated matrix of DH ³¡σj¢D´in terms of the basis BVDis: DH ³¡σj¢D´= AjCj 21(σj)D 22 +BjAjCj 31(σj)D 22 2Aj(σj)D 23(σj)D 22 AjCj 21(σj)D 33 AjCj 31(σj)D 33 +Bj2Aj(σj)D 23(σj)D 33 AjCj 21(σj)D 23 AjCj 31(σj)D 23 2Aj(σj)D 23(σj)D 23 +Bj , where Aj=θ(Tj)(q−2) ¯¯¯¡σj¢D¯¯¯q−4,(4.59) Bj=θ(Tj)¯¯¯¡σj¢D¯¯¯q−2,(4.60) Cj 21 = (σj)D 22 −(σj)D 11, Cj 31 = (σj)D 33 −(σj)D 11.
58 Chapter 4. Improved numerical solution The eigenvalues of DH ³¡σj¢D´are given by the expressions: ωj 1=ωj 2=Bj, ωj 3=Bj+2 3Ajσ2 eq,(4.61) where σ2 eq is the Von Mises equivalent stress defined by σ2 eq =3 2¡σj¢D:¡σj¢D. Moreover, the eigenvalues verify that ωj i≥0,i= 1,2,3. Proof. It is easy to prove that DH ³¡σj¢D´is given by the expression DH ³¡σj¢D´τ=θ(Tj)(q−2) ¯¯¯¡σj¢D¯¯¯q−4³¡σj¢D:τ´¡σj¢D+θ(Tj)¯¯¯¡σj¢D¯¯¯q−2τ= Aj³¡σj¢D:τ´¡σj¢D+Bjτ, for all τ∈VD. Then, in order to obtain its associated matrix in terms of the basis BVD, we calculate the image of each element τi∈BVD, i=1,2,3, and we write them in terms of the basis BVD: DH ³¡σj¢D´(τ1) =Aj¡(σj)D 22 −(σj)D 11¢¡σj¢D+Bjτ1=AjCj 21 ¡σj¢D+Bjτ1= AjCj 21(σj)D 11 −Bj0 0 0AjCj 21(σj)D 22 +BjAjCj 21(σj)D 23 0AjCj 21(σj)D 23 AjCj 21(σj)D 33 = ³AjCj 21(σj)D 22 +Bj´τ1+AjCj 21(σj)D 33τ2+AjCj 21(σj)D 23τ3, DH ³¡σj¢D´(τ2) =Aj¡(σj)D 33 −(σj)D 11¢¡σj¢D+Bjτ2=AjCj 31 ¡σj¢D+Bjτ2= AjCj 31(σj)D 11 −Bj0 0 0AjCj 31(σj)D 22 AjCj 31(σj)D 23 0AjCj 31(σj)D 23 AjCj 31(σj)D 33 +Bj = AjCj 31(σj)D 22τ1+³AjCj 31(σj)D 33 +Bj´τ2+AjCj 31(σj)D 23τ3, DH ³¡σj¢D´(τ3) =2Aj(σj)D 23 ¡σj¢D+Bjτ3= 2Aj(σj)D 23(σj)D 11 0 0 0 2Aj(σj)D 23(σj)D 22 2Aj(σj)D 23(σj)D 23 +Bj 0 2Aj(σj)D 23(σj)D 23 +Bj2Aj(σj)D 23(σj)D 33 = 2Aj(σj)D 23(σj)D 22τ1+ 2Aj(σj)D 23(σj)D 33τ2+¡2Aj(σj)D 23(σj)D 23 +Bj¢τ3. By using symbolic calculus, we obtain the eigenvalues of the associated matrix, that are given by expressions (4.61).
4.2. Algorithms for the numerical solution of the nonlinear constitutive law 59 Then, from Propositions 4.2.7 and 4.2.8, we can conclude the following result: Corollary 4.2.9. The optimal choice of parameters for the viscoplastic multiplier qjis γj p=qωj 1ωj 3, λj p=1 2γj p , where ωj 1and ωj 3are given by expression (4.61). Bermúdez-Moreno algorithm with variable parameters to solve the viscoplastic problem under plane strain assumption Summing up, the proposed algorithm to solve Problem (3.16)-(3.18) is: •Let (u0,σ0)and γ0 pbe given; we compute q0=∇Φγ0 p q(σD 0). •Then, for j≥0,(uj,σj,qj, γj p)known at time tj, we determine (uj+1,σj+1,qj+1, γj+1 p)at time tj+1, an approximated weak solution of Problem (3.16)-(3.18). To do so, we propose the following iterative algorithm: i) Given (uj+1 k−1,σj+1 k−1,qj+1 k−1),we compute (uj+1 k,σj+1 k)∈Uadh(tj+1)×Xh(tj+1)by solving the variational equality ZΩj+1 sV(Tj+1)"ε(uj+1 k) + ∆tγj psj+1 3tr(ε(uj+1 k))I#:ε(v)dx+ ZΩj+1 l Λ−1 lε(uj+1 k) : ε(v)dx =ZΩj+1 shV(Tj+1)∆tqj+1 k−1i:ε(v)dx− ZΩj+1 s£V(Tj+1)Fj¤:ε(v)dx +ZΩj+1 fj+1 ·vdx ∀v∈Uadh(tj+1), where Fj=−ε(uj)+Λs(Tj)σj−αs(Tj+1)(Tj+1 −Tj)¡1+∆tγj psj+1¢I+ ∆tγj psj+1 3µ1 sjtr(σj)−tr(ε(uj))¶Iin Ωj+1 s, and the updated stress tensor is given by σj+1 k= V(Tj+1)Ãε(uj+1 k)−∆tqj+1 k−1+∆tγj psj+1 3tr(ε(uj+1 k))I+Fj!in Ωj+1 s, Λ−1 lε(uj+1 k)in Ωj+1 l, and the discretized space Uadh(tj+1)is restricted to the plane strain case in the usual way.
60 Chapter 4. Improved numerical solution ii) The updated viscoplastic multiplier is given by qj+1 k= 2γj p"³σj+1 k´D+1 2γj p qj+1 k−1#Ã1−2 ηj+1 k!. •Once the convergence of the viscoplastic multiplier is achieved, the updated parameter γj+1 p is given by γj+1 p=qωj+1 1ωj+1 3, where ωj+1 1, ωj+1 3are given by (4.61). 4.3 Two new algorithms to solve the casting problem In this section, we are going to summarize two algorithms for the numerical resolution of the complete model (3.7)-(3.9). To do this, we combine the algorithm introduced in Section 4.1 to solve the contact condition with those introduced in Section 4.2 to deal with the constitutive law. As we will see in the next chapter, in the numerical simulations presented in this manuscript, we consider the two-dimensional problem under the assumption of plane strain on the symmetry plane [x1= 0]. Due to this, we are going to introduce the algorithms in the particular case of plane strain. The two algorithms are the following: •NM algorithm: In this case, we use the Bermúdez-Moreno algorithm together with Newton’s techniques, introduced in Sections 4.1 and 4.2.1. •NVPM algorithm: In order to overcome the poor convergence of the NM algorithm in the simulation of casting processes, we consider the following methodology: the Newton algorithm for the contact nonlinearity introduced in Section 4.1, together with the Bermúdez-Moreno algorithm with variable parameters for the viscoplastic nonlinearity presented in Section 4.2.2. Let us summarize these algorithms in the following sections. Notice that the discretized spaces are restricted to the plane strain case in the usual way. 4.3.1 NM algorithm 1. Let (u0,σ0)be given and let consider p0=−(σ0)nand q0=∇Φq¡σD 0¢. 2. Then, for j≥0,(uj,σj,qj, p j)known at time tj, we determine (uj+1,σj+1,qj+1, p j+1) at time tj+1, an approximated weak solution of Problem (3.7)-(3.9). To do so, we propose the following iterative algorithm: 2.1. Initialize uj+1 0=uj,qj+1 0=qj, p j+1 0=pj.
4.3. Two new algorithms to solve the casting problem 61 2.2. With qj+1 k−1, p j+1 k−1known, calculate (uj+1 k, p j+1 k)in two steps: 2.2.1. Solving the variational equality ZΩj+1 s¡Λs(Tj+1)¢−1ε(uj+1 k) : ε(v)dx +ZΩj+1 l Λ−1 lε(uj+1 k) : ε(v)dx+ 1 ²cZ(Γ+ C,k−1)j+1 ³uj+1 k´nvndγ =ZΩj+1 fj+1 ·vdx+ ZΩj+1 s¡Λs(Tj+1)¢−1³∆tqj+1 k−1−Fj´:ε(v)dx, ∀v∈Uh(tj+1),(4.62) where Fj= Λs(Tj)σj−ε(uj)−αs(Tj+1)(Tj+1 −Tj)Iin Ωj+1 s. 2.2.2. Updating the contact multiplier pj+1 k=1 ²c³uj+1 k´non (Γ+ C,k−1)j+1,(4.63) and the effective contact boundary (Γ+ C,k)j+1 ={C∈Sh;sj+1 k>0}, sj+1 k=³uj+1 k´n+λcpj+1 k.(4.64) 2.3. The updated stress tensor is defined by σj+1 k= ¡Λs(Tj+1)¢−1³ε(uj+1 k)−∆tqj+1 k−1+Fj´in Ωj+1 s, Λ−1 lε(uj+1 k)in Ωj+1 l. (4.65) 2.4. The updated viscoplastic multiplier is given by qj+1 k= 0,if κj+1 k=0, ³ηj+1 k−1´ λpóσj+1 k´D+³ηj+1 k−1´q−2 ηj+1 k κj+1 k! −³ηj+1 k−1´q−2 |κj+1 k|2³qj+1 k−1:κj+1 k´κj+1 k,if κj+1 k6=0, (4.66) in Ωj+1 s, where κj+1 k=³σj+1 k´D+λpqj+1 k−1, λp>1and ηj+1 kis the unique root of the equation ηq−1−ηq−2−λpθ0|κj+1 k|q−2= 0, in the interval [1,+∞). Figure 4.6 shows a flowchart of NM algorithm to solve aluminium casting problems. After initializing variables (Step 1 in Figure 4.6), a time step loop is carried out. In this loop, at the j+1 iteration, the displacement vector uj+1 is computed by constructing the fixed stiffness matrix (Step 2) and solving iteratively (Steps 5-7) equations (4.62)-(4.66), where the variable stiffness matrix is constructed at each iteration (Step 6).
68 Chapter 5. Numerical results Remark 5.1.1. Notice that, with stopping test (5.3), if the multiplier is close to zero we perform an absolute test with ˜ δ2and, in the other case, a relative test with ˜ δ. Figure 5.2 shows the displacements on the deformed configuration at two different time steps and the Von Mises norm of stresses at the last time step is represented in Figure 5.3. Figure 5.2: Displacements on the deformed configuration at two time steps. Figure 5.3: Von Mises norm of stresses at the last time step. In this test we compare the results obtained with the NM algorithm with those obtained with the FPM one. The L2relative error at the last time step obtained with both algorithms is 3.957 ×10−4min displacements and 7.452 ×10−3N/m2in stresses. Table 5.1 gathers the obtained results, showing that with NM algorithm the cpu-time decreases approximately 74% with respect to the FPM method and the number of iterations decreases 99%. In order to check the improvement of carrying out a matrix factorization adapted to the contact problem, the third column gathers the results for the NM algorithm without that factorization, denoted by NM∗.
5.2. Test 2: A problem with large gradients 69 Method FPM NM NM∗ cpu-time (s) 17480.61 4615.30 7315.46 iterations 258 3 3 Table 5.1: Cpu-time and number of iterations for ∆t= 0.1s. 5.2 Test 2: A problem with large gradients In this second test we want to solve an elasto-viscoplastic problem with displacements and stresses whose gradients have a magnitude similar to that of real aluminium casting. Nevertheless, it does not reproduce the diverse behaviours which appear in the real process. For this purpose, we use the NM algorithm introduced in Section 4.3.1, which involves an adimensionalization technique and a time step optimization, and the NVPM algorithm described in Section 4.3.2. As in the previous test, we compare the obtained results with those of FPM algorithm. Let (0,0.5s] be the time interval and Ωbe the cylindrical body of Test 1 (its axis is parallel to x1-direction and its section the square of dimensions 0.5m×0.5m in the plane x2x3). For the sake of simplicity, we consider Ωindependent of time. On ∂Ωwe distinguish four parts: ∂Ω = ¯ ΓD∪¯ ΓN∪¯ ΓC∪¯ Γ±, where ¯ ΓD=∂Ω∩[x2= 0],¯ ΓN=∂Ω∩([x2= 0.5] ∪[x3= 0.5]) , ¯ ΓC=∂Ω∩[x3= 0],¯ Γ±=∂Ω∩[x1=±1]. In this test, we consider that Ωis already a solidified block with constant temperature. Then, there are not thermal stresses and the material parameters corresponding to the constitutive law, independent of the temperature T, are E= 109N/m2, ν = 0.35, θ0= 1.953125 ×10−39m2/(sN), q = 6. These data are similar to those of an aluminium alloy. The problem to solve is div(σ) = 0in Ω, u=h(t) (0, x2,−x3)on ΓD, σn= 108gon ΓN, στ=0, un= 0 on Γ±, στ=0, σn≤0, un≤0, σnun= 0 on ΓC, ˙ ε(u) = Λ ˙ σ+θ0|σD|q−2σDin Ω, u(x, 0) = 0,σ(x, 0) = 0in Ω, (5.4) with h(t) = 1081 + ν Et+2 31040θ0t6and g(x, t) = ½tn,on ΓN∩[x2= 0.5], −tn,on ΓN∩[x3= 0.5].
70 Chapter 5. Numerical results Its solution is readily verifiable u(x, t) = h(t) (0, x2,−x3),σ(x, t) = 108 0 0 0 0t0 0 0 −t . ΓN ΓDΩ x2 ΓC x3 Figure 5.4: Computational domain. As we have mentioned, to reduce the size of the problem we consider the corresponding plane strain problem over the symmetry plane [x1= 0] (see Figure 5.4). To solve this problem we use a uniform spatial mesh with 800 elements and 441 nodes, corresponding to a discretization parameter ∆x= 0.025 m. For time discretization we use a small time step (∆t= 10−5) since the solution is nonlinear in time. We initialize the multipliers in this way q0=0, p0= 10−5, and we consider the following parameters: •for the FPM algorithm: γc= 0.5×1012, λc= 10−12, γp= 10−6, λp= 0.5×106, ϑ= 0.9; •for the NM algorithm: λc=λp= 1, ²c= 10−50, ϑ = 0.9. The test for convergence is performed on both multipliers. The stopping test on the contact multiplier pis detailed in (5.3).The stopping test on each component of the viscoplastic multiplier qis |(qk)il −(qk−1)il|<˜ δmax{˜ δ, |(qk)il|},∀K∈ Th,and |(qk)il −¡∇Φq(σD k)¢il |<˜ δmax{˜ δ, |¡∇Φq(σD k)¢il |},∀K∈ Th,(5.5) where ˜ δis a small parameter. Firstly, we have tested the Newton algorithm for the elasto-viscoplastic law considering a Dirichlet condition on ΓC. Nevertheless, convergence is not always achieved due to large gradient stresses. To overcome this difficulty, we use the adimensionalization technique introduced in Section 4.2.1 to reach the convergence of this algorithm.
5.2. Test 2: A problem with large gradients 71 Testing adimensionalization technique We consider that the adimensionalization parameter is ηδ= 10−8. From now on, we will denote by NMFthe Newton method together with this technique. Figure 5.5 shows the successive iterants calculated with and without adimensionalization. The cpu-time, the L2relative error at the last time step and the mean of the number of iterations are summarized in Table 5.2. Notice that with the NMFalgorithm the number of iterations decreases approximately 95% and the cpu-time 85% with respect to the FPM algorithm. Figure 5.5: Numerical approximation of viscoplastic multiplier. The top row of figures without adimensionalization and the second one with the proposed adimensionalization technique. The ∗ shows the successive iterants with the NMFalgorithm. Table 5.2: Problem (5.4) without contact with ∆t=˜ δ= 10−5. Method FPM NMF cpu-time (s) 18790.2 2744.0 iterations 57 3 u(m)-error 2.622 ×10−52.621 ×10−5 σ(N/m2)-error 1.245 ×10−51.244 ×10−5 Finally, adding the contact condition on ΓCto Problem (5.4), the reduction in cpu-time and in the number of iterations is analogous to that obtained previously (see Table 5.3).
72 Chapter 5. Numerical results Table 5.3: Problem (5.4) with ∆t=˜ δ= 10−5. Method FPM NMF cpu-time (s) 25256.8 5365.2 iterations 69 3 u(m)-error 2.918 ×10−52.917 ×10−5 σ(N/m2)-error 1.345 ×10−51.345 ×10−5 We notice that the NMFalgorithm needs less iterations for convergence, and also less cpu-time, than the FPM algorithm. Nevertheless, the FPM algorithm is more robust and it does not need the adimensionalization technique. Increasing the time step Since the time interval of interest in the aluminium casting is about 150s, if we want to employ the NMFalgorithm in casting simulation we must be able to increase the time step used in the previous simulations. If we do so, the FPM algorithm works well, but the NMFalgorithm does not always converge due to the large gradients. To solve this problem, we use the optimization technique on the time step, introduced in Section 4.2.1: given ∆t, if convergence is not achieved, we reduce the time step until the algorithm converges. Furthermore, to stabilize the Newton algorithm we employ an Armijo rule on the viscoplastic algorithm (see [20]). Summing up, we use the complete NM algorithm introduced in Section 4.3.1 to solve Problem 5.4. Table 5.4 shows the obtained results; the NM algorithm obtains smaller displacement and stress errors than the FPM algorithm; in this case, the time step is reduced from 0.0625 to 1.5625×10−2. Furthermore, the number of iterations with these techniques decreases 93% and the cpu-time 48%. Table 5.4: Problem (5.4) with ∆t= 0.0625 and ˜ δ= 10−3. Method FPM NM cpu-time (s) 6.58 3.39 iterations 254 18 u(m)-error 2.66 ×10−17.30 ×10−2 σ(N/m2)-error 1.63 ×10−12.70 ×10−2 Testing the NVPM algorithm From the above results we can state that NM method introduced in Section 4.3.1 is a good method to solve mechanical problems with elasto-viscoplastic laws of Maxwell-Norton type and a Signorini contact condition. Nevertheless, when using this method to solve the real casting
5.3. Numerical simulation of casting processes 73 problem, convergence of NM method slows down dramatically. Due to this, we propose to use the NVPM algorithm where the Newton algorithm in viscoplasticity is replaced by the Bermúdez- Moreno algorithm with variable parameters. The considered parameters for this algorithm are: λc= 1, ²c= 10−50, γ0 p= 0.1, ϑ = 0.9. Table 5.5 shows a comparison between the three mentioned algorithms. Although as expected the NM algorithm continues to obtain better results, the reduction in cpu-time (40%) and iterations (83%) of the NVPM algorithm with respect to the FPM one is considerable. Table 5.5: Problem (5.4) with ∆t= 0.0625 and ˜ δ= 10−3. Method FPM NM NVPM cpu-time (s) 6.58 3.39 3.92 iterations 254 18 43 u(m)-error 2.66 ×10−17.30 ×10−22.03 ×10−1 σ(N/m2)-error 1.63 ×10−12.70 ×10−28.40 ×10−2 5.3 Numerical simulation of casting processes In this section we are going to present the results obtained in the numerical simulation of a real casting process. Real data have been provided by ALCOA-INESPAL S.A. (A Coruña, Spain). We compare the efficiency between the FPM, NM and NVPM algorithms. The considered computational domain is shown in Figure 5.6. x3 x2 Gc G1 n(t) G2 n(t) Gsl(t) Gc Ws(t) Wl(t) Gsup(t) Gl(t) Gsym(t) Figure 5.6: Computational domain.
74 Chapter 5. Numerical results Mesh construction The mesh used for the simulation is built in a natural way. Firstly we make the mesh of the aluminium contained inside the bottom block before it begins to go down. After the start of the casting, at each time step we add the mesh of the liquid aluminium poured in that time interval. The mesh construction takes into account the peculiarities of casting processes: •The contact zone corresponds with the part of the slab rested on the mold. So, the mesh is constructed in such a way that the first nodes correspond to the contact nodes. •The slab deformation depends strongly on the thermal gradients, which are larger in the recently solidified zone, that varies with the time. So, the mesh is finest-grained where the thermal gradients are larger. •The computational domain grows with time. Then, the mesh is structured in layers to model the filling process and it is reconstructed at each time step adding the amount corresponding to the metal poured during that step. Figure 5.7 shows the mesh of the computational domain at the start and the end of the casting process. The initial mesh has 1320 elements and 732 nodes while the final mesh has 2520 elements and 1342 nodes. A)Initialmesh B)Meshafter145sofcasting Figure 5.7: Initial mesh and mesh after 145 s of casting The considered time interval is (0,145s] and in the numerical solution a time step ∆t= 1s is used. We must take into account that the NM algorithm has incorporated a time step optimization.
5.3. Numerical simulation of casting processes 75 Behaviour law data The parameters to characterize the thermo-elasto-viscoplastic law of aluminium (see (2.6)) have been introduced in the previous works [9, 10, 11] after carrying out a complete bibliographic search in the engineering literature (see [38, 57, 84]): •Elastic law: Values for Young’s modulus Eand Poisson’s coefficient νdepending on temperature are obtained from Table 5.6 by linear interpolation: Table 5.6: Elastic law parameters. T(◦C) E(109N/m2)ν 126.9 67 0.3134 326.9 58 0.3448 526.9 40 0.45 •Viscoplastic law: The viscoplastic parameter θ0is computed from the data of Table 5.7 as θ0=φ (√2k)q. Table 5.7: Viscoplastic law parameters. φ(s−1)k(109N/m2)G(Kcal/mol) q 4.0×1012 15.72 37.3 6 •Thermal law: We consider the density function introduced in [38], given in Kg/m3: ρ(T) = ½2700 −0.23Tif T≤Ts, 2360 if T=Tl, where Ts= 607◦Cand Tl= 649◦C. The density values in (Ts, Tl)are computed by linear interpolation. Algorithm parameters In the following, we detail the parameters associated to the algorithms used in this simulation: •for the FPM algorithm: γc= 0.5×106, λc= 10−6, γp= 1.25 ×10−12, λp= 0.4×1012, ϑ= 0.9;
76 Chapter 5. Numerical results •for the NM algorithm: λc= 1, ²c= 0.001, λp= 1, ϑ= 0.9, ηδ= 10−9; •for the NVPM algorithm: λc= 1, ²c= 0.001, γ0 p= 1.25 ×10−13, ϑ= 0.9, ηδ= 109; •for the stopping test in both multipliers, ˜ δ= 10−4. Finally, to impose the metallostatic pressure by using the fictitious domain method introduced in Section 3.5.1, we consider the following parameters in the liquid: ¯ λ= 1.2414 ×1011,¯µ= 1.3793 ×1010, α = 2, β = 1, ² = 10−3. Numerical results Figure 5.8 shows the isotherms obtained with the code developed in [63] after carry out a parameter adjustment, and the butt curl obtained at the last time step. Figure 5.8: Isotherms and butt curl deformation after 150s. The numerical simulation of the problem has been carried out by using the three algorithms. Although the NM algorithm has a very good convergence in the academic problems as we have seen in Sections 5.1 and 5.2, in real casting simulations, when a zone of the slab is recently solidified, sometimes the convergence is not achieved. Besides, the FPM algorithm has a good but slow convergence which depends strongly on the chosen parameters. Due to this, the NVPM algorithm incorporates the good convergence properties of the Newton algorithm in contact and
5.3. Numerical simulation of casting processes 77 improves the convergence of the fixed point algorithm in viscoplasticity, computing automatically the parameters. The cpu-time and the mean of the number of iterations obtained with the FPM and NVPM algorithms are summarized in Table 5.8. Data for the NM algorithm are not shown since, as we have said, convergence is not always achieved. Method FPM NVPM cpu-time (s) 2660.49 465.562 iterations 6593 89 Table 5.8: Cpu-time and number of iterations after 145sof casting. Notice that with the NVPM algorithm the cpu-time decreases approximately 83% with respect to the FPM algorithm and the number of iterations decreases 99%.
84 TABLE OF CONTENTS 10.3.3 Verification of assumption (H3) ........................ 133 10.3.4 Existence and uniqueness of solution for experiments I and II . . . . . . . . 134 10.4 Obtaining a new formula for the MOR . . . . . . . . . . . . . . . . . . . . . . . . . 135 10.4.1 Analytical solution. A new formula for MOR . . . . . . . . . . . . . . . . . 137 11 Computing the MOR 141 11.1 Numerical results for Experiment I . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 11.2 Numerical results for Experiment II . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 12 Conclusions 149 A Appendix: Modulus of rupture 151
Chapter 7 Introduction An increasing number of new challenges are constantly arising for engineers, comprising new applications, new materials or possibly both. Very often, the behaviour of a certain material under new conditions cannot be correctly predicted, and many trial tests are needed in order to determine whether it can be used or not. Therefore, an enormous effort is actually being made with the purpose of developing new prediction techniques that allow saving huge amounts of money and time in testing the materials. Strangely enough, there are very few studies about how to predict the behaviour of the different materials under simple conditions. Several factors, such as slight differences in the composition or in the microstructure, can dramatically affect the result of the considered tests, making useless previous tests to use on similar materials. This is particularly critical for the ceramic materials. In the common applications, materials are subjected to different forces, with complicated components of bending, torsion, compression and usually combinations of all of them. To predict the behaviour of a certain material, without performing expensive destructive tests, it is necessary to understand the mechanisms that affect its resistance to select the suitable failure criteria. The properties of the material are carefully measured in laboratory tests on a sample of material; these tests determine the value of their mechanical properties by measuring the applied forces and the corresponding relative changes of length. There exist many tests to determine material properties, classified into several groups attending to the applied forces, the sample conditions at the experiment instant, the environment conditions, etc. The most important experiments are: •Tensile or compression test: An axial load is applied to a sample in order to obtain a uniform stress distribution on the cross-section. •Torsion test: This test determine the behaviour of samples subjected to a pair of loads. •Bending or flexure test: This test determine the behaviour of materials subjected to simple loads, for example, concentrated or uniformly distributed loads. •Fatigue test: Used to study the behaviour of materials subjected to variable loads. •Impact test: This test is a particular case of previous tests; the load is applied suddenly. 85
86 Chapter 7. Introduction Figure 7.1: Three-point bending experiment. In Part II of this manuscript we will focus our study on bending tests. The relative simplicity of these tests has permitted its popularization for measuring the mechanical properties of brittle materials. One of the most common tests is the three-point bending test: a sample of a brittle material is placed between three cylinders, without additional support, while the upper cylinder applies an increasing gradually force until the beam breaks (see Figure 7.1). From this test a property known as modulus of rupture (MOR) is obtained, which represents the maximum surface stress of the bent beam at the instant of failure. In the engineering literature, the MOR is calculated by an explicit formula which involves the value of the load at failure, H, the distance, 2l, between the two lower cylinders or supports, and the second moment of inertia, I1, of the transversal section of the beam (see [25]). Nevertheless, discrepancies by no means negligible for MOR values have been noticed in experimental measurements depending on the specimen size, the ratio between the area of its section and the overall length or the relation between the distance of the two lower cylinders and the overall length of the beam. The fact that samples of same material can lead to different values for the MOR is something quite common but so far there are not many studies about this behaviour. The aim of Part II is to study the mathematical problem associated to three-point bending tests in order to increase the knowledge about the MOR for brittle materials. Specifically, we prove the existence of a unique solution to the elastic problem with frictionless unilateral contact that arises from the mechanical model of the three-point bending test. Furthermore, in order to justify the classic formula for the MOR, we are interested in finding its limit model by using asymptotic analysis techniques when the thickness of the sample beam goes to zero (see [2, 18, 36, 80, 81]). By using this method, classical theories like Bernouilli-Navier, Saint Venant, Timoshenko and Vlassov theories were studied and justified (see [18, 80]). Also beams with a unilateral contact with a rigid foundation have already been studied in [82, 83] but only when the beam is clamped at the ends. Despite the absence of Dirichlet conditions for the three-dimensional mechanical problem associated with the three-point bending test, the asymptotic method allow us to obtain the limit one-dimensional problems for bending and for traction and to identify the relationship between the failure load and the corresponding axial normal tension. All theoretical results are collected in [70]. Besides, we are interested in finding an effective method to compute the MOR from brittle materials. To do this, we use all the theoretical background presented here to compute various approximations of the MOR. The procedure is the following: Once we have carried out several laboratory tests with ceramic beams, we determine the maximum force that beams can support until breaking –the rupture load. By using the classic formula for the MOR, we calculate the first value for it –the theoretical MOR. Then, after making numerical simulations of the one and threedimensional models, considering the rupture load, we compute two new approximations for the MOR –the 1d and 3d numerical MOR. Finally, we obtain a fourth approach for the MOR by using a new expression obtained from the asymptotic analysis of our problem. All these methodologies are applied for cylindrical and rectangular beams made of porcelain, which is a brittle material used very frequently in engineering.
87 The new theoretical formula obtained in this manuscript leads us to a better approximation of the MOR for brittle materials. This expression takes into account not only the rupture load and the total length of the beam but also the distance between the two lower cylinders, the effect of the gravity and the distance between the ends of the beam and the lower cylinders (see [51]). The outline of Part II is as follows: In Chapter 8, we study the static behaviour of a three-dimensional elastic beam when is subjected to a three-point bending test. Under suitable compatibility conditions on the applied forces and on the geometry of the beam, we prove the existence of a unique solution for the associated contact elastic problem. These conditions of compatibility on the data come from the absence of a Dirichlet condition on the beam boundary. In Chapter 9, we study the asymptotic behaviour of this problem; in particular, we deduce the one-dimensional models associated to the displacement components, and we give the existence and uniqueness of solution for them. Moreover, we give an expression for the normal axial stress in the beam which is related to the MOR of brittle materials. In the final part of the chapter, we deal with the regularity of the solution for the bending problem and we prove some properties of the coincidence set. In Chapter 10, several examples for laboratory tests for cylindrical and rectangular beams made of porcelain are detailed. We show that the previous assumptions of compatibility are satisfied for these real experiments. Moreover, we present analytical solutions of the bending and axial displacements for these real experiments. Their solutions are obtained from their differential formulations given by means of the asymptotic analysis. Moreover, we introduce a new expression for the MOR. Finally, in Chapter 11 we gather all the approaches for the MOR of the porcelain by using the methodology described above. We apply this procedure for both cylindrical and rectangular beams. In Appendix A we offer a brief description of the classic expression for the MOR used in mechanics of materials given that this information is not usually included in the mathematical bibliography.
88 Chapter 7. Introduction
Chapter 8 Mathematical study of three-point bending tests In this chapter we introduce a mathematical model to simulate the deformation suffered by a beam subjected to the three-point bending test as well as its associated variational formulation. Moreover, we obtain a result of existence and uniqueness of solution under suitable compatibility conditions. In Section 8.1, we show the mathematical model associated with the three-point bending test: an elastic problem with a unilateral contact condition, without friction, formulated in displacements. In Section 8.2, we propose a weak formulation of the problem in the appropriate functional spaces. Despite not impose Dirichlet conditions, assuming a suitable condition of compatibility on the applied forces, we prove the existence of a solution in Section 8.3. Assuming additional conditions of symmetry in the geometry of the beam and in the applied forces, in Section 8.4 we demonstrate the uniqueness of the solution. 8.1 Mathematical model As usual in solid mechanics, Latin subscripts are understood to range over the integers {1,2,3} and Greek subscripts (other than ²) over the integers {1,2}. Moreover, Einstein summation over repeated subscripts is implied. Given two vectors u,v∈R3, with components (ui),(vi), respectively, u·v=uivirepresents their scalar product. We denote by R9 sthe space of second order symmetric tensors over R3, endowed with the usual scalar product σ:τ=σijτij,σ,τ∈R9 s. Let ωbe a bounded, open, connected Lipschitz subset of R2with area A. Given 0< ² ≤1, we define ω²=² ω and we consider the domain Ω²=ω²×(−L, L),L > 0, which corresponds to the reference configuration of the beam. Consequently, the beam length is 2Land the area of the 89
90 Chapter 8. Mathematical study of three-point bending tests cross section is A²=²2A. We refer the motion of the beam to a fixed Cartesian system 0x² 1x² 2x² 3 whose origin is situated at the center of gravity of the beam (see Figure 8.1). Figure 8.1: Sketch of the three-point bending test. The physical problem consists of determining the displacement vector field, u², and the stress tensor field, σ², that the beam ¯ Ω²suffers when it is subjected to a three-point bending test. 8.1.1 Boundary conditions Let Γ²be the boundary of the beam which is the union of the lateral boundary Γ² l=γ²×[−L, L], where γ²=∂ω², and the ends of the beam Γ² ±=ω²× {±L};Γ²is assumed to be Lipschitzcontinuous. An arbitrary point of ¯ Ω²is denoted by x²= (x² 1, x² 2, x² 3)and n²denotes the outward unit normal vector to Γ². Notice that n²(x²) = (n² 1, n² 2,0) on the lateral boundary Γ² l, where (n² 1(x² 1, x² 2), n² 2(x² 1, x² 2)) is the outward unit normal vector to ω²in (x² 1, x² 2)∈γ². Moreover, n²= (0,0,±1) on Γ² ±. Lateral boundary Γ² lis partitioned into three non-empty, open and disjointed parts Γ² C,Γ² N1 and Γ² N2satisfying ¯ Γl²=¯ Γ² C∪¯ Γ² N1∪¯ Γ² N2,where: •¯ Γ² Cis the region of the beam where, due to the presence of the two lower cylinders, the normal component of the displacement cannot be positive. To impose this constraint and take into account the reaction force exerted by the two lower cylinders when the contact is effective, we consider a Signorini condition (see [54]): σ² τ=0, σ² n≤0, u² n≤0, u² nσ² n= 0,on Γ² C(8.1) where u² nand σ² nare the normal displacement and the normal stress at each point of the boundary, respectively, and u² τand σ² τare the tangential displacement and the tangential
8.1. Mathematical model 91 component of σ²n², respectively. We suppose that the two-dimensional measure of Γ² Cis strictly positive. •Γ² N1corresponds to the region of the boundary in which the upper cylinder exerts a compression force with density h²known. •Γ² N2represents the remaining lateral boundary of the beam that is free of forces. The ends of the beam are also assumed to be free of forces. 8.1.2 Equilibrium equations Under the assumption of small displacements, the behaviour of the beam is governed by the usual equilibrium equations −div(σ²) = f²in Ω², where the divergence operator is given by the expression (div²(τ²))i=∂τ² ij ∂x² j , and f²denotes the volume density forces. In practice, we will only consider the gravitational forces. 8.1.3 Behaviour law Let u² idenote the icomponent of the displacement vector field u²; then the components of the infinitesimal strain field are given by ε² ij(u²) = 1 2Ã∂u² i ∂x² j +∂u² j ∂x² i!,(8.2) where ∂u² i ∂x² jdenotes the partial derivative of the component u² iwith respect to x² j. We assume that ε²(u²)is related to the stress tensor σ²through Hooke’s law (see [28] or [52]), σ²=Λ²ε²(u²).(8.3) In the above expression, Λ²is the fourth order linear elasticity tensor defined by Λ²τ=λ²tr(τ)I+ 2µ²τ,(8.4) where λ²and µ²are Lamé parameters of the material, related to Young’s modulus E²and Poisson’s ratio ν²by λ²=E²ν² (1 −2ν²)(1 + ν²), µ²=E² 2(1 + ν²).
92 Chapter 8. Mathematical study of three-point bending tests 8.1.4 Problem (P²) Summing up, the problem we must solve is the following: Problem (P²): Find the displacement vector field u(x)and the stress tensor field σ(x), at each point x∈Ω, satisfying: −div²(σ²) = f²in Ω²,(8.5) σ²n²=h²on Γ² N1,(8.6) σ²n²=0on Γ² N2∪Γ² ±,(8.7) σ² τ=0, σ² n≤0, u² n≤0, σ² nu² n= 0 on Γ² C,(8.8) σ²=Λ²ε²(u²) = λ²tr(ε²(u²))I+ 2µ²ε²(u²)in Ω².(8.9) 8.2 Weak formulation To carry out the mathematical analysis of Problem (P²), we introduce a weak formulation in a similar manner to those used in [40, 46, 54]. We consider the space of displacement vector fields H1(Ω²)=[H1(Ω²)]3which is a Hilbert space with the norm kv²k1,Ω²= 3 X i=1 kv² ik2 0,Ω²+ 3 X i,j=1 kε² ij(v²)k2 0,Ω² 1/2 , where k·k0,Ω²denotes the usual norm in L2(Ω²). The subset of kinematically admissible displacements is the convex set (see [54]) K(Ω²) = {v²∈H1(Ω²); v² n≤0on Γ² C}. Let us define the stress tensor space X(Ω²) = [L2(Ω²)]9 swhere [L2(Ω²)]9 s={τ²= (τ² ij)∈[L2(Ω²)]9;τ² ij =τ² ji}. X(Ω²)is a real Hilbert space endowed with the usual norm in [L2(Ω²)]9. Its subspace H(div,Ω²) = {τ²∈X(Ω²); div²(τ²)∈L2(Ω²)}, is also a Hilbert space with the norm kτ²kH(div,Ω²)=¡kτ²k2 0,Ω²+kdiv²(τ²)k2 0,Ω²¢1/2. We denote H1/2(Γ²) = [H1/2(Γ²)]3as the image set of the trace application from H1(Ω²). The following Green’s formula holds: ZΩ² τ²:ε²(v²)dx²+ZΩ² div²(τ²)·v²dx²=hτ²n²,v²iΓ², for all τ²∈H(div,Ω²)and v²∈H1(Ω²), where h·,·iΓ²represents the duality between H−1/2(Γ²) and H1/2(Γ²).
8.2. Weak formulation 93 8.2.1 Assumptions From now on we consider the following hypotheses. (H1²)The volume forces satisfy f²∈L2(Ω²)and the surface forces h²∈L2(Γ² N1). (H2²)The material is homogeneous and isotropic; so, the elasticity tensor Λ²defined in (8.9) is such that λ², µ²∈R, λ², µ²>0. (H3²)The applied forces verify the following compatibility condition ZΩ² f²·w²dx²+ZΓ² N1 h²·w²dΓ²≤0, for all w²∈K(Ω²)∩R(Ω²), where R(Ω²) = {w²∈H1(Ω²) : ε²(w²) = 0in [L2(Ω²)]9}, is the set of rigid displacements. Moreover, ZΩ² f²·w²dx²+ZΓ² N1 h²·w²dΓ²= 0, if and only if w²∈Rn(Ω²), where Rn(Ω²) = {w²∈R(Ω²); w² n= 0 a.e. on Γ² C}. (H4²)The set (K(Ω²)∩R(Ω²)) \Rn(Ω²)is non-empty. 8.2.2 Weak formulation of the problem To simplify the writing, the following notations are useful: •The external virtual work for the applied forces over the virtual displacement v²∈H1(Ω²) is defined as: F²(v²) = ZΩ² f²·v²dx²+ZΓ² N1 h²·v²dΓ².(8.10) •The internal virtual work is related to the bilinear form a²defined by a²(v²,w²) = ZΩ² Λ²ε²(v²) : ε²(w²)dx²,(8.11) for all v²,w²∈H1(Ω²).
100 Chapter 9. Asymptotic analysis of three-point bending tests 9.1 Rescaling functions and changing variables Following [80] and in order to simplify the calculations, we are going to introduce the mixed variational formulation of Problem (P²). Problem (MV P²): Find u²∈K(Ω²)and σ²∈X(Ω²)such that ZΩ²·1 + ν² E²σ²−ν² E²tr(σ²)I¸:τ²dx²=ZΩ² ε²(u²) : τ²dx²,∀τ²∈X(Ω²),(9.1) ZΩ² σ²:ε²(v²−u²)dx²≥F²(v²−u²),∀v²∈K(Ω²),(9.2) where E²and ν²are the Young’s modulus and the Poisson’s ratio of the material, respectively, related to the Lamé parameters by the usual relations. Remark 9.1.1. Under assumptions (H1²)-(H6²), there exists a unique solution of Problem (MV P²)in (H1 s(Ω²)∩K(Ω²)) ×X(Ω²)thanks to Theorem 8.4.1. Given that the dependence of the solution (u²,σ²)of Problem (MV P²)with respect to ²is very complex, we define an equivalent problem for which this dependence is explicit. To do so, we use a change of variable to a fixed domain and the subsequent rescaling of the displacements and stresses. The reference domain, independent of ², is Ω = ω×(−L, L),Ω = Ω1, γ=∂ω =γ1,Γ = ∂Ω = Γ1. We associate each point x∈Ωwith a point x²∈Ω²through the transformation (see [18]) Π²: Ω →Ω² x= (x1, x2, x3)→Π²(x) = x²= (²x1, ²x2, x3),(9.3) in such a way that Γ²= Π²(Γ),n²(x²) = n(x), where nis the outward unit normal vector to the surface Γ. We associate the functions v²∈H1(Ω²)with the scaled functions v(²)∈H1(Ω) through vα(²)(x) = ²v² α(x²), v3(²)(x) = v² 3(x²),(9.4) and the functions τ²∈X(Ω²)with the scaled functions τ(²)∈X(Ω) through ταβ(²)(x) = ²−2τ² αβ(x²), τ3β(²)(x) = ²−1τ² 3β(x²), τ33(²)(x) = τ² 33(x²).(9.5) Notice that (9.4)-(9.5) are not a priori assumptions on the order of magnitude of the displacements or stresses. It will be a result of the forthcoming analysis that the rescaled unknowns are actually bounded. In what follows, we consider the following additional assumptions on the data:
9.1. Rescaling functions and changing variables 101 (H7²)The Lamé parameters are independent of ²:λ²=λ,µ²=µ. Consequently, Young’s modulus and Poisson’s coefficient are independent of ². (H8²)There exist functions f∈[L2(Ω)]3and h∈[L2(ΓN1)]3independent of ²and such that f² α(x²) = ²fα(x), f² 3(x²) = f3(x),x∈Ω, h² α(x²) = ²2hα(x), h² 3(x²) = ²h3(x),x∈ΓN1. In addition, fand hverify: fα∈L2(Ω), f3∈H1((−L, L); L2(ω)), hα∈L2(ΓN1), h3∈H1((−δ, δ); L2(γN1)). Under the previous assumptions we obtain the following result: Theorem 9.1.2. Under assumptions (H1²)-(H8²), the scaled displacement and stress fields (u(²),σ(²)) of the solution (u²,σ²)∈(H1 s(Ω²)∩K(Ω²)) ×X(Ω²)of Problem (MV P²)is the unique solution of the problem: Problem (MV P (²)): Find u(²)∈H1 s(Ω) ∩K(Ω) and σ(²)∈X(Ω) such that a0(σ(²),τ) + ²2a2(σ(²),τ) + ²4a4(σ(²),τ)−b(τ,u(²)) = 0,∀τ∈X(Ω),(9.6) b(σ(²),v−u(²)) ≥F(v−u(²)),∀v∈H1 s(Ω) ∩K(Ω),(9.7) where a0(σ,τ) = ZΩ 1 Eσ33τ33dx, (9.8) a2(σ,τ) = ZΩ£2(1 + ν) Eσ3βτ3β−ν E(σ33τµµ +σµµτ33)¤dx, (9.9) a4(σ,τ) = ZΩ£1 + ν Eσαβ −ν Eσµµδαβ¤ταβdx, (9.10) b(τ,v) = ZΩ τijεij(v)dx, (9.11) F(v) = ZΩ fividx +ZΓN1 hividΓ.(9.12) Proof. Taking into account expression (9.4), we obtain the following relation between the linearized strain tensor ε²(v²)and the scaled linearized strain tensor ε(v): εαβ(v(²)) = ²2ε² αβ(v²), ε3β(v(²)) = ² ε² 3β(v²), ε33(v(²)) = ε² 33(v²). So, we have that ZΩ² f²·v²dx²+ZΓ² N1 h²·v²dΓ²=²2ÃZΩ f·v(²)dx +ZΓN1 h·v(²)dΓ!, ZΩ² σ²:ε²(v²)dx²=²2µZΩ σ(²) : ε(v(²))dx¶,
102 Chapter 9. Asymptotic analysis of three-point bending tests thanks to (9.3)-(9.5). Consequently, b(σ(²),v−u(²)) ≥F(v−u(²)),∀v∈H1 s(Ω) ∩K(Ω). Analogously, substituting expressions (9.4), (9.5) and (9.8)-(9.12) in Problem (MV P²)and applying the change of variable theorem, we obtain expression (9.6). 9.2 Asymptotic expansion method From now on in Part II, unless otherwise stated, we assume hypotheses (H1²)-(H8²). We suppose that the solution (u(²),σ(²)) ∈(H1 s(Ω) ∩K(Ω)) ×X(Ω) of Problem (MV P(²)) can be written as a formal expansion (u(²),σ(²)) = (u0,σ0) + ²2(u2,σ2) + ²4(u4,σ4) + o(²4),(9.13) where ui∈H1 s(Ω) ∩K(Ω) and σi∈X(Ω) are independent of ². Substituting expression (9.13) in Problem (MV P(²)) and identifying the coefficients of the same powers of ², we obtain formally that the terms of order zero must be a solution of the problem: Problem (MV P0): u0∈H1 s(Ω) ∩K(Ω),σ0∈X(Ω), a0(σ0,τ)−b(τ,u0) = 0,∀τ∈X(Ω),(9.14) b(σ0,v−u0)≥F(v−u0),∀v∈H1 s(Ω) ∩K(Ω). In order to characterize the element (u0,σ0), solution of Problem (MV P0), we consider the decomposition H1 s(Ω) = Ws2(Ω) ×H1 s(Ω), where Ws2(Ω) = {(v1, v2)∈[H1(Ω)]2; (v1, v2,0) ∈H1 s(Ω)}, H1 s(Ω) = {v∈H1(Ω); (0,0, v)∈H1 s(Ω)}. Taking into account that in the contact boundary ΓCthe outward normal vector is given by (n1, n2,0), we have that K(Ω) = K2(Ω) ×H1(Ω), where K2(Ω) = {(v1, v2)∈[H1(Ω)]2;vαnα≤0on ΓC}.
9.2. Asymptotic expansion method 103 Consequently, Problem (MV P 0)can be rewritten in the equivalent form: (u0 1, u0 2)∈Ws2(Ω) ∩K2(Ω), u0 3∈H1 s(Ω),σ0∈X(Ω),(9.15) ZΩ 1 Eσ0 33τ33dx −ZΩ ∂3u0 3τ33dx = 0,∀τ33 ∈L2(Ω),(9.16) ZΩ ε3β(u0)τ3βdx = 0,∀(τ3β)∈[L2(Ω)]2,(9.17) ZΩ εαβ(u0)ταβdx = 0,∀(ταβ)∈[L2(Ω)]4 s,(9.18) ZΩ (σ0 33∂3v3+σ0 3β∂βv3)dx =ZΩ f3v3dx +ZΓN1 h3v3dΓ,∀v3∈H1 s(Ω),(9.19) ZΩ [σ0 αβ∂α(vβ−u0 β) + σ0 3β∂3(vβ−u0 β)]dx ≥ZΩ fβ(vβ−u0 β)dx +ZΓN1 hβ(vβ−u0 β)dΓ,(9.20) ∀(vβ)∈Ws2(Ω) ∩K2(Ω). To characterize the solutions of Problem (9.15)-(9.20), we introduce the following notations: •The total force per unit of length is qi(x3) = Zω fi(x1, x2, x3)dx1dx2+χ[−δ,δ](x3)ZγN1(x3) hi(x1, x2, x3)dγ. (9.21) •The total moment about the origin of the cross-section per unit of length is defined as M(x3) = Zω (x2f3,−x1f3, x1f2−x2f1)dx1dx2+ χ[−δ,δ](x3)ZγN1(x3) (x2h3,−x1h3, x1h2−x2h1)dγ. (9.22) We also introduce the subspace of H1(Ω) corresponding to Bernouilli-Navier displacements (see [80]): VBN (Ω) = {v∈H1 s(Ω); vα∈H2(−L, L), v3=v3−xα∂3vα, v3∈H1(−L, L)}. Moreover, to take into account the hypotheses of symmetry, we consider the subspace of Hm(−L, L),Hm e(−L, L)and Hm o(−L, L),m= 1,2,corresponding to its even and odd functions, respectively. The following theorem characterizes the solutions of Problem (9.15)-(9.20): Theorem 9.2.1. Every solution of Problem (9.15)-(9.20) verifies: (i) The displacement field u0∈VBN (Ω);u0 3being given by u0 3(x1, x2, x3) = u0 3(x3)−x2∂3u0 2(x3),∀(x1, x2, x3)∈Ω,(9.23) with u0 3∈H1 o(−L, L). Moreover, u0 1is null and u0 2∈H2 e(−L, L).
104 Chapter 9. Asymptotic analysis of three-point bending tests (ii) The axial stress σ0 33 is given by σ0 33 =E∂3u0 3=E(∂3u0 3−x2∂33u0 2),(9.24) and the corresponding component of the bending moment is m0 2=Zω x2σ0 33dx1dx2=−EI1∂33u0 2. (iii) The stretch component u0 3is a solution of the axial problem: Problem (P0 3): Find u0 3∈H1 o(−L, L)such that ZL −L EA∂3u0 3∂3v0 3dx3=ZL −L q3v0 3dx3,(9.25) for all v0 3∈H1 o(−L, L). (iv) The displacement u0 2is a solution of the bending problem: Problem (P0 2): Find u0 2∈H2 e(−L, L),(0, u0 2)∈K2(Ω), such that ZL −L EI1∂33u0 2∂33(v0 2−u0 2)dx3≥ZL −L q2(v0 2−u0 2)dx3−ZL −L M1∂3(v0 2−u0 2)dx3,(9.26) for all v0 2∈H2 e(−L, L),(0, v0 2)∈K2(Ω). Proof. (i) By expression (9.18) we obtain that ∂1u0 1=∂2u0 2= 0, ∂1u0 2+∂2u0 1= 0. Then, we deduce u0 1(x1, x2, x3) = φ0 1(x2, x3), u0 2(x1, x2, x3) = φ0 2(x1, x3), ∂2φ0 1(x2, x3) + ∂1φ0 2(x1, x3) = 0. Given that the variables x1, x2and x3are independent, we find that ∂2φ0 1(x2, x3) = −∂1φ0 2(x1, x3) = z0(x3) and consequently, u0 1(x1, x2, x3) =z0 1(x3) + x2z0(x3), u0 2(x1, x2, x3) =z0 2(x3)−x1z0(x3). In addition, u0 1must be odd with respect to x1, so z0 1=z0≡0. Therefore, u0 1= 0 and u0 2=z0 2(x3). Moreover, (u0 1, u0 2)∈Ws2(Ω) and consequently, z0 2∈H1 e(−L, L). Furthermore, from (9.17), ε3β(u0) = 0, thus ∂1u0 3=−∂3u0 1= 0, ∂2u0 3=−∂3u0 2=−∂3z0 2.
9.2. Asymptotic expansion method 105 Integrating, we obtain that u0 3(x1, x2, x3) = u0 3(x3)−x2∂3z0 2(x3) = u0 3(x3)−x2∂3u0 2(x3).(9.27) Summing up, we have proved that: u0(x1, x2, x3) = (0, u0 2(x3), u0 3(x3)−x2∂3u0 2(x3)) ∈H1 s(Ω) ∩K(Ω).(9.28) In particular, ∂3u0 3=∂3u0 3−x2∂33u0 2∈L2(Ω), from which we infer that: ZΩ (∂3u0 3)2dx =AZL −L (∂3u0 3(x3))2dx3+I1ZL −L (∂33u0 2(x3))2dx3<∞; so u0 2∈H2 e(−L, L)and u0 3∈H1 o(−L, L). (ii) It is easy to prove (9.24) substituting expression (9.27) in equation (9.16). (iii) Substituting expression (9.24) in (9.19), we deduce that ZΩ [E(∂3u0 3−x2∂33u0 2)∂3v3+σ0 3β∂βv3]dx =ZΩ f3v3dx +ZΓN1 h3v3dΓ,∀v3∈H1 s(Ω). Choosing v3(x1, x2, x3) = v0 3(x3)∈H1 o(−L, L)as a test function, we have that ZΩ E(∂3u0 3−x2∂33u0 2)∂3v0 3dx =ZΩ f3v0 3dx +ZΓN1 h3v0 3dΓ. Taking into account (8.18), we conclude that ZL −L EA∂3u0 3∂3v0 3dx3=ZL −LµZω f3dx1dx2¶v0 3dx3+Zδ −δÃZγN1(x3) h3dγ!v0 3dx3, for all v0 3∈H1 o(−L, L). Therefore, from expression (9.21), u0 3is a solution of Problem (P0 3). (iv) By choosing (v1, v2)(x1, x2, x3) = (0, v0 2(x3)) ∈Ws2(Ω) ∩K2(Ω), with v0 2∈H2 e(−L, L), as a test function in inequality (9.20), we obtain that ZΩ σ0 32∂3(v0 2−u0 2)dx ≥ZΩ f2(v0 2−u0 2)dx +ZΓN1 h2(v0 2−u0 2)dΓ.(9.29) Moreover, choosing v3(x1, x2, x3) = x2∂3(v0 2−u0 2)(x3)∈H1 s(Ω) in expression (9.19), we deduce that ZΩ [σ0 33∂3(x2∂3(v0 2−u0 2)) + σ0 32∂2(x2∂3(v0 2−u0 2))]dx =ZΩ f3x2∂3(v0 2−u0 2)dx +ZΓN1 h3x2∂3(v0 2−u0 2)dΓ,
106 Chapter 9. Asymptotic analysis of three-point bending tests and, consequently, ZΩ [σ0 33x2∂33(v0 2−u0 2) + σ0 32∂3(v0 2−u0 2)]dx =ZΩ f3x2∂3(v0 2−u0 2)dx +ZΓN1 h3x2∂3(v0 2−u0 2)dΓ.(9.30) Now, taking into account expression (9.24), we have that ZΩ σ0 33x2∂33(v0 2−u0 2)dx =ZL −L ∂33(v0 2−u0 2)µZω x2E(∂3u0 3−x2∂33u0 2)dx1dx2¶dx3= −ZL −L EI1∂33u0 2∂33(v0 2−u0 2)dx3.(9.31) From (9.29)-(9.31) we infer that: ZΩ σ0 32∂3(v0 2−u0 2)dx =ZL −L EI1∂33u0 2∂33(v0 2−u0 2)dx3+ ZΩ x2f3∂3(v0 2−u0 2)dx +ZΓN1 x2h3∂3(v0 2−u0 2)dΓ, ZΩ σ0 32∂3(v0 2−u0 2)dx ≥ZΩ f2(v0 2−u0 2)dx +ZΓN1 h2(v0 2−u0 2)dΓ. Finally, from expressions (9.21) and (9.22), we deduce inequality (9.26) and then (iv). Theorem 9.2.2. Under assumptions (H1²), (H2²), (H7²)and (H8²), there exists a unique solution of Problem (P0 3). Proof. Since u0 3must be an odd function, the bilinear form associated to Problem (P0 3)is coercive on H1 o(−L, L)⊂ C([−L, L]) and the uniqueness is proved. Theorem 9.2.3. Under assumptions (H1²)-(H3²),(H7²)and (H8²), there exists a unique solution of Problem (P0 2). Proof. It can be proved by applying a well-known result of Lions and Stampacchia [58]. Then, we consider V=H2 e(−L, L)which is a dense subspace of H1 e(−L, L)⊂C0([−L, L]). We denote by p0the usual norm in H1(−L, L)which gives a pre-Hilbert structure to Vand by p1the seminorm in V: p1(v) = °°°° d2v dx2 3°°°°0,(−L,L) . The subspace Y={v∈V;p1(v) = 0}is finite-dimensional since Y={v∈V;v(x3) = c, c ∈R}.
9.3. Convergence of scaled displacements and stresses 107 Moreover, due to the compact injection V ,→H1 e(−L, L), we have that (see [74]) inf y∈Yp0(v−y)≤c1p1(v),∀v∈V. (9.32) If we denote by athe bilinear form in (9.26), we obtain that ais continuous on Vand semicoercive, i.e., a(v, v)≥c2p2 1(v),∀v∈V, c2>0.(9.33) In addition, the right-hand side of equation (9.26) can be identified as an element of V0 (F0, v) = ZL −L q2vdx3−ZL −L M1∂3vdx3, and, thanks to hypothesis (H3²)for ε= 1, it is easy to prove that (F0, y)<0,∀y∈Y, y 6= 0,(0, y)∈K2(Ω).(9.34) Relations (9.32)-(9.34) allow us to apply Theorem 5.1of Lions and Stampacchia [58] and obtain the existence of, at least, one solution of Problem (P0 2). Now, let us assume that uand ˜uare two different solutions of Problem (P0 2). Then, we can prove that there exists y∈Ysuch that ˜u=u+y. Indeed, given that uand ˜uare solutions, we can write a(u, ˜u−u)≥(F0,˜u−u) = (F0, y), a(˜u, ˜u−u)≤(F0,˜u−u) = (F0, y), that is, a(˜u−u, ˜u−u)≤0. Consequently, from (9.33), y∈Yand (F0, y) = 0; so, from (9.34), y= 0. In conclusion, there exists a unique solution of Problem (P0 2). 9.3 Convergence of scaled displacements and stresses In order to complete the asymptotic analysis of the problem, we are going to carry out a convergence analysis of the solution of Problem (MV P(ε)),(u(²),σ(²))ε>0, when ²goes to 0. For all v∈H1(Ω), we denote |v|1,Ω=kvk0,Ω+kε(v)k0,Ω, which is an equivalent norm to k·k1,Ωin H1(Ω); the first addend gives a pre-Hilbert structure to H1(Ω) and the second addend is a seminorm in H1(Ω). Definition 9.3.1. We define the subspace Y(Ω) = H1 s(Ω) ∩R(Ω) Y(Ω) = {w∈H1 s(Ω); kε(w)k0,Ω= 0}={w= (0, c, 0); c∈R}.
108 Chapter 9. Asymptotic analysis of three-point bending tests Lemma 9.3.2. Let Pbe the projection of H1 s(Ω) onto Y(Ω), with respect to the pre-Hilbert structure on H1 s(Ω) defined by k·k0,Ω. Then, Pv= (0,MΩ(v2),0) ,∀v∈H1 s(Ω), MΩ(v2) = 1 meas(Ω) ZΩ v2dx. Furthermore, inf w∈Y(Ω) kv−wk0,Ω=kv−Pvk0,Ω≤˜ Ckε(v)k0,Ω,∀v∈H1 s(Ω),˜ C > 0,(9.35) ˜ Cdepending only on Ω. Proof. The projection Pis well defined because Y(Ω) is a finite-dimensional space. In addition, given that Y(Ω) is a closed subspace, we have for each v∈H1 s(Ω) (see [26]) (v−Pv,w) = 0,∀w∈Y, that is, 0 = (v−Pv,(0, c, 0)) = cZΩ (v2−(Pv)2)dx, ∀c∈R, c 6= 0. Taking into account that (Pv)2is a constant, we can deduce (Pv)2=1 meas(Ω) ZΩ v2dx =MΩ(v2). Finally, MΩ(v1) = MΩ(v3)=0since v∈H1 s(Ω) and from Poincaré-Wirtinger inequality, we obtain (9.35). From now on, we will use an auxiliary stress tensor ˜σ(²)defined by ˜σ33(²) = σ33(²),˜σ3β(²) = ²σ3β(²),˜σαβ(²) = ²2σαβ(²).(9.36) Lemma 9.3.3. The following estimates hold: k˜ σ(²)k2 0,Ω≤C1³kfk0,Ω+khk0,ΓN1´ku(²)k1,Ω,(9.37) C2kε(u(²))k2 0,Ω≤F(u(²)),(9.38) kε(u(²))k2 0,Ω≤C3³kfk0,Ω+khk0,ΓN1´ku(²)k1,Ω.(9.39) Proof. If we consider v=0and v= 2u(²)in (9.7) and τ=σ(²)in (9.6) we obtain that a0(σ(²),σ(²)) + ²2a2(σ(²),σ(²)) + ²4a4(σ(²),σ(²)) = b(σ(²),u(²)) = F(u(²)).(9.40) From (9.8)-(9.12) and (9.40) and taking into account expressions (9.36), it is easy to prove that ZΩµ1 + ν E˜σij(²)−ν E˜σpp(²)δij¶˜σij(²)dx =ZΩ fiui(²)dx +ZΓN1 hiui(²)dΓ.
9.3. Convergence of scaled displacements and stresses 109 Thanks to the properties of the Lamé coefficients and the assumptions on the volume and surface forces, we have that C4k˜ σ(²)k2 0,Ω≤ZΩµ1 + ν E˜σij(²)−ν E˜σpp(²)δij¶˜σij(²)dx = ZΩ fiui(²)dx +ZΓN1 hiui(²)dΓ≤C5³kfk0,Ω+khk0,ΓN1´ku(²)k1,Ω, and consequently, C4k˜ σ(²)k2 0,Ω≤ZΩ f·u(²)dx +ZΓN1 h·u(²)dΓ = F(u(²)),(9.41) k˜ σ(²)k2 0,Ω≤C1³kfk0,Ω+khk0,ΓN1´ku(²)k1,Ω, where C1and C4are positive constants depending on Ω,Eand ν. So, we have proved estimate (9.37). From (9.6) and taking into account the definition of bilinear forms a0,a2and a4, we obtain that ¯¯¯¯ZΩ τ:ε(u(²))dx¯¯¯¯≤C6k˜ σ(²)k0,Ωkτk0,Ω,∀τ∈X(Ω), C6>0, since 0< ² ≤1. By choosing τ=ε(u(²)), we deduce that kε(u(²))k0,Ω≤C6k˜ σ(²)k0,Ω.(9.42) From inequalities (9.41) and (9.42), we have that C2kε(u(²))k2 0,Ω≤ZΩ f·u(²)dx +ZΓN1 h·u(²)dΓ = F(u(²)), C2>0, therefore, kε(u(²))k2 0,Ω≤C3³kfk0,Ω+khk0,ΓN1´ku(²)k1,Ω, and hence estimates (9.38) and (9.39) are verified. Theorem 9.3.4. There exist constants C7and C8, independent of ², such that for all 0< ² ≤1 the solution (u(²),σ(²)) ∈(H1 s(Ω) ∩K(Ω)) ×X(Ω) of Problem (MV P(ε)) verifies ku(²)k1,Ω≤C7,(9.43) kσ33(²)k0,Ω≤C8, ²kσ3β(²)k0,Ω≤C8, ²2kσαβ(²)k0,Ω≤C8.(9.44) Proof. Let us assume that (9.43) is false, so there exists a sequence Rn→+∞as n→+∞such that for all n > 0there exists ²nsatisfying ku(²n)k1,Ω=Rn.(9.45)
116 Chapter 9. Asymptotic analysis of three-point bending tests Multiplying by ²−2and taking limits, we see that ZΩ Ψαβςαβdx = 0.(9.69) Therefore, from (9.67)-(9.69) we obtain ZΩ Ψ33ε33(v−˜ u)dx ≥ZΩ f·(v−˜u)dx +ZΓN1 h·(v−˜u)dΓ,∀v∈VBN (Ω) ∩K(Ω). Taking into account the expression for Ψ33 obtained in (9.64), we can write: ZΩ [E(∂3˜u3−x2∂33 ˜u2) + νΨµµ]ε33(v−˜ u)dx ≥ZΩ f·(v−˜u)dx +ZΓN1 h·(v−˜u)dΓ, for all v∈VBN (Ω) ∩K(Ω). Then, ZΩ E∂3˜u3ε33(v−˜ u)dx −ZΩ Ex2∂33 ˜u2ε33(v−˜ u)dx +ZΩ νΨµµε33(v−˜ u)dx ≥ ZΩ f·(v−˜u)dx +ZΓN1 h·(v−˜u)dΓ,∀v∈VBN (Ω) ∩K(Ω).(9.70) In particular, considering in equation (9.70) v∈VBN (Ω) ∩K(Ω) such that v1= 0,v2depending only on x3and even with (0, v2)∈K2(Ω) and v3= ˜u3−x2∂3v2, we obtain ZΩ Ex2 2∂33 ˜u2∂33(v2−˜u2)dx ≥ZΩ f2(v2−˜u2)dx +ZΓN1 h2(v2−˜u2)dΓ− ZΩ x2f3∂3(v2−˜u2)dx −ZΓN1 x2h3∂3(v2−˜u2)dΓ + ZΩ νx2Ψµµ∂33(v2−˜u2)dx, thanks to (9.63). Taking into account (9.59), we conclude that ZL −L EI1∂33 ˜u2∂33(v2−˜u2)dx3≥ZL −L q2(v2−˜u2)dx3−ZL −L M1∂3(v2−˜u2)dx3, for all v2∈H2 e(−L, L),(0, v2)∈K2(Ω); so, ˜u2=u0 2thanks to Theorem 9.2.3. Similarly, choosing v∈VBN (Ω) ∩K(Ω) as a test function in (9.70), such that vα= ˜uαand v3= ˜u3±v3, with v3∈H1 o(−L, L), we find that ZL −L EA∂3˜u3∂3v3dx3+ZΩ νΨµµ∂3v3dx =ZL −L q3v3dx3, and from (9.59), ˜u3=u0 3thanks to Theorem 9.2.2, which proves (9.61) for the subsequence u(²). We can repeat this procedure for each subsequence of (u(²))²>0converging weakly in H1 s(Ω); then, the uniqueness of the limit implies that the entire sequence (u(²))²>0converges weakly to u0.
9.3. Convergence of scaled displacements and stresses 117 Remark 9.3.10. From previous Theorem and Corollary 9.3.8, εαβ(u(²)) * εαβ(u0) = 0,1 ²2εαβ(u(²)) * ςαβ,(9.71) ε3β(u(²)) * ε3β(u0) = 0,1 ²ε3β(u(²)) * ς3β,(9.72) ε33(u(²)) * ε33(u0) = ς33,(9.73) weakly in L2(Ω). Theorem 9.3.11. For each 0< ² ≤1, let (u(²),σ(²))²>0∈(H1 s(Ω) ∩K(Ω)) ×X(Ω) be the solution of problem (9.6)-(9.7). Let u0∈H1 s(Ω)∩K(Ω),σ0 33 ∈L2(Ω) be the first term of the asymptotic expansion (9.13). Then, when ²→0, the following convergences are obtained: ku(²)−u0k1,Ω→0,(9.74) kσ33(²)−σ0 33k0,Ω→0,(9.75) ²kσ3βk0,Ω→0,(9.76) ²2kσαβk0,Ω→0.(9.77) Proof. Convergences (9.75)-(9.77) are obtained by applying the same methodology given in Proof of Theorem 10.4 in [80]. In order to prove (9.74), since u(²)converges weakly in H1(Ω), there exists a subsequence, still denoted u(²), such that converges strongly in [L2(Ω)]3, ku(²)−u0k0,Ω→0.(9.78) From (9.6) and (9.14), we obtain the following identity: b(τ,u(²)−u0) = a0(σ(²)−σ0,τ) + ²2a2(σ(²),τ) + ²4a4(σ(²),τ).(9.79) Moreover, for each τ∈X(Ω) |a0(σ(²)−σ0,τ) + ²2a2(σ(²),τ) + ²4a4(σ(²),τ)| ≤ C12 £kσ33(²)−σ0 33k0,Ω+²kσ3β(²)k0,Ω+²2kσαβ(²)k0,Ω¤kτk0,Ω,(9.80) taking into account that 0≤² < 1and then ²2,²4can be bounded by ²and ²2, respectively, when necessary. Choosing τ=ε(u(²)−u0)in relations (9.79) and (9.80), we have that C17kε(u(²)−u0)k2 0,Ω≤a0(σ(²)−σ0,ε(u(²)−u0))+ ²2a2(σ(²),ε(u(²)−u0)) + ²4a4(σ(²),ε(u(²)−u0)) ≤ C12 ¡kσ33(²)−σ0 33k0,Ω+²kσ3β(²)k0,Ω+²2kσαβ(²)k0,Ω¢kε(u(²)−u0)k0,Ω,
118 Chapter 9. Asymptotic analysis of three-point bending tests therefore, kε(u(²)−u0)k0,Ω→0, and taking into account (9.78) we conclude (9.74) for the subsequence u(²). This procedure may be applied to any weakly convergent subsequence of (u(²))²>0in H1(Ω); so, the whole sequence converges strongly. 9.4 Returning to domain Ω² In previous section, we have proved that the first term of the asymptotic expansion (u0,σ0)is a good approximation of the scaled displacements and stresses (u(²),σ(²)) when the cross-section parameter ²is small. Then, by applying the inverse transformations of (9.3)-(9.5), we obtain an approximation of (u²,σ²), denoted by (u0²,σ0²), such that u0² α(x²) = 1 ²u0 α(x), u0² 3(x²) = u0 3(x), σ0² αβ(x²) = ²2σ0 αβ(x), σ0² 3β(x²) = ²σ0 3β(x), σ0² 33(x²) = σ0 33(x), The following theorem gives a characterization of this approximation: Theorem 9.4.1. The terms u0²,σ0² 33 are determined as follows: (i) The displacement field u0²belongs to the subspace VBN (Ω²), where u0² 1is null, u0² 2∈H2 e(−L, L) and u0² 3=u0² 3−x² 2∂3u0² 2, u0² 3∈H1 o(−L, L).(9.81) (ii) The stress component σ0² 33 is given by σ0² 33 =E(∂3u0² 3−x² 2∂33u0² 2).(9.82) (iii) The axial displacement u0² 3is the unique solution of the problem: Problem (P0² 3): Find u0² 3∈H1 o(−L, L)such that ZL −L EA²∂3u0² 3∂3v0 3dx3=ZL −L q² 3v0 3dx3,∀v0 3∈H1 o(−L, L),(9.83) (iv) The displacement u0² 2is the unique solution of the bending problem: Problem (P0² 2): Find u0² 2∈H2 e(−L, L),(0, u0² 2)∈K2(Ω²), such that ZL −L EI² 1∂33u0² 2∂33(v0 2−u0² 2)dx3≥ZL −L q² 2(v0 2−u0² 2)dx3−ZL −L M² 1∂3(v0 2−u0² 2)dx3,(9.84) for all v0 2∈H2 e(−L, L),(0, v0 2)∈K2(Ω²).
9.5. Some properties of the limit models 119 9.5 Some properties of the limit models In the axis of the beam, we distinguish: •Γ0 N1= (−δ, δ), the part of the axis corresponding to the beam slice in which some surface forces are applied. •Γ0 C= (−l−ˆ δ, −l)∪(l, l+ˆ δ), the parts of the beam axis under sections where their boundaries intersect with Γ² C. From now on and without loss of generality, we consider that the second component of the normal vector is negative, n² 2<0, on the contact boundary Γ² C. Theorem 9.5.1. The unique solution u0² 3∈H1 o(−L, L)of Problem (P0² 3)verifies the following differential formulation: EA²∆u0² 3=q² 3,(9.85) EA²∂3u0² 3(±L) = 0.(9.86) In addition, u0² 3∈H2 o(−L, L). Proof. Firstly, choosing an odd function φ∈ D(−L, L)⊂H1(−L, L), as test function in (9.83), and applying integration by parts, we obtain −ZL −L EA²∆u0² 3φdx3=ZL −L q² 3φdx3. So, expression (9.85) is obtained in the sense of distributions. Moreover, since q² 3is smooth, u0² 3∈H2 o(−L, L). Multiplying (9.85) by v∈H1 o(L, L)and applying again Green’s theorem, we have ZL −L q² 3vdx3=ZL −L EA²∂3u0² 3∂3vdx3−EA²∂3u0² 3(L)v(L) + EA²∂3u0² 3(−L)v(−L) = ZL −L q² 3vdx3−2EA²∂3u0² 3(L)v(L), taking into account expression (9.83), so we deduce (9.86). Theorem 9.5.2. Let u0² 2∈H2 e(−L, L)be the unique solution of Problem (P0² 2)and let denote the coincidence set I={x3∈¯ Γ0 C;u0² 2(x3)=0}. Then there exists a nonnegative Radon measure µ such that EI² 1∆2u0² 2=q² 2+∂3M² 1+µin (−L, L),(9.87) ∂3∆u0² 2(±L) = ∆u0² 2(±L) = 0,(9.88) with supp(µ)⊂I.
120 Chapter 9. Asymptotic analysis of three-point bending tests Proof. Firstly, let φ∈ D(J)⊂H2(−L, L)be an even function, with J= (−L, L)\¯ Γ0 C; taking v0 2=u0² 2±φ, as test function in (9.84), and applying integration by parts, we find that ZJ EI² 1∆u0² 2∆φdx3=ZJ q² 2φdx3+ZJ ∂3M² 1φdx3. So, we obtain in the sense of distributions EI² 1∆2u0² 2=q² 2+∂3M² 1in D0(J).(9.89) Therefore, multiplying (9.89) by a test function v, taking into account that H2 e(L, L)⊂ C1(−L, L) and since q² 2+∂3M1∈L2(−L, L), applying again integration by parts, we deduce that ZJ EI² 1∆2u0² 2vdx3=ZJ EI² 1∆u0² 2∆vdx3−EI² 1∂3∆u0² 2(−L)v(−L)+ EI² 1∂3∆u0² 2(L)v(L) + EI² 1∆u0² 2(−L)∂3v(−L)−EI² 1∆u0² 2(L)∂3v(L) = ZJ q² 2vdx3+ZJ ∂3M² 1vdx3, for all v∈H2 e(−L, L)such that v= 0 in ¯ Γ0 Cand hence, using variational formulation (9.84), we deduce (9.88). Next, let x0 3∈Γ0 C\I,x0 3>0, be arbitrary and Br(±x0 3) = Br(x0 3)∪Br(−x0 3). Because (0, u0² 2)∈K2(Ω²), there exist r > 0and φ∈ D(Br(±x0 3)) such that φ(x)>0in Br(±x0 3), u0² 2≥φin H2(Br(±x0 3)). Then, for all ξ∈ D(Br/2(±x0 3)), there exists ν > 0such that u0² 2+νξ ≥φ 2in H2(Br/2(±x0 3)). Choosing v0 2=u0² 2+νξ in (9.84) and dividing by ν, we have that ZBr/2(±x0 3) EI² 1∆u0² 2∆ξdx3≥ZBr/2(±x0 3) q² 2ξdx3+ZBr/2(±x0 3) M² 1∂3ξdx3. Analogously, applying the same procedure to −ξ, we obtain the equality ZBr/2(±x0 3) EI² 1∆u0² 2∆ξdx3=ZBr/2(±x0 3) q² 2ξdx3+ZBr/2(±x0 3) M² 1∂3ξdx3, for all ξ∈ D(Br/2(±x0 3)), and then, EI² 1∆2u0² 2=q² 2+∂3M² 1in D0(Γ0 C\I).(9.90) Now, we notice that from (9.84), ZL −L EI² 1∆u0² 2∆vdx3≥ZL −L q² 2vdx3−ZL −L M² 1∂3vdx3,
9.5. Some properties of the limit models 121 for all v∈H2 e(−L, L)with v≥0in Γ0 C; in particular, a(u0² 2, v)−(f² 0, v)≥0,∀v∈H2 e(−L, L)∩H1 0(−L, L), v ≥0in Γ0 C, where ais the bilinear form a(u, v) = ZL −L EI² 1∆u∆vdx3,∀v∈H2 e(−L, L), and f² 0∈[H2 e(−L, L)]0is given by (f² 0, v) = ZL −L q² 2vdx3+ZL −L ∂3M² 1vdx3,∀v∈H2 e(−L, L)∩H1 0(−L, L). Therefore, by Riesz-Schwartz theorem (see [76]), there exists a positive Radon measure µin (−L, L) such that a(u0² 2, v0 2)−(f² 0, v0 2) = ZL −L v0 2dµ, ∀v0 2∈H2 e(−L, L)∩H1 0(−L, L), v0 2≥0in Γ0 C. Finally, from (9.89) and (9.90), supp(µ)⊂Iand we deduce (9.87). Remark 9.5.3. The problem defined by equations (9.87) and (9.88) constitutes a generalization of the Euler-Bernouilli theory for the case of beams in contact with a rigid foundation, in which becomes into a singularity in the resultant of the applied loads. Other generalizations of Bernouilli- Euler models with singularities are shown in [21] and [80]. Proposition 9.5.4. The limit vertical displacement satisfies the following regularity properties: u0² 2∈ C2(−L, L)∩C3((−L, L)\I). Moreover, the lateral limits ∂3∆u0² 2(x± 3)exist for all x3∈Iand ∂3∆u0² 2(x+ 3)≥∂3∆u0² 2(x− 3),∀x3∈I. (9.91) Proof. From previous Theorem it is clear that ∆2u0² 2∈L2((−L, L)\I). Since Iis a closed interval, u0² 2∈H4((−L, L)\I), and then, u0² 2∈ C3((−L, L)\I). Moreover, from (9.87) and given that µis a positive Radon measure, and hence it can be considered as the derivative of an increasing function ϕ(see [76]), we can write EI² 1∂3∆u0² 2=G²+ϕin (−L, L),(9.92) thanks to (9.88), where G²is the continuous function in (−L, L)given by G²(x3) = Zx3 −L (q² 2(s) + ∂3M² 1(s))ds, (9.93) and ϕis defined by: ϕ(x3) = Zx3 −L dµ. (9.94)
122 Chapter 9. Asymptotic analysis of three-point bending tests By integrating (9.92) and thanks again to (9.88), we obtain EI² 1∆u0² 2(x3) = Zx3 −L (G²+ϕ)ds, x3∈(−L, L),(9.95) where the right-side term is a continuous function in (−L, L); therefore, ∆u0² 2is continuous in (−L, L)and, consequently, u0² 2∈ C2(−L, L). Finally, we know that ϕdefined by (9.94) is a function of bounded variation in all closed interval; hence the discontinuities of ϕare only jump discontinuities (see [3]) and hence, from (9.92), the lateral limits ∂3∆u0² 2(x± 3)exist and verify (9.91) thanks to ϕis growing. Proposition 9.5.5. i) If the axial forces F3and H3are null, the axial stress component σ0² 33 is null on int(I). ii) If I⊂(−L, L), the measure µverifies: µ(I) = −ZL −L q² 2dx3,(9.96) µ=−q² 2−∂3M² 1in D0(int(I)).(9.97) iii) The jump discontinuities of ∂3∆u0² 2lie in ∂I. Proof. i) It is easy to prove that u0² 3= 0. Then, let x3∈int(I); there exists α > 0such that u0² 2(y) = 0,∀y∈Bα(x3)⊂int(I), which implies that u0² 2(y) = ∂3u0² 2(y) = ∂33u0² 2(y) = ∂333u0² 2(y) = 0,∀y∈Bα(x3).(9.98) So, thanks to (9.82), σ0² 33 is null in Bα(x3), and therefore, σ0² 33 is null in int(I). ii) Multiplying expression (9.92) by ∂3v, with v∈H2 e(−L, L), such that v≥0in Γ0 Cand integrating in (−L, L), we find that: ZL −L EI² 1∂3∆u0² 2∂3vdx3=ZL −L G²∂3vdx3+ZL −L ϕ∂3vdx3. By applying Green’s identity and the even character of v, we can write −ZL −L EI² 1∆u0² 2∆vdx3=−ZL −L q² 2vdx3+µZL −L q² 2dx3¶v(L)+ ZL −L M² 1∂3vdx3+ZL −L ϕ∂3vdx3.(9.99)
9.5. Some properties of the limit models 123 Besides, from the variational formulation (9.84) we obtain that ZL −L EI² 1∆u0² 2∆vdx3≥ZL −L q² 2vdx3−ZL −L M² 1∂3vdx3,(9.100) By term-by-term addition of (9.99) and (9.100), we have that 0≥µZL −L q² 2dx3¶v(L) + ZL −L ϕ∂3vdx3,∀v∈H2 e(−L, L); v≥0in Γ0 C.(9.101) Let us choose vsuch that supp(v)∩supp(µ) = ∅and let us consider ¯v=v−v(L)∈ H2 e(−L, L)∩H1 0(−L, L); by using the definition of derivative in the sense of distributions, ZL −L ϕ∂3vdx3=ZL −L ϕ∂3¯vdx3=−ZL −L ¯vdµ =−ZL −L vdµ +v(L)µ(I) = v(L)µ(I).(9.102) So, substituting (9.102) in expression (9.101), we deduce that v(L)µZL −L q² 2dx3+µ(I)¶≤0,∀v∈H2 e(−L, L); v≥0in Γ0 C. Then, we conclude (9.96) considering test functions vsuch that v(L) = ±1. Finally, if v∈ D(int(I)) is an even function, from (9.99) we obtain that 0 = ZI EI² 1∆u0² 2∆vdx3=ZI q² 2vdx3+ZI ∂3M² 1vdx3+ZI vdµ, and so (9.97). iii) From the demonstration of Proposition 9.5.4, we know that ∂3∆u0² 2has only jump discontinuities that are in the coincidence set I. Due to (9.98), ∂3∆u0² 2is continuous in int(I) (identically zero) and, therefore, the discontinuities must appear on ∂I. Proposition 9.5.6. Let u0² 2be a smooth solution of the differential problem (9.87)-(9.88). If u0² 2≥0in Γ0 C, then u0² 2is the solution of Problem (P0² 2). Proof. Following the same reasoning that in Proposition 9.5.5, we can write: ZL −L EI² 1∂3∆u0² 2∂3vdx3=−ZL −L q² 2vdx3+v(L)ZL −L q² 2dx3+ZL −L M² 1∂3vdx3− −ZL −L vdµ −v(L)ZL −L q² 2dx3= −ZL −L q² 2vdx3+ZL −L M² 1∂3vdx3−ZL −L vdµ, for all v∈H2 e(−L, L)such that v≥0in Γ0 C.
124 Chapter 9. Asymptotic analysis of three-point bending tests Besides, by applying integration by parts to the first term in the previous expression and taking into account (9.88), we have that ZL −L EI² 1∆u0² 2∆vdx3=ZL −L q² 2vdx3−ZL −L M² 1∂3vdx3+ZL −L vdµ ≥ZL −L q² 2vdx3−ZL −L M² 1∂3vdx3, since µis a nonnegative measure and, therefore, the theorem is proved.
Chapter 10 An improved formula for the MOR As we have mentioned in the introduction of Part II, the relative simplicity of the three-point bending tests has permitted its popularization for measuring the MOR of brittle materials; nevertheless, the fact that samples of same material can lead to different values for the MOR is something quite common. Our guess is that this effect is due to the classic theoretical formula generally used in engineering to compute the MOR which, from our point of wiev, should be completed with some three-dimensional effects. By taking advantage of all the theoretical results collected in Chapters 8 and 9, in this chapter we are going to obtain a new formula for the MOR. This new expression takes into account not only the rupture load and the total length of the beam but also the distance between the two lower cylinders, the effect of the gravity and the distance between the ends of the beam and the lower cylinders. The outline of this chapter is as follows: In Section 10.1, we introduce the classic formula for the MOR of brittle materials. Section 10.2 is devoted to recall the mathematical model associated with this experiment: an elastic problem with a frictionless unilateral contact condition. We summarize the assumptions of compatibility on the forces needed to obtain the existence of a unique solution of the associated variational problem. Several examples for laboratory tests for cylindrical and rectangular beams are detailed in Section 10.3: The Experiments I and II. We show that the previous compatibility assumptions are satisfied for these real experiments. In Section 10.4, we obtain, from the differential formulation of the bending and axial displacements obtained in Section 9.5 by using the asymptotic analysis, analytical solutions for Experiments I and II. Moreover, we introduce the new theoretical expression for the MOR. 10.1 The theoretical MOR In the classical theory of beams, the MOR is given by an explicit expression which involves the value of the load of failure, H, the distance, 2l, between the two lower cylinders, and the second moment of inertia, I1, of the transversal section of the beam (see [5]). In particular, the theoretical MOR of a cylindrical beam with radius Rsubjected to a failure 125
132 Chapter 10. An improved formula for the MOR Deriving with respect to zin (10.18) and by using (10.16) we have that wθ,zz = 0. Again, deriving with respect to rin (10.18) and by using (10.21), we obtain wθ,zr =1 r2ˆ f0 1(θ).(10.23) Deriving with respect to zin (10.17) and by using (10.20) and (10.23) we can write wθ,z =1 rˆ f0 1(θ) + f0 2(θ),(10.24) and deriving now with respect to r wθ,zr =−1 r2ˆ f0 1(θ).(10.25) Expressions (10.23) and (10.25) lead to the condition ˆ f0 1(θ) = 0 ∀θ∈[0,2π)⇒ˆ f1(θ) = ˆ f1∈R∀θ∈[0,2π).(10.26) Therefore, from (10.21), (10.22), (10.24) and (10.26) we have the expressions wθ=f0 2(θ)z+f3(r, θ),(10.27) wz=ˆ f1−f2(θ)r. (10.28) Substituting (10.20), (10.27) and (10.28) in (10.15) and (10.17), we obtain 0 = (f00 2(θ) + f2(θ))z+f1(θ) + f3,θ,(10.29) 0 = f3,r −1 rf3(r, θ) + 1 rf0 1(θ).(10.30) Deriving with respect to rin (10.29) and with respect to θin (10.30), we deduce 0 = f3,rθ, 0 = −1 rf3,θ +1 rf00 1(θ), and then f3,θ =f00 1(θ). Therefore, we have f3(r, θ) = f0 1(θ) + f4(r).(10.31) In addition, substituting (10.31) in (10.30), we obtain that f4must verify the ordinary differential equation: f0 4(r) = 1 rf4(r), whose general solution is f4(r) = f4r, f4∈R; this leads to the following expression for f3: f3(r, θ) = f0 1(θ) + f4r. (10.32)
10.3. Verifying hypotheses for two laboratory experiences 133 Finally, substituting the previous expression in (10.29), we have 0 = (f00 2(θ) + f2(θ))z+ (f00 1(θ) + f1(θ)),∀(θ, z)∈[0,2π)×(−L, L),(10.33) so we obtain the homogeneous ordinary differential equations for f1and f2: f00 1(θ) + f1(θ) = 0, f00 2(θ) + f2(θ) = 0, whose general solutions are: f1(θ) = Asin θ+Bcos θ, A, B ∈R,(10.34) f2(θ) = Csin θ+Dcos θ, C, D ∈R.(10.35) In conclusion, from (10.20), (10.27), (10.28), (10.32) and (10.34)-(10.35) we have that wverifies (iii). 10.3.3 Verification of assumption (H3) In three-point bending tests, volume and surface forces are given by the following expressions in cartesian and cylindrical coordinates: f=−ρge2=−ρg(sin θer+ cos θeθ),(10.36) h=−he2=−h(sin θer+ cos θeθ), h > 0.(10.37) We are going to prove that these forces satisfy assumption (H3) in order to obtain the existence of a unique solution of Problem (V P). Lemma 10.3.2. In experiment I, volume and surface forces verify assumption (H3). Proof. Taking into account expression (10.9) and by using (10.36) and (10.37), we have ZΩ f·wdx+ZΓN1 h·wdΓ = −2R(πLRρg + 2δ2h)A, (10.38) for every w∈R(Ω) given by (iii) of Lemma 10.3.1. Then, hypothesis (H3) is verified if we prove that A≥0, for every w∈K(Ω) ∩R(Ω) and A= 0 if and only if wn= 0 on ΓC. Let us assume that w∈K(Ω) ∩R(Ω); since n=eron ΓC, we can write wn=w·n=Asin θ+Bcos θ+z(Csin θ+Dcos θ)≤0a.e. on ΓC. Consequently: 0≥ZΓC w·ndΓ = −4ˆ δcos Ã3π 2+ˆ δ 2!RA, (10.39)
134 Chapter 10. An improved formula for the MOR and so A≥0. In order to complete the proof of (H3), from (10.39), A= 0 if and only if ZΓC w·ndΓ = 0. Since wn≤0on ΓC, we can conclude that A= 0 if and only if wn= 0 on ΓC. Lemma 10.3.3. In experiment II, volume and surface forces verify assumption (H3). Proof. In this case we consider a rectangular domain, where x1∈(−a, a),x2∈(−b, b)and x3∈(−L, L). Let w∈K(Ω) ∩R(Ω), we have that ZΩ f·wdx+ZΓN1 h·wdΓ = −(8abLρg + 4aδh)c2,(10.40) for every w∈R(Ω) given by (ii) of Lemma 10.3.1. Then, hypothesis (H3) is verified if we prove that c2≥0, for every w∈K(Ω) ∩R(Ω) and c2= 0 if and only if w∈Rn(Ω). Let us assume that w∈K(Ω) ∩R(Ω); because w·n≤0a.e. on ΓC, we have that 0≥ZΓC w·ndΓ = −4aˆ δc2,(10.41) and so c2≥0. Moreover, c2= 0 if and only if ZΓC w·ndx= 0, and, since the sign of wnis constant on ΓC, we conclude that c2= 0 if and only if wn= 0 on ΓC. 10.3.4 Existence and uniqueness of solution for experiments I and II Proposition 10.3.4. There exists a unique weak solution u∈H1 s(Ω) of Problem (P), defined by (10.3)-(10.7), corresponding to experiments I and II for a porcelain beam. Proof. Assumptions (H1) and (H6) are satisfied thanks to expressions (10.36) and (10.37). As porcelain is a homogeneous and isotropic material, with positive and constant Lamé parameters, assumption (H2) is also verified. Lemmas 10.3.2 and 10.3.3 guarantee that assumption (H3) is satisfied for experiments I and II, respectively. Assumption (H5) is also verified due to the sample geometries. Furthermore, it is easy to verify assumption (H4). Indeed, if we consider w∈K(Ω) ∩R(Ω), w·n= (c1+d2x3−d3x2)n1+ (c2+d3x1−d1x3)n2≤0on ΓC.
10.4. Obtaining a new formula for the MOR 135 Finally, since n2is non null, choosing c1=d2=d3= 0 and c2, d16= 0, we obtain that there exists w∈(K(Ω) ∩R(Ω)) \Rn(Ω),w6=0, such that w·n= (c2−d1x3)n26= 0. So the result is obtained from Theorem 10.2.1. 10.4 Obtaining a new formula for the MOR In Section 9.4 we have obtained the one-dimensional variational problems associated to the bending and axial displacements in three-point bending tests by applying the well-known method of asymptotic expansions. We have proved that a good approximation of the solution of the threedimensional Problem (V P)is the limit displacement u0= (u0 1, u0 2, u0 3), where the transversal displacement, u0 1, is null, the bending displacement, u0 2, only depends on x3, and the axial displacement, u0 3, is given by: u0 3=u0 3−x2∂3u0 2. In particular, in the scenario of experiments I and II introduced in Section 10.3, the limit bending and axial displacements are the solution of the following variational problems: Problem (P0 2): Find u0 2∈H2 e(−L, L), such that u0 2≥0in Γ0 C= (−l−ˆ δ, −l)∪(l, l +ˆ δ)and ZL −L EI1∆u0 2∆(v0 2−u0 2)dx3≥ZL −L q2(v0 2−u0 2)dx3,(10.42) for all v0 2∈H2 e(−L, L)such that v0 2≥0in Γ0 C. Here, q2is the total vertical force per unit of length: q2(x3) = −ρgβ −χ(−δ,δ)(x3)αh, where: β=½πR2for experiment I, 4ab for experiment II, α=½2δR for experiment I, 2afor experiment II. Problem (P0 3): Find u0 3∈H1 o(−L, L), such that ZL −L ∂3u0 3∂3v0 3dx3= 0,∀v0 3∈H1 o(−L, L).(10.43) Moreover, in Section 9.5, we establish some properties of the bending and axial displacements which are collected in the following propositions (see [69]). Proposition 10.4.1. In the scenario of experiments I and II, the bending and axial problems (P0 2) and (P0 3)have unique solutions. Furthermore, u0 3is null and the axial normal stress is given by the expression σ0 33(x) = −Ex2∆u0 2(x3),∀x∈Ω.(10.44)
136 Chapter 10. An improved formula for the MOR Proposition 10.4.2. Let denote the coincidence set I={x3∈¯ Γ0 C;u0 2(x3) = 0}. i) The bending displacement u0 2satisfies the following differential problem: EI1∆2u0 2=q2+µin (−L, L),(10.45) ∂3∆u0 2(±L) = ∆u0 2(±L) = 0,(10.46) where µis a nonnegative Radon measure µsuch that supp(µ)⊂I. ii) u0 2verifies the regularity properties: u0 2∈ C2(−L, L)∩C3((−L, L)\I). iii) The lateral limits ∂3∆u0 2(x± 3)exist for all x3∈Iand ∂3∆u0 2(x+ 3)≥∂3∆u0 2(x− 3),∀x3∈I. iv) The measure µverifies: µ(I) = −ZL −L q2dx3. v) The reciprocal of i) is also true if u0 2≥0in Γ0 C. In the particular case of the experiments I and II, we can demonstrate that the effective support on each lower cylinder will be a point or an interval. Proposition 10.4.3. Let u0 2be the solution of Problem (P0 2). If x0, x1∈Iwith 0< x0< x1, therefore [−x1,−x0]∪[x0, x1]⊂I. Proof. The proof is very close to that given in [29]. From now on, we denote ˜ I= [−x1,−x0]∪ [x0, x1]. We define the following function: v=u0 2(1 −χ˜ I) which satisfies v∈H2 e(−L, L)as a straightforward consequence of the fact that u0 2∈H2 e(−L, L)⊂ C1(−L, L)and u0 2≥0on Γ0 C. Moreover, v≥0in Γ0 C. We consider w=u0 2−v∈H2 e(−L, L)that verifies w≥0on ˜ Iand w= 0 on (−L, L)\˜ I. Then, choosing v0 2=w+u0 2as test function in expression (10.42), we have that ZL −L EI1∆u0 2∆wdx3≥ZL −L q2wdx3, and from the definition of w, Z˜ I EI1∆u0 2∆u0 2dx3≥Z˜ I q2u0 2dx3.
10.4. Obtaining a new formula for the MOR 137 Now, choosing v0 2=vas test function in expression (10.42), we have that ZL −L EI1∆u0 2∆(v−u0 2)dx3≥ZL −L q2(v−u0 2)dx3, and taking into account the definition of v, Z˜ I EI1∆u0 2∆u0 2dx3≤Z˜ I q2u0 2dx3. Therefore, we deduce that 0≤Z˜ I EI1∆u0 2∆u0 2dx3=Z˜ I q2u0 2dx3≤0,(10.47) since q2is negative and u0 2positive. Relation (10.47) implies that ∆u0 2= 0 in ˜ Iand as u0 2= 0 in ∂˜ I, then u0 2= 0 on ˜ Iand we conclude the proof. Corollary 10.4.4. The coincidence set Iis either two isolated points or two closed intervals. 10.4.1 Analytical solution. A new formula for MOR It is easy to obtain an analytical solution of differential problem (10.45)-(10.46): Proposition 10.4.5. An analytical solution of differential problem (10.45)-(10.46) is given by: u0 2(x3) = 1 24EI1 −(αh +ρgβ)x4 3+b3x2 3+d3in (−δ, δ), −ρgβ(x4 3−l4) + b2(x3 3−l3) + c2(x2 3−l2) + d2(x3−l)in [δ, l), −ρgβ(x4 3−l4) + b1(x3 3−l3) + c1(x2 3−l2) + d1(x3−l)in [l, L), and by its even extension in [−L, −δ),where b1= 4ρgβL, b2=−4αhδ, b3= 6ρgβL(2l−L)−6αhδ(δ−2l), c1=−6ρgβL2, c2= 6ρgβL(2l−L) + 12αhδl, d1= 12ρgβLl2−4αhδ(δ2−3l2), d2=−4αhδ3, d3=ρgβl2¡l2−6L(2l−L)¢−αhδ(8l3+δ3−4δ2l). The Radon measure µis given by: µ= (αhδ +ρgβL)(δ{x3=−l}+δ{x3=l}),(10.48) and I= supp(µ) = {−l, l}. Moreover, the axial normal stress σ0 33 on each transversal section of the beam is expressed by: σ0 33(x) = x2 2I1 ρgβ(x2 3+L2−2Ll) + αh(x2 3+δ2−2δl),if x3∈(−δ, δ), ρgβ(x2 3+L2−2Ll)+2αhδ(x3−l),if x3∈(δ, l) ρgβ(L−x3)2,if x3∈(l, L). (10.49)
138 Chapter 10. An improved formula for the MOR Figure 10.2: Analytical solution of Problem (10.45)-(10.46) and its derivatives. Proof. It is easy to check that u0 2is a solution of equation (10.45) with µgiven by expression (10.48). So, the proof is a straightforward consequence from Proposition 10.4.2. Figures 10.2 and 10.3 show a graphical representation of the analytical solution of Problem (10.45)-(10.46) for the experiment II with the following data: 2L= 0.118m,2a= 2b= 0.0085m,2l= 0.06m, δ = 0.0025m, E= 6.9975 ×1010N/m2, ρ = 2340kg/m3, h = 5649412N/m2. These values will be used both in the experiments in laboratory and in the numerical simulations. Notice that the second derivative of u0 2is still continuous while in the third derivative there exist two discontinuities in the support points, ±l, which conform the coincidence set I. Corollary 10.4.6. i) If the following condition is satisfied: ρgβl ¡3(L−l)2−2l2¢−αhδ(3l2−δ2)<0,(10.50) then, there exists ˆ δ0>0such that for all ˆ δ < ˆ δ0the analytical solution given in Proposition 10.4.5 is the unique solution of variational problem (P0 2). ii) If condition (10.50) is satisfied and ˆ δverifies ρgβ h(2l+ˆ δ)((l+ˆ δ)2+l2+ 6L2)−4L((2l+ˆ δ)(l+ˆ δ)+4l2)i−4αhδ(3l2−δ2)<0.(10.51) then, u0 2is the unique solution of variational problem (P0 2).
10.4. Obtaining a new formula for the MOR 139 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 −0.5 0 0.5 1 1.5 2 2.5 3 3.5 x 107Maximum traction normal stress Figure 10.3: Axial normal stress obtained in the fiber x2=−b. Proof. i) Condition (10.50) guarantees that ∂3u0 2(l)>0, so there exists ˆ δ1>0such that u0 2 is increasing in (l, l +ˆ δ1). Moreover, since u0 2(l) = 0, we deduce that u0 2is positive in the interval (l, l +ˆ δ), with ˆ δ < ˆ δ0= min(ˆ δ1, L −l). A similar argument can be used at x3=−l. Therefore, u0 2is the solution of Problem (P0 2). ii) Condition (10.51) guarantees that the solution is positive at x3=l+ˆ δ. Thanks to condition (10.50) and taking into account that u0 2is a concave function in the interval (l, L), we can state that u0 2is positive in the interval (l, l +ˆ δ). From expression (10.49), we can deduce a new formula for the MOR. Recall that the MOR corresponds with the maximum surface stress of the bent beam at the instant of failure. As we will see, the maximum stress occurs in the opposite zone to the load application zone, which coincides with the results found in the classic literature. Proposition 10.4.7. If 2l > L, then the maximum fiber normal stresses are achieved at (0,±d, 0) and their values are: σmt =d 2I1 [αhδ (2l−δ)−ρgβL(L−2l)] ,(10.52) for the maximum in traction and σmc =−σmt, for the maximum in compression, where dis the farthest distance in the x2-direction from the axis of the beam. Proof. Let us suppose that x2is fixed and let us define f(x3) = σ0 33(x2, x3),
140 Chapter 10. An improved formula for the MOR with σ0 33 defined by (10.49). Then, df dx3 (x3) = x2 I1 x3(ρgβ +αh),if x3∈(−δ, δ), ρgβx3+αhδ, if x3∈[δ, l) −ρgβ(L−x3),if x3∈(l, L], In order to study the extrema of fwe have to analyze its critical points x3= 0,x3=±Land also the jump discontinuities of its derivative x3=±l. So, we can distinguish two situations: i) If x2<0, we have that x3= 0 is a local maximum that verifies: f(0) >0 = f(±L), if 2l > L. Furthermore, f(±l)≤0, so the maximum traction normal stress is achieved at (0,−d, 0). ii) Analogously if x2>0, we have that x3= 0 is a local minimum that verifies: f(0) <0 = f(±L), if 2l > L. Furthermore, f(±l)≥0, so the maximum compression normal stress is achieved at (0, d, 0). Corollary 10.4.8. The MOR for the experiment I is given by σ0 ft =H πR3(2l−δ)−F πR3(L−2l).(10.53) where Fand Hare the modulus of the total gravity forces and the load of failure, respectively. Analogously, the MOR for the experiment II is given by σ0 ft =3H 16ab2(2l−δ)−3F 16ab2(L−2l).(10.54) Proof. From expression (10.52) and taking into account that I1=πR4 4, I1=4 3ab3, for experiments I and II, respectively, we have that σmt =2 πR3[αhδ (2l−δ)−ρgβL(L−2l)] for experiment I, σmt =3 8ab2[αhδ (2l−δ)−ρgβL(L−2l)] for experiment II. Then, since the total gravity forces and the load of failure are given by F= 2Lρgβ, H = 2δαh, we conclude the proof. Notice that if the total gravity forces, F, can be neglected and the upper cylinder exerts really a punctual force, that is δ= 0, expressions (10.53) and (10.54) are in line with classic formulas (10.1) and (10.2).
Chapter 11 Computing the MOR In this chapter we are going to show the results from numerical simulations of experiments I and II described in Section 10.3. These experiments have been carried out in collaboration with the Institute of Ceramics of the Universidade de Santiago de Compostela. In this section, the methodology to compute the MOR is as follows: •Firstly, we have carried out various laboratory tests on porcelain beams in order to determine the experimental force of rupture, H. •Secondly, we calculate a theoretical value for the modulus of rupture, σft, by using the classic formulas of the MOR (10.1) and (10.2). •Thirdly, considering the experimental force of rupture, H, we compute numerically the corresponding three-dimensional stresses and so, we compute its 3d numerical MOR, σfn. •Fourthly, we compute numerically the displacements and stresses by using the one-dimensional model and, therefore, we compute its 1d numerical MOR, σ0 fn. •Finally, we also compute a new value for the MOR, the 1d analytical MOR, denoted by σ0 ft, by using expressions (10.53) and (10.54) obtained in Section 10.4. These all experiences are repeated for different sizes of samples and for different distances between the two lower cylinders. In order to carrying out the three-dimensional numerical simulations that we present below we use the finite element software package MSC.MARC1, that let us make simulations in solids mechanics when there exist contact conditions with a foundation or a deformable body. We show below figures corresponding to the computed stresses and displacements, and tables with computed values for the MOR, by using hexahedrical meshes with a characteristic size of e= 0.0025 m. Besides, for the one-dimensional numerical simulation, a MATLAB code has been developed which uses a generalized Newton technique to deal with the nonlinearity due to the contact conditions (see [7]). 1MSC.MARC is a registered trademark of MSC.Software Corporation. 141