Abstract
This work presents an analysis methodology based on the use of the Finite Element Method (FEM) nowadays considered one of the main numerical tools for solving Boundary Value Problems (BVPs). The proposed methodology, so-called cg-FEM (Cartesian grid FEM), has been implemented for fast and accurate numerical analysis of 2D linear elasticity problems. The traditional FEM uses geometry-conforming meshes; however, in cg-FEM the analysis mesh is not conformal to the geometry. This allows for defining very efficient mesh generation techniques and using a robust integration procedure, to accurately integrate the domain’s geometry. The hierarchical data structure used in cg-FEM together with the Cartesian meshes allow for trivial data sharing between similar entities. The cg-FEM methodology uses advanced recovery techniques to obtain an improved solution of the displacement and stress fields (for which a discretization error estimator in energy norm is available) that will be the output of the analysis. All this results in a substantial increase in accuracy and computational efficiency with respect to the standard FEM. cg-FEM has been applied in structural shape optimization showing robustness and computational efficiency in comparison with FEM solutions obtained with a commercial code, despite the fact that cg-FEM has been fully implemented in MATLAB.
1. Introduction
Researchers have devoted many efforts to solve Boundary Value Problems (BVPs) that in one way or another are relevant to different disciplines. The advent of computers allowed the researchers to develop methods to obtain numerical solutions of this type of problems. One of the most common methods used to solve BVPs is the Finite Element Method (FEM) which, nowadays, is broadly used for industrial applications. In the FEM the BVP is discretized subdividing (meshing) the problem domain into small subdomains (elements) of simple geometry. The original FEM strongly relies on a mesh that must conform the geometry of the domain to be analyzed and must be adapted to the solution of the problem, for example, using -adapted meshes where elements of smaller size are used in regions where the solution is less smooth. According to [1], a study at Sandia National Laboratories revealed that the process of creating an analysis-suitable geometry and the meshing of that geometry require about 80% of overall analysis time, whereas only 20% of overall time is devoted to the analysis itself. One way to decrease this 80% of overall time prior to the analysis is to make the geometry of the mesh used for the analysis independent of the domain to be analyzed. This approach can considerably reduce the time devoted to prepare an analysis-suitable model and mesh the domain and is especially useful in applications that would require continuous remeshings during the analysis: structural shape optimization problems, wear modelling, and so forth.
There have been several variations of the original FEM that followed this path to improve the efficiency of the FEM. Two of these improvements of the original FEM are the Extended Finite Element Method (XFEM) [2, 3] developed by Belytschko and his group at Northwestern University (USA) and the Generalized Finite Element Method (GFEM) [4, 5] developed by Melenk and Babuška at the University of Texas at Austin (USA). XFEM is mainly devoted to the analysis of inclusions or cracks. It uses the Partition of Unity Method (PUM) [6] to introduce enrichment functions to represent the displacement discontinuity between the crack faces, the singular fields around the crack tip, and the geometrical description of the crack by means of the Level Set Method (LSM) [7]. The method improves the accuracy of the results and is particularly interesting in crack growing problems as the mesh can remain unchanged when the crack evolves. GFEM [4] follows a similar approach also based in the PUM to include enrichment functions to describe the known behavior of the solution. In GFEM the mesh used for the analysis can be independent of the geometry. For example, a Cartesian grid is used in the GFEM III implementation described in [4]. Both XFEM and GFEM require the use of an integration mesh, purely for integration purposes, in the elements cut by the boundary to take into account the part of the element actually lying within the domain. Another approach to represent the geometry for nonconforming meshes is to use the LSM for boundary representation [8]. Other authors also combine the LSM for boundary representation with the XFEM to represent a solution with gradient discontinuities into an element containing more than one material [9, 10].
Another variation of the FEM, developed to improve its performance, consists of using an auxiliary domain which embeds the problem domain . In this approach, the discretized domain is instead of the problem domain. Generally is a domain with a simple geometry that can be easily meshed. This technique is therefore closely related to the GFEM. The analysis methodology presented in this paper is based on this idea. In our implementation for 2D problems is a square whose discretization into squared quadrilateral elements of uniform size is trivial. These techniques have been used both in the Finite Volume Method (FVM) and in the FEM. According to [11] the literature of these techniques starts in the 60s when VK Saul’ev published, in Russian, the paper. Solution of certain boundary-value problems on high-speed computers by the fictitious-domain method (Sibirsk. Mat. Z. 1963.4: 912-925). Subsequently it has been applied in various fields, such as acoustics [12–14], fluid dynamics and fluid-structure interaction [15, 16], biomedical problems [17, 18], convection-diffusion [19], and optimization [20–23]. The present work will only focus on the FEM framework.
These techniques have several names in the literature, such as Fictitious Domain [11, 12, 19, 20, 24], Implicit Meshing [25], Immersed FEM [26], Immersed Boundary Method [15, 27], and Fixed grid FEM [28, 29]. They have been described in [30] under the name Finite Elements in Ambient Space.
Since the mesh is not conforming with the geometry, these methods require the information of the problem domain to be available during the evaluation of element integrals. The accuracy of the results provided by these techniques depends on the accuracy of the integration process. Hence, the methodology proposed in this paper includes an efficient integration procedure which would be even able to consider the actual boundary providing the exact element integrals (up to the accuracy of the numerical integration and round-off errors).
One major difference of these methodologies with the standard FEM is the consideration of the Dirichlet boundary conditions as in the proposed methodology; in the general case, there are no nodes lying on these boundaries. A procedure based on the use of the Lagrange multipliers technique has been used to apply these boundary conditions.
All these previous approaches are mainly interested in decoupling the geometry representation from the FE mesh where the solution is interpolated. Generally, in these last techniques, the computational cost is concentrated along the elements falling on the boundary. The same occurs in the Boundary Element Method (BEM) where only the external boundary is considered for the analysis. In this framework, recent works of Simpson et al. [31] and Scott et al. [32] where the Isogeometric Analysis IGA [33] with NURBS and T-Splines, respectively, have been adapted to the BEM. Thus, the geometry and the BEM mesh are strongly coupled in this approach.
In this paper we will consider the FE solution of linear elasticity problems, where the Zienkiewicz and Zhu [34] (ZZ) error estimator in energy norm is commonly used to quantify the accuracy of the numerical solution. The information provided by the ZZ error estimator at element level can be used to improve the model by means of -adaptive procedures. The ZZ error estimator is widely used due to its robustness, accuracy, and easy implementation but also because it requires the evaluation of a recovered stress solution , more accurate than the raw FE solution , that can be used instead of . From all recovery techniques we can highlight those based on the Superconvergent Patch Recovery (SPR) technique [35] also developed by Zienkiewicz and Zhu. The publication of the original SPR technique was followed by several works aimed at improving its quality [36–38]. Ródenas and coworkers also proposed several improvements [39–42] of the original SPR in order to yield a highly accurate recovered field which locally fulfills the equilibrium and compatibility equations. These techniques allowed us to obtain the first procedures to get practical upper error bounds for FEM and XFEM based on recovery techniques [40, 42] instead of using the traditional residual-based error estimators [43].
Reference [25] indicates the following: “Unfortunately, for an implicit mesh it would be very difficult to implement such a superconvergent recovery scheme of the stress field for elements that intersect the boundary.” However in the XFEM framework, where the mesh is independent of the crack, efficient recovery techniques have been already proposed based on the Moving Least Squares (MLS) technique [44–47] and some on the SPR technique [42, 48], which introduce worthy improvements to the solution specially along the boundaries, even in elements trimmed by the crack. These SPR-based techniques have been adapted in this work to the context of Cartesian grids.
According to some authors [11, 25, 27, 28], the main drawback of the use of the techniques under the large umbrella of finite elements in ambient space is the low accuracy along the boundaries since they are not explicitly represented. In the proposed methodology, the recovery techniques developed by Ródenas and coworkers [39, 40, 42, 48] have been specially adapted both to be used in the Zienkiewicz and Zhu error estimator [34] that will guide the -adaptive refinement process to improve the quality of the solution and to neutralize the possible lack of accuracy along the boundaries in the cg-FEM framework providing an enhanced solution (for which an error estimator is available) that will be used instead of the FE solution.
In this contribution we have applied cg-FEM to a structural shape optimization analysis. This type of processes will benefit from the computational efficiency and accuracy of cg-FEM but also from the data structure that will allow for sharing information between different geometries analyzed during the process, further improving the behavior of the optimization process.
The paper is organized as follows. The problem to be solved is defined in Section 2. The main characteristics of the proposed cg-FEM methodology are presented in Sections 3 and 4, respectively, devoted to the use of Cartesian grids independent of the geometry for the FE analysis and to the -adaptive procedures specially adapted to the use of Cartesian grids. Numerical examples will be used in Section 5 to evaluate the performance of the cg-FEM methodology.
2. Problem Statement
Let us consider the 2D linear elastic problem. Let be the displacement field defined in , the solution of the boundary value problem: where and , with and , are two disjoint boundary parts where the Neumann and Dirichlet boundary conditions are applied. are the body loads, are the tractions along , and are the prescribed displacements along . The weak formulation of this problem reads as follow: find such that where is the standard test space for the linear elasticity problem and where is the Hook’s tensor and and are the stress and strain operators, respectively.
This problem will be solved numerically by means of the cg-FEM methodology described in this paper, which includes -adaptive refinement capabilities. The use of Cartesian grids and a hierarchical data structure [49] improves the manipulation information and reduces the computational cost of the analysis.
3. On the Use of Cartesian Grids for Efficient FE Structural Analysis
In the standard FEM the mesh must conform to the boundary of the domain. This requirement, together with the need to adequately refine the mesh in order to accurately represent the behavior of interest whilst maintaining the geometry of the elements as undistorted as possible, is the origin of the high cost (both in terms of computing time and analyst’s man-hours) of the process required to generate an adequate FE model. In cg-FEM the mesh does not need to conform the geometry. In other words, we will deal with two different domains, which is the domain to be meshed, a square in our case, and the domain of the problem to be solved . The only requirement between these domains is that .
3.1. Generation of the Analysis Mesh
cg-FEM is based on the use of a sequence of uniformly refined Cartesian meshes where hierarchical relations between the different mesh levels have been defined. This sequence of meshes used to discretize is called the Cartesian grid pile (see Figure 1(a)) which embeds all the problem domain and it is formed by bilinear (Q4) or biquadratic (Q8) squared elements of uniform size. A hierarchical data structure for -adaptive FE analysis based on element splitting was presented in [49]. This data structure took into account the hierarchical relations between the elements of different levels, obtained during the element splitting process, to speed up FE computations. The data structure has been adapted to the particular case of a sequence of meshes given by the Cartesian grid pile where all elements are geometrically similar to the element used in the coarsest level (level 0) of the Cartesian grid pile, called the reference element. One of the main benefits of the data structure is that, as described later in Section 3.3.1, in the linear elastic case all elements of the Cartesian grid pile will have the same stiffness matrix that will be evaluated only for the reference element and shared with the rest of elements in the pile, making the evaluation of element stiffness matrices trivial. This and other hierarchical relations considered in the data structure allow for a simplification of the mesh refinement process and the preevaluation of most of the information used by the FE code, remarkably influencing the efficiency of the code.
(a) Cartesian grid pile
| (b) Analysis mesh. nodes with multipoint constraints for continuity |
The code uses functions that directly provide the nodal coordinates, mesh topology, hierarchical relations, neighborhood relations, and so forth, in an efficient manner, when required. Therefore, there is no need to store this information in memory, considerably reducing memory usage.
The first step of the analysis consists of creating the analysis mesh used to obtain the FE solution of the problem. This mesh is built taking a set of nonoverlapped elements of different sizes (see Figure 1(b)) taken from the different levels of the Cartesian grid pile. A maximum difference of 1 refinement level is allowed between adjacent elements in the analysis mesh. Multipoint constraints (MPCs) [50, 51] are used to enforce continuity between adjacent elements of different levels.
3.2. Geometry-Mesh Intersection: Integration
The analysis mesh is formed by three element types (see Figure 2):(i)boundary elements: elements placed along the domain boundary, . Only a part of each of these elements remains into the problem domain . The evaluation of the stiffness matrix of these elements requires evaluating the intersection between the geometry and the sides of the elements. These intersection points will later be used in order to detect the element area situated into the domain,(ii)internal elements: elements fully located into the domain. These are treated as standard FE elements. These elements have the same stiffness matrix as the reference element as all of them are geometrically similar to this element, (iii)external elements: elements fully located outside of the domain and therefore not considered in the analysis.
The bibliography shows several methods to evaluate the intersection between the domain and the boundary elements and to perform the domain integrals in these elements [16, 28, 29, 52]. The process used in this work is similar to that shown in [29] and consists of three steps:(i)intersection of boundary with element edges: Figure 3(a) shows the intersections of the curves that define the boundary of the domain with the Cartesian Grid element edges. The colored zone represents the problem domain, (ii)addition of intermediate points: as shown in Figure 3(b) we will identify some extra points placed on the curves that define the boundary. The number of these extra points is related with the curvature of the boundary. This set of points together with the vertex nodes of the element are used to create a Delaunay tessellation that defines integration subdomains at each boundary element,(iii)selection of internal triangular subdomains: Figure 3(c) represents the integration subdomains selected to evaluate the domain integral.
(a) Geometry and problem domain, intersection with elements
(b) Intermediate point addition and Delaunay’s tessellation
(c) Integration subdomain generation
For each boundary element this process generates (a) a discretization of the element for domain integration and (b) a discretization of the boundary useful for boundary integration. Thereby, the integration in those elements will be performed using these integration subdomains by means of where is the domain of element , is the number of triangular integration subdomains in placed into the problem domain, and is each of the integration subdomains. The curved boundaries can be approximated with straight line segments in these triangular subdomains when Q4 elements are used. Higher order elements would require higher order approximations to the boundary to avoid modeling errors that would reduce the theoretical convergence rate of the FE analysis. Alternatively it is also possible to use the transfinite mapping techniques commonly used in -adaptive analysis [53] to consider the exact geometry to increase the accuracy of the results. The use of this mapping increases the computation cost per subdomain but reduces the number of triangular subdomains. Our numerical experience shows that in order to maintain the theoretical convergence rate for elements of order , at least -order polynomials must be used to approximate the boundary.
Boundary integration will be performed using a similar procedure, discretizing the boundary within element into integration subdomains. The parametric definition of the boundary provides the exact boundary representation of which is used for boundary integration instead of approximations like the straight line segments represented in Figure 3:
3.3. Element Data Sharing
The use of Cartesian grids together with the data structure used in the implementation of the cg-FEM methodology allows for simple and efficient information data sharing between elements that considerably reduces the total amount of calculations thus improving the computatcional efficiency of the FE code. This section describes, for the internal and boundary elements, how information is shared between elements during the analysis.
3.3.1. Internal Elements
Reference [49] showed that, in linear elasticity, the terms used to evaluate the stiffness matrix ( matrix of derivatives of shape functions and Jacobian matrix ) of geometrically similar elements are related by a constant value evaluated as a function of the scaling factor between the elements. In fact, the stiffness matrices of geometrically similar elements are simply related by a factor where is the number of spatial dimensions. Therefore the evaluation of the stiffness matrices of all the internal elements of the analysis mesh is trivial as, for constant material properties, all these elements share the same stiffness matrix, which will be evaluated only for the reference element (element used to define the coarsest level, level 0) and then shared with the rest of the internal elements through the hierarchical data structure. This implies a major increase in efficiency of the generation of the system of equations to be solved by the FE code.
3.3.2. Boundary Elements
The hierarchical data structure used enables the use on boundary elements of the so-called vertical data sharing and horizontal data sharing described next.
-Adaptive Refinement Process: Vertical Data Sharing. As previously explained, the evaluation of the stiffness matrix of all internal elements of the analysis mesh is trivial. As each boundary element is trimmed differently, each of these elements will require a particular evaluation of the element matrices, following the procedure exposed in Section 3.2. It could be said that the computational cost of the generation of the FE model for the analysis is a function of the number of boundary elements, that is, -dimensional.
In many -adaptive FE codes the previous meshes are discarded and completely new meshes are created as the -adaptive analysis evolves, thus preventing the reuse of information evaluated in previous meshes. In our case, the use of Cartesian grids together with the hierarchical data structure allows for reusing calculations performed in previous meshes. The hierarchical data structure provides the so-called vertical data sharing by means of which elements present in different meshes of the -adaptive process will not be reevaluated for the newer meshes.
As an example, Figure 4 shows two consecutive meshes obtained during the -adaptive analysis of a gravity dam. Note that elements colored in blue are present in both meshes. The vertical data sharing allows for the reuse in the finest mesh of the information evaluated in the coarsest mesh. Note that in the finest mesh new element matrices are only evaluated for white elements.
(a) Coarse mesh
(b) Fine mesh
Shape Optimization Problems: Horizontal Data Sharing. The use of optimization technologies for structural shape optimization of mechanical components has increased its use in industry during the last years. The main drawback of these iterative technologies is their high computational cost, which prevents their widespread use. Shape optimization algorithms are composed by two different levels. The higher level is in charge of generating the geometries to be considered during the iterative process. On the other hand, the lower level is in charge of analyzing each of the geometries proposed by the higher level. A numerical method is used to obtain an approximation to the exact behavior of each geometry. In our case we will use cg-FEM in this lower level because of the benefits in computational cost obtained when evaluating each of the different geometries but also because data can be shared between different geometries through the so-called horizontal data sharing to further improve the overall computational efficiency of the optimization process.
The optimization problem can be formally defined as follows: given a decision space (search space) , an objective space (objective values) , the objective function and a set of constraints :
In the particular case of structural optimization the objective function (OF) is, normally, the weight of the component, are the design variables (e.g., coordinates of control points) that define the geometry, are the constraints expressed, normally, in terms of displacements or stresses, and and define the side constrains.
There are different optimization algorithms that can be used in the higher level to create the different geometries to be analyzed during the optimization process. The benefits of the use of the cg-FEM methodology could be obtained with any optimization algorithm that would require the use of the FEM to analyze the different geometries during the iterative process. In our case, in the numerical examples, we have considered the use of the genetic algorithm (GA) proposed by Storn and Price [54]. More precisely, we use the Differential Evolution (DE) algorithm, version DE1.
With the traditional FEM it is almost impossible to enable the exchange of information between elements of different geometries because, in general, the elements of different meshes are different and completely unrelated as each geometry requires a different mesh conformal to the boundary. However if we use cg-FEM considering the same Cartesian grid pile for all the geometries to be analyzed, we will be able to relate elements used in different geometries making it possible to define a process for horizontal data sharing, that is, between elements of different geometries.
Note that the parametric definition of the boundary of the components to be analyzed can be subdivided into two parts:(i)the fixed part: this is the part of the boundary that remains fixed in all the geometries (such as the external boundary and the lower straight segment of the internal boundary of the gravity dams represented in Figure 5), (ii)the moving part: this is the part of the boundary that would be modified by the optimization algorithm (such as the curved part of the internal boundary of the gravity dams in Figure 5).
(a) Geometry i of the optimization process
(b) Geometry j of the optimization process
(c) Geometry i of the optimization process, an h-adaptive mesh
The horizontal data sharing consists in reusing the computations performed over the elements intersected by the fixed part of the boundary in the different geometries analyzed during the optimization process. Figure 5 shows an example of this horizontal data sharing. This figure shows two different geometries and analyzed during the iterative process. h-adaptive analysis is used to obtain an accurate solution for each geometry as the low accuracy results would negatively affect the performance of the optimization process. In this case green elements represented in Figures 5(a) and 5(b) for geometry are reused in geometry represented in Figure 5(c). Observe that the horizontal data sharing implies a significant reduction of calculations as the information required for most of the boundary elements used in geometry was already evaluated in geometry . The only element matrices evaluated for the analysis of geometry are those corresponding to the white elements.
4. h-Adaptive Refinement
The cost of FE analyses of complex structural components can be reduced by means of the use of -adaptive techniques. These techniques can provide an adequate sizing of the elements adapted to the characteristics of the problem. The use of these techniques can provide the required accuracy in the solution with optimized FE models where the number of elements has been minimized, thus reducing the computational cost of the analysis. In some cases, like in structural shape optimization problems, -adaptive analysis is a must because inaccurate FE results can negatively affect the behavior of the optimization algorithm [55], leading to nonoptimal solutions, reducing the convergence rate to the optimal solution or even preventing convergence.
cg-FEM implements two -refinement procedures. As in [56], the first one is based exclusively on geometrical criteria, whereas the second one considers the quality of the solution and is based on the minimization of the error in the energy norm.
4.1. Geometrical Refinement
The analysis starts adapting the Cartesian Grid dimensions domain () to the problem domain in order to ensure that is embedded into . A preliminary mesh (not used for FE analysis) of uniform element size defined by the user is created as the first step of the analysis process. This preliminary mesh is then intersected with the problem domain. The first analysis mesh is then created following a refinement process based on the geometry of the domain. This procedure consists in refining the boundary elements where the curvature of the boundary is too large with respect to the element size. A simple curvature indicator is defined in (7), where, as represented in Figure 6, the values of represent the distances between points over the boundary and a straight line segment of length defined by the intersection of the boundary with the element sides. is the number of those intersection points. The refinement process is repeated until the relative curvature indicator is small enough. The first analysis mesh is created as a result of this process. Figure 7 shows examples of preliminary meshes and first analysis meshes:
4.2. Solution-Based Refinement
After the FE solution of the first analysis mesh has been obtained, new meshes are created following a refinement procedure that takes into account the quality of the FE solution. This procedure aims to minimize the error in energy norm of the solution. In order to estimate the error in energy norm, we use the Zienkiewicz and Zhu (ZZ) error estimator (8), presented in [34], where is the FE stress field and is an improved stress field recovered from :
Particularizing (8) to each element domain, we would obtain the estimation of the error in energy norm at element level. With that information, and using a mesh optimization criterion based on the equidistribution of the error in the elements of the mesh to be created [57], we obtain the new levels (sizes) of the elements in each zone. Examples of analysis meshes obtained by this procedure are represented in Figure 8.
The numerical results presented in Section 5 will show that, as previously exposed, the use of the cg-FEM increases the computational efficiency of the FE analyses. On the other hand, some authors have reported that the main drawback of this technique is the low accuracy along the boundaries [25, 28]. Loon et al. [52] propose to slightly modify the shape of those elements over the boundary to fit them, generating a mesh conformal to the boundary. One of the benefits of the use of the ZZ error estimator is that the process requires the evaluation of a recovered solution . Thus, our approach to improve the solution, not only on boundary elements but on all the elements of the analysis mesh, consists in substituting the raw FE solution (obtained by cg-FEM) by the recovered solution obtained during the error estimation process to guide the adaptive process. The recovery technique used by cg-FEM is based on the Superconvergent Patch Recovery (SPR) scheme presented and developed in [35, 58–60] by Zienkiewicz and Zhu. The proposed technique is called SPR-CD and provides a recovered displacement field that is then used to obtain . A stress-based version of this technique was presented initially for the traditional FEM in [61], adapted for XFEM in [41, 42] and finally also adapted for cg-FEM in [62]. The displacement-based technique was first presented in [63]. The recovery process produces a recovered displacement field that satisfies the Dirichlet and Neumann boundary conditions and the internal equilibrium equation in each patch of elements attached to the vertex nodes. A continuous recovered field is then obtained using a partition of unity approach as described in the previous references. The recovered solution obtained with the SPR-CD technique is very accurate both along the boundary and in the interior of the domain.
The proposed SPR-CD recovery method is based on the minimization at each patch of elements (elements connected to a vertex node) of the difference between the recovered displacement field and the FE solution , according to (9): where with , and are the vectors of unknown coefficients to be evaluated. Minimizing (9), integrating numerically, and combining the two displacement components, we end up with the linear system of equations shown in (10): where is the determinant of the Jacobian matrix that provides the mapping between the derivatives of the physical and reference element and is the weight associate to each integration point.
The equilibrium and Dirichlet constrains are imposed via the Lagrange multipliers technique through matrix in the system of equations shown in (11) used to evaluate the unknown coefficients . Once have been evaluated, the recovered displacements can be evaluated at each patch . The recovered stresses will be evaluated at each patch from , yielding a solution pair , that is, both statically and kinematically admissible at each patch:
Once the solution is obtained at each patch, it is evaluated for the whole domain , (12), by using the Conjoint Polynomial Enhancement presented by Blacker and Belytschko in [37] to ensure continuity of both fields:
4.3. Efficiency of the Recovery Procedure
A maximum difference of one level is allowed by the -adaptive process between two adjacent elements. Because of this and the topological features of the Cartesian grid, only a reduced number of possible patch configurations can be obtained. Figure 9 represents the 19 possible configurations of internal patches for 2D. Note that internal patches are defined as the patches composed by internal elements only.
The polynomial coefficients used to describe the recovered displacement field according to (11) are obtained for a normalized coordinate system. Then, matrix in (11) will be exactly the same for all internal patches having the same configuration. This implies that we will only need to invert a maximum of 19 matrices to obtain the recovered field in all the internal patches. The first step of the recovery process for the internal patches consists in codifying the configuration of each patch and classifying them according to the configurations shown in Figure 9; then the values of are evaluated for each of the different patch configurations. Once the terms in (11) have been evaluated for each patch, the unknown coefficients are evaluated from . This procedure considerably reduces the computational cost associated to the evaluation of internal patches. The computational cost associated to the evaluation of internal patches is negligible with respect to the cost associated to the evaluation of patches that contain boundary elements as each of these patches will have different matrices. In practice, the computational cost of the recovery process is only depending on the number of patches along the boundary. That implies a -dimensional computational cost.
4.4. Error Estimation of the Recovered Solution: Stopping Criterion
With this novel recovery procedure we have an improved solution more accurate than the raw FE solution . The recovered solution can therefore be used as the standard output provided by cg-FEM instead of the FE solution. Then the accuracy of the recovered solution needs to be evaluated. This means that an error estimator for the recovered solution should be available. An error estimation of the recovered field was first introduced in the report [64]. This report shows that the error in energy norm of the recovered solution can be estimated from where is the exact error of the recovered displacement field, is the estimated error of the FE displacements, represents the lack of internal equilibrium of the recovered stresses, and represents the lack of equilibrium on Neumann boundaries, where is the matrix that projects the boundary stress field into tractions over the boundary. Reference [42] showed that as the recovery technique produces an almost exact satisfaction of the boundary equilibrium equation (), then, the second integral in (13) (boundary integral) is negligible when compared with the first integral (domain integral). Therefore, only the first integral in (13) will be considered for the error estimation of the recovered solution.
Figure 10 shows a scheme of the evolution of the exact errors in energy norm of the finite element and recovered solutions. Note that in SPR-based recovery techniques the recovered solution is in practice more accurate than the FE solution and has a higher convergence rate. Therefore, the number of degrees of freedom required by the recovered solution to obtain a prescribed accuracy level defined by blue horizontal line in Figure 10 is considerably smaller than that required by the FE solution (). We can thus define a highly efficient -adaptive refinement process that considers the accuracy of the recovered solution instead of the accuracy of the FE solution. The -adaptive refinement process of cg-FEM will, as usual, be guided by the error estimation of the FE solution using (8). The error reduction of the FE solution simultaneously produces an error reduction of the recovered solution that is estimated using (13). The refinement process will run until the estimated error of the recovered solution reaches the prescribed accuracy level. As a resume, the error of the FE solution is used to drive the -adaptive process whereas the error of the recovered solution is used to define the stopping criterion.
5. Numerical Examples
This section will present numerical examples to evaluate the improvements provided by cg-FEM.
5.1. Problem 1: Cylinder under Internal Pressure
Let us consider a thick wall cylinder subjected to an internal pressure and plain strain conditions, as described in Figure 11. The analytical solution of this problem is known and is also given in Figure 11. Only a quarter of the section has been considered in the model together with the appropriate symmetry conditions.
Figures 12 and 13 show the evolution of the errors in energy norm and the error convergence rate, respectively, for bilinear (Q4) and biquadratic (Q8) elements. In the graphs on the left of these figures we observe that the exact and estimated errors of the FE solution (blue and red lines) almost overlap. We can see that the cg-FEM results converge smoothly to the exact solution and that our recovery-based error estimator provides very accurate results. This graph also shows that the brown and black curves that represent the evolution of the exact and estimated errors of the recovered solution also overlap, clearly presenting a higher accuracy and convergence rate than the FE solution. The convergence rates of the exact error for the FE and recovered solutions are shown in the graphs on the right of the figures. The theoretical convergence rates of the exact error in energy norm of the FE solution with respect to the number of degrees of freedom are 0.5 and 1.0 for Q4 and Q8 elements. In all cases the recovered solution has presented higher convergence rates than the FE solution.
(a)
(b)
(a)
(b)
The accuracy of the error estimation techniques for the FE and recovered solutions can be clearly evaluated by means of the use of the effectivity index represented in Figure 14 for Q4 and Q8 elements and for the error estimators of the FE and recovered solutions. The effectivity index is defined as the ratio of the estimated to the exact errors in energy norm. Thus, values close to 1 would indicate accurate error estimators. This figure clearly shows that both error estimators are very accurate. Although the error estimator in energy norm of the FE solution is more accurate than that of the recovered solution, the error estimator of this last solution presents a very good behavior and can be used to effectively evaluate its accuracy.
(a)
(b)
Regarding the computational cost, we have used this problem to compare cg-FEM with the commercial code ANSYS 11 considering uniform meshes. This could not be taken as a strict comparison because the commercial code works with a more complex data structure. However, our purpose here is to show that the computational cost of the cg-FEM is comparable with that obtained by widely used commercial codes. In Figure 15 we show how the computational cost is always lower for the proposed code than that for the commercial code. On the right hand side we show the speedup obtained. It is worth to mention that the cg-FEM code is fully implemented in Matlab 2010b and no subroutine has been compiled. As represented in this figure the speedup obtained with our cg-FEM implementation is about 1.5. In addition, we also have to consider that the commercial code is postprocessing the solution with a nodal averaging technique, whereas cg-FEM uses an advanced recovery process and evaluates two error estimations, one for the FE solution and other one for the recovered solution. This is a numerical evidence that the cg-FEM methodology will be highly competitive with a standard FEM compilation.
(a)
(b)
5.2. Problem 2: Flywheel
This problem with a more complex geometry represents a flywheel under tangential tractions along the external surface and constrained displacements in the internal surface as shown in Figure 16. The material is aluminium with elastic modulus GPa and Poisson coefficient . Plane stress conditions are considered.
Figures 17 and 18 present the set of the three first -adapted meshes for the analysis of this component with Q4 and Q8 elements. The first mesh of each set has been obtained using the geometrical refinement criterion. The following meshes have been created with the -adaptivity process based on the error estimation in energy norm.
This problem has no analytical solution; however we can check the behavior of the cg-FEM code considering the convergence rate of the error of the FE solution [45–47]. We have to check that, in the asymptotic range, this convergence rate tends to the theoretical value ( for Q4 elements and for Q8). We can verify this in Figures 19 and 20 where the evolution of the estimated error and its convergence rate have been evaluated for the FE and recovered solutions. We can observe that the convergence rate of the FE results tends to theoretical values in both cases (blue lines). Furthermore, the estimated error of the recovered solution (red lines) presents a smooth behavior and a higher (and stable) convergence rate than the FE solution along the refinement process.
(a)
(b)
(a)
(b)
5.3. Structural Shape Optimization
5.3.1. Problem 3: Pipe Optimization
To check the performance of the cg-FEM combined with the optimization software, we have chosen an optimization problem with known analytical solution. Our objective in this problem is to minimize the cross-sectional area of a tube under internal pressure applied on the circular internal surface, with unknown external surface, where the Von Mises stresses must be below the yield stress .
As shown in Figure 21, two symmetry surfaces are considered. The external surface is given by a cubic spline defined by three points. The optimal analytical solution corresponds to a circular external surface. For and the optimal external radius is and the corresponding area is .
30 individuals per generation have been considered in the genetic algorithm used for the optimization. We have performed analyses with Q4 elements for different prescribed error levels to study the influence of the error in energy norm on the computational cost and the accuracy of the solution provided by the optimization process. We have considered two types of -adaptive analysis techniques to evaluate the numerical solution of each geometry.(i)Strategy a: we have used the FE solution as the output of the analysis together with a criterion to stop the adaptive process based on , see (8). (ii)Strategy b: we have used the recovered solution and a criterion to stop the adaptive process based on the error estimator for this solution ; see (13).
Figure 22 presents the results of this analysis. The graph on the left shows the evolution of the exact relative error in area (evaluated from the exact known analytical solution for this problem) of the best solutions provided by the process versus the number of generations. We can observe that the optimization process will provide accurate solutions only for the most restrictive error levels. Although we used the same values for both strategies, the best results are obtained for Strategy b ( curves). The reason for this could be that, given an FE solution and a recovered solution with similar global errors in energy norm, the maximum value of the Von Mises stress (used as a control parameter during the optimization process) is usually very accurate in the recovered solution as it will appear on the boundary where this solution has been forced to satisfy the boundary constraints. The graph on the right shows the evolution of the computational cost along the optimization process. It can be clearly observed that the use of Strategy b that takes into account the recovered solution to drive the optimization process results in a computing time considerably lower than the time required with Strategy a, driven by the FE solution. The speedup obtained with Strategy b is of about 3.3 for and 3.1 for . Thus, this speedup together with the considerable increase in the accuracy displayed in the graph on the left of the figure clearly shows that the use of Strategy b, that is, the use of the recovered solution to drive the optimization process, produces remarcable benefits with respect to the use of the FE solution.
(a)
(b)
6. Conclusions
In this work, we have presented a novel numerical technique, cg-FEM, to solve Boundary Value Problems. This technique, which has been implemented to solve 2D linear elasticity problems, is based on the use of -adaptive Finite Element Analysis and has three main characteristics:(1)use of Cartesian meshes, nonconformal to the boundary of the domain, to decrease mesh burden, focusing the computational cost of the mesh generation process on the domain’s boundary which must be intersected with the so-called boundary elements; (2)use of a hierarchical data structure that allows reusing previous computations, reducing to the minimum the total number of calculations; and(3)prominent use of a recovered solution obtained with the SPR-CD technique, for which an efficient error estimator in energy norm is available, with two main objectives: (a)to drive the -adaptive analysis based on recovery-based error estimation techniques,(b)to replace the raw FE solution , commonly used as the standard output of the analysis, in order to improve accuracy and reduce the size of the FE model required to reach the prescribed accuracy level.
The numerical results have shown that the FE solution of cg-FEM smoothly converges to the exact solution with (approximately) the theoretical convergence rate. Furthermore, the recovered solution, as a result of being based on the Superconvergent Patch Recovery technique, is more accurate and converges with a higher convergence rate than the FE solution. An error estimator for the recovered solution, whose accuracy is only supported by numerical tests, has been used. Further investigations to provide a theoretical basis for this error estimator are now in process.
The cg-FEM technology has allowed for decreasing the complexity of the problem to one less dimension in the vast majority of the routines. As a result basic numerical tests have shown that even its Matlab implementation can be computationally more efficient than commercial codes based on the traditional FEM. A 3D implementation of cg-FEM under development has already provided promising preliminary results.
The numerical results have shown how shape optimization processes can benefit from the use of cg-FEM, leading to a substantial reduction in computational cost that would help to increase the use of optimization processes in industrial applications.
Acknowledgments
This work has been developed within the framework of research project DPI2010-20542 of the Ministerio de Economía y Competitividad (Spain). The financial support of the FPU program (AP2008-01086), the funding from Universitat Politècnica de València, and Generalitat Valenciana (PROMETEO/2012/023) are also acknowledged. The authors also thank the support of the Framework Programme 7 Initial Training Network Funding under Grant no. 289361 “Integrating Numerical Simulation and Geometric Design Technology."