key: cord-0978582-p7q6nffs
authors: Reddy, J.N.; Nampally, Praneeth
title: A Dual Mesh Finite Domain Method for the Analysis of Functionally Graded Beams
date: 2020-06-27
journal: Compos Struct
DOI: 10.1016/j.compstruct.2020.112648
sha: ea252d16ada1a196b14c9a837d667b3e7189ccae
doc_id: 978582
cord_uid: p7q6nffs

A method that employs a dual mesh, one for primary variables and another for dual variables, for the numerical analysis of functionally graded beams is presented. The formulation makes use of the traditional finite element interpolation of the primary variables (primal mesh) and the concept of the finite volume method to satisfy the integral form afforganizationof the governing differential equations on a dual mesh. The method is used to analyze bending of straight, through-thickness functionally graded beams using the Euler–Bernoulli and the Timoshenko beam theories, in which the axial and bending deformations are coupled. Both the displacement and mixed models using the new method are developed accounting for the coupling. Numerical results are presented to illustrate the methodology and a comparison of the generalized displacements and forces/stresses computed with those of the corresponding finite element models. The influence of the coupling stiffness on the deflections is also brought out.

Laminated composite structures have sudden change in material properties from one layer to another layer (see Reddy [1] ). If two dissimilar materials are bonded together, there is a very high chance that debonding will occur at interface due to some extreme loading conditions, may it be mechanical or thermal load. Another problem in a laminated composite structure is the presence of residual stresses due to the difference in coefficients of thermal expansion between different material layers (see, e.g., [2] ). These problems can be resolved by gradually varying the volume fraction of the constituents rather than abruptly changing them over an interface.

A gradual variation of the material results in a very efficient structure tailored to suit the needs.

Structural elements made of functionally graded materials (see, e.g., [3] and [4] ) are a class of composite structures that have a gradual variation of material properties through the structural thickness (see, e.g., [5] - [12] ). Functionally graded materials (FGMs) are often used in thermal barrier structures for aerospace applications as well as in space structures. In fact, functionally graded material characteristics are present in most structures found in nature (e.g., sea shells, bones, etc), and perhaps a better understanding of the highly complex form of materials in nature will help us in synthesizing new materials (the science of so called "biomimetics"). A typical thermal barier plate or shell is made of a mixture of ceramic and metal for use in thermal environments. The ceramic constituent of the material provides high-temperature resistance due to its low thermal conductivity. The ductile metal constituent, on the other hand, prevents fracture due to high temperature gradient in a very short period of time. The gradation in properties of the material reduces thermal stresses, residual stresses, and stress concentration factors.

A two-constituent functionally graded through-thickness materials are charcaterized by a power-law variation, among other types of variations, of the modulus of elasticity while the Poisson ratio is kept constant. If the z-coordinate is taken along the thickness of the beam, plate, or shell, the modulus E(z) of an FGM along the thickness coordinate z is assumed to be represented by the simple power-law as (see Praveen and Reddy [7] and Reddy [1] )

where E 1 and E 2 are the material properties of the material 1 (top face) and material 2 (bottom face) of the beam, respectively, and n is known as the power-law index. Note that when n = 0, we obtain the single-material structure (with modulus E 1 ).

With the advent of computers a number of numerical methods have been developed during the last four decades. Among them the finite element method (FEM) and the finite volume method (FVM) have gained the most popularity, with the FEM dominating structural mechanics and coupled problems and the FVM dominating fluid dynamics. Typically, all approximate methods convert a differential equation described by the operator equation Au = f governing a variable u to a set of algebraic equations of the matrix form Ku = F, among the nodal values of the variable u and its dual variable F at a selected number of points (called nodes) in the domain and on its boundary. The duality pairs (u, F ) (or cause and effect; e.g., force and displacement; temperature and heat) exist in all fields of science and engineering. While a physical process may have more than one duality pair, the dualities are unique and no polygamy exists among the duality pairs (i.e., one variables is only dual to the other variable in the pair and not to a variable in another pair). The actual process that results in the final matrix equation Ku = F differs from one method to another. The FEM is based on the following three-fold idea [13] :

(1) the total domain Ω can be represented as a collection of a finite number of nonoverlapping but interconnected (at the boundaries of the) subdomains, called finite elements, Ω e ; the elements are of a particular geometry that allows the construction of approximation (or interpolation) functions;

(2) over each element Ω e , the dependent unknown u is interpolated through a set of points (nodes) of the element as u ≈ u j ψ j , u j being the value of u at the jth node and ψ j are suitable approximation functions, and the governing equation is converted to a set of algebraic equations K e u e = F e (called finite element model) using a method of approximation (e.g., weak-form Galerkin or Ritz, subdomain, least-squares, and so on); the element equations contain nodal variables from only the element under consideration; and

(3) the element equations from all elements are put together (element assembly) using balance and continuity conditions at element interfaces to obtain a global set of algebraic equations, Ku = F, which are then solved after applying the boundary conditions at the nodes.

We note that in the FEM, the domain (i.e., element) used for the approximation of dependent variable(s) and the satisfaction of the governing equation(s) in a weakform sense is the same. Of course, the mesh used for the approximation of the geometry can be different from that used for the solution variables, but in the isoparametric formulations, both meshes are the same.

In the FVM [14, 15, 16] one represents a given domain, much like in FEM, as a collection of non-overlapping domains, called control volumes. Then an integral (not a weighted-integral) statement of the governing equation, after invoking the Green-Gauss theorem to covert the domain integral to the boundary integral, is used over a typical control volume to derive the algebraic equations. In the FVM, at the centroid of each control volume lies a mesh point, and the derivatives of the dependent variables at the control volume interfaces are calculated in terms of the values of the dependent variables at the mesh points using Taylor's series approximations (i.e., "finite difference-like" approximations). Thus, there is no explicit interpolation (although there is an polynomial approximation implied by the truncated Taylor's series) of the dependent variables is employed in the FVM. The algebraic equations derived using a typical control volume involve mesh point values from the neighboring control volumes (a notable difference from FEM), naturally connecting the control volumes. The resulting algebraic equations resemble more like finite difference stencils, which are valid for a typical mesh point in the entire domain and include contributions from the neighboring mesh points to obtain the required algebraic equations of the entire mesh. Thus, in the FVM there is no formal assembly of control volumes is involved. The imposition of gradient type boundary conditions involves, sometimes, fictitious nodes from outside the domain, and there is no unique methodology followed (i.e., differs from an author to author) for the imposition of boundary conditions or the evaluation of integral expressions in the FVM. The major advantage of the FVM is the satisfaction of the global form of the governing equations exactly and thus resulting in a better accuracy for the secondary variables like fluxes and forces.

Recently, Reddy [17] introduced a numerical approach termed the dual mesh finite domain method (DMFDM) which uses ideas from both the FEM and FVM to solve second-order equations. In the DMFDM, the domain is represented with a primal mesh of finite elements, and a dual mesh is superimposed on the primal mesh such that the nodes of the primal mesh are at the center of the dual mesh of finite domains, expecpt for the nodes on the boundary. Then the governing equation is required to be satisfied in an integral sense (not a weighted-integral sense) over the finite (control) domain. The second-order terms in the differential equation are integrated by parts and expressed as quantities (dual variables) on the interfaces of the dual mesh. When the interfaces fall on the boundary, either the dual variables or their counterparts (i.e., primary variables) are known. Thus, the dual mesh finite domain method brings the best features of the FEM, namey, the interpolation of the variables and imposition of physical boundary conditions, and the satisfaction of the actual balance equations over the control domain of the FVM. The major merits of the DMFDM are that the method inherits the desirabe aspect of the FVM (in satisfying the global form of the governing equations exactly) and overcomes the disadvatage of the discontinuity of the secondary variables at the interfaces of the fnite elements by calculating them at the boundaries of the control domains, where they are continuous (i.e., uniquely defined). More details of the method are presented in the sequel.

The DMFDM method was introduced in [17] and illustrated with applications to heat transfer problems with a single unknown (also, see [18] ). In the present paper, the method is extended to multivariables problems involving functionally graded beams. In particular, the DMFDM is applied to the Euler-Bernoulli and Timoshenko beams with three dependent variables; displacement and mixed formulations involving second-order differential equations are developed and discretized using the DMFDM.

Following this introduction, a review of the governing equations of the Euler-Bernoulli and Timoshenko beam theories (see Reddy [19] ) as applied to functionally graded beams is presented in section 2. The usual equations in terms of the general-ized displacements as well as those derived in terms of the generalized displacements and bending moment (in place of the rotation) are presented. The finite element and dual mesh finite domain discretizations of these equations are presented in section 3. Several of these models are not readily found in the literature. Then, in section 4, numerical solutions are presented for pinned and clamped (at borth ends of the) beams, and the numerical solutions obtained with FEM and DMFDM models are compred with the exact solutions to assess the relative accuracy in predicting the generalized displacements and generalized forces. Finally, some remarks about DMFDM and its extensions are outlined in section 5. The linearized (i.e., small strains) equations of equlibrium of the EBT are (see Reddy [19] ):

where f (x) and q(x) are the axial and transverse distributed loads, respectively, on the beam, c f is the modulus of the foundation (if any) on which the beam rests, and N xx and M xx are the stress resultants defined by and expressed interms of the generalized displacements u and w, with θ x = −dw/dx as:

Here A xx , B xx , and D xx are the extensional, extensional-bending, and bending stiffness coefficients

We note that Eqs. 

To express the governing equations in terms of u, w, and M xx , we first rewrite Eqs.

(3a) and (3b) in the form (i.e., displacements in terms of stress resultants)

Equations (2a), (2b), and (7b) constitute the governing equations in terms of u, w, and M xx :

The equations of equilibrium of the TBT are

The stress resultants (N xx , M xx , Q x ) in the TBT can be expressed in terms of the generalized displacements (u, w, φ x ), where φ x denotes the rotation of the cross section about the y-axis, as

where S xz denotes the shear stiffness and K s is the shear correction factor, which is taken to be K s = 5/6. Substitution of Eqs. (10a) and (10b) into Eqs. (9a)-(9c) yield the governing equations in terms of the displacements: 

Then we rewrite Eqs. (10a) and (10b) for displacements in terms of stress resultants,

Differentiating Eq. (12c) with respect to x, and solving for dφ x /dx,

Equations (12a), (12b), and (14c) constitute the governing equations in terms of u, w, and M xx for the TBT:

3 Discretized Equations 3.1 Finite element models

The displacement finite element model of Eqs. (5a) and (5b) is based on their weak forms:

where Q i are the generalized forces at the element nodes (see Fig. (1)) defined by

Lagrange interpolation of u(x) and Hermite interpolation of w(x) are required with a minimum of linear for u(x) and Hermite cubic for w(x). The details of the finite element model are included in Appendix A. 

where

Here (Q e 1 , Q e 2 , Q e 3 ) and (Q e 4 , Q e 6 , Q e 6 ) denote the axial force, transverse shear force, and rotation at nodes 1 and 2, respectively (see Fig 2) . Clearly, the weak forms in Eqs. (19a)-(19c) admit Lagrange interpolations as low as linear for (u, w, M xx ). The details of the finite element model are presented in Appendix A. 

The displacement finite element model is based on the weak forms of Eqs. (11a)-(11c):

where Q i are the generalized forces at the element nodes (see Fig. ( 3)) defined by

The weak forms indicate that Lagrange interpolation of (u, w, φ x ) are required. The details of the finite element model are included in Appendix A. 

The DMFDM is best suited to discretize second-order equations. Therefore, we can only consider the mixed model of the EBT using Eqs. (8a)-(8c). In the DMFDM, we divide the domain Ω = (0, L) into a set of N finite elements separated by nodes, as shown in Fig. 4 , with each node having its own finite domain (a dual mesh). The nodes are numbered sequentially from the left to the right. We consider two adjacent elements connected at a typical interior node I and finite domain associated with that node (see Fig. 5 ).

1 I h - I x D Control domain 1 I - I 1 I + I h 1 2  1 N +  1 I x - ( ) I a x ( ) I b x x x L = ( ) N W

Finite element, Next, we derive the discretized equations associated with Eqs. (8a)-(8c). The complete steps of the DMFDM are presented by considing Eq. (8a) and then summarize results for the remaining equations and also for the other models described here.

The first step is to set up the integral statement of Eq. (8a) over the control domain:

where A and B refer to the left and right end locations of the control domain (associated with node I), which have the coordinates x (I)

a and x (I) b , respectively (see Fig.  6 for the nodal degrees of freedom). Unlike in a weighted-residual method (or weak form), we weaken the differentiability on u by carrying out the indicated integration and obtain 0 =

Here N Next, we use finite element approximations of u(x), w(x), and M xx (x) over a typical finite element, Ω (I) = (x I , x I+1 ). For example, u(x) is approximated using linear interpolation

where U I is the value of u at node I (i.e., U I ≈ u(x I )) and ψ 

Hence, we can calculate parts of N x

The integral of a function f (x) over the control domain (x

) can be evaluated using either exact integration or numerical integration (e.g., trapezoidal rule, onethird Simpson's rule, and so on).

Lastly, we write the discretized equations for the boundary nodes [see Fig. 8 and Appendix B]:

wherex is the local coordinate with origin at node 1 of element N and N are the boundary forces (at nodes 1 and N + 1, respectively), which are either specified or their duality counter parts, namely, the displacements U 1 and U (N +1) , are specified. This completes the discretization of Eq. (8a). The same procedure can be applied to Eqs. (8b) and (8c) to obtain the discretized equations for the interiror and boundary nodes. We have the following integral statements of Eqs. (8b) and (8c):

where V (I) 1

and V (I) 2

denote the secondary variables (shear forces acting upward positive) at the left and right interfaces of the control domain centered at node I,

and Θ (I) 1 and Θ (I) 2 denote the secondary variables (rotations in counterclock direction) at the left and right interfaces of the control domain centered at node I,

The discretized equations associated with Eq. (8b) at an interior node are:

where C (I) denotes the value of c f in element I. For the boundary nodes 1 and N + 1, we have

The discretized equations associated with Eq. (8c) are

for an interiror node. Here D I denotes the value of D xx in element I andB I denotes the value of B xx /D xx in element I. For the boundary nodes 1 and N + 1, we have

This completes the development of the discretized equations based on the DMFDM for the mixed formulation of the Euler-Bernoulli beam theory.

In order to derive the discretized equations associated with Eqs. (11a)-(11c), we follow the same procedure as described for the EBT. The integral statements over the Ith control domain centered around node I, occupying the domain between points A and B (see Fig. 9 for the nodal degrees of freedom) for each of these three equations are:

Here (N 

2 ) denote the axial forces, shear forces, and bending moments at the left and right interfaces, respectively, of the control domain centered at node I (see Fig. 9 ). Since the displacement model of the TBT suffers from shear locking, we evaluate the integral appearing in Eq. (45a) (i.e., the shear force) as a constant to avoid shear locking:

The discretized equations associated with Eqs. (11a)-(11c) for an interiror node I are

The discretized equations of the left boundary node are

For the node on the right boundary, we have

Lastly, we develop the DMFDM models of Eqs. (15a)-(15c). Without repeating the steps, which were amply illustrated in the preceding sections, we present the final discretized equations here. The discretized equations associated with Eq. (15a) at the Ith node, node 1, and node N + 1 are given by

The discretized equations associated with Eq. (15b) are [see Eqs. (37)-(39)]:

Finally, the discretized equations associated with Eq. (15c) at the Ith node, node 1, and node N + 1 are:

Here We investigate the effect of mesh and the power-law index n on the deflections and stresses. We consider a beam of length L = 100 in, b × h = 1 × 1 in 2 crosssectional dimensions, functionally graded through the height (h) (E 1 = 30 × 10 6 psi, E 2 = 3 × 10 6 psi, and ν = 0.3), and subjected to uniformly distributed transverse load of intensity q 0 lb/in. For the pinned-pinned and clamped-clamped boundary conditions, we can exploit the symmetry about x = L/2, and use the left half of the beam as the computational domain.

First, we consider the beam with both ends pinned. The boundary conditions on the primarly variables in various models for this problem are as follows: 

where ξ = x/L and

We note that for homogeneous beams u(x) = 0 everywhere. The EBT solutions are obtained from Eq. (55) by settingD xx = 0 and replacing φ x with −dw/dx. The bending stress, σ(x, z), is computed at x = L/2 (where the bending moment is the maximum) and z = h/2, h being the beam height. The stress is post-computed in the displacement models at the element center using the relation

On the other hand, the stress in the mixed models is computed using the bending moment M xx (x) according to the formula

where M xx (x) is the calculated from the finite element interpolation (i.e., σ(L/2, h/2) = M (L/2)h/2I, and M (L/2) is the nodal value). Tables 1 and 2 contain the normalized center deflection and stress, respectively, for homogeneous (D xx = D xx = EI andB xx = B xx = 0) pinned-pinned (P-P) beams for different number of elements in the half beam. The tabulated deflections and stresses are normalized as follows (with K s = 5/6):

where I is the moment of inertia. From the results it is clear that the mixed dual mesh finite domain models are the most accurate in predicting the displacements and stresses. In fact, all mixed models and the displacement model of the DMFDM give the exact stress for any number of elements. The shear locking is alleviated in the displacement finite element model and DMFDM of the TBT by the use of reduced integration. No such trick is used in the mixed FEM and mixed DMFDM. The DMFDM always gives the exact value of the moment at x = L/2. We also note that for this slender beam (L/h = 100), the effect of shear deformation is negligble and the EBT and TBT solutions forw are the same upto the fourth decimal point. Table 3 contains the normalized center deflection for functionally graded pinnedpinned (P-P) beams for different values of the power-law index n. All of the results are obtained using 16 elements in the half beam. All models predict solutions that match the exact solutions upto the fourth decimal point. The stresses in the FGM beams are exactly the same as those in the homogeneous beams, because the bending moment is independent of the stiffness coefficients. It is interesting to note that the effect of the power-law index n on the deflections is not monotonic. As n goes from zero to a value of about n = 2, the deflection decreases and then increases for n > 2, as shown in Fig. 10 . This is due to the fact that B xx is not a monotonically increasing or decreasing function of n (where M denotes the moduli ratio, M = E 1 /E 2 and B 0 = bh 2 ), as can be seen from Fig. 11 (see [20] ). Figure 12 shows the variation of the normalized center defelcetionw with the power-law index n. Next, we consider a beam clamped (C-C) at both ends. The boundary conditions on the primarly variables in various models for this problem are as follows (replace dw/dx with φ x for the TBT: The exact solutions for clamped-clamped beams are given by (ξ = x/L)

Tables 4 and 5 contain the normalized center deflection and stress, respectively, for the clamped-clamped (C-C) homogeneous beam. From the results it is clear that the mixed dual mesh finite domain models and finite element results are very close, if not identical. The displacement finite element model is the most accurate by virtue of the higher (Hermite cubic) approximation of the deflection. Table 6 contains the normalized center deflection for the functionally graded clamped-clamped beam for different values of n. All of the results were obtained with 16 elements in half beam. For the C-C beams the normalized delfections do not deviate significantly from each other (they differ only in the fourth or fifth decimal point). Figure 13 shows the variation of the normalized deflectionw versus n, which has the same form as that for the P-P beams. 

The dual mesh finite domain method (DMFDM) proposed by Reddy [17] makes use of two different meshes: one mesh of finite elements for interpolation of the variables and another mesh of control domains to satisfy the integral form of the governing equations. Thus, the DMFDM contains the merits of the finite element method (FEM) as well as the finite volume method. An advantage of the dual mesh finite domain method is that the secondary variables are computed at the finite domain interfaces, where the derivatives are continuous and thus they are uniquely determined.

One of the shortcomings of the finite volume and dual mesh finite domain method is that they cannot be used for solving differential equations of order higher than 2, unless they are recast as second-order differential equations. Of course, most problems of engineering are first-order when they are first derived using balance laws, and they become higher order only when some of the variable are eliminated using auxiliary (kinematic and/or constitutive) relations. In this paper, the DMFDM is extended to such multi-variable problems, using the equations governing the bending of functionally graded (FGM) beams. Mixed models (i.e., models that involve displacements and moments) of the Euler-Bernoulli and Timoshenko beams are formulatated using the FEM and DMFDM. The mixed models have the advantage of enabling stress calculation at the bounday nodes, instead of Gauss points in the interior of the domain. Numerical results indicate that the DMFDM gives very accurate results, especially for the bending moments, which are used to the compute stresses. A study of the FGM beams also revealed that the dimensionless bending deflections (w = wD xx /q 0 L 4 ) are not monotonic functions of the power-law index n because the coupling stiffness B xx is not a monotnically increasing or decreasing function of the modulus ratio.

Extensions of the DMFDM to multivariable problems in multidimensions, initialvalue problems, and nonlinear problems are awaited. In addition, application to laminated composite beams and plates and flows of viscous fluids, and a study of the mathematical aspects such as the existence, uniqueness, and error estimates of the DMFDM require attention.

Dedication. This paper is dedicated to the memory of Professor José Fernández Sáez of the University Carlos III of Madrid, Spain, who lost his life to COVID-19. Professor Fernández Sáez, affectionately known as "Pepe", was a scholar, passionate teacher, and kind human being. No words can entirely remove the pain of losing a friend like Pepe, but he lives in our memory as a great human being and an unconditional friend of his friends, students, colleagues, and collaborators.

The displacement finite element model of the EBT is of the form

and Q i are the generalized nodal forces, as shown in Fig. 1 ; ϕ I are the Hermite cubic functions and ψ i are the linear Lagrange interpolation functions; u denotes the vector of nodal displacements associated with the linear approximation of u(x); and ∆ denotes the nodal displacements (transverse deflection and rotation at each node) associated with the Hermite cubic interpolation of w(x).

The mixed finite element model of the EBT, based on the Lagrange interpolation of all variables, is given by

i (x b ) Q 6 (A4)

Here (ψ (1) i , ψ

i , ψ

i ) are the Lagrange interpolation functions used for (u, w, M xx ), respectively. In general they are different from each other, but here we took them to be the same for all variables.

The displacement finite element model of the TBT is of the form

i (x a )Q 3 + ψ

i (x b )Q 6 (A6)

The mixed finite element mode of the TBT has the same form as the mixed finite element model of the EBT, Eq. (A3), and the coefficients K αβ ij and F α i (α, β = 1, 2, 3) also remain the same as those in Eq. (A4), except for the following coefficient:

Mechanics of Laminated Composite Plates and Shells: Theory and Analysis

Hygro-thermomechanical modelling and analysis of multilayered plates with embedded functionally graded material layers

Enhanced Thermal Stress Resistance of Structural Ceramics with Thermal Conductivity Gradient

The Concept of FGM

Steady thermal stresses in a functionally gradient material plate

Thermomechanical Analysis of Functionally Graded Cylinders and Plates

Nonlinear Transient Thermoelastic Analysis of Functionally Graded Ceramic-Metal Plates

Analysis of Functionally Graded Plates

Exact solution for thermoelastic deformations of functionally graded thick rectangular plates

Thermal stresses in functionally graded beams

A nonlinear modified couple stress-based third-order theory of functionally graded plates

Nonlinear analysis of functionally graded microstructure-dependent beams

An Introduction to the Finite Element Method

Computational Fluid Dynamics, The Finite Volume Method

Computational Methods for Fluid Dynamics

Numerical Methods for Partial Differential Equations. Finite Difference and Finite Volume Methods

A dual mesh finite domain method for the numerical solution of differential equations

A dual mesh finite domain method for steady-state convection-diffusion problems

Energy Principles and Variational Methods in Applied Mechanics

An Introduction to Nonlinear Finite Element Analysis

Acknowledgements. The research reported herein, in parts, was supported by the Oscar S. Wyatt Jr. Endowed Chair and a project from Technical Data Analysis (TDA) to Texas A&M University.

 26) and (27), we can calculate integrals over a finite domain and values at the finite domain interfaces. The following formulas are used in the development of the DMFDM discretized equations, where the coefficients a and c are assumed to be element-wise constants: