Abstract
Direct modeling and simulation of engineering problems with various irregularities are computationally very inefficient and in some cases impossible, even in these days of massively parallel computational systems. As a result, in recent times, a number of schemes have been put forward to tract such problems in a computationally efficient manner. Needless to say, such schemes are still going through evolutionary stages. This paper addresses direct solution based on the selective use of different dimensional models at different regions of the problem domain. For the multidimensional approach, a higher-order transition element is developed to connect the different element types where two- and three-dimensional laminated elements based on higher-order subparametric concept are considered. Modeling simplicity and calculation efficiency of the multidimensional approach are shown for the analysis of cantilever plates with stepped section and patch-repaired plates.
1. Introduction
A finite element method is one of versatile numerical tools to solve some differential equations which model physical phenomena of various engineering problems. In last few decades, lower-order finite elements have been conventionally used to solve a wide range of practical problems. However, these elements are unable to provide accurate distributions of stress resultants in structures with free edges and stress singularities. Even though one uses a highly refined mesh to model problems having significant stress gradients, the accuracy of stresses obtained by using lower-order finite elements is rather poor [1]. The quest for robust finite elements has triggered researchers to develop higher-order finite elements. The latest developments in this field indicate that its futures lies in adaptive higher-order methods, which successfully respond to the increasing complexity of engineering simulations and satisfy the overall trend of simultaneous resolution of phenomena with multiple scales. The mathematical justification concerning the advantage of using a higher-order approximation was shown by Babuska et al. [2] who showed high accuracy, high convergence rate with coarser meshes, and improved performance in handling stress singularity problems. Additional advantages are that a few large elements selected to satisfy primarily the geometric domain, loading, and boundary requirements can cover an entire computational domain and ease effort of solution convergence. Thus, the user’s time is saved significantly. Past works of the finite element methods using the higher-order approximation were reviewed by Szabo et al. [3] and were recently carried out by Ahn et al. [4–8].
Meanwhile, some examples of irregularities in a structure, say, in plates, can appear as step change in thickness, presence of reentrant corner due to skewness, cutouts of various shape, cracks, and notches, free edges of laminates, local damage involving plasticity, strain-hardening, patch repair at damage locations, and so forth. The primary objective of multidimensional approach is to achieve computational efficiency by accurately capturing the local three-dimensional behavior using, as needed, three-dimensional models and to use lower dimensional model away from such locations. This requires special modeling as zone of transition between the two different dimensional models. In the case of composite materials, the concept of weight-averaged apparent constituent properties is used. As needed, orthotropic and transversely isotropic properties are considered. In the case of patch-repaired plates, the repair can be bonded or riveted and the patch material can be same as the parent material, or more conveniently, composite material bonded to the area with damage in the form of, say, holes, cracks, and so forth. Physically the patch repair work is yet another source of irregularity over and existing irregularities caused by the damage. Furthermore, the situation gets more complex due to the potential for interfacial failure in patch repair leading to debonding. Further consideration of the first model will, therefore, be focused with respect to patch-repaired problems. In addition, the scheme will consider layer by layer modeling of laminated composites in the context of both two and three dimensions.
For more efficient analysis in terms of solution accuracy and computational efficiency, techniques using sequential methods or methods based on a combination of different mathematical models have been proposed. For sequential methods [9–13], one of the main criticisms is that the equilibrium of forces along the boundaries between different models is not maintained while the displacement continuity is ensured. Also, another disadvantage of these methods is the problem of incorporation of nonlinearity into the analysis. Unlike the sequential methods, simultaneous mixed methods combine different mathematical models to analyze the entire computational domain, including the use of distinctly different levels of space- or function-refinements in specified local regions. Here different subregions with different mathematical models are explicitly accounted for and are, thus, easily amenable to nonlinear analysis. One simple scheme of simultaneous application of mixed methods to composite laminate analysis is the concept of selectively grouping the plies in the vicinity of the location where accurate stress values are desired [14–18]. Another scheme [19] to apply simultaneous mixed methods is to use multipoint constraint equations or Lagrange multipliers. Here the variational statement is supplemented with additional terms to enforce compatibility between adjacent subregions. Although the Lagrange multiplier approach can be used to couple subregions which use different mathematical models, the method to connect different mathematical models has rarely been used, because the implementation would be very cumbersome. The most popular approach for connecting different mathematical models (like simultaneous two-dimensional to three-dimensional modeling of plates and shells) is to implement special transition elements [20–23]. Those researches adopt finite element mesh using conventional lower-order functions. The performance of those finite elements has not been satisfactory in predicting smooth and accurate variations of stresses as well as displacements for plates with various irregularities. Therefore, this leads to a finite element formulation having superior continuity of shape functions.
Over the years, different -refined two-dimensional modeling schemes for laminates based on classical and first-order shear versions of equivalent single-layer theories based on assumptions regarding kinematics of deformations or stresses, higher-order theories assuming third- and higher-degree polynomial planar variations of displacements [24–26] and layerwise models in which the variation of displacements over the thickness of each layer is assumed, sometimes, leading to the zig-zag effect [27, 28]. Additional schemes have been put forward to optimize the computational process [29–31]. These schemes tried to provide transition between dissimilar finite element meshes in two and three dimensions. Modeling schemes based on -refinement were put forward by some researchers [31–35]. Contrary to -version of finite element analysis using low order basis functions of Lagrange type and requiring highly refined mesh adjacent to singular points like crack tips, in -version of finite element analysis, basis functions are of hierarchical Legendre type or integrals thereof with model refinement being done in the functional space by appending higher-order basis functions, leading to highly accurate solutions with fewer degrees of freedom [5–8, 32, 34, 36].
In this study, discussion for this multidimensional analysis will be confined to the higher-order approximation based on Lobatto shape functions with hierarchical properties. In many industrial applications the structures with laminated composite materials are composed of three-dimensional solid continuum with the thin shell-like portions connected to them. When the structures are modeled by function refined mesh technique to reduce negative effect from various irregularities, this multidimensional analysis proposed in this work can be considered. Also, this proposed method can be attempted to distribute limited computational resources in an optimal manner to achieve maximum solution accuracy with minimal solution cost, subject to certain problem-specific constraints.
2. Lobatto Shape Functions
One-dimensional shape functions are classified into two groups as shown in Figure 1. In the one-dimensional element on standard domain, two nodal modes are same to the linear Lagrange interpolation functions as follows:and additional modes often called internal shape functions, internal modes, or bubble modes follow the Lobatto shape functions like the following:where is the Legendre polynomials as follows:
In (3), Legendre polynomials take the values 1 and −1 at and −1, respectively. Therefore, they cannot qualify as hierarchical shape functions, since the values of all shape functions except nodal modes must be zero at two end points (). That is why (2), not only (3), is used to construct hierarchical shape functions. Also, the orthogonality property of Legendre polynomials in (3) implies where refers to the Kronecker delta. The use of shape functions with orthogonality can avoid round-off errors usually associated with polynomials of high degree and minimize coupling between hierarchical degrees of freedom. In addition, a more dominant diagonal form of the stiffness matrix can be obtained.
Two-dimensional shape functions can be built from the one-dimensional shape functions defined above. They are divided by three groups as shown in Figure 2. The four nodal modes are as follows:In any -level, edge modes on four sides for the two-dimensional standard element can be defined by Here the superscripts refer to edge index. Bubble modes to satisfy completeness of an element are valid for only and can be obtained from a combination of only internal modes on the one-dimensional element like this:Here, superscript 5 refers to bubble modes and subscripts and which are any positive integers must satisfy the inequality given byLike one-dimensional shape functions, two-dimensional shape functions are also hierarchical. That is, in present shape functions, all lower-order shape functions can be contained in the higher-order basis, while in some other -version finite element methods [37–39] all standard shape functions corresponding to each -level must be built each time. Zienkiewicz et al. [38] showed some advantages of the hierarchical shape functions. The concept leads to a stiffness matrix with a dominant diagonal. Also, the condition number of the stiffness matrix is improved by an order of magnitude.
3. Finite Element Formulations
3.1. Three-Dimensional Element with Layerwise Theory
For three-dimensional elements, displacement fields are given bywhere , , , and follow Einstein summation convention. is the number of variables on each edge on two-dimensional plane in , and 4, while in is the number of internal modes on two-dimensional plane. is the number of one-dimensional internal modes across thickness for displacement variation between bottom and top surfaces. Also, superscripts and are indexes for variables on bottom and top surfaces, respectively. The symbols like , , and refer to the variables corresponding to nonnodal modes, say internal modes on one-dimensional element and edge and bubble modes on a two-dimensional element. Stress-strain relationships are based on three-dimensional elasticity theory: Thus, all six strain components with respect to local axes (, , and ) are present as follows:Meanwhile, displacement filed of an element defined in (9) can be written by the following general form:Here, column vectors and refer to nodal and nonnodal variables, respectively. Difference of the two vectors is that a nodal variable itself has physical sense, while a nonnodal variable itself does not. It is important that accuracy of nodal variables can be improved by more nonnodal variables. The shape functions with respect to the nodal variables are given byAdditional shape functions related to nonnodal variables can be built by product of two- and one-dimensional modes such as , , , and shown in (9). Thus, the three-dimensional element equations can be derived using the principle of virtual work. More detailed process is omitted in this paper because it is similar with general formulation of conventional finite elements.
3.2. Two-Dimensional Shell Elements with Equivalent Single-Layer Theory
For two-dimensional modeling, an element is formulated by the dimensional reduction from three-dimensional solid to two-dimensional surface which satisfy the first-order shear deformation plate theory and plane stress theory. Each of in-plane displacement fields, and , is explained by translational and rotational variables. Out-of-plane displacement field is in accordance with that of the RP element. The form of the displacement fields is given byFor stress-strain relationships, the following stress-strain relationship based on plane stress assumption can be obtained like the following:Strains with respect to axes , , and have the form as follows:where superscripts , , and refer to membrane, bending, and transverse shear strains, respectively. Equation (14) denotes the assumption that the strains are continuous through the thickness. The stress resultants are given bywhere , , and are the membrane force, bending moment, and transverse force resultants, respectively. Also, refers to shear correction factor. For a homogeneous cross-section the shear correction factor is equal to 5/6. is the variable of distances away from the representative surface through thickness. and are the positions of bottom and top of the plate with respect to global -coordinate.
3.3. Transition Elements with Lobatto Shape Functions
To carry out connection of two different types of elements, the transition elements are developed. From the aforementioned (9) and (14), the assumed displacement fields of the -convergent transition element can be expressed asHere, and refer to indexes in the choice of nodal mode types corresponding to HP, RP, and 2DS elements for in-plane displacements. For defining nonnodal parts, the variables and are used for in-plane displacements. The form about the variables is as follows:Similarly, for out-of-plane displacements, the variables , , , and are utilized. Thus, the form is given byIn the present element, the validity of including either transversal normal strains or neglecting them primarily depends on its application. For instance, if the transition element is located in junction of a thick solid component connected to a moderately thick shell component, the transversal normal strains should be included in the finite element formulation. On the other hand, if the transition element is in the region where thin solid components are connected to thin shell components, transversal normal strains and stresses can be neglected in deriving the element properties. This requires selection of the proper stress and the strain components at the integration points. Figure 3 illustrates a typical mesh model using the -convergent transition element for a typical case where the three-dimensional elements are connected to a two-dimensional element by using a transition element. The proposed transition elements allow variable kinematic elements to be connected with satisfying interelement compatibilities. In addition, for the proposed transition element with hierarchical shape function, the internal modes follow the most economical mode type among types of vertex and edge modes applied to the element.
4. Numerical Examples Using Multidimensional Analysis
4.1. Stepped Plate Problem
For demonstrating the scheme, the simple problem of a 12 m long × 1 m wide stepped cantilever plate shown in Figure 4 is considered. The thicker part is 2 m long × 0.7 m thick and the remaining part is 0.1 m thick. Two loading cases, uniaxial tension (MPa) and lateral load (Pa), at the free end face are considered separately.
For convenience in selecting different trial finite element models, the problem domain is divided into four segments, as shown in Figure 5, and three different models are considered, as shown in Table 1. It may be noted that, in Segment 1, three-layered model is considered. The degrees of planar and along the thickness polynomial approximations are taken as 8 and 3, respectively.
4.1.1. In-Plane Loading
The center line variation of axial displacement is found to be identical with Models A and C. However, the results by purely 2D modeling (Model B) show discrepancy, as is evident from the plot of Figure 6, shown for . In the upper half of the thickness, the variation of axial displacement in the thicker part just at the step is shown in Figure 7, which tended to become uniform, as the left support was approached. So, Model B fails to represent the true behavior at this location and multidimensional Model C is as good as fully three-dimensional Model A. The same is evident from the normal stress () results in the thinner part at the step, as shown in Figure 8.
4.1.2. Out-of-Plane Loading
In this case, the transverse displacement profile for Cases A and C shows perfect agreement, whereas, Case B tends to underestimate it, especially in the thinner part. It is evident from Figure 9 that the variation of normal displacement at the step across the thickness of the thicker part shows full agreement for Models A and C, whereas, Model B results are grossly erroneous but tend to approach Models A and C results for sections away from the step. In the case of normal stresses too, similar relative performances were noticed. In Figure 10, the normal stress fringes obtained by Model C are shown, near the step region. The extremely large elastic stress values at the corners are indicative of the existence of stress singularity at the corner and in reality the material will undergo damage at this point by yielding and strain-hardening well before reaching the indicated stress levels.
4.2. Patch-Repaired Plate with 15 mm Edge Crack
A uniaxial stretched single-edge-cracked aluminum plate is shown in Figure 11, as repaired by a bonded patch. The size of the parent plate is 200 mm long × 40 mm wide × 1.5 mm thick. The glass-epoxy patch dimensions are 26 mm × 26 mm × 1.4 mm, and the adhesive film is 0.2 mm thick. The material properties are (a) for aluminum GPa, ; (b) for adhesive film GPa, GPa, , , GPa, GPa; and (c) for adhesive film GPa, . The applied remote stress to the aluminum plate is MPa. Taking advantage of one-way symmetry, only a symmetric half the plate is analyzed using a multidimensional model comprising a 4 × 5 mesh, as shown in Figure 12. It may be noted that in the mesh the gray-shaded region depicts three-dimensional elements encircling the crack only and the rest is comprised of two-dimensional elements covering patched and the unpatched regions with transition elements located in-between. In the patched region, the aluminum plate, composite patch, and adhesive film are treated as the three layers of the element. Due to the presence of one-sided patch, the plate undergoes flexure and hence the stress intensity factors vary along the thickness. Consequently, the average value is considered as the representative value to compare it with published values [39] based on conventional (-version) 10,920 three-dimensional solid elements and 40,176 degrees of freedom. Average stress intensity factors as the current results are calculated by -levels of 7 in planar direction and 3 in the thickness direction with 20 elements and 4,868 degrees of freedom based on strain energy release rate technique and is found to be 334.35 MPa which agrees well with published value, 341.84 35 MPa. Also, out of the various possibilities, the present mesh was found to be most efficient.
5. Conclusions
Multidimensional analyses with example problems were presented. In the case of multidimensional modeling three-dimensional and two-dimensional -version finite element formulations with Lobatto shape functions were presented. In order to provide deformation continuity between the two models when multidimensional modeling is undertaken, suitable -version transition elements were developed. These elements were successfully applied to singular problems, namely, one of which is a cantilever stepped plate undergoing extension and flexure and the other is a cracked panel with bonded composite patch repair, leading to accurate results based on significantly simpler models and drastically reduced computational effort as compared to conventional approaches.
Future works would consider that the application of the proposed method will be investigated on problems composed of other composite materials like functionally graded materials which have a continuous variation of material properties from one surface to another. While present work was implemented to show mechanical response on geometrical irregularities and material irregularities of laminated systems, future works on functionally graded materials will consider mechanical response on high temperature environments including thermal shock.
Competing Interests
The author declares that there is no conflict of interests regarding the publication of the paper.
Acknowledgments
This work was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2015R1D1A1A01060909).