Line 319: | Line 319: | ||
|} | |} | ||
| colspan='2' style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(10) | | colspan='2' style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(10) | ||
+ | |} | ||
+ | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" | ||
+ | |- | ||
+ | | | ||
+ | {| style="text-align: center; margin:auto;width: 100%;" | ||
|- | |- | ||
| colspan='2' style="text-align: center;white-space: nowrap;"|<math>{\sigma }_{x}=\frac{{q}_{v}}{2\pi }\left[ \frac{\pi }{2}-\frac{xyz}{B\left( {x}^{2}+{z}^{2}\right) }- {tan}^{-1}\frac{zB}{xy}+\left( 1-2\nu \right) \left( {tan}^{-1}\frac{y}{x}- {tan}^{-1}\frac{yB}{xz}\right) \right]</math> | | colspan='2' style="text-align: center;white-space: nowrap;"|<math>{\sigma }_{x}=\frac{{q}_{v}}{2\pi }\left[ \frac{\pi }{2}-\frac{xyz}{B\left( {x}^{2}+{z}^{2}\right) }- {tan}^{-1}\frac{zB}{xy}+\left( 1-2\nu \right) \left( {tan}^{-1}\frac{y}{x}- {tan}^{-1}\frac{yB}{xz}\right) \right]</math> | ||
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(11) | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(11) | ||
+ | |} | ||
+ | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" | ||
+ | |- | ||
+ | | | ||
+ | {| style="text-align: center; margin:auto;width: 100%;" | ||
|- | |- | ||
| colspan='2' style="text-align: center;white-space: nowrap;"|<math>{\sigma }_{y}=\frac{{q}_{v}}{2\pi }\left[ \frac{\pi }{2}-\frac{xyz}{B\left( {y}^{2}+{z}^{2}\right) }- {tan}^{-1}\frac{zB}{xy}+\left( 1-2\nu \right) \left( {tan}^{-1}\frac{x}{y}- {tan}^{-1}\frac{xB}{yz}\right) \right]</math> | | colspan='2' style="text-align: center;white-space: nowrap;"|<math>{\sigma }_{y}=\frac{{q}_{v}}{2\pi }\left[ \frac{\pi }{2}-\frac{xyz}{B\left( {y}^{2}+{z}^{2}\right) }- {tan}^{-1}\frac{zB}{xy}+\left( 1-2\nu \right) \left( {tan}^{-1}\frac{x}{y}- {tan}^{-1}\frac{xB}{yz}\right) \right]</math> | ||
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(12) | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(12) | ||
+ | |} | ||
+ | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" | ||
+ | |- | ||
+ | | | ||
+ | {| style="text-align: center; margin:auto;width: 100%;" | ||
+ | |- | ||
|- | |- | ||
| colspan='2' style="text-align: center;white-space: nowrap;"|<math>B=\sqrt{{x}^{2}+{y}^{2}+{z}^{2}}</math> | | colspan='2' style="text-align: center;white-space: nowrap;"|<math>B=\sqrt{{x}^{2}+{y}^{2}+{z}^{2}}</math> |
Many of the engineering problems are analyzed using numerical methods such as the finite element (FEM) whose results provide a basis to make basic decisions regarding the design of many important works. It is commonly accepted that FEM computations are reliable; however, the results may be affected by the configuration of the finite element mesh to simulate the medium to be analyzed, this is particularly true when the internal and external boundaries are time dependent, as is the case of soil consolidation. Accordingly, a thorough investigation was carried out with the main purpose of eliminating this shortcoming. The main steps to carried out the development of the innovative geometric procedure to automatically refine finite element tetrahedra-type (3D) are described. This geometric algorithm is based on the theory of fractals and is a generalization of the algorithm for triangular element finite element meshes (2D) [1,2]. This paper presents the fundaments of this new algorithm and shows its great approximation using 3D close form solutions, and its versatility to adapt the original Finite Element Mesh when the load boundary conditions are modified (Neumann conditions).
Keywords: Finite element analysis, remeshing algorithm, adaptive method, civil engineering, geotechnics, numerical analysis
It is acknowledged that the natural basic element to discretize complex 3D continuum media are the tetrahedral-shaped elements. Having this in mind the goal of this paper is to present the development of an algorithm for 3D finite element meshes that draws principles of the theory of fractals.
The finite element method is a powerful, practical and versatile numerical analysis tool for solving differential equations through continuum discretization, but its efficiency and effectiveness can be notably hampered if local or total mesh refinements during the analyses are made manually, which can likely lead to continuum discretization blunders. As it is well known, the accuracy of finite element calculations depends on how well the mesh mimics the continuum, particularly where abrupt changes in the material properties originally existed or developed due to local material plastization of one or a number of fine elements. Likewise, the mesh has to be fine enough to accurately describe the geometry of the problem constituents (i.e., sharp corners of soil-structure problems) and capture abrupt changes in local stress and/or strain gradients stemming from material heterogeneities (i.e., soil-concrete interfaces of embedded foundations).
Finite element approximations improve either by reducing the size of finite elements (refining specific zones of the mesh by automatic remeshing) or increasing the order of the polynomial shape function. (Evidently, a combination of these two approaches is also a viable alternative). The former is used in this paper, by means of the adaptive algorithm presented herein. At present time, several approaches exist to generate and refine locally 3D and 2D finite element meshes [3,4]. Furthermore, some mesh generation techniques such as the Delaunay triangulation [5,6], the advancing front approach [7], mesh generation using contours [8], the coring technique [9], the “quadtree” technique [10], among others, are available for automatically generate gradation triangular-element meshes [11]. Moreover, a number of studies focused on modelling contact regions [12], also to avoid large aspect ratios in order to maintain an adequate approximation level [13] and to compute accurate solutions when large mesh distortions develop during computations [14]. In addition, algorithms for refinement/derefinement of nested meshes using fractal concepts and iterated function systems (IFS) have been proposed [15-18]. Evolution problems have been addressed elsewhere [19]. In addition, the skeleton based refinement (SBR) has been used to refine unstructured meshes for 2D and 3D problems [20]. Dunyach et al. [21] proposed a curvature-adaptive isotropic remeshing technique for real-time mesh deformation. You et al. [22] developed an adaptive remeshing procedure for finite element analysis of heterogeneous materials. Adaptive finite element techniques have been advanced and used for more than twenty years. Adaptive procedures try to automatically refine, coarsen, or relocate a mesh and/or adjust the basis to achieve a solution having a specified accuracy in an optimal fashion. Most known methods of adaptive mesh refinement include: a) r-method (refines locations of nodes), b) h-method (refines an element into multiple elements, and c) p-method (refines polynomial order of element shape functions). Adaptive remeshing implies the adjustment of mesh control parameters based on error estimates or attributes, deleting an existing mesh, and regenerating a new mesh. Adaptive remeshing is a geometry/topology based process and is commonly performed as a pre/post-processing function.
Broadly speaking, any remeshing procedure must include the following features:
(i) A mesh-refining criterion related to an error estimator of the calculation accuracy, or to the local distortion mesh or to the spatial discretization of the boundaries.
(ii) A mesh generator assuming element size distribution inside the domain or around its boundaries depending on the chosen error estimator.
(iii) A transfer procedure that carries out a mapping of history dependent variables from the old to the new mesh for each remeshing step.
The 3D algorithm included herein is a generalization of the 2D algorithm [1], where the triangular elements (2D) are expanded to tetrahedral elements (3D). Such changes imply from the way in which the elements interact, since now there are two types of neighboring elements (face and edge), to the number of elements obtained after partitioning.
The algorithm integrates a number of features arranged in a “software” with modular structure, which makes it a versatile tool capable of 1) Maintaining all problem initial symmetries and 2) simulating nonlinear analysis, using the proper material constitutive model.
To the authors´ knowledge, this algorithm is one of the few available refining procedures considering innovative concepts educed from fractal theory. As has been shown earlier for the 2D case [1,2], the algorithm provides almost identical results regardless of mesh initial geometrical characteristics and of its roughness.
To demonstrate the algorithm’s reliability, the procedure is assessed comparing the stresses it yields with those computed with the closed form solutions of a flexible shallow square slab footing on the surface of a soft clayey homogeneous deposit, considering first only vertical loading and then under vertical and horizontal loading. From the results included herein it becomes clear the importance of adapting the mesh to the ongoing load boundary conditions.
It is worth mentioning that currently there are commercial and non-commercial software very useful to build finite element meshes in 2D and 3D of excellent quality, e.g., Onelab [23], Agros2D [24], Fenics Project [25], but they do not automatically refine the region exceeding the stress-strain behavior thresholds and hence, the user must indicate where the mesh should be refined according to his experience. If the remeshing is carried out manually, mishaps in defining the topology of the small elements are likely to occur, or the refining is not sufficient to obtain reliable results.
As it is known, tetrahedral elements are used to discretize media that can later be analyzed by some numerical method like the Finite Element. As will be shown below, the tetrahedral element allows a region to be refined locally without inducing an important propagation of the refinement outside of the zone being refined. Accordingly, in this work the tetrahedral element is used to refine finite element meshes where needed at different local zones.
In this section some of the geometric characteristics of the tetrahedron are presented as background information to make easier the understanding of developing the parameter quality indicator for qualifying the degeneration of the tetrahedron. This parameter is used to define the skewedness of the tetrahedral elements that integrate the refined finite element mesh.
Let , , , and be points in that are not contained in one plane. Denoting by T the tetrahedron with vertices and (Figure 1), forming the simplest closed convex polyhedron, which has 4 triangular faces and 6 edges. The volume of can be calculated by the following expression:
|
(1) | |
|
(2) |
Defining:
|
(3) |
where is the radius of the circumscribed ball into , where is the surface of the boundary of .
Figure 1. Tetrahedron |
The radius of the circumscribed ball about T can be computed as [26,27]:
|
(4) |
where
|
(5) |
In the above formula and are the Euclidean lengths of opposite edges of for :
, , , , ,
The dihedral angles of a tetrahedron are the six angles between each pair of faces of . They are defined as the complementary angles of outward unit normal to those facets and can be calculated by means of the inner product [27]:
where and are the outward unit normal of particular faces.
It is common in FE analysis to qualitatively distinguish between so-called “well-shaped” (i.e., close to regular) tetrahedra and “badly-shaped” ones (i.e., close to degenerate). Some classifications of badly-shaped tetrahedra are given in [27]. The degree of degeneration of a tetrahedron T is often measured in terms of the quality indicator [28]:
|
(6) |
Tetrahedra with quality indicator near 1 are almost regular, whereas those with near 0 are nearly degenerate.
The method of subdivision of tetrahedral elements proposed in this work consists of three categories: the first establishes the sectioning of a central element, the second consists in the sectioning of a neighbor-face element and the third induces a partition of a neighbor-edge element. Different tetrahedral sectioning techniques can be found elsewhere [11,27,29]. The proposal in this paper makes a combination of three forms reported in the literature [27]: red, green and red-green, for the central elements, neighbor-face and neighbor edge respectively as seen later.
The procedure proposed herein is purely geometric (branded ARFEM3D-Automatic Refinement of Finite Element Meshes-3D; the software may be obtained directly from the senior author) and was inspired by fractal technology, i.e., from a given fractal figure one is capable of subdividing it in as many as desired self-similar figures having much smaller dimensions.
The algorithm uses the Sierpinski’s idea that consists in dividing subsequently a tetrahedral element following specific rules [1,2], as graphically indicated in Figure 2. Herein are presented the key aspects of a new procedure to refine finite element tetrahedral meshes as problem-solving computations proceed. The main difference between the refinement process proposed herein and the Sierpinski tetrahedron is that any tetrahedral element that does not satisfy the threshold criterion is not extracted but refined until the threshold criterion is fulfilled.
Figure 2. Sierpinski tetrahedron generation |
In FE analysis and computation, one needs sequences (infinite or finite) of partitions that have certain properties. One can define various kinds of “well-shapedness”, usually called regularity, in the sense that certain properties of the tetrahedral elements are supposed to hold uniformly over all partitions of the family [27].
Consider the section mesh of Figure 3(a), conformed by three tetrahedral elements. Now assume that the gray-shaded element (nodes 1, 2, 3 and 4) does not fulfill the (i.e., octahedral stress) threshold criterion and hence it should be refined. To this end, it is precise to consider the contiguous tetrahedral elements that delimit the zone to be refined. In Figure 3(a), the area to be refined is defined by numbering the elements and nodal points. It is noted that there are two types of element neighboring the central one: the first is of the "face" type, since it shares one of its four faces with the central element (element in Figure 3, conformed by nodes 1, 2, 4 and 5), and the second is of the "edge" type, because it only touches the central element at one of its six edges (element of Figure 3(b), formed by nodes 2, 4, 5 and 6).
Figure 3. Region to be refined: central element Ec and its neighbors (face type) and (edge type) |
Therefore, when the central element is fragmented, the subdivision of neighboring elements will be induced, according to its type. Refinement of the central element is carried out bisecting its four edges generating eight tetrahedrons (), as shown in Figure 4(b), wherein the generated elements have been contracted and colored to show in detail the fragmentation of the central element . Note that the unrefined central element is now integrated by eight self-similar smaller elements such that . In doing this partition one introduces six nodal points (between and ), (between and ), (between and ), (between and ), (between and ) and (between and ), which are located at mid-distances among existing nodal points , , and , respectively. These new nodal points violate the connectivity rule of the finite element. To overcome this problem, one has to refine locally the mesh by subdividing neighbor elements ( and ) as indicated in the inset of Figure 5, which generates four new tetrahedral elements from and , and two tetrahedral elements from and , which implies that and . It is important to note that refinement of the initial central triangle (Figure 4) compels to refine also all adjacent tetrahedral elements (Figure 5), but the algorithm is properly designed to preclude the refinement be extended to nearby tetrahedral elements. This is obviously advantageous because when the refinement of a large mesh is required at different zones, only the elements at each zone will be involved; hence, the remeshing task is carried out efficiently. However, each of the zones where the integrating elements were refined can be refined further within the constrained zone. In addition, notice that the algorithm produces a kind of element-refining implosion and all new tetrahedral elements are self-similar. At all times it is verified that the volumes before and after the partition do not change, hence the following identities must be fulfilled.
Central element volume:
|
(7) |
Neighbor element volume:
|
(8) |
Neighbor element volume:
|
(9) |
Figure 4. (a) Central element after first iteration. (b) Shrink elements generated |
Figure 5. Central and adjacent elements after first iteration. in blue and in magenta. Central element in gray |
So far, we have considered that there is only one central element. However, there may be situations when it is advantageous to consider two or more adjacent tetrahedral elements as the central tetrahedron from which refinement starts. To accommodate these situations, we introduced the following rules:
Rule 1: If two elements (red and gray contours, in Figure 6(a) share one of their faces, both are considered central elements and remeshing will proceed as if each element was independent. Notice that only side-to-side or face-to-face adjacent elements to are subdivided. Figure 6(b) shows graphically such procedure.
Figure 6. Remeshing considering two adjacent central elements |
Rule 2: If two or more elements (i.e., red and blue tetrahedra, in Figure 7(a) are adjacent to other common neighbor element, all three are considered central as depicted in Figure 7(b) and mesh refinement follows Rule 1. By deducting reasoning, if elements are adjacent each of them can be considered central and hence the refinement Rule 1 applies likewise.
Figure 7. Remeshing considering three adjacent central elements |
It is important to mention that the partitioning of central elements may lead to adjacent tetrahedral elements rather skewed, as portrayed in Figure 9. To avoid the formation of much skewed tetrahedral elements within the mesh; the algorithm includes a tool to analyze all the tetrahedra according to their side ratios. This technique is illustrated in Figure 10. Considering as initial mesh the one shown in Figure 8 and assuming the white element has to be refined, the resulting mesh would be the one indicated in Figure 10 It is seen that sharper elements result from the fractioning the elements next to the central selected element.
Figure 8. Tetrahedral element to be refined |
Figure 9. Elongated elements due to refinement of central element |
The algorithm avoids the formation of highly distorted elements by performing a previous analysis of the degree of degeneration () of all tetrahedral elements. If the specified shape ratio and the element exceed specified thresholds, then these elements are considered central elements as indicated in Figure 10. The procedure follows until none of the new elements violates the element ratio, as graphically depicted in Figure 10 for ratios equal to 0.5. It is worth mentioning that these features improve appreciably poorly designed meshes making sure the results are meaningful. When defining the threshold ratio, it is important to have in mind that the ratio of equilateral tetrahedra is one, and all internal angles are equal among theme, and when bisected, their ratio will be greater or equal to 0.5.
Figure 10. Automatic searching of elongated elements and refinement towards their elimination |
Before any finite element computations are carried out, the algorithm is applied to smooth the initial mesh to eliminate all distorted (sharp) elements as indicated by the block diagram included as inset in Figure 11. The full smooth-refining algorithm depicted in Figure 11 was coded yielding the ARFEM3D software briefly commented next.
Figure 11. Flow chart of the automatic remeshing procedure |
After the initial mesh is smoothed, the automatic refining procedure begins by forming a file containing all the initial characteristics of the mesh and characteristics of the problem itself (i.e., boundary conditions, material properties, present and future boundary conditions, etc.).
Having defined the initial conditions of the problem, a loop where the boundary conditions are modified, is initiated (From to NBoundCond). For each boundary condition, a stress-strain finite element analysis is carried out and at the end of each cycle of this loop, material properties can be modified to account, if required, for nonlinear materials characteristics. Once the new stress state of the full mesh is known, a search of elements to be refined begins according to the threshold initially specified. The element-refining cycle is carried out until none of the elements surpasses the threshold specified.
It is worth stressing the fact that the finite element refining proceeds whenever the threshold is exceeded; to make element-mesh discretization adequate so the problem at hand is properly modeled. This element refinement leads to a smoothed mesh that minimizes the inclusion of spurious stiffnesses and numerical inaccuracies in the finite element computations. As the block diagram of Figure 11 indicates, these tasks are performed within every cycle until the number of iterations defined are completed. Accordingly, every time a new mesh is generated, the new stress-strain state is computed. Once the remeshing cycle is completed, the boundary conditions if required are modified and the above-mentioned procedure is performed all over again.
In this chapter, the algorithm developed and shown in this paper is compared in terms of vertical stresses and vertical displacements along of a number of profiles to the results computed with the corresponding close form solutions. These comparisons show that whenever the problem under analysis is not properly discretized (although at deep scrutineers seemed otherwise) the results more often than not will be misleading and hence the designs might leave much to be desired.
Stresses in an infinite homogeneous ground mass generated by the action of a vertical uniformly distributed load () on a rectangular flexible foundation on the ground surface (Figure 12) can be calculated under a corner of the loaded surface using the following expressions [30,31]:
|
(10) |
|
Published on 12/01/21
Accepted on 16/11/20
Submitted on 02/03/20
Volume 37, Issue 1, 2021
DOI: 10.23967/j.rimni.2020.11.001
Licence: CC BY-NC-SA license
Are you one of the authors of this document?