Finite element modelling of exponentially graded composites with microstructure

Abstract

Composite material properties are dependent on their microstructure. To adequately model these microstructural effects, a formulation of elasticity that accounts for microstructural effects must be considered. Size-dependent behaviour is an inherent property of such materials, resulting in a need for non-classical continuum theories for adequate characterization. A three-dimensional model for an exponentially graded composite with microstructural effects is developed under the framework of Cosserat elasticity. The mixed boundary value problem is formulated and solved through the finite element method, and implemented in the commercial software Abaqus through two user developed elements. Recovery of the classical limit for the Dirichlet problem is discussed as an initial model validation until complete experimental measurements can be conducted.

1. Introduction

Composites are materials in which two or more constituent materials are combined to form a continuum. The resulting mixture possesses material properties different to the individual components, with many composites being tailored for specific applications. Some examples of composites include concrete, rubber, and polycrystalline aggregates [1], with the majority of applications appearing in the automotive, aviation, and construction sectors.

Due to the difference in the properties of the composites constituents, the nature in which they are combined significantly influences the mechanical behaviour of the new material. Functional grading has become a topic of great interest due to the manufacturers ability to design the mechanical properties of the material to suit the desired application. Functional grading involves creating a gradient in the homogenized material properties to achieve the desired behaviour as a function of position. The mixture of the constituents can be altered using various methods, most notably using volume fraction, to achieve the gradient in properties that result in the desired behaviour.

Common functional grades include stepwise, power law, polynomial, exponential, and periodic. Examples of functionally graded composite materials include high performance sports equipment, human bone, and thermal shielding. Composite materials that do not have a functional grading, but rather are a perfect mixture of the constituents, can be considered as functionally graded materials with constant grading. Due to the rising interest in these materials and their applicability to future structural components, there is a fundamental need to understand their behaviour.

Many composites hold a microstructure that is statistically uniform (fibre and matrix), resulting in a macroscopically homogeneous material [2]. Such composite materials typically fail under significant loading conditions through the separation of the fibres and matrix—known as delamination or debonding [3,4]. The delamination is a result of significant shear stresses occurring at the interface where the constituents with differing properties are bonded [5]. Modelling of the debonding processes in such materials is challenging due to interface conditions, typically involving high stress concentrations.

Recently, functionally graded materials (FGMs) have been used in an attempt to solve the issue of debonding due to stress concentrations. By varying their elastic properties with position, FGMs create a smooth gradient that eliminates high stress concentrations and interface conditions [5]. The elimination of these sharp interfaces, where failure is typically initiated, creates a material with desirable properties [4,6] that can withstand extreme conditions (high temperature gradients and stress concentrations) [3,7,8]. FGMs are classified as microscopically heterogeneous materials due to the gradients that dictate their macroscopic properties [4,9]. In practice, the desired properties are obtained by spatially varying the volume fractions of the constituents in a preferred direction.

Since composite materials are heavily dependent on their microstructure, and several experiments have shown an inability of classical theory to capture these effects [10–19], a formulation of elasticity that accounts for microstructural effects must be considered. Adequate characterization of such materials necessitates non-classical continuum theories since size-dependent behaviour is inherent to them.

To attempt to supplement classical elasticity theory and account for experiments in which microstructural effects prevail, the micropolar (Cosserat) theory of elasticity has experienced a re-emergence. Mora and Wass [20] state that Cosserat elasticity is a good foundation for characterizing materials that exhibit these effects, since it includes higher order effects that are not present in classical continuum theories. These effects consider force-like quantities that develop in bodies when strain gradients exist, and consider rotations on the microstructure scale as kinematically independent from displacement. By incorporating these effects, Cosserat elasticity inherently includes a length scale [21] in its formulation, which can account for the observed size effects. This enriched model considers 6-degrees of freedom, the three classical displacements, as well as three microrotations (pointwise orientation). The spatial gradients of these rotations result in internal curvature measures, referred to as twist, which form work conjugates [22] to the resulting internal bending moments.

Furthermore, equilibrium of a material point includes a transfer of moments in addition to axial forces, which gives rise to the aforementioned kinematically independent microrotations. Considering moment equilibrium for a body, the addition of microrotations gives rise to stress-like quantities [23], known as couple stress, which consist of a couple per unit area. The inclusion of couple-stress terms in the moment equilibrium equations result in the stress tensor being asymmetric, contrasting classical stress tensors which can be simplified using Voigt (symmetric) notation.

This paper considers formulating the three-dimensional governing equations of Cosserat elasticity for an exponentially graded material with microstructure. The associated variational problem is constructed, and solved using the finite element method. Exponential grading is considered to demonstrate the formulation and solution procedure; however, any grading method can be used in conjunction with the presented procedure. The solution is then implemented in the commercial software Abaqus [24] by means of a user defined element (UEL). Experimental data containing measurements of Cosserat elastic constants are lacking, and as such complete experimental validation remains elusive until adequate measurements are made. Until such a time when adequate measurements are made, experimental validation is limited to recovery of the classical limit.

2. Preliminaries

Throughout what follows, tensor and index notation are utilized to describe the necessary functions. The tensors used cover three dimensions with indices assuming values 1, 2, and 3. Repeated indices denote a summation from 1 to 3 as per the Einstein summation convention, unless stated otherwise. ${\mathcal{M}}_{\mathrm{mxn}}$ denotes the space of $(m\times n)$ matrices, ${\mathcal{I}}_n$ denotes the identity matrix in ${\mathcal{M}}_{\mathrm{nxn}}$ , and $(\dots ){,}_{\alpha }$ denotes differentiation with respect to the coordinate $x_{\alpha }$ .

Consider a general Cosserat elastic body that occupies a domain $\mathrm{\Omega }$ in ${\mathbb{R}}^3$ , bounded by surface $\partial \mathrm{\Omega }$ . If the body is subject to the action of an external force $\overrightarrow{\mathbf{X}}=(X_1,X_2,X_3{)}^T$ and an external moment $\overrightarrow{\mathbf{Y}}=(Y_1,Y_2,Y_3{)}^T$ , the deformation of the body can be characterized by a displacement vector:

$$\begin{array}{c}\overrightarrow{\mathit{u}}(\overrightarrow{\mathit{x}})=(u_1(\overrightarrow{\mathit{x}}),u_2(\overrightarrow{\mathit{x}}),u_3(\overrightarrow{\mathit{x}}))\hfill \end{array},$$

(1)

and a microrotation vector, containing kinematically independent microrotations about each axis:

$$\overrightarrow{\mathit{\phi }}(\overrightarrow{\mathit{x}})=({\phi }_1(\overrightarrow{\mathit{x}}),{\phi }_2(\overrightarrow{\mathit{x}}),{\phi }_3(\overrightarrow{\mathit{x}})),$$

(2)

where $\overrightarrow{\mathit{x}}=(x_1,x_2,x_3{)}^T$ is a point in ${\mathbb{R}}^3$ . In some cases, it is convenient to group the displacements and microrotations into a vector containing all degrees of freedom, such that

$$\overrightarrow{\mathit{U}}(\overrightarrow{\mathit{x}})={\left(\overrightarrow{\mathit{u}}(\overrightarrow{\mathit{x}}),\overrightarrow{\mathit{\phi }}(\overrightarrow{\mathit{x}})\right)}^T.$$

(3)

One can then extract the states of strain, twist, stress, and couple-stress [25], defined, respectively, as

$$\begin{array}{c}\underset{¯}{\mathit{\gamma }}=\left[\begin{array}{ccc}{\gamma }_{11}\hfill & {\gamma }_{12}\hfill & {\gamma }_{13}\hfill \\ \multicolumn{1}{c}{{\gamma }_{21}}& {\gamma }_{22}\hfill & {\gamma }_{23}\hfill \\ \multicolumn{1}{c}{{\gamma }_{31}}& {\gamma }_{32}\hfill & {\gamma }_{33}\hfill \end{array}\right],\underset{¯}{\mathit{\varkappa }}=\left[\begin{array}{ccc}{\varkappa }_{11}\hfill & {\varkappa }_{12}\hfill & {\varkappa }_{13}\hfill \\ \multicolumn{1}{c}{{\varkappa }_{21}}& {\varkappa }_{22}\hfill & {\varkappa }_{23}\hfill \\ \multicolumn{1}{c}{{\varkappa }_{31}}& {\varkappa }_{32}\hfill & {\varkappa }_{33}\hfill \end{array}\right]\hfill \\ \underset{¯}{\mathbf{\sigma }}=\left[\begin{array}{ccc}{\sigma }_{11}\hfill & {\sigma }_{12}\hfill & {\sigma }_{13}\hfill \\ \multicolumn{1}{c}{{\sigma }_{21}}& {\sigma }_{22}\hfill & {\sigma }_{23}\hfill \\ \multicolumn{1}{c}{{\sigma }_{31}}& {\sigma }_{32}\hfill & {\sigma }_{33}\hfill \end{array}\right],\underset{¯}{\mathit{\mu }}=\left[\begin{array}{ccc}{\mu }_{11}\hfill & {\mu }_{12}\hfill & {\mu }_{13}\hfill \\ \multicolumn{1}{c}{{\mu }_{21}}& {\mu }_{22}\hfill & {\mu }_{23}\hfill \\ \multicolumn{1}{c}{{\mu }_{31}}& {\mu }_{32}\hfill & {\mu }_{33}\hfill \end{array}\right]\hfill \end{array}$$

(4)

Next, the equilibrium equations are formulated in terms of the defined tensors. Considering force and moment equilibrium, one obtains

$$\{\begin{array}{c}{\sigma }_{ji}{,}_j+X_i=0,\\ {{ϵ}_{ij}}_k{\sigma }_{jk}+{\mu }_{ji,j}+Y_i=0,\end{array}$$

(5) (6)

where ${ϵ}_{\mathrm{ijk}}$ is the alternating tensor.

The presence of the couple stresses in the formulation of micropolar media introduces asymmetry to the stress and couple-stress tensors. Furthermore, the asymmetry of the tensors in the governing equations for a micropolar solid results in a more complicated formulation. It is also important to note that in the absence of couple stresses, the above equations reduce to the classical formulation.

From these governing equations, one can extract basic relations between the quantities of interest. Namely, one can derive relations between the deformation tensors (strain, twist) and the corresponding degrees of freedom (displacements and microrotations in ${\mathbb{R}}^3$ ). The relations are

$$\{\begin{array}{c}{\gamma }_{ij}=u_{j,i}-{ϵ}_{ijk\phi k,}\\ {\varkappa }_{ij}={\phi }_{j,i}\end{array}$$

(7) (8)

The constitutive law for a general, anisotropic, and homogeneous micropolar elastic body proposed by Eringen [26,27] is given by

$$\{\begin{array}{c}{\sigma }_{ij}=A_{ijkl\gamma kl}+B_{ijkl\varkappa kl,}\\ {\mu }_{ij}=B_{klij\gamma kl}+C_{ijkl\varkappa kl,}\end{array}$$

(9) (10)

where $A_{\mathrm{ijkl}}=A_{\mathrm{klij}}$ and $C_{\mathrm{ijkl}}=C_{\mathrm{klij}}$ for material characteristic constants $A_{\mathrm{ijkl}}$ , $B_{\mathrm{ijkl}}$ , and $C_{\mathrm{ijkl}}$ .

Considering an isotropic material, $A_{\mathrm{ijkl}}$ , $C_{\mathrm{ijkl}}$ have the same values in every coordinate system and can be re-written as a general isotropic tensor of rank 4 such that

$$\{\begin{array}{c}{A}_{ijkl}=(\mu +\alpha ){\delta }_{ik}{\delta }_{jl}+(\mu +\alpha ){\delta }_{il}{\delta }_{jk}+\lambda {\delta }_{ij}{\delta }_{kl}\\ C_{ijkl}=\left(\gamma +\epsilon \right){\delta }_{ik}{\delta }_{jl}+\left(\gamma +\epsilon \right){\delta }_{il}{\delta }_{jk}+\beta {\delta }_{ij}{\delta }_{kl}\end{array}$$

(11) (12)

where $\mu$ , $\lambda$ , $\alpha$ , $\beta$ , $\gamma$ , and $\epsilon$ are the elastic micropolar constants.

Applying the isotropic tensors to the stresses yields the final form for a linear isotropic micropolar material,

$${\sigma }_{\mathrm{ij}}=(\mu +\alpha ){\gamma }_{\mathrm{ij}}+(\mu -\alpha ){\gamma }_{\mathrm{ji}}+\lambda {\delta }_{\mathrm{ij}}{\gamma }_{\mathrm{kk}},$$

(13)

$${\mu }_{\mathrm{ij}}=(\gamma +\epsilon ){\varkappa }_{\mathrm{ij}}+(\gamma -\epsilon ){\varkappa }_{\mathrm{ji}}+\beta {\delta }_{\mathrm{ij}}{\varkappa }_{\mathrm{kk}},$$

(14)

In an FGM, the elastic constants, and consequently, the associated engineering constants, are functions of position, with the constants in an exponentially graded material written as

$$\mu =\mu (\overrightarrow{\mathit{x}})={\mu }_0{e}^{\overrightarrow{\mathit{G}}·\overrightarrow{x}},$$

(15)

$$\lambda =\lambda (\overrightarrow{\mathit{x}})={\lambda }_0{e}^{\overrightarrow{\mathit{G}}·\overrightarrow{x}},$$

(16)

$$\alpha =\alpha (\overrightarrow{\mathit{x}})={\alpha }_0{e}^{\overrightarrow{\mathit{G}}·\overrightarrow{x}},$$

(17)

$$\beta =\beta (\overrightarrow{\mathit{x}})={\beta }_0{e}^{\overrightarrow{\mathit{G}}·\overrightarrow{x}},$$

(18)

$$\gamma =\gamma (\overrightarrow{\mathit{x}})={\gamma }_0{e}^{\overrightarrow{\mathit{G}}·\overrightarrow{x}},$$

(19)

$$\epsilon =\epsilon (\overrightarrow{\mathit{x}})={\epsilon }_0{e}^{\overrightarrow{\mathit{G}}·\overrightarrow{x}},$$

(20)

where $\overrightarrow{\mathit{G}}=(G_1,G_2,G_3)$ is the chosen grading vector. Substituting the above definitions into the constitutive relations, and ultimately the equilibrium equations, yields a system of partial differential equations in terms of the 6-degrees of freedom. The general system can be written as a mixed boundary value problem, such that

$$\begin{array}{c}\underset{¯}{\mathit{L}}\overrightarrow{\mathit{U}}=\underset{¯}{0}\mathrm{in}\mathrm{\Omega },\hfill \\ \overrightarrow{\mathit{U}}=\overrightarrow{\tilde{\mathit{U}}}\mathrm{on}\partial {\mathrm{\Omega }}_1,\hfill \\ \overrightarrow{\mathit{t}}=\overrightarrow{\tilde{\mathit{t}}}\mathrm{on}\partial {\mathrm{\Omega }}_2,\hfill \end{array}$$

(21)

where $\overrightarrow{\tilde{\mathit{U}}}$ is the generalized displacement vector, prescribed on the Dirichlet part of the boundary $\partial {\mathrm{\Omega }}_1$ , and $\overrightarrow{\tilde{\mathit{t}}}$ is the generalized traction vector prescribed on the Neumann part of the boundary $\partial {\mathrm{\Omega }}_2$ , with $\partial \mathrm{\Omega }=\mathrm{\Gamma }=\partial {\mathrm{\Omega }}_1\cup \partial {\mathrm{\Omega }}_2$ .

The equivalent variational formulation of the equilibrium equations is presented, such that

$$\begin{array}{l}{\int }_{\mathrm{\Omega }}\left({\overrightarrow{\overline{\mathit{\gamma }}}}^T\overrightarrow{\mathit{\sigma }}+{\overrightarrow{\overline{\mathit{\varkappa }}}}^T\overrightarrow{\mu }\right)d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}({\overrightarrow{\overline{\mathit{u}}}}^T{\overrightarrow{\mathit{p}}}_v+{\overrightarrow{\overline{\mathit{\phi }}}}^T{\overrightarrow{\mathit{m}}}_v)d\mathrm{\Omega }\\ +{\int }_{\mathrm{\Gamma }}({\overrightarrow{\overline{\mathit{u}}}}^T{\overrightarrow{\mathit{p}}}_s+{\overrightarrow{\overline{\mathit{\phi }}}}^T{\overrightarrow{\mathit{m}}}_s)d\mathrm{\Gamma }=\overrightarrow{0},\end{array}$$

(22)

where ${\overrightarrow{\mathit{p}}}_s$ and ${\overrightarrow{\mathit{m}}}_s$ are prescribed stresses and couple stresses defined on the surface $\mathrm{\Gamma }$ , and ${\overrightarrow{\mathit{p}}}_v$ and ${\overrightarrow{\mathit{m}}}_{\upsilon }$ are body stresses and couple stresses defined in the interior domain $\mathrm{\Omega }$ . Equation (22) follows from the principle of stationary total potential energy, which states that the function that extremizes the total potential energy will be the one that satisfies equilibrium. That is to say, of the kinematically admissible solution states, the one that is required must satisfy the equilibrium equations and boundary conditions. As such, the required state must minimize the potential energy functional at stable equilibrium. Minimizing the functional consists of two parts. First, the first variation of the functional must vanish, thereby obtaining an extrememum of the functional (equation (22)). Next, the second variation must be shown to be non-negative, establishing the obtained extremum as a minimum [28]. The proof of existence and uniqueness of the solution to the mixed boundary value problem can be found in Barrage et al. [28]. The last two terms in equation (22) constitute the work done by body forces and body couples, and by external forces and couples, respectively. In what follows, the finite element method will be used to find the solution, with the proposed implementation focusing on the Dirichlet problem.

3. Finite element discretization and implementation into commercial software

3.1. Finite element formulation

To simplify the expressions for use in the finite element method, some grouping of vectors is introduced.

First, the degrees of freedom are grouped such that $\overrightarrow{\mathit{U}}=(u_x,u_y,u_z,{\phi }_x,{\phi }_y,{\phi }_z{)}^T$ . Next, the strains are grouped such that

$$\overrightarrow{\mathcal{E}}=\left[\begin{array}{c}\overrightarrow{\mathit{\gamma }}\hfill \\ \multicolumn{1}{c}{\overrightarrow{\mathit{\varkappa }}}\end{array}\right],$$

(23)

Following the grouping of the strains, the corresponding stresses are grouped such that

$$\overrightarrow{\mathit{\Sigma }}=\left[\begin{array}{c}\overrightarrow{\mathbf{\sigma }}\hfill \\ \multicolumn{1}{c}{\overrightarrow{\mathit{\mu }}}\end{array}\right].$$

(24)

Considering the constitutive laws for a linear elastic Cosserat material with central symmetry, defined in equations (13) and (14), the expressions

$$\begin{array}{c}{\sigma }_{\mathrm{ji}}=a_{\mathrm{jikl}}{\gamma }_{\mathrm{kl}}\hfill \\ \Rightarrow \overrightarrow{\mathbf{\sigma }}=\mathit{A}\overrightarrow{\mathit{\gamma }},\hfill \end{array}$$

(25)

$$\begin{array}{c}{\mu }_{\mathrm{ji}}=b_{\mathrm{jikl}}{\varkappa }_{\mathrm{kl}}\hfill \\ \Rightarrow \overrightarrow{\mathit{\mu }}=\mathit{B}\overrightarrow{\mathit{\varkappa }},\hfill \end{array}$$

(26)

can be grouped such that

$$\overrightarrow{\mathit{\Sigma }}=\left[\begin{array}{cc}\underset{¯}{\mathit{A}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \underset{¯}{\mathit{B}}\hfill \end{array}\right]\overrightarrow{\mathcal{E}}=\underset{¯}{\mathit{C}}\overrightarrow{\mathcal{E}}.$$

(27)

Finally, the body and surface stresses and couple stresses are grouped such that ${\overrightarrow{\mathit{T}}}_v=({\overrightarrow{\mathit{p}}}_v,{\overrightarrow{\mathit{m}}}_v{)}^T$ and ${\overrightarrow{\mathit{T}}}_{\mathbf{s}}=({\overrightarrow{\mathit{p}}}_s,{\overrightarrow{\mathit{m}}}_s{)}^T$ , respectively. With these groupings, the variational formulation from equation (14) reduces to

$${\int }_{\mathrm{\Omega }}{\overrightarrow{\overline{\mathcal{E}}}}^T\overrightarrow{\mathit{\Sigma }}d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}{\overrightarrow{\overline{\mathit{U}}}}^T{\overrightarrow{\mathit{T}}}_{v}d\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}{\overrightarrow{\overline{\mathit{U}}}}^T{\overrightarrow{\mathit{T}}}_{s}d\mathrm{\Gamma }=\overrightarrow{0}.$$

(28)

Due to the varying material properties within a functionally graded continuum, careful consideration is required in selecting the appropriate elements. Martínez-Pañedas [29] work on functionally graded finite elements suggests a need for quadratic elements to avoid introducing error related to non-linearly varying elastic properties. As such, quadratic, isoparameteric elements are considered. To properly discretize both simple and complex geometries, two elements are constructed: A quadratic, 20-noded isoparametric quadrilateral, and a quadratic, 15-noded isoparametric wedge. In what follows, the standard isoparametric approximations for the nodal degrees of freedom are used, such that

$$\overrightarrow{\mathit{U}}={\underset{¯}{\mathit{N}}}^T{\overrightarrow{\mathit{d}}}^e,$$

(29)

$$\overrightarrow{\overline{\mathit{U}}}={\underset{¯}{\mathit{N}}}^T{\overrightarrow{\overline{\mathit{d}}}}^e,$$

(30)

where $\mathit{N}$ has dimensions $6n\times 6$ , with $n$ corresponding to the number of nodes in the element, matrices of shape functions, and ${\overrightarrow{\mathit{d}}}^e={[u_{x_1},\cdots ,{\phi }_{z_1},\cdots ,u_{x_n},\cdots ,{\phi }_{z_n}]}^T$ is a $6n\times 1$ vector containing the element’s nodal degrees of freedom. The matrix $\underset{¯}{\mathit{N}}$ can be represented in block form, such that

$$\underset{¯}{\mathit{N}}={\left[\begin{array}{ccccccc}\underset{¯}{{\mathit{N}}_1}\hfill & \underset{¯}{0}\hfill & \underset{¯}{{\mathit{N}}_2}\hfill & \underset{¯}{0}\hfill & \cdots \hfill & \underset{¯}{{\mathit{N}}_{\mathit{n}}}\hfill & \underset{¯}{0}\hfill \\ \multicolumn{1}{c}{\underset{¯}{0}}& \underset{¯}{{\mathit{N}}_1}\hfill & \underset{¯}{0}\hfill & \underset{¯}{{\mathit{N}}_2}\hfill & \cdots \hfill & \underset{¯}{0}\hfill & \underset{¯}{{\mathit{N}}_{\mathit{n}}}\hfill \end{array}\right]}^T,$$

(31)

and

$$\underset{¯}{{\mathit{N}}_{\mathit{i}}}=\left[\begin{array}{ccc}{N}_i\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& N_i\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & N_i\hfill \end{array}\right],$$

(32)

represent matrices containing the element shape function $N_i$ at node $i$ . Next, the deformation-degree of freedom relations (equations (7) and (8)) are used to construct the deformation matrix $\underset{¯}{\mathit{D}}$ ( $18\times 6n$ ), mapping the vector ${\overrightarrow{\mathit{d}}}^e$ to the deformations $\overrightarrow{\mathcal{E}}$ . The deformation matrix $\underset{¯}{\mathit{D}}$ can be represented in block form, such that

$$\underset{¯}{\mathit{D}}=\left[\begin{array}{cccc}\underset{¯}{{\mathit{D}}_{\mathit{1}}}\hfill & \underset{¯}{{\mathit{D}}_{\mathit{2}}}\hfill & \cdots \hfill & \underset{¯}{{\mathit{D}}_{\mathit{n}}}\hfill \end{array}\right],$$

(33)

and

$$\underset{¯}{{\mathit{D}}_{\mathit{i}}}=\left[\begin{array}{cccccc}\frac{\partial N_i}{\partial x}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \frac{\partial N_i}{\partial y}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{\partial N_i}{\partial z}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \frac{\partial N_i}{\partial x}\hfill & 0\hfill & 0\hfill & 0\hfill & -N_i\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{\partial N_i}{\partial x}\hfill & 0\hfill & N_i\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{\partial N_i}{\partial y}\hfill & -N_i\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{\partial N_i}{\partial y}}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & N_i\hfill \\ \multicolumn{1}{c}{\frac{\partial N_i}{\partial z}}& 0\hfill & 0\hfill & 0\hfill & -N_i\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \frac{\partial N_i}{\partial z}\hfill & 0\hfill & N_i\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \frac{\partial N_i}{\partial x}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \frac{\partial N_i}{\partial y}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & \frac{\partial N_i}{\partial z}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \frac{\partial N_i}{\partial x}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & \frac{\partial N_i}{\partial x}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & \frac{\partial N_i}{\partial y}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \frac{\partial N_i}{\partial y}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \frac{\partial N_i}{\partial z}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \frac{\partial N_i}{\partial z}\hfill & 0\hfill \end{array}\right],$$

(34)

Furthermore, the chain rule used to relate differentiation in the global system to differentiation in the local system is given by

$$\left[\begin{array}{c}\frac{\partial N_i}{\partial x}\hfill \\ \multicolumn{1}{c}{\frac{\partial N_i}{\partial y}}\\ \multicolumn{1}{c}{\frac{\partial N_i}{\partial z}}\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial \xi }{\partial x}\hfill & \frac{\partial \eta }{\partial x}\hfill & \frac{\partial \zeta }{\partial x}\hfill \\ \multicolumn{1}{c}{\frac{\partial \xi }{\partial y}}& \frac{\partial \eta }{\partial y}\hfill & \frac{\partial \zeta }{\partial y}\hfill \\ \multicolumn{1}{c}{\frac{\partial \xi }{\partial z}}& \frac{\partial \eta }{\partial z}\hfill & \frac{\partial \zeta }{\partial z}\hfill \end{array}\right]\left[\begin{array}{c}\frac{\partial N_i}{\partial \xi }\hfill \\ \multicolumn{1}{c}{\frac{\partial N_i}{\partial \eta }}\\ \multicolumn{1}{c}{\frac{\partial N_i}{\partial \zeta }}\end{array}\right]={\underset{¯}{\mathit{J}}}^{-1}\left[\begin{array}{c}\frac{\partial N_i}{\partial \xi }\hfill \\ \multicolumn{1}{c}{\frac{\partial N_i}{\partial \eta }}\\ \multicolumn{1}{c}{\frac{\partial N_i}{\partial \zeta }}\end{array}\right],$$

(35)

where the Jacobian matrix of transformation $\underset{¯}{\mathit{J}}$ is defined as

$$\underset{¯}{\mathit{J}}=\left[\begin{array}{ccc}\frac{\partial x}{\partial \xi }\hfill & \frac{\partial y}{\partial \xi }\hfill & \frac{\partial z}{\partial \xi }\hfill \\ \multicolumn{1}{c}{\frac{\partial x}{\partial \eta }}& \frac{\partial y}{\partial \eta }\hfill & \frac{\partial z}{\partial \eta }\hfill \\ \multicolumn{1}{c}{\frac{\partial x}{\partial \zeta }}& \frac{\partial y}{\partial \zeta }\hfill & \frac{\partial z}{\partial \zeta }\hfill \end{array}\right],$$

(36)

Substituting the definitions into the variational formulation yields the familiar finite element formulation

$${\underset{¯}{\mathit{K}}}^e{\overrightarrow{\mathit{d}}}^e={\underset{¯}{\mathit{F}}}^{e,\mathrm{\Omega }}-{\underset{¯}{\mathit{F}}}^{e,\mathrm{\Gamma }}={\underset{¯}{\mathit{F}}}^e,$$

(37)

where

$${\underset{¯}{\mathit{K}}}^e=\left[{\int }_{\mathrm{\Omega }}{\underset{¯}{\mathit{D}}}^T\underset{¯}{\mathit{C}}\underset{¯}{\mathit{D}}d\mathrm{\Omega }\right],$$

(38)

$${\underset{¯}{\mathit{F}}}^{e,\mathrm{\Omega }}=\left[{\int }_{\mathrm{\Omega }}\underset{¯}{\mathit{N}}{\overrightarrow{\mathit{T}}}_{v}d\mathrm{\Omega }\right],$$

(39)

$${\underset{¯}{\mathit{F}}}^{e,\mathrm{\Gamma }}=\left[{\int }_{\mathrm{\Gamma }}\underset{¯}{\mathit{N}}{\overrightarrow{\mathit{T}}}_{s}d\mathrm{\Gamma }\right],$$

(40)

and where ${\underset{¯}{\mathit{K}}}^e$ is the element stiffness matrix, ${\underset{¯}{\mathit{F}}}^{e,\mathrm{\Omega }}$ represents the body force vector, and ${\underset{¯}{\mathit{F}}}^{e,\mathrm{\Gamma }}$ represents the external force vector, which vanishes in the considered Dirichlet problem.

3.2. Finite element implementation

User-defined material behaviour is typically implemented in Abaqus [24] in the form of a user-defined material subroutine (UMAT). The material behaviour is then incorporated into an element definition selected from the Abaqus element library. Such elements are not suitable for Cosserat analyses, as the 3 additional kinematically independent degrees of freedom (microrotations) are not accounted for in the available library of elements. To establish a dedicated Cosserat element as a foundation for future works, this work instead constructs a user element using the UEL subroutine.

The UEL subroutine defines all the relevant quantities for element computations, and can include a constitutive update to avoid having to define a separate user material. Visualization of the results is an area of concern with user defined elements, as Abaqus is incapable of natively displaying contour plots for quantities of interest, other than the three classical nodal displacements. Abaqus employs Voigt notation to define stresses, which renders its native output database structure unsuitable for the storage of results of elements with asymmetric stress states.

The issue of contour plot visualization is circumvented by Shi et al. [30] through the use of a ghost mesh approach, where the results from the user material simulation are overlaid onto the results form another simulation in which the material had negligible elastic properties. While this approach allows for native visualization of user elements that hold the same degrees of freedom as the available elements in the Abaqus library, it does not solve the issue for elements with additional degrees of freedom. Another approach by Eremeyev et al. [31] overcomes the issue of asymmetric stress storage by generating three output databases, storing 6 of 18 stress/couple-stress components in each, and keeping track of the variables manually. While effective, this approach can present some confusion to the user. This work instead bypasses the Abaqus output database altogether, writing the desired output to a series of comma separated value (CSV) files for post-processing.

A post-processing code is developed in Python, which reads the input file used for simulation, as well as the generated CSV output, and uses object oriented programming to replicate the output database hierarchical structure that Abaqus employs. The processed data is written to a series of VTU (virtual toolkit) files which can be viewed in open source visualization software such as Paraview. The approach used to write the VTU files is based on the approach by Liu et al. [32].

3.2.1. 20-noded quadrilateral

The 20-noded isoparametric quadratic quadrilateral element is designed with the nodal coordinates shown in Figure 1. The chosen node numbering follows the convention employed by Abaqus to facilitate integration into the software.

Figure 1.

Nodal locations in local coordinate system for quadrilateral element.

The associated shape functions are presented for reproducibility (following Abaqus’s convention) and take the form

$$N_i=\frac{1}{8}(1+{\xi }_i\xi )(1+{\eta }_i\eta )(1+{\zeta }_i\zeta )({\xi }_i\xi +{\eta }_i\eta +{\zeta }_i\zeta -2),i=1,\dots ,8,$$

(41)

$$N_i=\frac{1}{4}(1-{\xi }^2)(1+{\eta }_i\eta )(1+{\zeta }_i\zeta ),i=10,12,14,16,$$

(42)

$$N_i=\frac{1}{4}(1+{\xi }_i\xi )(1-{\eta }^2)(1+{\zeta }_i\zeta ),i=9,11,13,15,$$

(43)

$$N_i=\frac{1}{4}(1+{\xi }_i\xi )(1+{\eta }_i\eta )(1-{\zeta }^2),i=17,18,19,20,$$

(44)

in the natural coordinate system, where $({\xi }_i,{\eta }_i,{\zeta }_i)$ represents the local coordinates of node $i$ . An 8 × 8 × 8 Gauss quadrature scheme is used for numerical integration.

3.2.2. 15-noded wedge

To construct a 15-noded isoparametric quadratic wedge element, in a similar manner to the quadrilateral element, the node numbering convention used by Abaqus is followed, and the nodal coordinates can be seen in Figure 2. The associated shape functions are presented for reproducibility (following Abaqus’s convention) and take the form

$$N_1=\frac{1}{2}(1-\xi -\eta )(\zeta -1)(2\xi +2\eta +\zeta ),$$

(45)

$$N_2=\frac{1}{2}\xi (1-\zeta )(2\xi -\zeta -2),$$

(46)

$$N_3=\frac{1}{2}\eta (1-\zeta )(2\eta -\zeta -2),$$

(47)

$$N_4=\frac{1}{2}(1-\xi -\eta )(1+\zeta )(-2\xi -2\eta +\zeta ),$$

(48)

$$N_5=\frac{1}{2}\xi (1+\zeta )(2\xi +\zeta -2),$$

(49)

$$N_6=\frac{1}{2}\eta (1+\zeta )(2\eta +\zeta -2),$$

(50)

$$N_7=2\xi (1-\xi -\eta )(1-\zeta ),$$

(51)

$$N_8=2\xi \eta (1-\zeta ),$$

(52)

$$N_9=2\eta (1-\xi -\eta )(1-\zeta ),$$

(53)

$$N_{10}=2\xi (1-\xi -\eta )(1+\zeta ),$$

(54)

$$N_{11}=2\xi \eta (1+\zeta ),$$

(55)

$$N_{12}=2\eta (1-\xi -\eta )(1+\zeta ),$$

(56)

$$N_{13}=(1-\xi -\eta )(1-{\zeta }^2),$$

(57)

$$N_{14}=\xi (1-{\zeta }^2),$$

(58)

$$N_{15}=\eta (1-{\zeta }^2),$$

(59)

Figure 2.

Nodal locations in local coordinate system for wedge element.

To numerically integrate the relevant quantities, a nine-point formula for integration over triangular domains is used, and extended to three dimensions using the eight-point rule for the quadrilateral. The approach for quadrature over triangles is outlined in Cowper et al. [33].

4. Numerical examples—classical limit

Typically, verification of solutions to finite element analyses requires comparison with experimental tests. While measurement tools are becoming more refined as the field of nano fabrication expands, testing and measurement of microrotations are currently difficult due to the size at which measurements must take place. Until such a time where these measurements become more commonplace, verification of the designed user elements, and by extension the solution to the formulated boundary value problems, involves recovery of the classical limit. Such a verification only provides a foundational step to establishing the model as an analysis tool, as a full verification can only be performed once the aforementioned measurements can take place. Nevertheless, recovery of the classical limit is an important first step.

To verify the recovery of the classical limit, the Cosserat elastic constants ( $\alpha ,\beta ,\gamma ,\epsilon$ ) are set to zero, which recovers the classical equilibrium equations. The elements are tested and compared to Abaqus equivalents under uniaxial compression, shear, and torsion loading to ensure recovery of the classical limit, with the results for torsion being presented as a representative case. The designed quadrilateral and wedge elements can be compared to their classical counterparts in Abaqus, labelled C3D20, C3D15, a 20-noded quadrilateral element and 15-noded wedge element, respectively. The material analysed represents a typical polymer with a Young’s modulus of $1250\mathrm{MPa}$ and Poisson’s ratio of 0.25. These properties translate to values of $\lambda =\mu =500\mathrm{MPa}$ , respectively.

In all comparisons, the individual displacement components, as well as stress components, are compared to ensure the classical limit is recovered. To ensure an adequate comparison, the same boundary conditions are employed in both models, with the addition of setting the 3 microrotational degrees of freedom to 0 in the user element. Throughout what follows, the Dirichlet problem is solved, with prescribed displacements imposing the required loading scenario. All displayed units are measured in $\mathrm{N},\mathrm{mm},$ and $\mathrm{MPa}$ . The problem is stated such that

$$\begin{array}{c}\underset{¯}{\mathit{L}}\overrightarrow{\mathit{U}}=\underset{¯}{0}\mathrm{in}\mathrm{\Omega },\hfill \\ \overrightarrow{\mathit{U}}=\overrightarrow{\tilde{\mathit{U}}}\mathrm{on}\partial {\mathrm{\Omega }}_1,\hfill \end{array}$$

(60)

4.1. Quadrilateral element assembly

The solution to the Dirichlet problem for a quadrilateral element assembly comprising 27 quadrilateral elements (representing a $6\times 6\times 6$ $\mathrm{mm}$ cube) is solved in the case of torsion, with the torsion loading taking the form of a prescribed displacement couple of $0.1\mathrm{mm}$ acting along the positive and negative x-axis shown in Figure 3, with fixed conditions (all displacements and rotations set to zero) prescribed on the bottom surface, and denoted by the orange pins.

Figure 3.

Boundary conditions and loading for 27 element quadrilateral element assembly (user/C3D20) in torsion.

Figures 5–7 (in Appendix 1) compare the displacement components, normal stresses, and shear stresses of the user element assembly to an equivalent C3D20 element assembly at the final analysis increment. The numerical values corresponding to the compared solutions displayed in Figures 5–7 (in Appendix 1) are identical up to Single Precision order. As such, all displacements are accurate to nine significant digits.

4.2. Wedge element assembly (cylinder)

The solution to the Dirichlet problem for a cylindrical wedge element assembly comprising of 36 wedge elements (representing a $2$ $\mathrm{mm}$ diameter cylinder with a length of $6$ $\mathrm{mm}$ ) is solved in the case of torsion, with the torsion loading taking the form of two prescribed displacement couples of $0.1\mathrm{mm}$ acting along the positive and negative x and y axes shown in Figure 4, with fixed conditions (all displacements and rotations set to 0) prescribed on the bottom surface, and denoted by the orange pins. Figures 8–10 (in Appendix 1) compare the displacement components, normal stresses, and shear stresses of the user element assembly to an equivalent C3D15 element assembly at the final analysis increment.

Figure 4.

Boundary conditions and loading for 36 element cylindrical wedge element assembly (user/C3D15) in torsion.

The numerical values corresponding to the compared solutions displayed in Figures 8–10 (in Appendix 1) are identical up to single precision order (accurate to nine significant digits), completing the recovery of the classical limit in the case of torsion for assemblies of the designed elements. While not presented in this work, the cases of uniaxial compression loading and shear loading also result in recovery of the classical limit, maintaining single precision accuracy.

5. Numerical example—conceptual case

To demonstrate a conceptual application of the model, a case in which fictitious material parameters are used is considered. The Dirichlet boundary value problem in equation (60), in the case of torsion loading (Figure 4) for the cylindrical wedge assembly is once again considered. The material elastic constants considered are ${\lambda }_0={\mu }_0=500\mathrm{MPa}$ in all cases, ${\alpha }_0=500\mathrm{MPa},$ ${\beta }_0={\gamma }_0={\epsilon }_0=500\mathrm{N}$ in all cases with Cosserat effects, and a chosen grading vector of $\overrightarrow{\mathit{G}}=(1,1,1)$ in all cases where exponential grading is considered. These values are chosen arbitrarily as no robust measurements of the parameters have been made experimentally, and only serve to display a proof of concept.

Figure 11 (in Appendix 1) displays a comparison of the magnitude of the stress vector for all four considered cases. A notable localization of stresses occurs in the case with only Cosserat effects at the loading nodes, with a more significant localization occurring at the loading nodes in the direction of grading in the case with both exponential grading and Cosserat effects.

The localized rotations at the load application nodes induced by the torsional loading give rise to significant differences between the classical predictions and Cosserat predictions. The most significant differences between classical and Cosserat predictions are therefore likely to arise in load cases where the loading induces more significant nodal rotations; this corresponds to the cases of torsional and shear loading, as opposed to compression loading.

A comparison of the numerical values for the magnitude of the stress vector highlights the differences between the material models, with the introduction of Cosserat effects only presenting a 270% increase in maximum magnitude (redistribution and localization). Similarly, exponential grading only offers a 45% increase in maximum stress magnitude from the classical case, as the stresses are localized in the direction of the grading. Introducing both effects results in a similar increase from the classical case, with the peak magnitudes being 24, 90, 35, and 110 $\mathrm{MPa}$ for the classical, Cosserat only, exponential grading only, and Cosserat with exponential grading cases, respectively.

6. Discussion

A three-dimensional, linear elastic, isotropic model for a functionally graded composite material is developed, and presented as a mixed boundary value problem. The model considers microstructural effects through Cosserat elasticity. The associated variational (weak form) formulation is presented and used to discretize the problem through the finite element method. Two quadratic elements (20-noded quadrilateral, 15-noded wedge) are developed and implemented in Abaqus through the UEL subroutine.

Preliminary validation of the elements through recovery of the classical limit in the Dirichlet problem is conducted in the cases of uniaxial compression, shear, and torsion loading, with the torsion loading case presented as a representative case. All displacements are accurate up to single precision order (nine significant digits).

Future availability of experimental measurements of microrotations allows for complete validation of the model, and its use as an analysis tool. Exponential grading is presented in this work; however, any chosen grading scheme can be used with the presented method. For alternative grading schemes, careful consideration must be put into the selected quadrature scheme to ensure an accurate response.

A demonstration of the applicability of the model is presented through the use of fictitious elastic parameters, demonstrating notable differences in the responses when compared to the classical case. The differences arise in the form of stress localization, in particular at the load application nodes, and in the direction of functional grading.