Numerical verification of elastic asymptotic solution in cracks in adhesive interfaces and bridged cracks

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Numerical verification of elastic asymptotic solution in cracks in adhesive interfaces and bridged cracks Autor: Eduardo Gómez Pérez Tutor: Vladislav Mantič Leščišin - Sara Jiménez Alfaro Dpto. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Numerical verification of elastic asymptotic solution in cracks in adhesive interfaces and bridged cracks Autor: Eduardo Gómez Pérez Tutor: Vladislav Mantič Leščišin - Sara Jiménez Alfaro Dpto. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024 Trabajo Fin de Grado: Numerical verification of elastic asymptotic solution in cracks in adhesive interfaces and bridged cracks Autor: Eduardo Gómez Pérez Tutor: Vladislav Mantič Leščišin - Sara Jiménez Alfaro El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos Amí familia por darme la posibilidad de estudiar y su apoyo en estos cuatro largos años. Eduardo Gómez Pérez Sevilla, 2024 I Resumen L a industria aeronáutica recientemente ha comenzado a implementar el uso de materiales compuestos en el desarrollo de las distintas aeronaves de nueva generación que se pretenden lanzar al mercado. El estudio del principal elemento que permite unir distintas láminas de material compuesto (el material adhesivo) es fundamental para entender y predecir el comportamiento de estas estructuras. Pero incluso manteniendo la presencia de materiales homogéneos como el tradicional aluminio, la industria está especialmente enfocada en la reducción de pesos. Esta reducción de pesos puede hacerse sustituyendo los tradicionales remaches por componentes con interfaz adhesiva. En este contexto, se desarrollará el análisis por el Método de los Elementos Finitos de cuatro distintos problemas de fractura con adhesivo, verificando las soluciones analíticas propuestas anteriormente mediante el software de Elementos Finitos FEniCS. III 2Capítulo 1. Introduction composed of steel and aluminum. Consequently, we can assert that the industry is notably seeking any way to reduce weight, and that this is one of the trends that will shape the future. Therefore, any contribution that can be made in terms of weight reduction, such as the study of adhesive interfaces, is positive for the industry. As mentioned above, even if we maintain the presence of homogeneous materials such as aluminum, the removal of the hundreds of rivets needed in each aircraft can represent a significant weight reduction, in addition to the elastic advantage of avoiding stress concentrations that occur at these holes through which we introduce the rivets. The primary alternative to riveting is the use of adhesive materials (adhesive interfaces). These adhesive interfaces, like continuous materials, can experience the presence of cracks generated, for example, by cyclic fatigue loads. Therefore, it is crucial that we understand the elastic behavior of these adhesives: more specifically, it is important that we understand how the elastic problem behaves in the case of a crack growing in this adhesive, as well as studying how this adhesive may be reinforcing a crack that has already appeared in an initially homogeneous material. Additionally, there is a clear trend towards reducing the number of experiments for the verification of new materials and components in aircraft design: an experimental verification can be costly as well as sophisticated. In this context, virtual certification has become one of the major objectives in recent years. 1.2 Objectives Numerical simulations are becoming more and more important for the industry. Hence, it is very important to develop numerical methodologies that provide reliable results with low computational complexity. Among others, one of the most costly simulations is the numerical study of adhesive joints around the crack tip region, which require very refined meshes. In 2024, a special singular finite element was proposed at University of Seville to reduce the computational complexity of these simulations in antiplane elasticity problem. This finite element was developed by Mantič, Vázquez-Sánchez, Romero-Laborda, Muñoz-Reja, Jiménez-Alfaro, Távara [6]. This element is defined from the analytical solution developed in 2020 at University of Seville by Jiménez-Alfaro, Mantič [7] and specifically adapted in 2023 for mode III cracks by Jiménez-Alfaro, Mantič [8], being numerically verified in Jiménez-Alfaro, Mantič [2]. In the future, this element is intended to be developed for plane elasticity, and thus the development of analytical solutions and their corresponding verification is essential. In the context of fracture mechanics, elastic solutions of bridged cracks were developed by Sara Jiménez and Vladislav Mantič for different modes of fracture (this previous work can be visualized in the original article by Jiménez-Alfaro, Mantič [2]). In addition to these bridged cracks (which we will describe later as the Dirichlet-Robin problem), solutions were also obtained for the problem of a crack propagating in an adhesive interface (the Robin-Neumann problem). Notice that these problems are radically different: the Dirichlet-Robin problem (or D-R) corresponds to what we know as bridged cracks, that is, where an adhesive element reinforces a specific crack that has formed in a homogeneous material. The Robin-Neumann problem (or R-N) corresponds to a crack propagating within an adhesive interface. The region that the crack has already traversed is completely separated, without any reinforcement or adhesive element, while in the region yet to be traversed, we observe the typical behavior of a continuous distribution of springs. The main objective of this project is the verification of the analytical solutions developed for the R-N and D-R problem for different modes of fracture. To that aim, the numerical solution of a benchmark problem is obtained using the Finite Element Method. If simulations and analytical solutions were to converge, the analytical solutions could be confirmed as correct, enhancing our understanding of the elastic behaviour of these two distinct situations (the bridged cracks problem, or D-R, and the R-N problem, where a crack is propagating in an adhesive interface). 1.3 Fracture modes 3 This numerical simulation involves simulating the area nearby the crack tip with the Finite Element Method (FEM) by implementing the solution in an open-source software. The simulations will be run in FEniCS, Finite Element Computational Software. FEniCS is an open-source computing platform for solving partial differential equations with the finite element method. Scientific models can be translated into efficient finite element codes programmed in Python and C++ interfaces. A more specific development of the FEniCS software can be found on its website, with tutorials such as Dokken [3] or Bleyer [4]. The initial phase of this study involved acquiring proficiency in programming with FEniCS and the comprehension of the application of the Finite Element Method to derive elastic solutions. Due to constraints on document length, detailed exposition of this knowledge will be omitted herein. Now, the fundamental elements of the elastic problem to be solved will be explained, focusing on the fracture problem, as well as an initial approach to the implementation of the finite element method. 1.3 Fracture modes Classical elastical fracture mechanics literature has detailed crack propagation under three distinct scenarios, commonly referred to as the three modes of fracture. Figure 1.2: Three modes of fracture, Jiménez-Alfaro, Mantič [2]. Mode I, depicted in Figure 1.2(a), is also known as symmetric or opening mode. Displacements occur perpendicular to the fracture surface. Mode II, illustrated in Figure 1.2(b), is recognized as shear mode. Displacements occur in the same direction as crack propagation. Lastly, Mode III, commonly referred to as the anti-plane problem of fracture, involves displacements perpendicular to the domain. Notice that, for the two scenarios (the R-N and the D-R problem), the three Modes of fracture must be studied. This means six problems must be solved in order to fully understand the behaviour of these elastic problems. Unfortunately, we will only study four of the six explained problems: the Robin-Neumann problem, corresponding to the crack propagating in an adhesive, has been extensively studied in its three fracture Modes. However, the Dirichlet-Robin problem, corresponding to what we know as a bridged crack, has only been studied for Mode III of fracture, as the analytical solution for the D-R problem corresponding to Modes II and I has not been fully developed yet. Lastly, a more comprehensive description of the analytical solution of these different scenarios can be found in the already mentioned work developed by Jiménez-Alfaro, Mantič [2]. Now we will provide a brief summary of the elastic problem to be solved, following the approach of the authors in the aforementioned work. 4Capítulo 1. Introduction 1.4 Elastic problem A first approach to the elastic problem, in general terms, can be found in París [5]. The elastic problem can be defined in terms of stress and strain tensors that are defined in a certain domain, Ω . Notice that these equations describe the elastic problem for an isotropic and linear-elastic material with indicial notation: σi j,j+fi=0in Ω(1.1) εi j =1 2(∂ui ∂xj +∂uj ∂xi )in Ω(1.2) σi j =2Gεi j +λεkkδi j in Ω(1.3) εi j =1+ν Eσi j −ν Eσkkδi j in Ω(1.4) u=¯u on ΓD(1.5) t=σn=¯ t on ΓN(1.6) where G is the shear modulus, λ is the Lamé Constant, E is the Young’s modulus, ν is the Poisson’s ratio and f are the body forces. u represents the displacements field. σ stands for the stress tensor, while ε represents the strain tensor. (1.1) represents the equilibrium equation. (1.2) and (1.3) are the Lamé and Hooke equations. (1.4) represents the compatibility equation. A displacement boundary condition can be applied as described in the case of the Dirichlet boundary condition, as shown in (1.5), on the Dirichlet boundaries of the domain. In the case of a Neumann-type boundary condition, we must apply (1.6) on the Neumann boundary. A more specific boundary condition, the Robin boundary condition, can relate stress and displacements with the spring equation. This will be detailed later. The elastic problem, to be fully described, must be accompanied by boundary conditions on displacements and stresses. This work will be focused in the plane and anti-plane elastic problems, which can be reduced to simpler problems: In the scenario of a plane elastic problem, instead of the usual 15 unknowns (6 components of the stress tensor, 6 components of the strain tensor and the 3 displacements) we can reduce the problem to 8 unknowns: 2 displacements, 3 components of the strain tensor and 3 components of the stress tensor. If we assume that body forces are zero, the problem can be reduced to ∇4φ(x,y) = 0 , where φ is called Airy function. The Airy function is called a biharmonic function. When this function is obtained, the stress tensor and the displacements field can be obtained. Additionally, we must differentiate between two distinct problems: the case of plane strain and the plane stress case. The problem of plane strain refers to a scenario where a structural component experiences deformations predominantly in two dimensions while remaining geometrically thin in the third dimension. This typically occurs in structures such as plates or shells undergoing bending or stretching, where one dimension is much smaller compared to the other two. On the other hand, the problem of plane stress involves analyzing a structural element where stresses predominantly act within a single plane, with negligible stress components in the perpendicular direction. In this project, we will work with plane strain unless otherwise specified. In the scenario of an anti-plane elastic problem, only one displacement must be considered (only the displacement perpendicular to the domain will be different from zero, as the in-plane displacements are equally zero). If we assume that body forces are zero, the problem can be reduced to a Laplace equation: ∇2Uz(x,y) = 0 . Only σxz and σyz components of the stress tensor will be non-zero. 1.5 Finite Element Method 5 The elastic problem, defined in these terms, will be solved by applying the Finite Element Method: a method that allows obtaining an approximate solution of the actual displacements over a certain domain, through the definition of the partial differential equations that govern the problem. The domain is discretized through nodes and elements. An element is an abstract entity that connects two or more nodes. The different nodes and their elements form what we know as the mesh of the domain. The partial differential equations will be applied to these nodes, and the solution will be obtained. Thus, we can obtain information about the elastic behaviour of the problem by implementing this method. 1.5 Finite Element Method The original formulation of the Virtual Work Theorem can be expressed as ZΩ σ:δε dV =ZΩ f·δudV +ZΓN ¯ t·δudS,∀δu∈V(1.7) where σ stands for the stress tensor, δε stands for the strain tensor of a virtual displacements solution, f represents body forces, g represents an applied load and δu represents the required displacements. Vis the space of all admissible displacements for the problem. In this equation, we can define: σ=λtr(ε)I+2Gε,(1.8) λ∗=2λG λ+2G(1.9) notice that constant λ can also be defined in terms of plane stress if needed. If we want to work under plane stress condition, λ must be replaced by λ∗ . (1.8) is the Lamé equation. I stands for an identity matrix, G is the shear modulus and tr stands for the trace of the tensor. ε stands for the strain tensor. The Finite Element Method (FEM) is a powerful computational technique used to solve complex problems in engineering and physics, particularly those involving structural analysis, heat transfer, and fluid dynamics. At its core, FEM involves discretizing a large, complex domain into smaller, manageable subdomains known as elements, which are connected at discrete points called nodes. The primary objective is to approximate the unknown field variables, such as displacements, within each element using interpolating functions, often polynomials, that are defined in terms of nodal values. These interpolating functions, also known as shape functions, are constructed to satisfy the governing equations and boundary conditions of the problem within each element. For instance, consider the approximation of displacement uwithin an element: ue(x) = n ∑ i=1 Ni(x)ui(1.10) where Ni(x) are the shape functions, ui are the nodal displacements, and n is the number of nodes in the element. The formulation leads to a system of equations derived from the principle of virtual work or from Galerkin’s method. The stiffness matrix Kand the force vector Fare assembled as follows: K= m ∑ e=1 Ke,F= m ∑ e=1 Fe(1.11) 6Capítulo 1. Introduction where Ke and Fe are the element stiffness matrices and force vectors, respectively, and m is the number of elements. The global system of equations can be written as: Ku =F(1.12) Solving this system provides the nodal displacements u , which in turn give an approximate solution to the original problem. This discretization allows for handling complex geometries and material properties, making FEM an indispensable tool in modern engineering analysis. 2 Problem definition 2.1 Domain of work and boundary conditions T he defined domain for our problem contains various boundaries, with the crack tip located in the origin of the Cartesian coordinates. This domain will be very useful for us to verify our solution: this is a benchmark problem, where the analytical solution is imposed on contourns 1,2 and 3, and the Robin, Neumann and Dirichlet boundary conditions are imposed on ΓR , ΓN and ΓD . If the analytical solution is really representing the Dirichlet-Robin or Robin-Neumann problem, the numerical solution should be the same as the analytical solution in the domain and the Robin, Dirichlet and Neumann contours. For this reason, the error between the analytical and numerical solutions will also be measured. For now, we will focus solely on the top of the domain, as the problem exhibits symmetry with respect to the x-axis. Figure 2.1.1: Domain for the benchmark problem, Robin-Neumann. 7 8Capítulo 2. Problem definition Figure 2.1.2: Domain for the benchmark problem, Dirichlet-Robin. In Figures 2.1.1 and 2.1.2 we can visualize the domain of our problem, as described in Jiménez- Alfaro, Mantič [2]. Boundaries labeled as Γ1 D , Γ2 D , and Γ3 D will be subjected to different nonhomogeneus Dirichlet boundary conditions, corresponding to the application of the analytical solution that is being verified. The boundary identified as ΓR will correspond to a Robin boundary condition. Finally, the boundary marked as ΓD or the boundary labeled as ΓN will be subject to different homogeneous Dirichlet or Neumann condition, depending on the specific simulation. The Dirichlet boundary condition can be expressed as a function of (x,y). Displacements are imposed for these parts of the boundary. The Robin boundary condition, where the springs are located, has been described in the variational formulation of the problem. The Neumann boundary condition can be described as stress-free condition. 2.2 Implementation of the Virtual Work Theorem We will now focus on how to apply the general Virtual Work Theorem equation in our problem: for every boundary, a different boundary condition will be applied as it was shown in Figures 2.1.1 and 2.1.2. In terms of mathematical expressions, this can be visualized: −∇(∇(u)) = f in Ω(2.1) u=uDon ΓD(2.2) −∂u ∂n=g on ΓN(2.3) −∂u ∂n=κu on ΓR(2.4) where κ stand for the stiffness of the spring, Ω stands for the domain and Γ stands for the different boundaries: ΓN represents the Neumann boundary, ΓR represents the Robin boundary and ΓD represents the Dirichlet boundary. u represents the displacements field. Notice that (2.3) is derived from (1.6). (2.4) represents a general Robin boundary condition. (2.2) is analogue to equation (1.5). 2.2 Implementation of the Virtual Work Theorem 9 We now have to multiply by a test function, v, and integrate by parts. −ZΩ ∇(∇(u))vdx =ZΩ ∇(v)∇(u)dx −ZδΩ ∂u ∂nvds (2.5) The Dirichlet boundary vanishes from the boundary integral since v=0 (Dirichlet boundaries will be subjected to a specific boundary condition described by a function. Therefore, the virtual displacement function must be equal to zero in the Dirichlet boundaries). Contributions from Neumann boundaries and Robin boundaries are described inserting the boundary conditions: −ZδΩ ∂u ∂nvds = n ∑ i=1ZΓN gids + n ∑ i=1ZΓR κiuds (2.6) F(u,v) = ZΩ ∇(u)∇(v)ds + n ∑ i=1ZΓN givds + n ∑ i=1ZΓR κiuds −ZΩ f vdx =0(2.7) Notice that, in (2.6), since we can have multiple Neumann and Robin boundaries, we must evaluate the integral on each of the different boundaries i and sum their contributions. We now separate terms: a(u,v) will contain all members that are a function of the test function, v , and the desired displacements function, u , also known as the trial function. L(v) will only be a function of the test function, v. Therefore, to obtain a solution, we must force a(u,v) = L(v): a(u,v) = ZΩ ∇(v)(∇(u))dx + n ∑ i=1ZΓR κiuvds (2.8) L(v) = ZΩ f vdx − n ∑ i=1ZΓN givds (2.9) 2.2.1 Anti-plane elasticity implementation We previously identified Mode III fracture as the anti-plane problem of fracture, where displacements are perpendicular to the domain. The implementation of the Virtual Work Theorem varies slightly between Mode III and Modes II and I. In Mode III fracture, only one displacement needs to be evaluated in each simulation: the displacement perpendicular to the domain. In the case of Modes II and I, two displacements need to be evaluated. Therefore, the displacement in Mode III can be described with a scalar, whereas the displacements in Modes II and I have to be described with a vector. In the case of Mode III, we can define the implementation as follows: a(u,v) = ZΩ ∇(u)∇(u)dx +ZΓR κ·u·vds (2.10) L(v) = ZΩ f·vdx −ZΓN g·vds (2.11) In a Python - FEniCS environment, this implementation can be achieved with the following lines. Notice that scalars g , f and κ must have been defined before, as well as the vectors to express integrals at subdomains dx and dx . The Robin term must be adapted to our specific simulation, where Rstands for the label of the Robin boundary: 10 Capítulo 2. Problem definition dx = Measure(’dx’, domain=mesh, subdomain_data = Boundary_markers) ds = Measure(’dx’, domain=mesh, subdomain_data = Boundary_markers) #In this implementation, r stands for Kappa a = dot(grad(u),grad(v))*dx + r*u*v*ds(R) L = f*v*dx - g*v*ds This implementation was achieved following Dokken [3]. 2.2.2 Plane elasticity implementation For Modes II and I of fracture, it is easier to describe the variational formulation in terms of stress and deformation. We firstly describe stress and deformation as shown: σ=λtr(ε)I+2Gε,(2.12) ε=sym(∇(u)) (2.13) Notice that sym stands for the symmetric part of the gradient of function u . Scalars κ,g and f must now be described as vectors, since this problem cannot be reduced to a single displacement. Now, this expressions can be introduced into the general formulation of the Virtual Work Theorem, as described in (1.7). a(u,v)and L(v)can now be specified: a(u,v) = ZΩ σ:εdx +ZΓR κ·u·vds (2.14) L(v) = ZΩ f·vdx −ZΓN g·vds (2.15) In a Python - FEniCS environment, this implementation can be achieved with the following lines. Once again, R stands for the label of the Robin boundary: dx = Measure(’dx’, domain=mesh, subdomain_data = Boundary_markers) ds = Measure(’dx’, domain=mesh, subdomain_data = Boundary_markers) def eps(u): return sym(grad(u)) def sigma(u): return lmbda*tr(eps(u))*Identity(2) + 2.0*mu*eps(u) #In this implementation, r stands for Kappa a = inner(sigma(u), eps(v))*dx + (r[0]*u[0]*v[0] + r[1]*u[1]*v[1])*ds(R) L = dot(f, v)*dx - inner(g,v)*ds This implementation was achieved following Bleyer [4]. 2.3 Transformation of coordinate system Notice that the analytical solutions that we want to verify are described in the polar coordinates system (r,θ). Unfortunately, FEniCS software requires our expressions to be described in the (x,y) coordinates system in order to properly implement the simulations. Therefore, a brief example of how to transform a vector from the (r,θ) to the (x,y) coordinates system is shown: 2.4 Dirichlet-Robin or Robin-Neumann 11 x=r·cos(θ)(2.16) y=r·sin(θ)(2.17) r=px2+y2(2.18) θ=arctany x(2.19) L(r,θ→x,y)=cosθ−sinθ sinθcosθ(2.20) u(r,θ)=ur uθ(2.21) u(x,y)=L(r,θ→x,y)·u(r,θ)(2.22) However, if our boundary conditions are expressed in the Cartesian coordinate system (x, y), the calculation of stresses will be referred to this same reference system. Nevertheless, since the problem can also be described in the reference system (r, θ ) and this system provides us with more information about the elastic behavior of the crack (as we can beforehand determine some components of the stress tensor), we must transform the stress tensor from the Cartesian coordinate system to the polar coordinate system. Now, in order to calculate stresses, we will follow the usual implementation described in París [5]: L(x,y→r,θ)=cosθsinθ −sinθcosθ(2.23) σ1=σxx σxy σxy σyy(2.24) σ2=σrr σrθ σrθσθθ (2.25) σ2=Lσ1LT(2.26) Therefore, all the important stresses will be obtained as shown in the (r,θ) coordinates system. 2.4 Dirichlet-Robin or Robin-Neumann Considering this domain and the elastic problem, various issues regarding bridged cracks and adhesive interfaces can be addressed. Specifically, two distinct problems will be examined in this study, following the approach outlined in Jiménez-Alfaro, Mantič [2] work. The first problem under investigation will involve the application of Dirichlet-Robin boundary conditions to the boundaries lying on y=0 . Conversely, for the second problem, Robin-Neumann boundary conditions will be employed for the boundaries located at y=0 . These represent two physically distinct scenarios: 1) The Dirichlet-Robin problem can be associated with a bridged crack. In this case, the Robin region, corresponding to the damaged area, i.e. a crack, is reinforced by particles, fibers or inclusions, while the Dirichlet region, representing the undamaged portion, is subject to perfect joining conditions. Refer to Figure 2.2 for a schematic representation. 18 Capítulo 3. Anti-plane elasticity, Mode III Figure 3.4: Computed error. Notice that the maximum error doesn’t exceed 0.5%, which is considered low. Only two terms of the analytical solution were used for this analysis. This means the analytical analysis was very robust and precise, and there’s no need to implement more terms. The implementation can be visualized in the codes provided at the programming appendix. 3.2 Neumann-Robin problem Once again, boundaries labeled as Γ1 D , Γ2 D , and Γ3 D in Figure 2.1.1 will be subjected to a specific boundary condition obtained by Jiménez-Alfaro, Mantič [2]. Notice that there is a slight variation with what we have previously considered: for this scenario, the scenario of Mode III of fracture, we are not evaluating the Robin-Neumann problem, where the crack is propagating in the positive X direction. We are evaluating the Neumann-Robin problem, where the crack is propagating in the negative X direction. This is because the analytical expressions do represent the Neumann-Robin problem, not the Robin-Neumann problem. Therefore, this variation must be taking into account. Again, this solution is an infinite series of analytical terms: each term has a first dominant term plus multiple shadow terms. Uz,j(r,θ) = U0 z,j(r,θ)+U1 z,j(r,θ)+... +UN z,j(r,θ)(3.8) The three first terms of the displacements’boundary condition are displayed: U0 z,0(r,θ) = 1(3.9) U1 z,0(r,θ) = r π·(θ·sinθ−logr·cosθ)(3.10) U2 z,0(r,θ) = r2cos(2θ) 4π2[log2r−logr−θ2]+ r2θsin(2θ) 4π2[1−2logr](3.11) Boundary labeled as Γ0 Nwill be subjected to: 3.2 Neumann-Robin problem 19 ∂u ∂θ =0for Γ0 N(3.12) Lastly, the boundary related to Robin’s condition will be subjected to: 1 r ∂u ∂θ +κu=0for Γ0 R(3.13) For this simulation, none of the boundaries will be subjected to the homogeneous Dirichlet boundary condition. Notice the behaviour around the crack tip: displacements have a greater gradient, resulting in a highly stressed area, with stress levels increasing infinitely nearby the crack tip. As in the previously evaluated Dirichlet Robin problem, the singularities are fully visible. Notice how the Neumann boundary condition can be visualized: we can identify the Neumann boundary as the boundary where stress-free condition is visible, that is, y=0 and x>0 , where both stress components are null. Simulations are run with E=1,G=1and ν=0. Figure 3.5: Displacements, Uz. 20 Capítulo 3. Anti-plane elasticity, Mode III Figure 3.6: σrz. Figure 3.7: σθz. In order to verify the analytical solution, relative error has also been computed for this simulation. In this situation, we will compute the relative error with respect to the displacements at the crack tip. Throughout the work we will call this normalised error. e=Uz,FEM −Uz,Analytical Uz,cracktip·100 (3.14) 3.2 Neumann-Robin problem 21 where Uz,FEM represents the displacements obtained by the Finite Element Method meanwhile Uz,Analytical stands for the displacements obtained by the analytical approach. Figure 3.8: Normalised error. Notice that the maximum error is below 0.2%, which is considered very low. Only two terms of the analytical solution were used for this analysis. This means the analytical analysis was very robust and precise and the Finite Element Method converges to the proposed analytical solution, and there’s no need to implement more terms of the infinite analytical solution. The implementation can be visualized in the codes provided at the programming appendix. 4 Robin-Neumann problem, Mode II B oundaries labeled as Γ1 D , Γ2 D , and Γ3 D in Figure 2.1.1 will be subjected to a specific boundary condition obtained by Jiménez-Alfaro [9]. Again, this solution is an infinite series of analytical terms, where each term of the series is defined by a first dominant term plus multiple shadow terms. As explained below, in this case we have chosen to use a semicircular mesh, so these three contours are merged into a single contour (refer to Figure 4 for a schematic representation): the contour of the semicircular mesh with R=0.01. Figure 4: Schematic for the semicircular mesh. Notice that, in this fracture mode, the problem can not be reduced to one displacement. Two displacements must be considered. Furthermore, as mentioned in the chapter 1 for the variational formulation (2.2), the constants κ,g and f which refer to the spring stiffness, the body forces and the loads applied on domain must be specified as vectors, not scalars. The original solution proposed can be expressed in the Cartesian coordinate system. Therefore, the analytical solution can be written as a vector u= (ux(x,y),uy(x,y)) (4.1) 23 24 Capítulo 4. Robin-Neumann problem, Mode II A first example can be given, with the main term and one shadow term: ux(r,θ) = 1 8π(6r(π−θ)sin(θ)+ 4rlog(r)cos(θ)+4π)(4.2) uy(r,θ) = 1 4π(r(sin(θ)+(π−θ)cos(θ))) (4.3) r=px2+y2(4.4) θ=arctany x(4.5) In this case, the expressions used for the boundary condition will not always be made explicit since we are going to evaluate the problem with multiple shadow terms, which can be large and costly to read. This is because the Python code developed for this problem is connected to another code where solutions are automatically generated and copied to the problem we want to solve. This code was developed for this specific finite element, tailored for this solution. However, this code is beyond the scope of our work, although we obviously use these analytical expressions when imposing the boundary condition. Simulations are run with E=1,G=1 and ν =0 (except for the last scenario, which will be specified). The implementation can be visualized in the codes provided at the programming appendix. 4.1 Impact of adding new terms 4.1.1 Analysis with one shadow term We will focus on the analysis of the problem with two terms, that is, with a main term and a shadow term. A semicircular mesh will be used, since the rectangular mesh exhibited issues at the corners of the domain: the rectangular mesh has corners and these corners can lead to some stress peaks to the fact that there is a strong change in the boundary conditions and the distance to the crack of the tip. In Mode I and Mode II, these peaks were too huge and therefore a semicircular mesh has been chosen. In the semicircular mesh we also find peaks, and they are also located in the corners of the mesh (y=0,xmax) or (y=0,xmin) , but they are less large. The radius of the semicircle, as mentioned before, will be R=0.01. Therefore, we can determine xmax =0.01 and xmin =−0.01. Stresses are shown in these images. Notice the singularities near the crack tip (x = 0, y = 0), where the gradient of the displacements field, u, is higher: 4.1 Impact of adding new terms 25 Figure 4.1: σrr. Figure 4.2: σrθ. 26 Capítulo 4. Robin-Neumann problem, Mode II Figure 4.3: σθθ . In the case of the σθθ stresses (Figure 4.3), we obtain a result that we already expected: due to the geometry of the problem, the stress tensor has components that can be known in advance. An expected result, which serves as a check, is that the σθθ stresses are zero for y=0. On the other hand, the Neumann condition has been imposed in weak form. This means that in the σrθgraph (Figure 4.2), on the Neumann boundary (y=0, x>0), the stress is not zero pointwise, but on average. Stress singularities must be observed: there should be singularities in the three stress components. It is a logarithmic singularity, since the analytical expression has logarithms. These singularities explain why we need to refine mesh in the vicinity of the crack tip and increase the computational complexity to obtain the numerical solution in this problem in the vicinity of the crack tip. Additionally, to verify the implementation of the Finite Element Method, the normalised error of the applied boundary condition has been calculated. This time, two different conditions must be checked (the boundary condition applied for the x displacement and the boundary condition applied for the y displacement). Firstly, displacements can be visualized: 4.1 Impact of adding new terms 27 Figure 4.4: Displacements. For Mode II, we will be normalising the error with respect to the displacements at the crack tip. Utip =qU2 x+U2 y(4.6) ex=Ux,FEM −Ux,Analytical Utip ·100 (4.7) ey=Uy,FEM −Uy,Analytical Utip ·100 (4.8) 34 Capítulo 4. Robin-Neumann problem, Mode II Figure 4.15: σθθ . For the sake of brevity, the conclusions that have already been remarked on the behaviour of the stresses will not be repeated. We will mainly focus on the improvement of the computational error. Figure 4.16: Displacements. 4.1 Impact of adding new terms 35 Figure 4.17: Error in the x-displacement. Figure 4.18: Error in the y-displacement. Again, from a visual point of view, we do not observe large differences in the displacements, nor in the calculation of the normalised error. It should be checked whether there has been an improvement in the calculation of the displacements, The error for the X displacement is 0.2801498626080477%. The error for the Y displacement is 0.10835584730544187%. The error in the X displacement has slightly reduced from the previous simulation, with the error in the Y displacement remaining with almost no variation. We can 36 Capítulo 4. Robin-Neumann problem, Mode II conclude a proper implementation can be made with only four terms: there’s no need to add the last two terms that we added for this simulation, since the implementation with one dominant term and three shadow terms was already highly precise. Notice that this precision would be enough for the chosen domain: maybe the implementation of a high number of shadow terms in a bigger domain can reduce the error further. It is interesting to verify whether these results agree with those previously obtained by Jiménez- Alfaro [9], that is, whether it is logical that the analysis with more shadow terms does not significantly reduce the error. In this case, we see that for the domain with radius R=0.01 in the aforementioned work, no significant improvements in error calculation are observed: Figure 4.19.1: Urdisplacements obtained for different shadow terms by Jiménez-Alfaro [9]. 4.2 Impact of mesh density 37 Figure 4.19.2: Uθdisplacements obtained for different shadow terms by Jiménez-Alfaro [9]. 4.2 Impact of mesh density The Finite Element Method is affected by the size of the elements that are creating the desired mesh: generally, a denser mesh will result in higher accuracy at the expense of increased computational cost. Initially, our semicircular domain will be meshed by a mesh with elements of size 0.005. We will later simulate with element sizes of 0.001 and 0.0001, and compare to the previous simulation. The size of the element may not give much information about the density of the mesh: in the first case, the mesh shall be composed of 3339 nodes. In the latter case, the mesh shall be composed of 10505 nodes. Finally, for the last case, the mesh is composed of 70960 nodes. Utip =qU2 x+U2 y(4.9) ex=Ux,FEM −Ux,Analytical Utip ·100 (4.10) ey=Uy,FEM −Uy,Analytical Utip ·100 (4.11) 38 Capítulo 4. Robin-Neumann problem, Mode II Figure 4.20: Impact of mesh density. The impact of reducing the element size from 0.005 to 0.001 was noticeable, with improvements in the simulation up to 5%. However, the reduction from 0.001 to 0.0001 did not yield such significant changes. The variation is imperceptible. Overall, we improved the X-displacement error by 5% while the Y-displacement error, interestingly, worsened by 4%. This implies that we have improved our simulation by lowering the element size. 4.3 Impact of spring stiffness Now, an analysis will be conducted in which we modify the stiffness of the springs. Previously, a unitary stiffness was chosen. Now, we will vary this value to 5 and 10, checking for possible variations in the calculation of the error compared to the unit case. Utip =qU2 x+U2 y(4.12) ex=Ux,FEM −Ux,Analytical Utip ·100 (4.13) ey=Uy,FEM −Uy,Analytical Utip ·100 (4.14) Figure 4.21: Impact of spring stiffness, r stands for κ. As expected, the error increases when the spring stiffness increases. This conclusion is in good agreement with the conclusions obtained in Jiménez Alfaro [9], where the error in the Robin boundary condition of the analytical solution was obtained. When the spring stiffness increases, this error is also bigger, and therefore we need to increase the number of shadow terms to obtain a more accurate solution. Consequently, the numerical solution obtained in the benchmark is farther from the exact solution of the Neumann-Robin problem. The Finite Element Solution is slightly different than the analytical solution. The error is still low, but not as low as in the simulations with κ= 1. These results align with what is expected by Jiménez-Alfaro [9]. This study demonstrates that by increasing the stiffness of the springs, either we add more shadow terms or either we accept the error to increase. In this simulations, the same number of shadow terms has been maintained, hence the increase in error. 4.4 Impact of Poisson’s coefficient 39 4.4 Impact of Poisson’s coefficient Previous simulations were run with the hypothesis of a null Poisson’s coefficient. We will now analyze the effect of the Poisson’s coefficient rising to 0.5. Also, we must vary the variational formulation as this problem must be described in terms of plain stress: λ will be expressed in terms of plain stress. The implementation was, again, highly robust: Utip =qU2 x+U2 y(4.15) ex=Ux,FEM −Ux,Analytical Utip ·100 (4.16) ey=Uy,FEM −Uy,Analytical Utip ·100 (4.17) With the error in the X displacement being of 0.11158% and the error in the Y displacement being of 0.5935%. This means a higher error than in the previous implementation, with null Poisson’s coefficient, but low enough to verify the analytical solutions. Therefore, the implementation of the Poisson’s coefficient different from zero did not improve the simulations, but it didn’t notably worsen them either. The Finite Element Method converged with lower precision than the null Poisson’s coefficient scenario, but still with a high precision. Stresses and displacements are shown: Figure 4.22: σrr,ν=0.5. 40 Capítulo 4. Robin-Neumann problem, Mode II Figure 4.23: σrθ,ν=0.5. Figure 4.24: σθθ ,ν=0.5. We do not obtain relevant information from σrr , but we do obtain information from σθθ stresses that are zero for y=0 : this is related to the antisymmetry of the problem. σrθ stress is, on average, zero for the Neumann boundary. 4.5 Presence of singularities 41 Figure 4.25: Displacements, ν=0.5. 4.5 Presence of singularities It can be observed that on the boundary where y=0 , singularities appear to be present at the origin (0,0) in the already mentioned work by Jiménez-Alfaro [9]. It is crucial to determine whether these singularities are real or merely elevated displacement gradients (but not singularities). Therefore, we employed additional software, Paraview, to visualize the stresses σxx and σyy along this line, y=0. This simulation is for null Poisson’s coefficient and 4 shadow terms. Figure 4.26: Stress tensor y=0. Notice that, in Figure 4.25, Stress0 stands for σxx and Stress1 stands for σyy . We do observe a singularity in the tip of the crack: FEniCS is calculating the average stress in every element, and 42 Capítulo 4. Robin-Neumann problem, Mode II it might not look like a singularity because of this reason, but a singularity can be observed if we evaluate the stress components in (0,0). 4.6 Conclusions After all these different analyses under different conditions and for different problems, we can draw conclusions about the elastic behaviour of the solid and the convergence of the analytical and finite element solutions. First and foremost, as expected, we can state that in general terms, reducing the mesh size and increasing the number of terms in the boundary condition improves the accuracy of the solution. There may be simulations in which, occasionally, the implementation of a shadow term does not notably improve the accuracy, but in general terms the improvement is clear. Secondly, as we increase the stiffness constant of the spring, we deviate more and more from the spring condition (since we bring a lot of stiffness to the problem). This makes the accuracy of the simulation worse. We can say that the solutions are sufficiently accurate, but worse than with the unit spring constant. We would need to add even more shadow terms to avoid reducing accuracy if we increase the stiffness. In the case of Poisson’s coefficient, the solution is precise enough to claim the convergence of the analytical solution and the Finite Element Method solution. But the solution was worsened when we evaluated the problem with Poisson’s coefficient set to 0.5. This means a worse simulation than the previous one, with null Poisson’s coefficient. Overall, all simulations converged and we can ensure that analytical solutions and Finite Element Method’s solutions are precise enough to describe the elastic behaviour of the different situations. They are in good agreement with the described analytical solutions, concluding that these solutions are verified for mode II of fracture. 5 Robin-Neumann problem, Mode I O nce again, boundaries labeled as Γ1 D , Γ2 D , and Γ3 D in Figure 2.1.1 will merge to a single boundary, and will be subjected to a specific boundary condition obtained by Jiménez-Alfaro [9]. Again, this solution is an infinite series of analytical terms. Each of this terms can be defined by a first dominant term plus multiple shadow terms. As it happened for Mode II of fracture, we have chosen to use a semicircular mesh, with radius 0.01. Refer to figure 4 for a schematic representation. Notice that, again, in this fracture mode (Mode I), the problem can not be reduced to one displacement. Two displacements must be considered. It is also recalled that the constants κ , g and f, as explained in Mode II, must be specified as vectors. The original solution proposed by Jiménez-Alfaro [9] can be expressed in the Cartesian coordinate system. Therefore, the analytical solution can be written as a vector: u= (ux(x,y),uy(x,y)) (5.1) A first example can be given, with the main term and one shadow term: ux(r,θ) = 1 4π(r((−π+θ)cos(θ)+(1+2log(r))sin(θ))) (5.2) uy(r,θ) = 1 4π(2π+2rcos(θ)log(r)+ r(π−θ)sin(θ)) (5.3) r=px2+y2(5.4) θ=arctany x(5.5) Once again, for this scenario the expressions used for the boundary condition will not be made explicit since we are going to evaluate the problem with multiple shadow terms, which can be large and costly to read. The expressions are imported from an original Python code that generates these expressions. Simulations are run with E=1 , G=1 and ν =0 (except for the last scenario, which will be specified). The implementation can be visualized in the codes provided at the programming appendix. 5.1 Impact of adding new terms 5.1.1 Analysis with one shadow term Now we will focus on the analysis of the problem with two terms, that is, with a dominant term and a shadow term. This time, we will use a semicircular mesh, since the rectangular mesh exhibited 43 50 Capítulo 5. Robin-Neumann problem, Mode I Figure 5.10: Displacements. Figure 5.11: Error in the x-displacement. 5.1 Impact of adding new terms 51 Figure 5.12: Error in the y-displacement. Finally, the maximum error for the x boundary condition in displacement is 1.3124200460535968%. The maximum error for the y boundary condition in displacement is 0.04791299727270815%. This means that our latest simulation has a similar error to the two previous simulations. The error is small enough to claim that the analytical solutions have been developed correctly, but the effect of adding different shadow terms does not seem to be relevant, as the decrease in error is not truly remarkable. 5.1.4 Analysis with five shadow terms As in Mode II, our final analysis will be comprising six terms: the dominant term and five shadow terms. This mean our analysis should be precise and robust. Displacements and stresses are also displayed. 52 Capítulo 5. Robin-Neumann problem, Mode I Figure 5.13: σrr. Figure 5.14: σrθ. 5.1 Impact of adding new terms 53 Figure 5.15: σθθ . Once again, similar conclusions to the ones described before can be made for this simulation, in terms of the stress tensor. For the sake of brevity, we will not repeat these conclusions. Figure 5.16: Displacements. 54 Capítulo 5. Robin-Neumann problem, Mode I Figure 5.17: Error in the x-displacement. Figure 5.18: Error in the y-displacement. The error for the X displacement remains located at 1.372428625709376%. The error for the Y displacement is situated at 0.046284048088260754%. Simulating with such a high number of terms means a high computational cost but only slightly varies the precision from the previous simulations. Our simulation was slightly better in the Y displacement while it was slightly worse for the X displacement. Overall, we did improve the simulation by adding new shadow terms, but the variation is not remarkable, while the computational cost exponentially increased. This is consistent with observations made by Jiménez-Alfaro [9], where it is described that the 5.1 Impact of adding new terms 55 displacement solution for R=0.01 is identical for different number of shadow terms, with almost imperceptible variations. Figure 5.19.1: Urdisplacements obtained for different shadow terms by Jiménez-Alfaro [9]. Figure 5.19.2: Uθdisplacements obtained for different shadow terms by Jiménez-Alfaro [9]. 56 Capítulo 5. Robin-Neumann problem, Mode I 5.2 Impact of mesh density The size of the element will have a major impact in the application of the Finite Element Method: generally, a denser mesh will result in higher accuracy at the expense of increased computational cost, while a coarser mesh will also be less accurate but much faster. The simulations will be run with three shadow terms. Initially, our semicircular domain will be meshed by a mesh with elements of size 0.005. We will later simulate with element sizes of 0.001 and 0.0001, and compare with the previous simulation to obtain information of the behaviour of the elastic problem. Once again, this may not give information about the mesh: it might be better to specify the number of nodes evaluated. In the first case, the mesh shall be composed of 3339 nodes. In the latter case, the mesh shall be composed of 10505 nodes. Finally, for the last case, the mesh is composed of 70960 nodes. Figure 5.20: Effect of mesh density. The impact of reducing the mesh size from 0.005 to 0.001 has been significant: we have seen an improvement in accuracy of up to 36%. However, the reduction from 0.001 to 0.0001 did not yield such significant changes, with variations being up to 2%. Therefore, we can state that with the second simulation, the accuracy is already sufficiently high, with the third simulation (the most computationally expensive) having low impact. 5.3 Impact of spring stiffness Now, an analysis will be conducted in which we modify the stiffness of the springs. Previously, a unitary stiffness was chosen. Now, we will vary this value to 5 and 10, checking for possible variations in the calculation of the error compared to the unit case. We expect computed error to increase, as increasing the stiffness of the springs means we deviate from the original spring behaviour. Figure 5.21: Effect of spring stiffness, r stands for κ. As expected, the error increases markedly as the behaviour of the springs becomes stiffer. In this case, the error for κ=10 is 13%, which is much higher than the errors we have had in previous simulations. Notice that this simulation was run with a refined mesh, in order to be precise in our analysis. This means that there is no precise convergence between the finite element method and the analytical solutions for this problem. 5.4 Impact of Poisson’s coefficient 57 Once again, these results align with what is expected by Jiménez-Alfaro [9]. If we increase the stiffness of the springs, we have to increase the shadow terms or accept a higher error. As these simulations were run with constant number of shadow terms, a higher error is expected. 5.4 Impact of Poisson’s coefficient A new analysis will be made, where the Poisson’s coefficient is set to 0.5, in order to verify the behaviour of the elastic solid. We will evaluate stress, displacements and the computed error. As mentioned for mode II of fracture, this analysis must be done under plain stress conditions. Figure 5.22: σrr,ν=0.5. 58 Capítulo 5. Robin-Neumann problem, Mode I Figure 5.23: σrθ,ν=0.5. Figure 5.24: σθθ ,ν=0.5. In the case of the σθθ stresses, we obtain a result that we already expected: due to the geometry of the problem, the stress tensor has components that can be known in advance. An expected result, is that the σθθ stresses are close to zero for y=0 and x<0. On the other hand, the Neumann condition has been tacitly imposed in weak form for the Neumann boundary (x>0,y=0). More information on the elastic behaviour of the problem and specially the convergence of the analytical solution and the Finite Element Method can be obtained by evaluating the normalised 5.4 Impact of Poisson’s coefficient 59 error. Figure 5.25: Displacements, ν=0.5. Utip =qU2 x+U2 y(5.9) ex=Ux,FEM −Ux,Analytical Utip ·100 (5.10) ey=Uy,FEM −Uy,Analytical Utip ·100 (5.11) 8 Programming appendix 8.1 Mesh Generation 8.1.1 Semicircular mesh generation import gmsh from math import pi, cos, sin gmsh.initialize() gdim = 2 R = 0.01 m0 = 0.0001 num_points = 1000 angle_increment = pi / (num_points - 1) points = [] for i in range(num_points): angle = i * angle_increment x = R * cos(angle) y = R * sin(angle) points.append(gmsh.model.geo.addPoint(x, y, 0, m0)) extra_points = [] num_extra_points = 10 x_increment = (R - 0.0001) / (num_extra_points + 1) for i in range(1, num_extra_points + 1): x = -0.0001 * i extra_points.append(gmsh.model.geo.addPoint(x, 0, 0, m0)) for i in range(1, num_extra_points + 1): x = 0.0001 * i extra_points.append(gmsh.model.geo.addPoint(x, 0, 0, m0)) 67 68 Capítulo 8. Programming appendix lines = [] for i in range(len(points) - 1): lines.append(gmsh.model.geo.addLine(points[i], points[i + 1])) p_center = gmsh.model.geo.addPoint(0, 0, 0, m0) lines.append(gmsh.model.geo.addLine(points[-1], extra_points[ num_extra_points - 1])) for i in range(num_extra_points - 1, 0, -1): lines.append(gmsh.model.geo.addLine(extra_points[i], extra_points[i - 1])) lines.append(gmsh.model.geo.addLine(extra_points[0], p_center)) lines.append(gmsh.model.geo.addLine(p_center, extra_points[ num_extra_points])) for i in range(num_extra_points, len(extra_points) - 1): lines.append(gmsh.model.geo.addLine(extra_points[i], extra_points[i + 1])) lines.append(gmsh.model.geo.addLine(extra_points[-1], points[0])) semicircle_loop = gmsh.model.geo.addCurveLoop(lines) surface = gmsh.model.geo.addPlaneSurface([semicircle_loop]) gmsh.model.geo.synchronize() domain_id = 1 gmsh.model.addPhysicalGroup(gdim, [surface], domain_id) gmsh.model.setPhysicalName(gdim, domain_id, "domain") boundary_id = 2 gmsh.model.addPhysicalGroup(gdim - 1, lines, boundary_id) gmsh.model.setPhysicalName(gdim - 1, boundary_id, "boundary") gmsh.model.mesh.generate(gdim) gmsh.write("malla_semicircular.msh") gmsh.fltk.run() gmsh.finalize() 8.1.2 Rectangular mesh generation import gmsh 8.2 Anti-plane elasticity problems 69 from math import* gmsh.initialize() gdim = 2 L, W, e = 0.2, 0.2, 0 m0 = 0.0005 p1 = gmsh.model.geo.addPoint(0,0,0,m0) p2 = gmsh.model.geo.addPoint(L/2,0,0,m0) p3 = gmsh.model.geo.addPoint(L/2,W,0,m0) p4 = gmsh.model.geo.addPoint(-L/2,W,0,m0) p5 = gmsh.model.geo.addPoint(-L/2,e,0,m0) l1 = gmsh.model.geo.addLine(p1,p2) l2 = gmsh.model.geo.addLine(p2,p3) l3 = gmsh.model.geo.addLine(p3,p4) l4 = gmsh.model.geo.addLine(p4,p5) l5 = gmsh.model.geo.addLine(p5,p1) cloop1 = gmsh.model.geo.addCurveLoop([l1,l2,l3,l4,l5]) surface1 = gmsh.model.geo.addPlaneSurface([cloop1],1) gmsh.model.geo.synchronize() ps1 = gmsh.model.addPhysicalGroup(gdim, [surface1], gdim) gmsh.model.setPhysicalName(gdim, ps1, "domain") pl1 = gmsh.model.addPhysicalGroup(gdim-1, [l1,l2,l3,l4,l5],1) gmsh.model.setPhysicalName(gdim-1, pl1 , "up_dofs") gmsh.model.mesh.generate(gdim) gmsh.write("malla_modo3.msh") print(gdim) gmsh.fltk.run() gmsh.finalize() 8.2 Anti-plane elasticity problems 8.2.1 Dirichlet-Robin from fenics import* import numpy as np from msh2xdmf_cambios import* 70 Capítulo 8. Programming appendix import matplotlib.pyplot as plt from matplotlib import cm from matplotlib.ticker import LinearLocator from mpl_toolkits.mplot3d import Axes3D from dolfin import* from matplotlib import rc msh2xdmf(’malla_modo3.msh’,2,directory=’.’) mesh, mesh_function, association_table = import_mesh( prefix=’malla_modo3’, dim=2, directory=’.’) print(mesh) V = FunctionSpace(mesh, ’CG’,1) Boundary_markers = MeshFunction(’size_t’, mesh, mesh.geometry().dim()-1) class Bound1(SubDomain): def inside(self,x,on_boundary): return x[0]<0 and x[1]==0 class Bound2(SubDomain): def inside(self,x,on_boundary): return x[0]>=0 and x[1]==0 boundary1 = Bound1() boundary2 = Bound2() boundary1.mark(Boundary_markers, 1) boundary2.mark(Boundary_markers, 2) class Bound3(SubDomain): def inside(self,x,on_boundary): return near (x[0], -0.1) and on_boundary class Bound4(SubDomain): def inside(self,x,on_boundary): return near (x[1], 0.2) and on_boundary class Bound5(SubDomain): def inside(self,x,on_boundary): return near (x[0], 0.1) and on_boundary boundary3 = Bound3() boundary4 = Bound4() boundary5 = Bound5() 8.2 Anti-plane elasticity problems 71 boundary3.mark(Boundary_markers,3) boundary4.mark(Boundary_markers,4) boundary5.mark(Boundary_markers,5) u_D1 = Expression(’pow(x[0]*x[0]+x[1]*x[1],0.25)*sin(atan2(x[1],x[0])/2) -2/(3*pi)*pow(x[0]*x[0]+x[1]*x[1],0.75)*(atan2(x[1],x[0])*cos(3* atan2(x[1],x[0])/2)+std::log(sqrt(x[0]*x[0]+x[1]*x[1]))*sin(3*atan2( x[1],x[0])/2))’,degree=3) u_D2 = Constant(0.0) r = Constant(1.0) g = Constant(0.0) f = Constant(0.0) bcd0 = DirichletBC(V, u_D1, boundary3, method=’topological’) bcd2 = DirichletBC(V, u_D1, boundary4, method=’topological’) bcd4 = DirichletBC(V, u_D2, boundary2, method=’topological’) bcd5 = DirichletBC(V, u_D1, boundary5, method=’topological’) bcs = [bcd0, bcd2, bcd4, bcd5] ds = Measure(’ds’, domain = mesh, subdomain_data = Boundary_markers) dx = Measure(’dx’, domain = mesh, subdomain_data = Boundary_markers) u = TrialFunction(V) v = TestFunction(V) a = dot(grad(u),grad(v))*dx + r*u*v*ds(1) L = f*v*dx + g*v*ds A = assemble(a) b = assemble(L) u = Function(V, name=’Displacement’) solve (a==L,u,bcs) G=1 gradu = nabla_grad(u) V = FunctionSpace(mesh, ’CG’, 1) uxp = project(gradu[0],V) ux = uxp.compute_vertex_values(mesh) uyp = project(gradu[1],V) uy = uyp.compute_vertex_values(mesh) x = mesh.coordinates()[:,0] y = mesh.coordinates()[:,1] 72 Capítulo 8. Programming appendix sigmarz = G * np.true_divide(np.multiply(ux,x),np.sqrt(np.multiply(x,x)+ np.multiply(y,y))) + G * np.true_divide(np.multiply(uy,y),np.sqrt(np. multiply(x,x)+np.multiply(y,y))) sigmathetaz = np.true_divide(G,np.sqrt(np.multiply(x,x)+np.multiply(y,y) )) * (np.multiply(ux,-y) + np.multiply(uy,x)) fig = plt.figure(1) ax = fig.add_subplot(111, projection=’3d’) surf_sigmathetaz = ax.plot_trisurf(x, y, sigmathetaz, cmap=cm.coolwarm, linewidth=0, antialiased=False, vmin=np.nanmin( sigmathetaz), vmax=np.nanmax(sigmathetaz)) ax.set_xlabel(’X-axis’) ax.set_ylabel(’Y-axis’) ax.set_zlabel(’Sigmathetaz’) fig.colorbar(surf_sigmathetaz, shrink=0.5, aspect=5) fig = plt.figure(2) ax = fig.add_subplot(111, projection=’3d’) surf_sigmarz = ax.plot_trisurf(x, y, sigmarz, cmap=cm.coolwarm, linewidth=0, antialiased=False, vmin=np.nanmin( sigmarz), vmax=np.nanmax(sigmarz)) ax.set_xlabel(’X-axis’) ax.set_ylabel(’Y-axis’) ax.set_zlabel(’Sigmarz’) fig.colorbar(surf_sigmarz, shrink=0.5, aspect=5) vtkfile = File(’results/displacements.pvd’) vtkfile << u intervalo_x = 0.001 intervalo_y = 0.001 Xp = np.arange(-0.1, 0.1, intervalo_x) Yp = np.arange(0, 0.2, intervalo_y) X, Y = np.meshgrid(Xp, Yp) Z = np.zeros_like(X) for i in range(X.shape[0]): for j in range(X.shape[1]): u.set_allow_extrapolation(True) Z[i, j] = u(X[i, j], Y[i, j]) fig = plt.figure(3) ax = fig.add_subplot(111, projection=’3d’) 8.2 Anti-plane elasticity problems 73 surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm, linewidth=0, antialiased=False) ax.set_xlabel(’X-axis’) ax.set_ylabel(’Y-axis’) ax.set_zlabel(’U_z’) fig.colorbar(surf, shrink=0.5, aspect=5) u_tip = u(-0.1,0.2) def bcd(x0,x1): pi = np.pi return np.power(x0*x0 + x1*x1, 0.25)*np.sin(np.arctan2(x1,x0)/2) -2/(3*pi)*np.power(x0*x0 + x1*x1, 0.75)*(np.arctan2(x1,x0)*np.cos (3*np.arctan2(x1,x0)/2) + np.log(np.sqrt(x0*x0 + x1*x1))*np.sin (3*np.arctan2(x1,x0)/2)) Zbc = bcd(X,Y) fig = plt.figure(4) ax = fig.add_subplot(111, projection=’3d’) surf = ax.plot_surface(X, Y, abs(Zbc-Z)/abs(u_tip)*100, cmap=’viridis’) ax.set_xlabel(’X-axis’) ax.set_ylabel(’Y-axis’) ax.set_zlabel(’Normalised error W’) fig.colorbar(surf, shrink=0.5, aspect=5) print(np.max(abs(Zbc-Z))/u_tip*100) plt.show() 8.2.2 Robin-Neumann from fenics import* import numpy as np from msh2xdmf_cambios import* import matplotlib.pyplot as plt from matplotlib import cm from matplotlib.ticker import LinearLocator from mpl_toolkits.mplot3d import Axes3D from dolfin import* msh2xdmf(’malla_modo3.msh’,2,directory=’.’) mesh, mesh_function, association_table = import_mesh( 74 Capítulo 8. Programming appendix prefix=’malla_modo3’, dim=2, directory=’.’) print(mesh) V = FunctionSpace(mesh, ’CG’,1) Boundary_markers = MeshFunction(’size_t’, mesh, mesh.geometry().dim()-1) class Bound1(SubDomain): def inside(self,x,on_boundary): return x[0]<0 and x[1]<0.01 class Bound2(SubDomain): def inside(self,x,on_boundary): return x[0]>=0 and x[1]==0 boundary1 = Bound1() boundary2 = Bound2() boundary1.mark(Boundary_markers, 1) boundary2.mark(Boundary_markers, 2) class Bound3(SubDomain): def inside(self,x,on_boundary): return near (x[0], -0.1) and on_boundary class Bound4(SubDomain): def inside(self,x,on_boundary): return near (x[1], 0.2) and on_boundary class Bound5(SubDomain): def inside(self,x,on_boundary): return near (x[0], 0.1) and on_boundary boundary3 = Bound3() boundary4 = Bound4() boundary5 = Bound5() boundary3.mark(Boundary_markers,3) boundary4.mark(Boundary_markers,4) boundary5.mark(Boundary_markers,5) u_D1 = Expression(’1+pow(x[0]*x[0]+x[1]*x[1],0.5)/pi*(atan2(x[1],x[0])* sin(atan2(x[1],x[0]))-std::log(pow(x[0]*x[0]+x[1]*x[1],0.5))*cos( atan2(x[1],x[0])))’,degree=3) u_D2 = Constant(0.0) r = Constant(1.0) 8.2 Anti-plane elasticity problems 75 g = Constant(0.0) f = Constant(0.0) bcd0 = DirichletBC(V, u_D1, boundary3, method=’topological’) bcd2 = DirichletBC(V, u_D1, boundary4, method=’topological’) bcd5 = DirichletBC(V, u_D1, boundary5, method=’topological’) bcs = [bcd0, bcd2, bcd5] ds = Measure(’ds’, domain = mesh, subdomain_data = Boundary_markers) dx = Measure(’dx’, domain = mesh, subdomain_data = Boundary_markers) u = TrialFunction(V) v = TestFunction(V) a = dot(grad(u),grad(v))*dx + r*u*v*ds(1) L = f*v*dx + g*v*ds A = assemble(a) b = assemble(L) u = Function(V, name=’Displacement’) solve (a==L,u,bcs) G=1 gradu = nabla_grad(u) V = FunctionSpace(mesh, ’CG’, 1) uxp = project(gradu[0],V) ux = uxp.compute_vertex_values(mesh) uyp = project(gradu[1],V) uy = uyp.compute_vertex_values(mesh) x = mesh.coordinates()[:,0] y = mesh.coordinates()[:,1] sigmarz = G * np.true_divide(np.multiply(ux,x),np.sqrt(np.multiply(x,x)+ np.multiply(y,y))) + G * np.true_divide(np.multiply(uy,y),np.sqrt(np. multiply(x,x)+np.multiply(y,y))) sigmathetaz = np.true_divide(G,np.sqrt(np.multiply(x,x)+np.multiply(y,y) )) * (np.multiply(ux,-y) + np.multiply(uy,x)) fig = plt.figure(1) ax = fig.add_subplot(111, projection=’3d’) surf = ax.plot_trisurf(x, y, sigmathetaz, cmap=cm.coolwarm, linewidth=0, antialiased=False, vmin=np.nanmin( sigmathetaz), vmax=np.nanmax(sigmathetaz)) ax.set_xlabel(’X-axis’) ax.set_ylabel(’Y-axis’) 82 Capítulo 8. Programming appendix fig = plt.figure(figsize=(14, 7)) ax1 = fig.add_subplot(121, projection=’3d’) surf1 = ax1.plot_surface(Xp, Yp, ux_values, cmap=’viridis’) fig.colorbar(surf1, ax=ax1, shrink=0.5, aspect=5) ax1.set_title(’U X displacements’) ax1.set_xlabel(’x’) ax1.set_ylabel(’y’) ax1.set_zlabel(’ux’) ax2 = fig.add_subplot(122, projection=’3d’) surf2 = ax2.plot_surface(Xp, Yp, uy_values, cmap=’plasma’) fig.colorbar(surf2, ax=ax2, shrink=0.5, aspect=5) ax2.set_title(’U Y displacements’) ax2.set_xlabel(’x’) ax2.set_ylabel(’y’) ax2.set_zlabel(’uy’) plt.tight_layout() uaux = u(0,0) ux_tip = uaux[0] uy_tip = uaux[1] u_tip = np.sqrt(ux_tip*ux_tip + uy_tip*uy_tip) def bcd_x(x0,x1): pi = np.pi r = np.sqrt(x0**2 + x1**2) theta = np.arctan2(x1, x0) return 1/(8*pi)*(6*r*(pi-theta)*np.sin(theta)+4*r*np.log(r)*np.cos( theta)+4*pi) def bcd_y(x0,x1): pi = np.pi r = np.sqrt(x0**2 + x1**2) theta = np.arctan2(x1, x0) return 1/(4*pi)*(r*(sin(theta)+(pi-theta)*cos(theta))) Zbc_x = np.zeros_like(Xp) Zbc_y = np.zeros_like(Yp) aux_x = np.zeros_like(Xp) aux_y = np.zeros_like(Yp) Zx = np.zeros_like(Xp) Zy = np.zeros_like(Yp) for i in range(Xp.shape[0]): for j in range(Xp.shape[1]): 8.3 Robin-Neumann, Mode II 83 u_val = u(Xp[i, j], Yp[i, j]) Zx[i, j] = u_val[0] Zy[i, j] = u_val[1] Zbc_x[i, j] = bcd_x(Xp[i, j], Yp[i, j]) Zbc_y[i, j] = bcd_y(Xp[i, j], Yp[i, j]) aux_x[i, j] = Zbc_x[i, j] - Zx[i, j] aux_y[i, j] = Zbc_y[i, j] - Zy[i, j] fig = plt.figure(7) ax = fig.add_subplot(111, projection=’3d’) surf = ax.plot_surface(Xp, Yp, abs(aux_y)/u_tip*100, cmap=’viridis’) ax.set_xlabel(’X-axis’) ax.set_ylabel(’Y-axis’) ax.set_zlabel(’normalised error U_Y’) fig.colorbar(surf, shrink=0.5, aspect=5) fig = plt.figure(6) ax = fig.add_subplot(111, projection=’3d’) surf = ax.plot_surface(Xp, Yp, abs(aux_x)/u_tip*100, cmap=’viridis’) ax.set_xlabel(’X-axis’) ax.set_ylabel(’Y-axis’) ax.set_zlabel(’normalised error U_X’) fig.colorbar(surf, shrink=0.5, aspect=5) print(np.max(abs(aux_x))/u_tip*100) print(np.max(abs(aux_y))/u_tip*100) def eps(v): return sym(grad(v)) def sigma(v): return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) T = TensorFunctionSpace(mesh,’DG’,2) stress = Function(T,name="Stress") stress.assign(project(sigma(u),T)) vtkfile = File(’requestSara/tensiones.pvd’) vtkfile << stress plt.show() 84 Capítulo 8. Programming appendix 8.4 Robin-neumann, Mode I For the sake of brevity, the script will only be displayed with one shadow term. More terms can be added in the same expression. from fenics import* import numpy as np import sympy from msh2xdmf_cambios import* import matplotlib.pyplot as plt from matplotlib import cm from matplotlib.ticker import LinearLocator from mpl_toolkits.mplot3d import Axes3D import dolfin from trans_r_theta_M1 import calculate_displacements from math import* msh2xdmf(’malla_semicircular.msh’,2,directory=’.’) mesh, mesh_function, association_table = import_mesh( prefix=’malla_semicircular’, dim=2, directory=’.’) V = VectorFunctionSpace(mesh, ’CG’,1) Boundary_markers = MeshFunction(’size_t’, mesh, mesh.geometry().dim()-1) class Bound1(SubDomain): def inside(self,x,on_boundary): return x[0]<0 and x[1]==0 class Bound2(SubDomain): def inside(self,x,on_boundary): return x[0]>=0 and x[1]==0 boundary1 = Bound1() boundary2 = Bound2() boundary1.mark(Boundary_markers, 1) boundary2.mark(Boundary_markers, 2) R2 = 0.01 class Bound3(SubDomain): def inside(self,x,on_boundary): return near(sqrt(x[0]**2+x[1]**2),R2,0.0001) boundary3 = Bound3() 8.4 Robin-neumann, Mode I 85 boundary3.mark(Boundary_markers,3) disp_x0 = Expression(’1/(4*pi)*(r*((-pi+theta)*cos(theta)+(1+2*std::log( r))*sin(theta)))’,r = Expression(’sqrt(x[0]*x[0]+x[1]*x[1])’,degree =3),theta = Expression(’atan2(x[1],x[0])’,degree=3),degree=3) disp_x1 = Expression(’1/(4*pi)*(2*pi+2*r*cos(theta)*std::log(r)+r*(pitheta)*sin(theta))’, r = Expression(’sqrt(x[0]*x[0]+x[1]*x[1])’, degree=3),theta = Expression(’atan2(x[1],x[0])’,degree=3),degree=3) disp_null = Constant(0) bcd0 = DirichletBC(V.sub(0), disp_x0, boundary3, method=’topological’) bcd1 = DirichletBC(V.sub(1), disp_x1, boundary3, method=’topological’) bcs = [bcd0, bcd1] ds = Measure(’ds’, domain = mesh, subdomain_data = Boundary_markers) dx = Measure(’dx’, domain = mesh, subdomain_data = Boundary_markers) u = TrialFunction(V) v = TestFunction(V) lmbda_or = Constant(1) mu = Constant(1) lmbda = (2*lmbda_or*mu)/(lmbda_or + mu*2) def eps(u): return sym(grad(u)) def sigma(u): return lmbda*tr(eps(u))*Identity(2)+2.0*mu*eps(u) r = Constant((0.0, 1.0)) g = Constant((0.0, 0.0)) f = Constant((0.0, 0.0)) a = inner(sigma(u), eps(v))*dx + (r[0]*u[0]*v[0] + r[1]*v[1]*u[1])*ds(2) L = dot(f, v)*dx + inner(g,v)*ds(1) u = Function(V, name=’Displacement’) solve (a==L,u,bcs) vtkfile = File(’results/displacements.pvd’) vtkfile << u Vp = 1 file = File(str(’./RobinNeumannEdu’+str(Vp)+’/RobinNeumannEdu’+str(Vp)+’ .pvd’)) file << Boundary_markers 86 Capítulo 8. Programming appendix File("domains.pvd")<<Boundary_markers T = TensorFunctionSpace(mesh, ’CG’, 1) stress = Function(T, name="Stress") stress.assign(project(sigma(u),T)) vtkfile = File(’results/stress.pvd’) vtkfile << stress E = Constant(1) G = Constant(1) nu = Constant(0) lmbda_or = Constant(1) E_val = E.values() nu_val = nu.values() G_val = G.values() lmbda_val_or = lmbda_or.values() mu_val = 1 lmbda_val = (2*lmbda_val_or*G_val)/(lmbda_val_or + 2*G_val) gradu = eps(u) V = FunctionSpace(mesh, ’DG’, 1) uxp = project(gradu[0,0],V) ux = uxp.compute_vertex_values(mesh) uyp = project(gradu[1,1],V) uy = uyp.compute_vertex_values(mesh) uxyp = project(gradu[0,1],V) uxy = uyp.compute_vertex_values(mesh) uyxp = project(gradu[1,0],V) uyx = uyp.compute_vertex_values(mesh) x = mesh.coordinates()[:,0] y = mesh.coordinates()[:,1] sigmaxx = np.multiply(lmbda_val, ux+uy) + 2*np.multiply(mu_val,ux) sigmayy = np.multiply(lmbda_val, ux+uy) + 2*np.multiply(mu_val,uy) sigmaxy = mu_val * (uxy + uyx) ntol = 1E-15 sigmarr = np.zeros(len(x)) sigmartheta = np.zeros(len(x)) sigmathetatheta = np.zeros(len(x)) for i in range(len(x)): a = x[i] / (ntol + np.sqrt(np.multiply(x[i], x[i]) + np.multiply(y [i], y[i]))) b = y[i] / (ntol + np.sqrt(np.multiply(x[i], x[i]) + np.multiply(y [i], y[i]))) 8.4 Robin-neumann, Mode I 87 L = np.array([[a, b], [-b, a]]) Lt = np.array([[a, -b], [b, a]]) tens = np.array([[sigmaxx[i], sigmaxy[i]], [sigmaxy[i], sigmayy[i ]]]) sigma = np.dot(np.dot(L, tens), Lt) sigmarr[i] = sigma[0, 0] sigmartheta[i] = sigma[0, 1] sigmathetatheta[i] = sigma[1, 1] fig = plt.figure(1) ax = fig.add_subplot(111, projection=’3d’) surf = ax.plot_trisurf(x, y, sigmarr, cmap=cm.coolwarm, linewidth=0, antialiased=False, vmin=np.nanmin( sigmarr), vmax=np.nanmax(sigmarr)) ax.set_xlabel(’X-axis’) ax.set_ylabel(’Y-axis’) ax.set_zlabel(’Sigmarr’) fig.colorbar(surf, shrink=0.5, aspect=5) fig = plt.figure(2) ax = fig.add_subplot(111, projection=’3d’) surf = ax.plot_trisurf(x, y, sigmartheta, cmap=cm.coolwarm, linewidth=0, antialiased=False, vmin=np.nanmin( sigmartheta), vmax=np.nanmax(sigmartheta)) ax.set_xlabel(’X-axis’) ax.set_ylabel(’Y-axis’) ax.set_zlabel(’Sigmartheta’) fig.colorbar(surf, shrink=0.5, aspect=5) fig = plt.figure(3) ax = fig.add_subplot(111, projection=’3d’) surf = ax.plot_trisurf(x, y, sigmathetatheta, cmap=cm.coolwarm, linewidth=0, antialiased=False, vmin=np.nanmin( sigmathetatheta), vmax=np.nanmax( sigmathetatheta)) ax.set_xlabel(’X-axis’) ax.set_ylabel(’Y-axis’) ax.set_zlabel(’Sigmathetatheta’) fig.colorbar(surf, shrink=0.5, aspect=5) theta = np.linspace(0, np.pi, 1000) r = np.linspace(0.0001, 0.01, 1000) theta, r = np.meshgrid(theta, r) 88 Capítulo 8. Programming appendix Xp = r * np.cos(theta) Yp = r * np.sin(theta) ux_values = np.zeros_like(Xp) uy_values = np.zeros_like(Yp) for i in range(Xp.shape[0]): for j in range(Xp.shape[1]): u.set_allow_extrapolation(True) ux, uy = u(Xp[i, j], Yp[i, j]) ux_values[i, j] = ux uy_values[i, j] = uy fig = plt.figure(figsize=(14, 7)) ax1 = fig.add_subplot(121, projection=’3d’) surf1 = ax1.plot_surface(Xp, Yp, ux_values, cmap=’plasma’) fig.colorbar(surf1, ax=ax1, shrink=0.5, aspect=5) ax1.set_title(’U X displacements’) ax1.set_xlabel(’x’) ax1.set_ylabel(’y’) ax1.set_zlabel(’ux’) ax2 = fig.add_subplot(122, projection=’3d’) surf2 = ax2.plot_surface(Xp, Yp, uy_values, cmap=’plasma’) fig.colorbar(surf2, ax=ax2, shrink=0.5, aspect=5) ax2.set_title(’U Y displacements’) ax2.set_xlabel(’x’) ax2.set_ylabel(’y’) ax2.set_zlabel(’uy’) plt.tight_layout() uaux = u(0,0) ux_tip = uaux[0] uy_tip = uaux[1] u_tip = np.sqrt(ux_tip*ux_tip + uy_tip*uy_tip) def bcd_x(x0,x1): pi = np.pi r = np.sqrt(x0**2 + x1**2) theta = np.arctan2(x1, x0) return 1/(4*pi)*(r*((-pi+theta)*np.cos(theta)+(1+2*np.log(r))*np.sin (theta))) def bcd_y(x0,x1): pi = np.pi r = np.sqrt(x0**2 + x1**2) theta = np.arctan2(x1, x0) 8.4 Robin-neumann, Mode I 89 return 1/(4*pi)*(2*pi + 2*r*np.cos(theta)*np.log(r) + r*(pi-theta)* sin(theta)) Zbc_x = np.zeros_like(Xp) Zbc_y = np.zeros_like(Yp) aux_x = np.zeros_like(Xp) aux_y = np.zeros_like(Yp) Zx = np.zeros_like(Xp) Zy = np.zeros_like(Yp) for i in range(Xp.shape[0]): for j in range(Xp.shape[1]): u_val = u(Xp[i, j], Yp[i, j]) Zx[i, j] = u_val[0] Zy[i, j] = u_val[1] Zbc_x[i, j] = bcd_x(Xp[i, j], Yp[i, j]) Zbc_y[i, j] = bcd_y(Xp[i, j], Yp[i, j]) aux_x[i, j] = Zbc_x[i, j] - Zx[i, j] aux_y[i, j] = Zbc_y[i, j] - Zy[i, j] fig = plt.figure(7) ax = fig.add_subplot(111, projection=’3d’) surf = ax.plot_surface(Xp, Yp, abs(aux_y)/u_tip*100, cmap=’viridis’) ax.set_xlabel(’X-axis’) ax.set_ylabel(’Y-axis’) ax.set_zlabel(’Normalised error U_Y’) fig.colorbar(surf, shrink=0.5, aspect=5) fig = plt.figure(6) ax = fig.add_subplot(111, projection=’3d’) surf = ax.plot_surface(Xp, Yp, abs(aux_x)/u_tip*100, cmap=’viridis’) ax.set_xlabel(’X-axis’) ax.set_ylabel(’Y-axis’) ax.set_zlabel(’Normalised error U_X’) fig.colorbar(surf, shrink=0.5, aspect=5) print(np.max(abs(aux_x))/u_tip*100) print(np.max(abs(aux_y))/u_tip*100) plt.show() def eps(v): return sym(grad(v)) def sigma(v): 90 Capítulo 8. Programming appendix return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) T = TensorFunctionSpace(mesh,’CG’,2) stress = Function(T,name="Stress") stress.assign(project(sigma(u),T)) vtkfile = File(’requestSara/tensiones.pvd’) vtkfile << stress