All papers are copyrighted by either a publisher, the research group of G. H. Paulino, or the authors. Publications available here are for educational or academic use only. All rights of reproduction or distribution in any form are reserved.
 2009  2008  2007  2006  2005  2004  2003  2002  2001  2000 and prior
Publications In Press and Online  
H. Chi, Y. Zhang, T. L. E. Tang, L. Mirabella, L. Dalloro, L. Song, and G. H. Paulino. "Universal Machine Learning for Topology Optimization." Computer Methods in Applied Mechanics and Engineering. Accepted. View/Hide Abstract
In this work, we put forward a general machine learningbased topology optimization framework which greatly accelerates the design process of largescale problems without any sacrifice in accuracy. The proposed framework has three novel features. First, a novel online training concept is proposed using the history data of topology optimization. Thus, the training is done during rather than before the topology optimization. Second, the framework adopts a tailored twoscale topology optimization formulation and introduces a localized online training strategy. The localized online training strategy can improve both the scalability and accuracy of the proposed framework. Third, the framework incorporates an online updating scheme which continuously improves the prediction accuracy of the machine learning models by providing them with new data generated from physical simulations Through numerical investigations and design examples, we demonstrate that the proposed framework is highly scalable and can efficiently handle design problem with a wide range of discretization levels, different load and boundary conditions, and various design considerations (e.g., the presence of nondesignable regions).


W. Li, P. Suryanarayana, and G. H. Paulino. "Accelerated fixedpoint formulation of topology optimization: Application to compliance minimization problems." Mechanics Research Communications. In Press. https://doi.org/10.1016/j.mechrescom.2019.103469. View/Hide Abstract
We present a simple, effective, and scalable approach for significantly accelerating the convergence in Topology Optimization simulations. Specifically, treating the design process as a fixedpoint iteration, we propose employing a recently developed acceleration technique in which Anderson extrapolation is applied periodically, with simple weighted relaxation used for the remaining steps. Through selected examples in compliance minimization, we show that the proposed approach is able to accelerate the overall simulation several fold, while maintaining the quality of the solution.


H. Chi, A. Pereira, I. F. M. Menezes, and G. H. Paulino. "Virtual element method (VEM)based topology optimization: an integrated framework." Structural and Multidisciplinary Optimization. Available online. DOI 10.1007/s0015801902268w. Link to PDF View/Hide Abstract
We present a virtual element method (VEM)based topology optimization framework using polyhedral elements, which allows for easy handling of nonCartesian design domains in three dimensions. We take full advantage of the VEM properties by creating a unified approach in which the VEM is employed in both the structural and the optimization phases of the framework. In the structural problem, the VEM is adopted to solve the threedimensional elasticity equation. Compared to the finite element method (FEM), the VEM does not require numerical integration and is less sensitive to degenerated elements (e.g., ones with skinny faces or small edges). In the optimization problem, we introduce a continuous approximation of material densities using VEM basis functions. As compared to the standard elementwise constant one, the continuous approximation enriches geometrical representations of structural topologies. Through two numerical examples with exact solutions, we verify the convergence and accuracy of both the VEM approximations of the displacement and material density fields. We also present several design examples involving nonCartesian domains, demonstrating the main features of the proposed VEMbased topology optimization framework. The source code for a MATLAB implementation of the proposed work, named PolyTop3D, will be made available in the (electronic) Supplementary Material accompanying this publication.


J. Chun, G. H. Paulino, and J. Song. "Reliabilitybased topology optimization by ground structure method employing a discrete filtering technique." Structural and Multidisciplinary Optimization. Available online. DOI 10.1007/s00158019022551 Link to PDF View/Hide Abstract
When conventional filtering schemes are used in reliabilitybased topology optimization (RBTO), identified solutions may violate probabilistic constraints and/or global equilibrium. In order to address this issue, this paper proposes to incorporate a discrete filtering technique termed the discrete filtering method (Ramos Jr. and Paulino 2016) into RBTO using the elastic formulation of the ground structure method. The discrete filtering method allows the optimizer to achieve more physically realizable truss designs in which thin bars are eliminated while ensuring global equilibrium. The method uses a potentialenergybased approach with Tikhonov regularization to solve the singular system of equations that may result from imposing the discrete filter. Combining this method with RBTO allows us to use the reliabilitybased truss sizing optimization for the purpose of topology optimization under uncertainties. Furthermore, a singleloop approach is adopted to enhance the computational efficiency of the proposed RBTO method. Numerical examples of two and threedimensional engineering designs demonstrate useful features of the proposed method and illustrate the influence of the discrete filter and parameter uncertainties on the optimization results. In order to check if the optimal topologies obtained by the proposed approach satisfy the constraints on the failure probabilities, structural reliability analysis is also performed using the firstorder reliability method and Monte Carlo simulations.


T. Zegard, S. R. M. Almeida, I. A. Torres, and G. H. Paulino. "Topology optimization of tall buildings under multiple loads." Submitted. View/Hide Abstract
Structural systems for supertall buildings have lateral resistance at their core; for these structures are governed by lateral loads, typically wind. The analysis of these structures has to be threedimensional in order to correctly model and capture the interactions that occur between the building's planes. The direction of these lateral loads is often unknown, and therefore an analysis (and optimization) considering multiple load directions is required. This, together with the dead and live loads, comprise the load combinations which must be considered to obtain meaningful results. For the optimized designs to be practical, manufacturing constraints such as symmetry and repetition are enforced. This manuscript explores the usage of densitybased topology optimization, considering multiple load cases, symmetry, and pattern repetition, to obtain optimal layouts of material, which can be interpreted as an optimized structural system for supertall buildings.


Back to Top 
2020  
O. GiraldoLondoño, Lucia Mirabella, Livio Dalloro, and G. H. Paulino. "Multimaterial thermomechanical topology optimization with applications to additive manufacturing: Design of main composite part and its support structure." Computer Methods in Applied Mechanics and Engineering. Vol. 363, No. 112812, pp. 127, 2020. Link to PDF View/Hide Abstract
This paper presents a densitybased topology optimization formulation for the design of multimaterial thermoelastic structures. The formulation is written in the form of a multiobjective topology optimization problem that considers two competing objective functions, one related to mechanical performance (mean compliance) and one related to thermal performance (either thermal compliance or temperature variance). To solve the optimization problem, we present an efficient design variable update scheme, which we have derived in the context of the Zhang–Paulino–Ramos (ZPR) update scheme by Zhang et al. (2018). The new update scheme has the ability to solve nonselfadjoint topology optimization problems with an arbitrary number of volume constraints, which can be imposed either to a subset of the candidate materials, or to subregions of the design domain, or to a combination of both. We present several examples that explore the ability of the formulation to obtain candidate Pareto fronts and to design support structures for additive manufacturing. Enabled by the ZPR update scheme, we are able to control the complexity and the length scale of the support structures by means of regional volume constraints.


K. Park, H. Chi, and G. H. Paulino. "Numerical recipes for elastodynamic virtual element methods with explicit time integration." International Journal for Numerical Methods in Engineering. Vol. 121, No. 1, pp. 131, 2020. Link to PDF View/Hide Abstract
We present a general framework to solve elastodynamic problems by means of the virtual element method (VEM) with explicit time integration. In particular, the VEM is extended to analyze nearly incompressible solids using the B‐bar method. We show that, to establish a B‐bar formulation in the VEM setting, one simply needs to modify the stability term to stabilize only the deviatoric part of the stiffness matrix, which requires no additional computational effort. Convergence of the numerical solution is addressed in relation to stability, mass lumping scheme, element size, and distortion of arbitrary elements, either convex or nonconvex. For the estimation of the critical time step, two approaches are presented, ie, the maximum eigenvalue of a system of mass and stiffness matrices and an effective element length. Computational results demonstrate that small edges on convex polygonal elements do not significantly affect the critical time step, whereas convergence of the VEM solution is observed regardless of the stability term and the element shape in both two and three dimensions. This extensive investigation provides numerical recipes for elastodynamic VEMs with explicit time integration and related problems.


Back to Top 
2019  
K. Liu, T. Tachi, and G. H. Paulino. "Invariant and smooth limit of discrete geometry folded from bistable origami leading to multistable metasurfaces." Nature Communications. 10 No. 4238, 2019. Link to PDF  Link to SI Movie 1  Movie 2  Movie 3 View/Hide Abstract
Origami offers an avenue to program threedimensional shapes via scaleindependent and nondestructive fabrication. While such programming has focused on the geometry of a tessellation in a single transient state, here we provide a complete description of folding smooth saddle shapes from concentrically pleated squares. When the offset between square creases of the pattern is uniform, it is known as the pleated hyperbolic paraboloid (hypar) origami. Despite its popularity, much remains unknown about the mechanism that produces such aesthetic shapes. We show that the mathematical limit of the elegant shape folded from concentrically pleated squares, with either uniform or nonuniform (e.g. functionally graded, random) offsets, is invariantly a hyperbolic paraboloid. Using our theoretical model, which connects geometry to mechanics, we prove that a folded hypar origami exhibits bistability between two symmetric configurations. Further, we tessellate the hypar origami and harness its bistability to encode multistable metasurfaces with programmable nonEuclidean geometries.


K. Park, H. Chi, and G. H. Paulino. "On nonconvex meshes for elastodynamics using virtual element methods with explicit time integration." Computer Methods in Applied Mechanics and Engineering. Vol. 356, pp. 669684, 2019. Link to PDF View/Hide Abstract
While the literature on numerical methods (e.g. finite elements and, to a certain extent, virtual elements) concentrates on convex elements, there is a need to probe beyond this limiting constraint so that the field can be further explored and developed. Thus, in this paper, we employ the virtual element method for nonconvex discretizations of elastodynamic problems in 2D and 3D using an explicit time integration scheme. In the formulation, a diagonal matrixbased stabilization scheme is proposed to improve performance and accuracy. To address efficiency, a critical time step is approximated and verified using the maximum eigenvalue of the local (rather than global) system. The computational results demonstrate that the virtual element method is able to consistently handle general nonconvex elements and even nonsimply connected elements, which can lead to convenient modeling of arbitrarilyshaped inclusions in composites.


K. Liu, T. Zegard, P. P. Pratapa, and G. H. Paulino. "Unraveling tensegrity tessellations for metamaterials with tunable stiffness and bandgaps." Journal of the Mechanics and Physics of Solids. Vol. 131, pp. 147166, 2019. Link to PDF  Supplementary Material  Movie View/Hide Abstract
Tensegrity structures resemble biological tissues: A structural system that holds an internal balance of prestress. Owing to the presence of prestress, biological tissues can dramatically change their properties, making tensegrity a promising platform for tunable and functional metamaterials. However, tensegrity metamaterials require harmony between form and force in an infinitely–periodic scale, which makes the design of such systems challenging. In order to explore the full potential of tensegrity metamaterials, a systematic design approach is required. In this work, we propose an automated design framework that provides access to unlimited tensegrity metamaterial designs. The framework generates tensegrity metamaterials by tessellating blocks with designated geometries that are aware of the system periodicity. In particular, our formulation allows creation of Class1 (i.e., floating struts) tensegrity metamaterials. We show that tensegrity metamaterials offer tunable effective elastic moduli, Poisson’s ratio, and phononic bandgaps by properly changing their prestress levels, which provide a new dimension of programmability beyond geometry.


O. GiraldoLondoño, G. H. Paulino, and W. G. Buttlar. "Fractional calculus derivation of a ratedependent PPRbased cohesive fracture model: Theory, implementation, and numerical results." International Journal of Fracture. Vol. 216, pp. 129, 2019. Link to PDF View/Hide Abstract
Ratedependent fracture has been extensively studied using cohesive zone models (CZMs). Some of them use classical viscoelastic material models based on springs and dashpots. However, such viscoelastic models, characterized by relaxation functions with exponential decay, are inadequate to simulate fracture for a wide range of loading rates. To improve the accuracy of existing models, this work presents a mixedmode ratedependent CZM that combines the features of the Park–Paulino–Roesler (PPR) cohesive model and a fractional viscoelastic model. This type of viscoelastic model uses differential operators of noninteger order, leading to powerlawtype relaxation functions with algebraic decay. We derive the model in the context of damage mechanics, such that undamaged viscoelastic tractions obtained from a fractional viscoelastic model are scaled using two damage parameters. We obtain these parameters from the PPR cohesive model and enforce them to increase monotonically during the entire loading history, which avoids artificial selfhealing. We present three examples, two used for validation purposes and one to elucidate the physical meaning of the fractional differential operators. We show that the model is able to predict ratedependent fracture process of rubberlike materials for a wide range of loading rates and that it can capture ratedependent mixedmode fracture processes accurately. Results from the last example indicate that the order of the fractional differential operators acts as a memorylike parameter that allows for the fracture modeling of long and shortterm memory processes. The ability of fractional viscoelastic models to model this type of process suggests that relaxation functions with algebraic decay lead to accurate fracture modeling of materials for a wide range of loading rates.


T. Zhao, A. S. Ramos Jr., G. H. Paulino. "Material Nonlinear Topology Optimization Considering the von Mises Criterion through an Asymptotic Approach: Max Strain Energy and Max Load Factor Formulations." International Journal of Numerical Methods in Engineering. Vol. 118, pp.804828, 2019. Link to PDF View/Hide Abstract
This paper addresses material nonlinear topology optimization considering the vonMises criterion bymeans of an asymptotic analysis using a fictitious nonlinear elastic model. In this context,we consider the topology optimization problem subjected to prescribed energy, which leads to robust convergence in nonlinear problems. Two nested formulations are considered. In the first, the objective is to maximize the strain energy of the system in equilibrium, and in the second, the objective is to maximize the load factor. In both cases, a volume constraint is imposed. The sensitivity analysis is quite effective and efficient in the sense that there is no extra adjoint equation. In addition, the nonlinear structural equilibrium problem is solved using direct minimization of the structural strain energy using Newton's method with an inexact linesearch strategy. Four numerical examples demonstrate the features of the proposed material nonlinear topology optimization framework for approximating standard von Mises plasticity.


K. Liu and G. H. Paulino. "Tensegrity topology optimization by force maximization on arbitrary ground structures. " Structural and Multidisciplinary Optimization. Vol. 59, pp. 20412062, 2019. Link to PDF View/Hide Abstract
This paper presents an optimization approach for design of tensegrity structures based on graph theory. The formulation obtains tensegrities from ground structures, through force maximization using mixed integer linear programming. The method seeks a topology of the tensegrity that is within a given geometry, which provides insight into the tensegrity design from a geometric point of view. Although not explicitly enforced, the tensegrities obtained using this approach tend to be both stable and symmetric. Borrowing ideas from computer graphics, we allow “restriction zones” (i.e., passive regions in which no geometric entity should intersect) to be specified in the underlying ground structure. Such feature allows the design of tensegrities for actual engineering applications, such as robotics, in which the volume of the payload needs to be protected. To demonstrate the effectiveness of our proposed design method, we show that it is effective at extracting both wellknown tensegrities and new tensegrities from the ground structure network, some of which are prototyped with the aid of additive manufacturing.


P. P. Pratapa, K. Liu, G. H. Paulino. "Geometric mechanics of origami patterns exhibiting Poisson’s ratio switch by breaking mountain and valley assignment." Physical Review Letters. Vol. 122, pp. 1555011 to 6, 2019. Link to PDF Link to SI Movie 1  Movie 2  Movie 3  Movie 4 View/Hide Abstract
Exploring the configurational space of specific origami patterns [e.g., Miuraori (flat surface with parallelogram crease patterns), eggbox] has led to notable advances in science and technology. To augment the origami design space, we present a pattern, named "Morph," which combines the features of its parent patterns. We introduce a fourvertex origami cell that morphs continuously between a Miura mode and an eggbox mode, forming an homotopy class of configurations. This is achieved by changing the mountain and valley assignment of one of the creases, leading to a smooth switch through a wide range of negative and positive Poisson’s ratios. We present elegant analytical expressions of Poisson’s ratios for both inplane stretching and outofplane bending and find that they are equal in magnitude and opposite in sign. Further, we show that by combining compatible unit cells in each of the aforementioned modes through kinematic bifurcation, we can create hybrid origami patterns that display unique properties, such as topological mode locking and tunable switching of Poisson’s ratio.


H. Chi, D. L. Ramos, A. S. Ramos Jr., and G. H. Paulino. "On structural topology optimization considering material nonlinearity: Plane strain versus plane stress solutions." Advances in Engineering Software. Vol. 131, pp. 217231, 2019. Link to PDF View/Hide Abstract
We present a structural topology optimization framework considering material nonlinearity by means of a tailored hyperelastic formulation. The nonlinearity is incorporated through a hyperelastic constitutive model, which is capable of capturing a range of nonlinear material behavior under both plane strain and plane stress conditions. We explore both standard (i.e. quadrilateral) and polygonal finite elements in the solution process, and achieve smooth convergence in both the optimization process and the solution of nonlinear state equations. Numerical examples are presented, which demonstrate that the topology optimization framework can effectively capture the influence of various material behaviors, load levels and loading conditions (i.e. plane stress versus plane strain) on the optimal topologies.


J. Vila, G. H. Paulino, and M. Ruzzene. "Role of nonlinearities in topological protection: Using magnetically coupled fidget spinners." Physical Review B. Vol. 99, pp. 1251161 to 13, 2019. Link to PDF Movie 1  Movie 2  Movie 3  Movie 4 View/Hide Abstract
We investigate and experimentally observe the existence of topologically protected interface modes in a onedimensional mechanical lattice and we report on the effect of nonlinearities on topological protection. The lattice consists of a onedimensional array of spinners with nearestneighbor coupling resulting from magnetic interactions. The distance between the spinners is spatially modulated to obtain a diatomic configuration, and to produce a nontrivial interface by breaking spatial inversion symmetry. For small amplitudes of motion, the interactions are approximately linear, and the system supports topologically protected interface modes at frequencies inside the bulk band gap of the lattice. Nonlinearities induced by increasing amplitude of motion cause the interface modes to shift and merge with the bulk bands. The resulting edgetobulk transition causes the extinction of the topologically protected interface mode and extends it to the entire length of the chain. Such transition is predicted by analytical calculations and verified by experimental observations. The paper thus investigates topologically protected interface modes obtained through broken spatial inversion symmetry, and documents the lack of robustness in the presence of nonlinearities.


H. Chi, L. Beirao da Veiga, and G. H. Paulino. "A simple and effective gradient recovery scheme and a posteriori error estimator for the Virtual Element Method (VEM)." Computer Methods in Applied Mechanics and Engineering. Vol 347, pp. 2158, 2019. Link to PDF View/Hide Abstract
This paper introduces a general recoverybased a posteriori error estimation framework for the Virtual Element Method (VEM) of arbitrary order on general polygonal/polyhedral meshes. The framework consists of a gradient recovery scheme and a posteriori error estimator based on the recovered displacement gradient. A skeletal error, which accurately mimics the behavior of the error of the displacement gradient by only sampling the displacement gradient on the mesh skeleton, is introduced. Through numerical studies on various polygonal/polyhedral meshes, we demonstrate that the proposed gradient recovery scheme can produce considerably more accurate displacement gradient than the original VEM solutions, and that the a posteriori error estimator is able to accurately capture both local and global errors without the knowledge of exact solutions. We also demonstrate the use of the VEM skeletal error estimators to guide adaptive mesh refinement.


J. Chun, J. Song, and G. H. Paulino. "Systemreliabilitybased design and topology optimization of structures under constraints on firstpassage probability." Structural Safety. Vol. 76, pp. 8194, 2019. Link to PDF View/Hide Abstract
For the purpose of reliability assessment of a structure subject to stochastic excitations, the probability of the
occurrence of at least one failure event over a time interval, i.e. the firstpassage probability, often needs to be evaluated. In this paper, a new method is proposed to incorporate constraints on the firstpassage probability
into reliabilitybased optimization of structural design or topology. For efficient evaluations of firstpassage
probability during the optimization, the failure event is described as a series system event consisting of instantaneous failure events defined at discrete time points. The probability of the series system event is then
computed by use of a system reliability analysis method termed as the sequential compounding method. The
adjoint sensitivity formulation is derived for calculating the parameter sensitivity of the firstpassage probability
to facilitate the use of efficient gradientbased optimization algorithms. The proposed method is successfully demonstrated by numerical examples of a space truss and building structures subjected to stochastic earthquake ground motions.


E. T. Filipov, G. H. Paulino, and T. Tachi. "Deployable sandwich surfaces with high outofplane stiffness." ASCE Journal of Structural Engineering. Vol. 145, pp. 040182441 to 15, 2019. Link to PDF View/Hide Abstract
This paper presents a set of deployable origami tube structures that can create smooth functional surfaces while simultaneously maintaining a high outofplane stiffness both during and after deployment. First, a generalized geometric definition for these tubes is presented such that they can globally have straight, curved, or segmented profiles, while the tubes can locally have skewed and reconfigurable cross sections. Multiple tubes can be stacked to form continuous and smooth assemblies in order to enable applications, including driving surfaces, roofs, walls, and structural hulls. Threepoint bending analyses and physical prototypes were used to explore how the orthogonal stiffness of the tubular structures depends on the geometric design parameters. The coupled tube structures typically had their highest outofplane stiffness when near to a fully deployed state. Tubes with skewed cross sections and more longitudinal variation (i.e., that had more zigzags) typically had a higher stiffness during deployment than tubes that were generally straight.


Back to Top 
2018  
S. A. Nauroze, L. S. Novelino, M. M. Tentzeris, and G. H. Paulino. "Continuousrange tunable multilayer frequencyselective surfaces using origami and inkjet printing
." Proceedings of the National Academy of Sciences. Vol. 115, No. 52, pp. 1321013215, 2018. Link to PDF  Link to SI  Movie 1  Movie 2  Movie 3 View/Hide Abstract
The tremendous increase in the number of components in typical electrical and communication modules requires lowcost, flexible and multifunctional sensing, energy harvesting, and communication modules that can readily reconfigure, depending on changes in their environment. Current subtractive manufacturingbased reconfigurable systems offer limited flexibility (limited finite number of discrete reconfiguration states) and have high fabrication cost and time requirements. Thus, this paper introduces an approach to solve the problem by combining additive manufacturing and origami principles to realize tunable electrical components that can be reconfigured over continuousstate ranges from folded (compact) to unfolded (large surface) configurations. Special “bridgelike” structures are introduced along the traces that increase their flexibility, thereby avoiding breakage during folding. These techniques allow creating truly flexible conductive traces that can maintain high conductivity even for large bending angles, further enhancing the states of reconfigurability. To demonstrate the idea, a MiuraOri pattern is used to fabricate spatial filters—frequencyselective surfaces (FSSs) with dipole resonant elements placed along the fold lines. The electrical length of the dipole elements in these structures changes when the MiuraOri is folded, which facilitates tunable frequency response for the proposed shapereconfigurable FSS structure. Higherorder spatial filters are realized by creating multilayer MiuraFSS configurations, which further increase the overall bandwidth of the structure. Such multilayer MiuraFSS structures feature the unprecedented capability of onthefly reconfigurability to different specifications (multiple bands, broadband/narrowband bandwidth, wide angle of incidence rejection), requiring neither specialized substrates nor highly complex electronics, holding frames, or fabrication processes.


E. D. Sanders, A. Pereira, M. A. Aguiló, and G. H. Paulino. "PolyMat: an efficient Matlab code for multimaterial topology optimization." Structural and Multidisciplinary Optimization. Vol. 58, pp. 27272759, 2018. Link to PDF  Download PolyMat View/Hide Abstract
We present a Matlab implementation of topology optimization for compliance minimization on unstructured polygonal finite element meshes that efficiently accommodates many materials and many volume constraints. Leveraging the modular structure of the educational code, PolyTop, we extend it to the multimaterial version, PolyMat, with only a few modifications. First, a design variable for each candidate material is defined in each finite element. Next, we couple a Discrete Material Optimization interpolation with the existing penalization and introduce a new parameter such that we can employ continuation and smoothly transition from a convex problem without any penalization to a nonconvex problem in which material mixing and intermediate densities are penalized. Mixing that remains due to the density filter operation is eliminated via continuation on the filter radius. To accommodate flexibility in the volume constraint definition, the constraint function is modified to compute multiple volume constraints and the design variable update is modified in accordance with the ZhangPaulinoRamos Jr. (ZPR) update scheme, which updates the design variables associated with each constraint independently. The formulation allows for volume constraints controlling any subset of the design variables, i.e., they can be defined globally or locally for any subset of the candidate materials. Borrowing ideas for mesh generation on complex domains from PolyMesher, we determine which design variables are associated with each local constraint of arbitrary geometry. A number of examples are presented to demonstrate the many material capability, the flexibility of the volume constraint definition, the ease with which we can accommodate passive regions, and how we may use local constraints to break symmetries or achieve
graded geometries.


L. A. M. Camones, E. A. Vargas Jr., R. Q. Velloso, and G. H. Paulino. "Simulation of hydraulic fracturing processes in rocks by coupling the lattice
Boltzmann model and the ParkPaulinoRoesler potentialbased cohesive
zone model. " International Journal of Rock Mechanics and Mining Sciences. Vol. 112, pp. 339353, 2018. Link to PDF View/Hide Abstract
This paper contributes a scheme for twodimensional (2D) numerical simulation of hydraulic fracturing processes
in geological materials, which consists of a numerical coupling between the finite element method (FEM)
and the lattice Boltzmann (LB) model. The ParkPaulinoRoesler potentialbased cohesive zone model (PPR) is
used to simulate the fracture propagation by means of interface finite elements, such that cohesive forces act on
the fracture surface, capturing the softening process. The PPR model is used because it is a generalized fracture
model that can represent the fracture process for mode I, mode II and mixed mode III, and can be applied to
various materials, including heterogeneous materials, such as rock. The FEM and LB are coupled in an iterative
process. The paper describes implementation details including procedures for coupling both methods. Examples
of hydraulic fracturing process, modeled with the proposed FEMLB coupling demonstrate the potential of this
numerical procedure to model hydraulic fracturing processes in geomaterials of complex geometries.


Z. Zhao, X. Kuang, J. Wu, Q. Zhang, G. H. Paulino, H. J. Qi, and D. Fang. "3D printing of complex origami assemblages for reconfigurable structures." Soft Matter. Vol. 14, pp. 80518059, 2018. Link to PDF Link to ESI, Movie 1 (MP4), Movie 2 (MP4), Movie 3 (MP4) View/Hide Abstract


J. C. R. Albino, C. A. Almeida, I. F. M. Menezes, and G. H. Paulino. "Corotational 3D beam element for nonlinear dynamic analysis of risers manufactured with functionally graded materials (FGMs). " Engineering Structures. Vol. 173, pp. 283299, 2018. Link to PDF View/Hide Abstract
This paper considers a corotational beam formulation for beams, which is used for the finite element analysis of flexible risers and pipelines made of functionally graded materials. The influence of material gradation is addressed using an exponential variation of properties throughout the thickness of the pipe. Space discretization of the equilibrium equations is derived based on the Euler–Bernoulli assumptions considering twonode Hermitian beam elements which are referred to a corotation coordinate system attached to the element local frame of coordinates. The geometric nonlinear effects of the beam are considered under large displacement and rotations, but under smallstrain conditions. The deflections of the riser result from forces caused by selfweight, buoyancy, sea currents, waves, the action of floaters, seabedstructure interactions, and ship’s motion. We provide numerical examples and compare our results with the ones available in the literature. In addition, applications related to practical offshore engineering situations are considered to highlight the behavior of functionally graded materials (FGMs) as compared to homogeneous risers.


E. D. Sanders, M. A. Aguiló, and G. H. Paulino. "Multimaterial continuum topology optimization with arbitrary volume and mass constraints." Computer Methods in Applied Mechanics and Engineering. Vol. 340, pp. 798823, 2018. Link to PDF View/Hide Abstract
A framework is presented for multimaterial compliance minimization in the context of continuum based topology optimization. We adopt the common approach of ﬁnding an optimal shape by solving a series of explicit convex (linear) approximations to the volume constrained compliance minimization problem. The dual objective associated with the linearized subproblems is a separable function of the Lagrange multipliers and thus, the update of each design variable is dependent only on the Lagrange multiplier of its associated volume constraint. By tailoring the ZPR design variable update scheme to the continuum setting, each volume constraint is updated independently. This formulation leads to a setting in which sufﬁciently general volume/mass constraints can be speciﬁed, i.e., each volume/mass constraint can control either all or a subset of the candidate materials and can control either the entire domain (global constraints) or a subregion of the domain (local constraints). Material interpolation schemes are investigated and coupled with the presented approach. The key ideas presented herein are demonstrated through representative examples in 2D and 3D.


P. P. Pratapa, P. Suryanarayana, and G. H. Paulino. "Bloch wave framework for structures with nonlocal interactions: Application to the design of origami acoustic metamaterials. " Journal of the Mechanics and Physics of Solids. Vol. 118, pp. 115132, 2018. Link to PDF View/Hide Abstract
We present a generalized Bloch wave framework for the dynamic analysis of structures with nonlocal interactions and apply it to the design of origami acoustic metamaterials. Specifically, we first discretize the origami structures using a customized structural barandhinge model that minimizes the degrees of freedom in the associated unit cell, while being sufficiently accurate to capture the behavior of interest. Next, observing that this discretization results in nonlocal structural interactions—the stiffness matrix has nonzeros between nodes that are not nearest neighbors due to the coupled deformations arising during folding or bending—we generalize the standard Bloch wave approach used in structural analysis to enable the study of such systems. Utilizing this framework, choosing the geometry of the unit cell as well as the folded state of the structure as design variables, we design tunable and programmable Miuraori and eggbox strips, sheets, and composites that are large band, low frequency acoustic switches. In doing so, we find that the number of bandgaps in the sheets is significantly smaller than their strip counterparts and also occur at relatively higher frequencies, a limitation which is overcome by considering composite structures that have individual panels made of different materials. Overall, we have found origami structures to be ideal candidates as acoustic metamaterials for noise control, vibration isolation, impact absorption, and wave guides.


X. S. Zhang, G. H. Paulino, and A. S. Ramos Jr. "Multimaterial topology optimization with multiple volume constraints: Combining the ZPR update with a ground structure algorithm to select a single material per overlapping set. " International Journal of Numerical Methods in Engineering. Vol. 114, pp. 1053–1073, 2018. Link to PDF View/Hide Abstract
Multimaterial topology optimization often leads to members containing composite materials. However, in some instances, designers might be interested in using only one material for each member. Therefore, we propose an algorithm that selects a single preferred material from multiple candidate materials. We develop the algorithm, based on the evaluation of both the strain energy and the crosssectional area of each member, to control the material profile (i.e., number of materials) in each subdomain of the final design. This algorithm actively and iteratively selects materials to ensure a single material is used for each member. In this work, we adopt a multimaterial formulation that handles an arbitrary number of volume constraints and candidate materials. To efficiently handle such volume constraints, we employ the ZPR (ZhangPaulinoRamos) design update scheme for multimaterial optimization, which is based upon the separability of the dual objective function of the convex subproblem with respect to Lagrange multipliers. We provide an alternative derivation of this update scheme based on the KarushKuhnTucker (KKT) conditions. Through numerical examples, we demonstrate that the proposed algorithm, which can be readily implemented in multimaterial optimization, along with the ZPR update scheme, is robust and effective for selecting a single preferred material among multiple candidate materials.


Y. Jiang, T. Zegard, W. F. Baker, and G .H. Paulino. "Formfinding of gridshells using the ground structure and potential energy methods: a comparative study and assessment. " Journal of Structural and Multidisciplinary Optimization. Vol. 57, pp. 11871211, 2018. Link to PDF View/Hide Abstract
The structural performance of a gridshell
depends directly on the geometry of the design. Formfinding
methods, which are typically based on the search for
bendingfree configurations, aid in achieving structurally
efficient geometries. This manuscript proposes two formfinding
methods for gridshells: one method is the potential
energy method, which finds the form in equilibrium by
minimizing the total potential energy in the system; the
second method is based on an augmented version of the
ground structure method, in which the load application
points become variables of the topology optimization problem.
The proposed methods, together with the wellknown
force density method, are evaluated and compared using
numerical examples. The advantages and drawbacks of the
methods are reviewed, compared and highlighted.


O. GiraldoLondoño, D. W. Spring, G. H. Paulino, and W. G. Buttlar "An efficient mixedmode ratedependent cohesive fracture model using sigmoidal functions. " Engineering Fracture Mechanics. Vol. 192, pp. 307327, 2018. Link to PDF View/Hide Abstract
Ratedependent fracture processes can be investigated by means of cohesive zone models (CZMs). For instance, one approach enhances existing CZMs with phenomenological expressions used to represent the fracture energy, cohesive strength, and/or maximum crack opening as a function of the crack opening rate. Another approach assumes a viscoelastic CZM in front of the crack tip. Although computationally less expensive, the former approach misses most of the physics driving the ratedependent fracture process. The latter approach better represents the physics driving the ratedependent fracture process, yet it is computationally more expensive. This work presents a methodology for studying mixedmode ratedependent fracture that is both efficient and approximates the viscoelastic material behavior in front of the crack tip. In this mixedmode approach, we approximate the viscoelastic behavior in front of the crack tip using two ratedependent springs. One spring acts in the normal direction to the crack plane, while the other acts in the tangential direction. In order to mimic a viscoelastic CZM, we assume that the stiffness of each spring is a function of the crack opening rate and enforce that their tractions are continuous with respect to changes in the crack opening rates. To account for damage, we scale the tractions from the ratedependent springs using two damage parameters extracted from the ParkPaulinoRoesler (PPR) cohesive fracture model. The ratedependent model is implemented as a user defined element (UEL) subroutine in Abaqus. While attaining a high level of accuracy, the present approach allows for significant savings in computational cost when compared with a CZM based on fractional viscoelastic theory.


D. W. Spring and G. H. Paulino. "Achieving pervasive fracture and fragmentation in threedimensions: an unstructuringbased approach. " International Journal of Fracture. Vol. 210, pp. 113136, 2018. Link to PDF View/Hide Abstract
There are few methods capable of capturing the full spectrum of pervasive fracture behavior in threedimensions. Throughout pervasive fracture simulations, many cracks initiate, propagate, branch and coalesce simultaneously. Because of the cohesive element method framework, this behavior can be captured in a regularized manner. However, since the cohesive element method is only able to propagate cracks along element facets, a poorly designed discretization of the problem domain may introduce artifacts into the simulated results. To reduce the influence of the discretization, geometrically and constitutively unstructured means can be used. In this paper, we present and investigate the use of threedimensional nodal perturbation to introduce geometric randomness into a finite element mesh. We also discuss the use of statistical methods for introducing randomness in heterogeneous constitutive relations. The geometrically unstructured method of nodal perturbation is then combined with a random heterogeneous constitutive relation in three numerical examples. The examples are chosen in order to represent some of the significant influencing factors on pervasive fracture and fragmentation; including surface features, loading conditions, and material gradation. Finally, some concluding remarks and potential extensions are discussed.


J. Park, A. Sutradhar, J. J. Shah, and G. H. Paulino. "Design of complex bone internal structure using topology optimization with perimeter control. " Computers in Biology and Medicine. Vol. 94, pp. 7484, 2018. Link to PDF View/Hide Abstract
Large facial bone loss usually requires patientspecific bone implants to restore the structural integrity and functionality that also affects the appearance of each patient. Titanium alloys (e.g., Ti6Al4V) are typically used in the interfacial porous coatings between the implant and the surrounding bone to promote stability. There exists a property mismatch between the two that in general leads to complications such as stressshielding. This biomechanical discrepancy is a hurdle in the design of bone replacements. To alleviate the mismatch, the internal structure of the bone replacements should match that of the bone. Topology optimization has proven to be a good technique for designing bone replacements. However, the complex internal structure of the bone is difficult to mimic using conventional topology optimization methods without additional restrictions. In this work, the complex bone internal structure is recovered using a perimeter control based topology optimization approach. By restricting the solution space by means of the perimeter, the intricate design complexity of bones can be achieved. Three different bone regions with wellknown physiological loadings are selected to illustrate the method. Additionally, we found that the target perimeter value and the pattern of the initial distribution play a vital role in obtaining the natural curvatures in the bone internal structures as well as avoiding excessive island patterns..


X. S. Zhang, G. H. Paulino, and A. S. Ramos Jr. "Multimaterial topology optimization with multiple volume constraints: A general approach applied to ground structures with material nonlinearity. " Journal of Structural and Multidisciplinary Optimization. Vol. 57, pp. 161182, 2018. Link to PDF View/Hide Abstract
Multimaterial topology optimization is a practical tool that allows for improved structural designs. However, most studies are presented in the context of continuum topology optimization  few studies focus on truss topology optimization. Moreover, most work in this field has been restricted to linear material behavior with limited volume constraint settings for multiple materials. To address these issues, we propose an efficient multimaterial topology optimization formulation considering material nonlinearity. The proposed formulation handles an arbitrary number of candidate materials with flexible material properties, features freely specified material layers, and includes a generalized volume constraint setting. To efficiently handle such arbitrary volume constraints, we derive a design update scheme that performs robust updates of the design variables associated with each volume constraint independently. The derivation is based on the separable feature of the dual problem of the convex approximated primal subproblem with respect to the Lagrange multipliers, and thus the update of design variables in each volume constraint only depends on the corresponding Lagrange multiplier. Through examples in 2D and 3D, using combinations of Ogdenbased, bilinear, and linear materials, we demonstrate that the proposed multimaterial topology optimization framework with the presented update scheme leads to a design tool that not only finds the optimal topology but also selects the proper type and amount of material. The design update scheme is named ZPR (phonetically, zipper), after the initials of the authors' last names (ZhangPaulinoRamos Jr.).


Back to Top 
2017  
K. Liu and G. H. Paulino."Nonlinear mechanics of nonrigid origami: an efficient computational approach. " Proceedings of the Royal Society A. Vol. 473, pp. 128, 2017. Link to PDF View/Hide Abstract
Origamiinspired designs possess attractive
applications to science and engineering (e.g.
deployable, selfassembling, adaptable systems).
The special geometric arrangement of panels and
creases gives rise to unique mechanical properties of
origami, such as reconfigurability, making origami
designs well suited for tunable structures. Although
often being ignored, origami structures exhibit
additional soft modes beyond rigid folding due to the
flexibility of thin sheets that further influence their
behaviour. Actual behaviour of origami structures
usually involves significant geometric nonlinearity,
which amplifies the influence of additional soft
modes. To investigate the nonlinear mechanics
of origami structures with deformable panels,
we present a structural engineering approach for
simulating the nonlinear response of nonrigid
origami structures. In this paper, we propose a fully
nonlinear, displacementbased implicit formulation
for performing static/quasistatic analyses of nonrigid
origami structures based on ‘barandhinge’
models. The formulation itself leads to an efficient
and robust numerical implementation. Agreement
between real models and numerical simulations
demonstrates the ability of the proposed approach to
capture key features of origami behaviour.


E. D. Sanders, A. S. Ramos Jr., and G. H. Paulino."A maximum filter for the ground structure method: An optimization tool to harness multiple structural designs. " Engineering Structures. Vol. 151, pp. 235252, 2017. Link to PDF View/Hide Abstract
The ground structure method seeks to approximate Michell optimal solutions for realworld design problems
requiring truss solutions. The single solution extracted from the ground structure is typically too
complex to realize directly in practice and is instead used to inform designer intuition about how the structure behaves. Additionally, a postprocessing step required to filter out unnecessary truss members in the final design often leads to structures that no longer satisfy global equilibrium. Here, a maximum filter is proposed that, in addition to guaranteeing structures that satisfy global equilibrium, leads to several design perspectives for a single problem and allows for increased user control over the complexity of the final design. Rather than applying a static filter in each optimization iteration, the maximum filter employs an interval reducing method (e.g., bisection)to find the maximum allowable filter value that can be imposed in a given optimization iteration such that the design space is reduced while preserving global equilibrium and limiting local increases in the objective function. Minimization of potential energy with Tikhonov regularization is adopted to solve the singular system of equilibrium equations resulting from the filtered designs. In addition to reducing the order of the state problem, the maximum filter reduces the order of the optimization problem to increase computational efficiency. Numerical examples are presented to demonstrate the capabilities of the maximum filter, including a problem with multiple load cases, and its use as an endfilter in the traditional plastic and nested elastic approaches of the ground structure method.


X. S. Zhang, E. de Sturler, and G. H. Paulino. "Stochastic sampling for deterministic structural topology optimization with many load cases: Densitybased and ground structure approaches. " Computer Methods in Applied Mechanics and Engineering. Vol. 325, pp. 463487, 2017. Link to PDF View/Hide Abstract
We propose an efficient probabilistic method to solve a fully deterministic problem  we present a randomized optimization approach that drastically reduces the enormous computational cost of optimizing designs under many load cases for both continuum and truss topology optimization. Practical structural designs by deterministic topology optimization typically involve many load cases, possibly hundreds or more. The optimal design minimizes a, possibly weighted, average of the compliance under each load case (or some other objective). This means that, in each optimization step, a large finite element problem must be solved for each load case, leading to an enormous computational effort. On the contrary, the proposed randomized optimization method with stochastic sampling requires the solution of only a few (e.g., 5 or 6) finite element problems (large linear systems) per optimization step. Based on simulated annealing, we introduce a damping scheme for the randomized approach. Through numerical examples in two and three dimensions, we demonstrate that the randomization algorithm drastically reduces computational cost to obtain similar final topologies and results (e.g., compliance) to those of standard algorithms. The results indicate that the damping scheme is effective and leads to rapid convergence of the proposed algorithm.


E. T. Filipov, K. Liu, T. Tachi, M. Schenk, and G. H. Paulino. "Bar and hinge models for scalable analysis of origami. " International Journal of Solids and Structures. Vol. 124, pp. 2645, 2017. Link to PDF View/Hide Abstract
Thin sheets assembled into three dimensional folding origami can have various applications from reconfigurable architectural structures to metamaterials with tunable properties. Simulating the elastic stiffness and estimating deformed shapes of these systems is important for conceptualizing and designing practical engineering structures. In this paper, we improve, verify, and test a simplified bar and hinge model that can simulate essential behaviors of origami. The model simulates three distinct behaviors: stretching and shearing of thin sheet panels; bending of the initially flat panels; and bending along prescribed fold lines. The model is simple and efficient, yet it can provide realistic representation of stiffness characteristics and deformed shapes of origami structures. The simplicity of this model makes it well suited for the growing community of origami engineers, and its efficiency makes it suitable for design problems such as optimization and parametrization of geometric origami variations.


K. Liu, J. Wu, G. H. Paulino, and H. J. Qi. "Programmable Deployment of Tensegrity Structures by StimulusResponsive Polymers. " Nature Scientific Reports. Vol. 7, Article No. 3511, pp. 18. 2017. Link to PDF View/Hide Abstract
Tensegrity structures with detached struts are naturally suitable for deployable applications, both in terrestrial and outerspace structures, as well as morphing devices. Composed of discontinuous struts and continuous cables, such systems are only structurally stable when selfstress is induced; otherwise, they lose the original geometrical configuration (while keeping the topology) and thus can be tightly packed. We exploit this feature by using stimulus responsive polymers to introduce a paradigm for creating actively deployable 3D structures with complex shapes. The shapechange of 3D printed smart materials adds an active dimension to the configurational space of some structural components. Then we achieve dramatic global volume expansion by amplifying componentwise deformations to global configurational change via the inherent deployability of tensegrity. Through modular design, we can generate active tensegrities that are relatively stiff yet resilient with various complexities. Such unique properties enable structural systems that can achieve gigantic shape change, making them ideal as a platform for super lightweight structures, shapechanging soft robots, morphing antenna and RF devices, and biomedical devices.


X. S. Zhang, A. S. Ramos Jr., and G. H. Paulino. "Material nonlinear topology optimization using the ground structure method with a discrete filtering scheme. " Journal of Structural and Multidisciplinary Optimization. Vol. 55, No. 6, pp. 20452072. 2017. Link to PDF View/Hide Abstract
Topology optimization of truss lattices, using the
ground structure method, is a practical engineering tool that
allows for improved structural designs. However, in general,
the final topology consists of a large number of undesirable
thin bars that may add artificial stiffness and degenerate the
condition of the system of equations, sometimes even leading
to an invalid structural system. Moreover, most work
in this field has been restricted to linear material behavior,
yet real materials generally display nonlinear behavior.
To address these issues, we present an efficient filtering
scheme, with reducedorder modeling, and demonstrate its
application to two and threedimensional topology optimization
of truss networks considering multiple load cases
and nonlinear constitutive behavior. The proposed scheme
accounts for proper load levels during the optimization
process, yielding the displacement field without artificial
stiffness by simply using the truss members that actually exist in the structure (spurious members are removed), and
improving convergence performance. The nonlinear solution
scheme is based on a NewtonRaphson approach with
line search, which is essential for convergence. In addition,
the use of reducedorder information significantly reduces
the size of the structural and optimization problems within
a few iterations, leading to drastically improved computational
performance. For instance, the application of our
method to a problem with approximately 1 million design
variables shows that the proposed filter algorithm, while
offering almost the same optimized structure, is more than
40 times faster than the standard ground structure method.


B. C. Hill, O. GiraldoLondoño, G. H. Paulino, and W. G. Buttlar. "Inverse estimation of cohesive fracture properties of asphalt mixtures using an optimization approach." Experimental Mechanics. Vol. 57, No. 4, pp. 637648. 2017. Link to PDF View/Hide Abstract
Tensile cracking in asphalt pavements due to vehicular
and thermal loads has become an experimental and numerical
research focus in the asphalt materials community.
Previous studies have used the discrete element method
(DEM) to study asphalt concrete fracture. These studies used
trialanderror to obtain local fracture properties such that the
DEMmodels approximate the experimental loadcrack mouth
opening displacement response. In the current study, we identify
the cohesive fracture properties of asphalt mixtures via a
nonlinear optimization method. The method encompasses a
comparative investigation of displacement fields obtained
using both digital image correlation (DIC) and heterogeneous
DEMfracture simulations. The proposed method is applied to
two standard fracture test geometries: the singleedge notched
beam test, SE(B), under threepoint bending, and the diskshaped
compact tension test, DC(T). For each test, the
Subset Splitting DIC algorithm is used to determine the displacement
field in a predefined region near the notch tip.
Then, a given number of DEM simulations are performed on
the same specimen. The DEM is used to simulate the fracture
of asphalt concrete with a linear softening cohesive contact
model, where fracturerelated properties (e.g., maximum tensile
force and maximum crack opening) are varied within a
predefined range. The difference between DIC and DEM displacement
fields for each set of fracture parameters is then
computed and converted to a continuous function via multivariate
Lagrange interpolation. Finally, we use a Newtonlike
optimization technique to minimize Lagrange multinomials,
yielding a set of fracture parameters that minimizes the difference
between the DEM and DIC displacement fields. The
optimized set of fracture parameters from this nonlinear optimization
procedure led to DEM results which are consistent
with the experimental results for both SE(B) and DC(T)
geometries.


H. Chi, L. Beirao da Veiga, and G. H. Paulino. "Some basic formulations of the Virtual Element Method (VEM) for finite deformations." Computer Methods in Applied Mechanics and Engineering. Vol 318, pp. 148192. 2017. Link to PDF View/Hide Abstract
We present a general virtual element method (VEM) framework for finite elasticity, which emphasizes two issues: elementlevel volume change (volume average of the determinant of the deformation gradient) and stabilization. To address the former issue, we provide exact evaluation of the average volume change in both 2D and 3D on properly constructed local displacement spaces. For the later issue, we provide a new stabilization scheme that is based on the trace of the material tangent modulus tensor, which captures highly heterogeneous and localized deformations. Two VEM formulations are presented: a twofield mixed and an equivalent displacementbased, which is free of volumetric locking. Convergence and accuracy of the VEM are verified by means on numerical examples, and engineering applications are demonstrated.


Back to Top 
2016  
G. H. Paulino and E. C. N. Silva. "Multiscale, Multifunctional and Functionally Graded Materials
(MM&FGM) Editorial." Mechanics Research Communications. Special Issue on Multiscale, Multifunctional and Functionally Graded Materials. Vol. 78, pp. 1. 2016. Link to PDF 

D. W. Spring, O. GiraldoLondoño, and G. H. Paulino. "A Study on the Thermodynamic Consistency of the ParkPaulinoRoesler (PPR) Cohesive Fracture Model." Mechanics Research Communications. Special Issue on Multiscale, Multifunctional and Functionally Graded Materials. Vol. 78, pp. 100109. 2016. Link to PDF View/Hide Abstract
Although the ParkPaulinoRoesler (PPR) potenitalbased cohesive zone fracture model was not derived based on a thermodynamics consistency principle, we investigate the thermodynamic consistency of the PPR model under conditions of loading, unloading and reloading. First, we present a general anisotropic Helmholtz free energy function. Then, we reformulate the PPR model into the anisotropic Helmholtz form, and investigate the consistency of the model and the various unloading/reloading relations, which have been proposed for use with the model. By recasting the PPR model into the Helmholtz form, we illustrate that the PPR cohesive potential, while not designed with thermodynamic consistency in mind, is thermodynamically consistent under the pure loading conditions for which it was designed (as expected). We also demonstrate that the unloading/reloading relations, which are commonly used with the PPR model, are not thermodynamically consistent; however, through our investigation, we develop a new coupled unloading/reloading relation, which maintains the thermodynamic consistency of the PPR cohesive model. The considerations addressed in this paper are aimed at achieving a better understanding of the PPR model and other models of similar nature.


K. Park, H. Choi, and G. H. Paulino. "Assessment of cohesive tractionseparation relationships in ABAQUS: A comparative study." Mechanics Research Communications. Special Issue on Multiscale, Multifunctional and Functionally Graded Materials. Vol 78, pp. 7178. 2016. Link to PDF View/Hide Abstract
The definition of a tractionseparation relationship is essential in cohesive zone models because it describes the nonlinear fracture process zone. A few models are investigated in this paper and a comparative study is conducted. Among various tractionseparation relationships, the one in Abaqus is assessed by evaluating the cohesive traction and its tangent stiffness according to a given separation path. The results demonstrate that the tractionseparation relationship in Abaqus can lead to nonphysical responses because of a pathological positive tangent stiffness under softening condition. This is reflected in cohesive tractions that increase and decrease repeatedly while the cohesive separation monotonically increases. Thus, together with supporting information, this paper conveys the message that a tractionseparation relationship should be developed and selected with great caution, especially under mixedmode conditions.


H. Chi and G. H. Paulino. "On variational formulations with rigidbody constraints for finite elasticity: Applications to 2D and 3D finite element simulations." Mechanics Research Communications. Special Issue on Multiscale, Multifunctional and Functionally Graded materials. Vol 78, pp. 1526. 2016. Link to PDF View/Hide Abstract
We investigate a recently proposed variational principle with rigidbody constraints and present an
extension of its implementation in three dimensional finite elasticity problems. Through numerical examples,
we illustrate that the proposed variational principle with rigidbody constraints is applicable to
both single field and mixed finite elements of arbitrary order and geometry, e.g. triangular/tetrahedral
and quadrilateral/hexagonal elements, in two and three dimensions. Moreover, we demonstrate that, as
compared to the commonly adopted approach of discretizing the rigid domains with standard finite elements,
the proposed formulation requires neither discretization nor numerical integration in the interior
of each rigid domain. As a comparative result, the variational formulation may reduce the total number
of degrees of freedom of the resulting finite element system and provide improved accuracy.


M. Zhou, G. Allaire, G. Cheng, J. Du, M. Gilbert, X. Guo, J. Guest, R. Haftka, A. Kim, T. Lewinski, K. Maute, J. Norato, N. Olhoff, G. H. Paulino, T. Sokol, M. Wang, R. J. Yang, and B. D. Youn. "Special issue dedicated to Founding Editor George Rozvany." Journal of Structural and Multidisciplinary Optimization. Vol. 54, pp. 11071111. 2016. Link to PDF View/Hide Abstract
This special issue is a collective effort from the Guest Editors
and contributors to pay tribute to Professor George Rozvany,
our Founding Editor, who passed away in July 2015. As the
editorial team mourned George’s passing, they collectively decided
to organize a special issue to honor him. George founded
this journal with Dr. Jaroslaw Sobieski in 1989, and was running
the journal as EditorinChief with highest dedication until
his death. He had a highly productive and long research life.


S. L. Vatanabe, T. N. Lippi, C. R. de Lima, G. H. Paulino, and E. C. N. Silva. "Topology optimization with manufacturing constraints: A unified projectionbased approach." Advances in Engineering Software. Vol. 100, pp. 97112. 2016. Link to PDF View/Hide Abstract
Despite being an effective and a general method to obtain optimal solutions, topology optimization generates
solutions with complex geometries, which are neither costeffective nor practical from a manufacturing
(industrial) perspective. Manufacturing constraint techniques based on a unified projectionbased
approach are presented herein to properly restrict the range of solutions to the optimization problem. The
traditional stiffness maximization problem is considered in conjunction with a novel projection scheme
for implementing constraints. Essentially, the present technique considers a domain of design variables
projected in a pseudodensity domain to find the solution. The relation between both domains is defined
by the projection function and variable mappings according to each constraint of interest. The following
constraints have been implemented: minimum member size, minimum hole size, symmetry, pattern repetition,
extrusion, turning, casting, forging and rolling. These constraints illustrate the ability of the projection
scheme to efficiently control the optimization solution (i.e. without adding a large computational
cost). Illustrative examples are provided in order to explore the manufacturing constraints in conjunction
with the unified projectionbased approach.


A. Pereira, C. Talischi, G. H. Paulino, I. F. M. Menezes, M. S. Carvalho. "Fluid flow topology optimization in PolyTop: stability and computational implementation." Journal of Structural and Multidisciplinary Optimization. Vol. 54, pp. 13451364. 2016. Link to PDF  Download PolyTopFluid View/Hide Abstract
We present a Matlab implementation of topology
optimization for fluid flow problems in the educational
computer code PolyTop (Talischi et al. 2012b). The underlying
formulation is the wellestablished porosity approach
of Borrvall and Petersson (2003), wherein a dissipative
term is introduced to impede the flow in the solid (nonfluid)
regions. Polygonal finite elements are used to obtain
a stable loworder discretization of the governing Stokes
equations for incompressible viscous flow. As a result, the
same mesh represents the design field as well as the velocity
and pressure fields that characterize its response. Owing
to the modular structure of PolyTop, incorporating new
physics, in this case modeling fluid flow, involves changes
that are limited mainly to the analysis routine. We provide
several numerical examples to illustrate the capabilities and
use of the code. To illustrate the modularity of the present
approach, we extend the implementation to accommodate
alternative formulations and cost functions. These include topology optimization formulations where both viscosity
and inverse permeability are functions of the design; and
flow control where the velocity at a certain location in the
domain is maximized in a prescribed direction.


X. S. Zhang, S. Maheshwari, A. S. Ramos Jr., and G. H. Paulino. "Macroelement and macropatch approaches to structural topology optimization using the ground structure method." ASCE Journal of Structural Engineering. Vol 142, No. 11, pp. 040160901 to 14. 2016. Link to PDF View/Hide Abstract
Topology optimization can be divided into continuum and discrete types, the latter being the emphasis of the present work. In the field of discrete structural topology optimization of trusses, the generation of an initial ground structure is crucial. Thus this paper examines the generation of ground structures for generic structural domains in two and three dimensions. It compares two methods of discretization, Voronoibased and structured quadrilateral discretizations, and proposes two simple and effective ground structure generation approaches: the macroelement and macropatch approaches. Both can be implemented with either types of discretization. This work presents several features of these approaches, including efficient generation of initial ground structures, a reduction in matrix bandwidth for the global stiffness matrix, finer control of bar connectivity, and a reduction in the number of overlapped bars. Generic examples and practical structural engineering designs are presented, they display the features of the proposed approaches, and highlight the comparison with results from either the literature, the traditional ground structure generation, or the continuum optimization method.


K. Liu, G. H. Paulino, and P. Gardoni. "Reliabilitybased topology optimization using a new method for sensitivity approximation  application to ground structures." Journal of Structural and Multidisciplinary Optimization. Vol 54, pp. 553571. 2016. Link to PDF View/Hide Abstract
This paper proposes an efficient gradientbased optimization approach for reliabilitybased topology optimization of structures under uncertainties. Our objective is to find the optimized topology of structures with minimum weight which also satisfy certain reliability requirements. In the literature, those problems are primarily performed with approaches that use a firstorder reliability method (FORM) to estimate the gradient of the probability of failure. However, these approaches may lead to deficient or even invalid results because the gradient of probabilistic constraints, calculated by first order approximation, might not be sufficiently accurate. To overcome this issue, a newly developed segmental multipoint linearization (SML) method is employed in the optimization approach for a more accurate estimation of the gradient of failure probability. Meanwhile, this implementation also improves the approximation of the probability evaluation at no extra cost. In general, adoption of the SML method leads to a more accurate and robust approach. Numerical examples show that the new approach, based on the SML method, is numerically stable and usually provides optimized structures that have more of the desired features than conventional FORMbased approaches. The present approach typically does not lead to a fully stressed design, and thus this feature will be verified by numerical examples.


A. Alhadeff, S. E. Leon, W. Celes, and G. H. Paulino. "Massively parallel adaptive mesh refinement and coarsening for dynamic fracture simulations." Engineering with Computers. Vol 32, pp. 533552. 2016. Link to PDF View/Hide Abstract
We use the graphical processing unit (GPU) to
perform dynamic fracture simulation using adaptively
refined and coarsened finite elements and the interelement
cohesive zone model. Due to the limited memory available
on the GPU, we created a specialized data structure for
efficient representation of the evolving mesh given. To
achieve maximum efficiency, we perform finite element
calculation on a nodal basis (i.e., by launching one thread
per node and collecting contributions from neighboring
elements) rather than by launching threads per element,
which requires expensive graph coloring schemes to avoid
concurrency issues. These developments made possible the
parallel adaptive mesh refinement and coarsening schemes
to systematically change the topology of the mesh. We
investigate aspects of the parallel implementation through
microbranching examples, which has been explored
experimentally and numerically in the literature. First, we
use a reducedscale version of the experimental specimen
to demonstrate the impact of variation in floating point
operations on the final fracture pattern. Interestingly, the
parallel approach adds some randomness into the finite
element simulation on the structured mesh in a similar way
as would be expected from a random mesh. Next, we take
advantage of the speedup of the implementation over a
similar serial implementation to simulate a specimen
whose size matches that of the actual experiment. At this
scale, we are able to make more direct comparisons to the
original experiment and find excellent agreement with
those results.


A. Sutradhar, J. Park, D. Carrau, T. H. Nguyen, M. J. Miller, and G. H. Paulino. "Designing patientspecific 3D printed craniofacial implants using a novel topology optimization method." Medical & Biological Engineering & Computing. Vol 54, No. 7, pp. 11231135. 2016. Link to PDF View/Hide Abstract
Large craniofacial defects require efficient
bone replacements which should not only provide good
aesthetics but also possess stable structural function. The
proposed work uses a novel multiresolution topology optimization
method to achieve the task. Using a compliance
minimization objective, patientspecific bone replacement
shapes can be designed for different clinical cases that
ensure revival of efficient load transfer mechanisms in the
midface. In this work, four clinical cases are introduced
and their respective patientspecific designs are obtained
using the proposed method. The optimized designs are then
virtually inserted into the defect to visually inspect the viability
of the design . Further, once the design is verified by
the reconstructive surgeon, prototypes are fabricated using
a 3D printer for validation. The robustness of the designs
are mechanically tested by subjecting them to a physiological
loading condition which mimics the masticatory activity.
The fullfield strain result through 3D image correlation
and the finite element analysis implies that the solution
can survive the maximum mastication of 120 lb. Also, the
designs have the potential to restore the buttress system
and provide the structural integrity. Using the topology
optimization framework in designing the bone replacement shapes would deliver surgeons new alternatives for rather
complicated midface reconstruction.


K. Liu, G. H. Paulino, and P. Gardoni." Segmental multipoint linearization for parameter sensitivity approximation in reliability analysis. Structural Safety. Vol. 62, pp. 101115. 2016. Link to PDF View/Hide Abstract
This paper proposes an efficient method named segmental multipoint linearization for accurate approximation
of the sensitivity of failure probability with respect to design parameters. The method is established
based on a piecewise linear fitting of the limit state surface and the analytical integral of the
gradient of the failure probability with respect to parameters in the limit state function. The proposed
method presents an attractive ratio of accuracy to computational cost. The general framework is scalable
such that, by adjusting its complexity, different levels of accuracy can be achieved. The method is especially
suitable to be implemented in gradientbased routines for solving reliabilitybased design optimization
problems where accuracy of the parameter sensitivity is essential for convergence to an
optimal solution.


A. S. Ramos Jr. and G. H. Paulino. "Filtering structures out of ground structures  a discrete filtering tool for structural design optimization."Journal of Structural and Multidisciplinary Optimization. Vol. 54, pp. 95116. 2016. Link to PDF View/Hide Abstract
This paper presents an efficient scheme to filter structures out of ground structures, which is implemented using a nested elastic formulation for compliance minimization. The approach uses physical variables and allows control of the minimum ratio between the minimum and maximum areas in the final topology. It leads to a singular problem which is solved using a Tikhonov regularization on the structural problem (rather than on the optimization problem). The filter allows a multiple choice solution in which the user can control the global equilibrium residual in the final structural topology and limit variations of the objective function between consecutive iterations (e.g. compliance). As a result, an unambiguous discrete solution is obtained where all the bars that belong to the topology have welldefined finite areas. This filter feature with explicit control on the member areas allows the user (e.g. engineer or architect) to play with different alternatives prior to selecting a specific structural configuration. Several examples are provided to illustrate the properties of the present approach and the fact that the technique does not always lead to a fully stressed design. The method is efficient in the sense that the finite element solution is computed on the filtered structure (reduced order model) rather than on the full ground structure.


H. Chi, C. Talischi, O. LopezPamies, and G. H. Paulino "A paradigm for higherorder polygonal elements in finite elasticity using a gradient correction scheme." Computer Methods in Applied Mechanics and Engineering. Vol. 306, pp. 216251. 2016. Link to PDF View/Hide Abstract
Recent studies have demonstrated that polygonal elements possess great potential in the study of nonlinear elastic materials
under finite deformations. On the one hand, these elements are well suited to model complex microstructures (e.g. particulate
microstructures and microstructures involving different length scales) and incorporating periodic boundary conditions. On the
other hand, polygonal elements are found to be more tolerant to large localized deformations than the standard finite elements,
and to produce more accurate results in bending and shear. With mixed formulations, lower order mixed polygonal elements are
also shown to be numerically stable on Voronoitype meshes without any additional stabilization treatment. However, polygonal
elements generally suffer from persistent consistency errors under mesh refinement with the commonly used numerical integration
schemes. As a result, nonconvergent finite element results typically occur, which severely limit their applications. In this work,
a general gradient correction scheme is adopted that restores the polynomial consistency by adding a minimal perturbation to the
gradient of the displacement field. With the correction scheme, the recovery of optimal convergence for solutions of displacementbased
and mixed formulations with both linear and quadratic displacement interpolants is confirmed by numerical studies of
several boundary value problems in finite elasticity. In addition, for mixed polygonal elements, the various choices of the pressure
field approximations are discussed, and their performance on stability and accuracy are numerically investigated. We present
applications of those elements in physicallybased examples including a study of filled elastomers with interphasial effect and a
qualitative comparison with cavitation experiments for fiber reinforced elastomers.


E. T. Filipov, J. Chun , G. H. Paulino, and J. Song. "Polygonal multiresolution topology optimization (PolyMTOP) for structural dynamics." Journal of Structural and Multidisciplinary Optimization. Vol. 53, pp. 673694. 2016. Link to PDF View/Hide Abstract
We use versatile polygonal elements along with a multiresolution scheme for topology optimization to achieve computationally efficient and high resolution designs for structural dynamics problems. The multiresolution scheme uses a coarse finite element mesh to perform the analysis, a fine design variable mesh for the optimization and a fine density variable mesh to represent the material distribution. The finite element discretization employs a conforming finite element mesh. The design variable and density discretizations employ either matching or nonmatching grids to provide a finer discretization for the density and design variables. Examples are shown for the optimization of structural eigenfrequencies and forced vibration problems.


T. Zegard and G. H. Paulino. "Bridging topology optimization and additive manufacturing."Journal of Structural and Multidisciplinary Optimization. Vol. 53, pp. 175192. 2016. Link to PDF  Download TOPSlicer v0.8.4 View/Hide Abstract
Topology optimization is a technique that allows for increasingly efficient designs with minimal a priori decisions. Because of the complexity and intricacy of the solutions obtained, these techniques were often bounded to research and theoretical studies. Additive manufacturing, a rapidly evolving field, fills the gap between topology optimization and application. Additive manufacturing has minimal limitations on the shape and complexity of the design, and is currently evolving towards new materials, higher precision and larger build sizes. Two topology optimization methods are addressed: the ground structure method and densitybased topology optimization. The results obtained from these topology optimization methods require some degree of postprocessing before they can be manufactured. A simple procedure is described by which output suitable for additive manufacturing can be generated. In this process, some inherent issues of the optimization technique may be magnified resulting in an unfeasible or bad product. In addition, this work aims to address some of these issues and propose methodologies by which they may be alleviated. The proposed framework has applications in a number of fields, with specific examples given from the fields of health, architecture and engineering. In addition, the generated output allows for simple communication, editing, and combination of the results into more complex designs. For the specific case of threedimensional densitybased topology optimization, a tool suitable for result inspection and generation of additive manufacturing output is also provided.


H. Chi, O. LopezPamies, and G. H. Paulino. "A variational formulation with rigidbody constraints for finite elasticity: theory, finite element implementation, and applications." Computational Mechanics. Vol. 57, pp. 325338. 2016. Link to PDF View/Hide Abstract
This paper presents a new variational principle
in finite elastostatics applicable to arbitrary elastic solids that
may contain constitutively rigid spatial domains (e.g., rigid
inclusions). The basic idea consists in describing the constitutive
rigid behavior of a given spatial domain as a set
of kinematic constraints over the boundary of the domain.
From a computational perspective, the proposed formulation
is shown to reduce to a set of algebraic constraints that can
be implemented efficiently in terms of both singlefield and
mixed finite elements of arbitrary order. For demonstration
purposes, applications of the proposed rigidbodyconstraint
formulation are illustrated within the context of elastomers,
reinforced with periodic and random distributions of rigid
filler particles, undergoing finite deformations.


E. T. Filipov, G. H. Paulino, and T. Tachi. "Origami tubes with reconfigurable polygonal crosssections." Proceedings of the Royal Society A. Vol. 472, pp. 123. 2016. Link to PDF View/Hide Abstract
Thin sheets can be assembled into origami tubes to
create a variety of deployable, reconfigurable and
mechanistically unique threedimensional structures.
We introduce and explore origami tubes with
polygonal, translational symmetric crosssections that
can reconfigure into numerous geometries. The
tubular structures satisfy the mathematical definitions
for flat and rigid foldability, meaning that they can
fully unfold from a flattened state with deformations
occurring only at the fold lines. The tubes do not need
to be straight and can be constructed to follow a nonlinear
curved line when deployed. The crosssection and
kinematics of the tubular structures can be reprogrammed
by changing the direction of folding at some folds. We
discuss the variety of tubular structures that can be
conceived and we show limitations that govern the
geometric design. We quantify the global stiffness of
the origami tubes through eigenvalue and structural
analyses and highlight the mechanical characteristics
of these systems. The twoscale nature of this work
indicates that, from a local viewpoint, the crosssections
of the polygonal tubes are reconfigurable
while, from a global viewpoint, deployable tubes of
desired shapes are achieved. This class of tubes has
potential applications ranging from pipes and microrobotics
to deployable architecture in buildings.


A. N. Amirkhanian, D. W. Spring, J. R. Roesler, and G. H. Paulino. "Forward and inverse analysis of concrete fracture using the diskshaped compact tension test." Journal of Testing and Evaluation. Vol. 44, No. 1, pp. 625634. 2016. Link to PDF View/Hide Abstract
A new concrete fracture geometry is presented, which can quantify multiple fracture
properties from a single specimen test. The diskshaped compact tension (DCT) geometry
allows specimens to be fabricated from laboratory cylinders or field cores. The DCT fracture
test characterizes the concrete’s critical stress intensity factor, K_{IC}, critical cracktip opening
displacement, CTOD_{c}, and initial fracture energy, G_{f}, as well as the specimendependent
total fracture energy, G_{F}. The DCTbased fracture properties have the same experimental
variation as the singleedge notched beam test. The experimentally derived fracture
parameters were implemented into a cohesive zone model, which enabled estimation of
concrete tensile strength from fieldextracted cores.


Back to Top 
2015  
T. Zegard and G. H. Paulino. "GRAND3 — Ground structure based topology optimization for arbitrary 3D domains using MATLAB." Journal of Structural and Multidisciplinary Optimization. Vol. 52, pp. 11611184. 2015. Link to PDF  Download GRAND3 v1.0 View/Hide Abstract
Since its introduction, the ground structure method has been used in the derivation of closed–form analytical solutions for optimal structures, as well as providing information on the optimal load–paths. Despite its long history, the method has seen little use in three– dimensional problems or in problems with non–orthogonal domains, mainly due to computational implementation difficulties. This work presents a methodology for ground structure based topology optimization in arbitrary three– dimensional (3D) domains. The proposed approach is able to address concave domains and with the possibility of holes. In addition, an easy–to–use implementation of the proposed algorithm for the optimization of least–weight trusses is described in detail. The method is verified against three–dimensional closed–form solutions available in the literature. By means of examples, various features of the 3D ground structure approach are assessed, including the ability of the method to provide solutions with different levels of detail. The source code for a MATLAB implementation of the method, named GRAND3 — GRound structure ANalysis and Design in 3D, is available in the (electronic) Supplementary Material accompanying this publication.


L. S. Duarte, W. Celes, A. Pereira, I. F. M. Menezes, and G. H. Paulino. "PolyTop++: an efficient alternative for serial and parallel topology optimization on CPUs & GPUs." Journal of Structural and Multidisciplinary Optimization. Vol. 52, pp. 845859. 2015. Link to PDF View/Hide Abstract
This paper presents the PolyTop++, an efficient
and modular framework for parallel structural topology
optimization using polygonal meshes. It consists of a C++
and CUDA (a parallel computing model for GPUs) alternative
implementations of the PolyTop code by Talischi
et al. (Struct Multidiscip Optim 45(3):329–357 2012b).
PolyTop++ was designed to support both CPU and GPU
parallel solutions. The software takes advantage of the C++
programming language and the CUDA model to design
algorithms with efficient memory management, capable of
solving largescale problems, and uses its objectoriented
flexibility in order to provide a modular scheme. We
describe our implementation of different solvers for the
finite element analysis, including both direct and iterative
solvers, and an iterative ‘matrixfree’ solver; these were all
implemented in serial and parallel modes, including a GPU
version. Finally, we present numerical results for problems
with about 40 million degrees of freedom both in 2D and 3D.


R. Hashemi, D. W. Spring, and G. H. Paulino. "On small deformation interfacial debonding in composite materials containing multicoated particles." Journal of Composite Materials. Vol. 49, No. 27, pp. 34393455. 2015. Link to PDF View/Hide Abstract
This paper presents an integrated theoretical and computational investigation into the macroscopic behavior of composite materials containing multiphase reinforcing particles with simultaneous nonlinear debonding along the microconstituent interfaces. The interfacial debonding is characterized by the nonlinear ParkPaulinoRoesler potentialbased cohesive zone model. The extended MoriTanaka method is employed as the basis for the theoretical model, which enables micromechanical formulations for composite materials with high particle volume fractions. The computational analysis is performed using a threedimensional finite elementbased cohesive zone model with intrinsic cohesive elements. To place the generality and robustness of the proposed technique in perspective, we consider several examples of composite materials with single or double separation along the interfaces of coated particles. The effect of many microstructural parameters, such as the geometry of the microstructure, the location of debonding, the material properties of the coating layer (i.e. homogenous and functionally graded coatings), and the fracture parameters are comprehensively investigated by both theoretical and numerical approaches. We verify that both theoretical and numerical results agree well with one another in estimating the macroscopic constitutive relationship of corresponding composite materials. The strong dependence of the overall response of composite materials on their microstructure is well recognized for all hardening, snapback and softening stages.


L. L. Beghini, A. Beghini, W. F. Baker, and G. H. Paulino. "Integrated discrete/continuum topology optimization framework for stiffness or global stability of highrise buildings." Journal of Structural Engineering. Vol. 141, No. 8, pp. 040142071 to 10. 2015. Link to PDF View/Hide Abstract
This paper describes an integrated topology optimization framework using discrete and continuum elements for buckling and
stiffness optimization of highrise buildings. The discrete (beam/truss) elements are optimized based on their crosssectional areas, whereas
the continuum (polygonal) elements are concurrently optimized using the commonly known density method. Emphasis is placed on
linearized buckling and stability as objectives. Several practical examples are given to establish benchmarks and illustrate the proposed
methodology for highrise building design.


E. T. Filipov, T. Tachi, and G. H. Paulino. "Origami tubes assembled into stiff, yet reconfigurable structures and metamaterials." Proceedings of the National Academy of Sciences. Vol. 112, No. 40. pp. 12321–12326. 2015. Link to PDF Link to Supporting Information View/Hide Abstract
Thin sheets have long been known to experience an increase in stiffness when they are bent, buckled, or assembled into smaller interlocking structures. We introduce a unique orientation for coupling rigidly foldable origami tubes in a “zipper” fashion that substantially increases the system stiffness and permits only one flexible deformation mode through which the structure can deploy. The flexible deployment of the tubular structures is permitted by localized bending of the origami along prescribed fold lines. All other deformation modes, such as global bending and twisting of the structural system, are substantially stiffer because the tubular assemblages are overconstrained and the thin sheets become engaged in tension and compression. The zippercoupled tubes yield an unusually large eigenvalue bandgap that represents the unique difference in stiffness between deformation modes. Furthermore, we couple compatible origami tubes into a variety of cellular assemblages that can enhance mechanical characteristics and geometric versatility, leading to a potential design paradigm for structures and metamaterials that can be deployed, stiffened, and tuned. The enhanced mechanical properties, versatility, and adaptivity of these thin sheet systems can provide practical solutions of varying geometric scales in science and engineering.


M. Eidini and G. H. Paulino. "Unraveling metamaterial properties in zigzagbase folded sheets." Science Advances. Vol. 1, No. 8, e1500224. 2015. Link to PDF Link to Supporting Information View/Hide Abstract
Creating complex spatial objects from a flat sheet of material using origami folding techniques has attracted attention
in science and engineering. In the present work, we use the geometric properties of partially folded zigzag
strips to better describe the kinematics of known zigzag/herringbonebase folded sheet metamaterials such as
Miuraori. Inspired by the kinematics of a one–degree of freedom zigzag strip, we introduce a class of cellular folded
mechanical metamaterials comprising different scales of zigzag strips. This class of patterns combines origami
folding techniques with kirigami. Using analytical and numerical models, we study the key mechanical properties
of the folded materials. We show that our class of patterns, by expanding on the design space of Miuraori, is
appropriate for a wide range of applications from mechanical metamaterials to deployable structures at small
and large scales. We further show that, depending on the geometry, these materials exhibit either negative or
positive inplane Poisson’s ratios. By introducing a class of zigzagbase materials in the current study, we unify
the concept of inplane Poisson’s ratio for similar materials in the literature and extend it to the class of zigzagbase
folded sheet materials.


A. Alhadeff, W. Celes, and G. H. Paulino. "Mapping cohesive fracture and fragmentation simulations to graphics processor units." International Journal for Numerical Methods in Engineering. Vol. 103, No. 12, pp. 859893. 2015. Link to PDF View/Hide Abstract
A graphics processor units(GPU)based computational framework is presented to deal with dynamic failure
events simulated by means of cohesive zone elements. The work is divided into two parts. In the first part,
we deal with preprocessing of the information and verify the effectiveness of dynamic insertion of cohesive
elements in large meshes in parallel. To this effect, we employ a novel and simplified topological data
structure specialized for meshes with triangles, designed to run efficiently and minimize memory occupancy
on the GPU. In the second part, we present a parallel explicit dynamics code that implements an extrinsic
cohesive zone formulation where the elements are inserted ‘onthefly’, when needed and where needed.
The main challenge for implementing a GPUbased computational framework using an extrinsic cohesive
zone formulation resides on being able to dynamically adapt the mesh, in a consistent way, by inserting
cohesive elements on fractured facets. In order to handle that, we extend the conventional data structure
used in the finite element method (based on element incidence) and store, for each element, references to the
adjacent elements. This additional information suffices to consistently insert cohesive elements by duplicating
nodes when needed. Currently, our data structure is specialized for triangular meshes, but an extension
to tetrahedral meshes is feasible. The data structure is effective when used in conjunction with algorithms
to traverse nodes and elements. Results from parallel simulations show an increase in performance when
adopting strategies such as distributing different jobs among threads for the same element and launching
many threads per element. To avoid concurrency on accessing shared entities, we employ graph coloring. In
a preprocessing phase, each node of the dual graph (bulk elements of the mesh as graph nodes) is assigned
a color different from the colors assigned to adjacent nodes. In that fashion, elements of the same color
can be processed in parallel without concurrency. All the procedures needed for the insertion of cohesive
elements along fracture facets and for computing nodal properties are performed by threads assigned to triangles,
invoking one kernel per color. Computations on existing cohesive elements are also performed based
on adjacent bulk elements. Experiments show that GPU speedup increases with the number of nodes and
bulk elements.


D. W. Spring and G. H. Paulino. "Computational homogenization of the debonding of particle reinforced composites: The role of interphases in interfaces." Computational Materials Science. Vol. 109, pp. 209–224. 2015. Link to PDF View/Hide Abstract
There are four primary factors which influence the macroscopic constitutive response of particle reinforced composites: component properties, component concentrations, interphases, and interfacial debonding. Interphases are often a byproduct of surface treatments applied to the particles to control agglomeration. Alternatively, in polymer based materials such as carbonblack reinforced rubber, an interphase or ``bound rubber'' phase often occurs at the particlematrix interface. This interphasial region has been known to exist for many decades, but is often omitted in computational investigations of such composites. In this paper, we present an investigation into the influence of interphases on the large deformation response of particle reinforced composites. In addition, since particles tend to debond from the matrix at large deformations, we investigate the influence of interfacial debonding on the macroscopic constitutive response. The investigation considers two different microstructures; both a simplified single particle model, and a more complex polydisperse representative unit cell. Cohesive elements, which follow the ParkPaulinoRoesler tractionseparation relation, are inserted between each particle and its corresponding interphase to account for debonding. To account for friction, we present a new, coupled cohesivefriction model and detail its formulation and implementation. For each microstructure, we discuss the influence of the interphase thickness and stiffness on the global constitutive response in both uniaxial tension and simple shear. To validate the computational framework, comparisons are made with experimental results available in the literature.


G. H. Paulino and A. L. Gain. "Bridging art and engineering using Escherbased virtual elements." Journal of Structural and Multidisciplinary Optimization. Vol. 51, pp. 867883. 2015. Link to PDF View/Hide Abstract
The geometric shape of an element plays a key role in computational methods. Triangular and quadrilateral shaped elements are utilized by standard finite element methods. The pioneering work of Wachspress laid the foundation for polygonal interpolants which introduced polygonal elements. Tessellations may be considered as the next stage of element shape evolution. In this work, we investigate the topology optimization of tessellations as a means to coalesce art and engineering. We mainly focus onM.C. Escher’s tessellations using recognizable figures. To solve the state equation, we utilize a Mimetic Finite Difference inspired approach, known as the Virtual Element Method. In this approach, the stiffness matrix is constructed in such a way that the displacement patch test is passed exactly in order to ensure optimum numerical convergence rates. Prior to exploring the artistic aspects of topology optimization designs, numerical verification studies such as the displacement patch test and shear loaded cantilever beam bending problem are conducted to demonstrate the accuracy of the present approach in twodimensions.


A. L. Gain, G. H. Paulino, L. S. Duarte, and I. F. M. Menezes. "Topology optimization using polytopes." Computer Methods in Applied Mechanics and Engineering. Vol. 293, pp. 411430. 2015. Link to PDF View/Hide Abstract
Meshing complex engineering domains is a challenging task. Arbitrary polyhedral meshes
can provide the much needed flexibility in automated discretization of such domains. The geometric
property of polyhedral meshes such as its unstructured nature and the connectivity
of faces between elements makes them specially attractive for topology optimization applications.
Numerical anomalies in designs such as the single node connections and checkerboard
pattern can be naturally circumvented with polyhedrons. In the current work, we solve
the governing threedimensional elasticity state equation using the Virtual Element Method
(VEM) approach. The main characteristic difference between VEM and standard finite element
methods (FEM) is that in VEM the canonical basis functions are not constructed
explicitly. Rather the stiffness matrix is computed directly utilizing a projection map which
extracts the linear component of the deformation. Such a construction guarantees the satisfaction
of the patch test (used by engineers as an indicator of optimal convergence of
numerical solutions under mesh refinement). Finally, the computations reduce to the evaluation
of matrices which contain purely geometric surface facet quantities. The present work
focuses on the firstorder VEM in which the degrees of freedom are associated with the vertices.
The features of the current optimization approach are demonstrated using numerical
examples for compliance minimization and compliant mechanism problems.


T. Goudarzi, D. W. Spring, G. H. Paulino, and O. LopezPamies. "Filled elastomers: A theory of filler reinforcement based on hydrodynamic and interphasial effects." Journal of the Mechanics and Physics of Solids. Vol. 80, pp. 3767. 2015. Link to PDF View/Hide Abstract
Experimental evidence has by now established that i) the hydrodynamic effect and ii) the presence of stiff interphases (commonly referred to as bound rubber) “bonding” the underlying elastomer to the fillers are the dominant microscopic mechanisms typically responsible for the enhanced macroscopic stiffness of filled elastomers. Yet, because of the technical difficulties of dealing with these finescale effects within the realm of finite deformations, the theoretical reproduction of the macroscopic mechanical response of filled elastomers has remained an open problem.
The object of this paper is to put forward a microscopic field theory with the capability to describe, explain, and predict the macroscopic response of filled elastomers under arbitrarily large nonlinear elastic deformations directly in terms of: i) the nonlinear elastic properties of the elastomeric matrix, ii) the concentration of filler particles, and iii) the thickness and stiffness of the surrounding interphases. Attention is restricted to the prominent case of isotropic incompressible elastomers filled with a random and isotropic distribution of comparatively rigid fillers. The central idea of the theory rests on the construction of a homogenization solution for the fundamental problem of a Gaussian elastomer filled with a dilute concentration of rigid spherical particles bonded through Gaussian interphases of constant thickness, and on the extension of this solution to nonGaussian elastomers filled with finite concentrations of particles and interphases by means of a combination of iterative and variational techniques. For demonstration purposes, the theory is compared with full 3D finiteelement simulations of the large deformation response of Gaussian and nonGaussian elastomers reinforced by isotropic distributions of rigid spherical particles bonded through interphases of various finite sizes and stiffnesses, as well as with experimental data available from the literature. Good agreement is found in all of these comparisons. The implications of this agreement are discussed. 

J. Chun, J. Song, and G. H. Paulino. "Parameter sensitivity of system reliability using sequential compounding method." Structural Safety. Vol. 55, pp. 2636. 2015. Link to PDF View/Hide Abstract
Computation of sensitivities of the ‘system’ failure probability with respect to various parameters is essential in reliability based design optimization (RBDO) and uncertainty/risk management of a complex engineering system. The system failure event is defined as a logical function of multiple component events representing failure modes, locations or time points. Recently, the sequential compounding method (SCM) was developed for efficient calculations of the probabilities of largesize, general system events for a wide range of correlation properties. To facilitate the use of SCM in RBDO and uncertainty/risk management under a constraint on the system failure probability, a method, termed as Chun–Song–Paulino (CSP) method, is developed in this paper to compute parameter sensitivities of system failure probability using SCM. For a parallel or series system, the derivative of the system failure probability with respect to the reliability index is analytically derived at the last step of the sequential compounding. For a general system, the sensitivity of the probability of the set involving the component of interest and the sensitivity of the system failure probability with respect to the supercomponent representing the set are computed respectively using the CSP method and combined by the chainrule. The CSP method is illustrated by numerical examples, and successfully tested by examples covering a wide range of system event types, reliability indices, number of components, and correlation properties. The method is also applied to compute the sensitivity of the firstpassage probability of a building structure under stochastic excitations, modeled by use of finite elements.


C. Talischi, A. Pereira, I.F. Menezes, and G. H. Paulino. "Gradient correction for polygonal and polyhedral finite elements." International Journal for Numerical Methods in Engineering. Vol. 102, No. 34, pp. 728747. 2015. Special Issue: Tribute to Ted Belytschko. Link to PDF View/Hide Abstract
Previous studies have shown that the commonly used quadrature schemes for polygonal and polyhedral finite elements lead to consistency errors that persist under mesh refinement and subsequently render the approximations nonconvergent. In this work, we consider minimal perturbations to the gradient field at the element level in order to restore polynomial consistency and recover optimal convergence rates when the weak form integrals are evaluated using quadrature. For finite elements of arbitrary order, we state the accuracy requirements on the underlying volumetric and boundary quadrature rules and discuss the properties of the resulting corrected gradient operator. We compare the proposed approach with the pseudoderivative method developed by Belytschko and coworkers and, for linear elliptic problems, with our previous remedy that involves splitting of polynomial and nonpolynomial of elemental energy bilinear form. We present several numerical results for linear and nonlinear elliptic problems in two and three dimensions that not only confirm the recovery of optimal convergence rates but also suggest that the global error levels are close to those of approximations obtained from exact evaluation of the weak form integrals.


M. Alfano, G. Lubineau, and G. H. Paulino. "Special issue computational and experimental mechanics of advanced materials preface." International Journal of Solids and Structures. Vol. 55, pp. 11, 2015. Link to PDF 

M. Alfano, G. Lubineau, and G. H. Paulino. "Global sensitivity analysis in the identification of cohesive models using fullfield kinematic data." International Journal of Solids and Structures. Vol. 55, pp. 6678, 2015. Link to PDF View/Hide Abstract
Failure of adhesive bonded structures often occurs concurrent with the formation of a nonnegligible fracture process zone in front of a macroscopic crack. For this reason, the analysis of damage and fracture is effectively carried out using the cohesive zone model (CZM). The crucial aspect of the CZM approach is the precise determination of the tractionseparation relation. Yet it is usually determined empirically, by using calibration procedures combining experimental data, such as loaddisplacement or crack length data, with finite element simulation of fracture. Thanks to the recent progress in image processing, and the availability of lowcost CCD cameras, it is nowadays relatively easy to access surface displacements across the fracture process zone using for instance Digital Image Correlation (DIC). The rich information provided by correlation techniques prompted the development of versatile inverse parameter identification procedures combining finite element (FE) simulations and full field kinematic data. The focus of the present paper is to assess the effectiveness of these methods in the identification of cohesive zone models. In particular, the analysis is developed in the framework of the variance based global sensitivity analysis. The sensitivity of kinematic data to the sought cohesive properties is explored through the computation of the socalled Sobol sensitivity indexes. The results show that the global sensitivity analysis can help to ascertain the most influential cohesive parameters which need to be incorporated in the identification process. In addition, it is shown that suitable displacement sampling in time and space can lead to optimized measurements for identification purposes.


A. S. Ramos Jr. and G. H. Paulino. "Convex topology optimization for hyperelastic trusses based on the groundstructure approach." Journal of Structural and Multidisciplinary Optimization. Vol. 51, No. 2, pp. 287304. 2015. Link to PDF View/Hide Abstract
Most papers in the literature, which deal with
topology optimization of trusses using the ground structure
approach, are constrained to linear behavior. Here
we address the problem considering material nonlinear
behavior. More specifically, we concentrate on hyperelastic
models, such as the ones by Hencky, SaintVenant, Neo
Hookean and Ogden. A unified approach is adopted using the
total potential energy concept, i.e., the total potential is used
both in the objective function of the optimization problem and
also to obtain the equilibrium solution. We proof that the
optimization formulation is convex provided that the specific
strain energy is strictly convex. Some representative examples
are given to demonstrate the features of each model. We conclude
by exploring the role of nonlinearities in the overall
topology design problem.


C.Talischi and G. H. Paulino. "A closer look at consistent operator splitting and its extensions for topology optimization." Computer Methods in Applied Mechanics and Engineering. Vol. 283, pp. 573598. 2015. Link to PDF View/Hide Abstract
In this work, we explore the use of operator splitting algorithms for solving regularized structural topology optimization problems. The context is a classical structural design problem (e.g., compliance minimization and compliant mechanism design), parametrized by means of density functions, whose illposedness is addressed by introducing a Tikhonov regularization term. The proposed forward–backward splitting algorithm treats the constituent terms of the cost functional separately, which allows for suitable approximations of the structural objective. We will show that one such approximation, inspired by the reciprocal expansions underlying the optimality criteria method, improves the convergence characteristics and leads to an update scheme resembling the heuristic sensitivity filtering method. We also discuss a twometric variant of the splitting algorithm that removes the computational overhead associated with bound constraints on the density field without compromising convergence and quality of optimal solutions. We present several numerical results and investigate the influence of various algorithmic parameters.


H. Chi, C. Talischi, O. LopezPamies, and G. H. Paulino. "Polygonal finite elements for finite elasticity." International Journal of Numerical Methods in Engineering. Vol. 101, pp. 305328. 2015. Link to PDF View/Hide Abstract
Nonlinear elastic materials are of great engineering interest, but challenging to model with standard finite elements. The challenges arise because nonlinear elastic materials are characterized by nonconvex storedenergy functions as a result of their ability to undergo large reversible deformations, are incompressible or nearly incompressible, and often times possess complex microstructures. In this work, we propose and explore an alternative approach to model finite elasticity problems in two dimensions by using polygonal discretizations. We present both lower order displacementbased and mixed polygonal finite element approximations, the latter of which consist of a piecewise constant pressure field and a linearlycomplete displacement field at the element level. Through numerical studies, the mixed polygonal finite elements are shown to be stable and convergent. For demonstration purposes, we deploy the proposed polygonal discretization to study the nonlinear elastic response of rubber filled with random and periodic distributions of rigid particles, as well as the development of cavitation instabilities in elastomers containing vacuous defects. These physicallybased examples illustrate the potential of polygonal finite elements in studying and modeling nonlinear elastic materials with complex microstructures under finite deformations.


Back to Top 
2014  
S. E. Leon, D. Spring and G.H. Paulino. "Reduction in mesh bias for dynamic fracture using adaptive splitting of polygonal finite elements." International Journal of Numerical Methods in Engineering. Vol. 100, pp. 555576. 2014. Link to PDF View/Hide Abstract
We present a method to reduce mesh bias in dynamic fracture simulations using the finite element method with adaptive insertion of extrinsic cohesive zone elements along element boundaries. The geometry of the domain discretization is important in this setting because cracks are only allowed to propagate along element facets and can potentially bias the crack paths. To reduce mesh bias, we consider unstructured polygonal finite elements in this work. The meshes are generated with centroidal Voronoi tessellations to ensure element quality. However, the possible crack directions at each node are limited, making this discretization a poor candidate for dynamic fracture simulation. To overcome this problem, and significantly improve crack patterns, we propose adaptive element splitting, whereby the number of potential crack directions is increased at each crack tip. Thus, the crack is allowed to propagate through the polygonal element. Geometric studies illustrate the benefits of polygonal element discretizations employed with element splitting over other structured and unstructured discretizations for crack propagation applications. Numerical examples are performed and demonstrate good agreement with previous experimental and numerical results in the literature.


A. L. Gain, C. Talischi, and G. H. Paulino. "On the Virtual Element Method for threedimensional linear elasticity problems on arbitrary polyhedral meshes." Computer Methods in Applied Mechanics and Engineering. Vol. 282, pp. 132160. 2014. Link to PDF View/Hide Abstract
We explore the recentlyproposed Virtual Element Method (VEM) for the numerical solution of boundary value problems on arbitrary polyhedral meshes. More specifically, we focus on the linear elasticity equations in threedimensions and elaborate upon the key concepts underlying the firstorder VEM. While the point of departure is a conforming Galerkin framework, the distinguishing feature of VEM is that it does not require an explicit computation of the trial and test spaces, thereby circumventing a barrier to standard finite element discretizations on arbitrary grids. At the heart of the method is a particular kinematic decomposition of element deformation states which, in turn, leads to a corresponding decomposition of strain energy. By capturing the energy of linear deformations exactly, one can guarantee satisfaction of the patch test and optimal convergence of numerical solutions. The decomposition itself is enabled by local projection maps that appropriately extract the rigid body motion and constant strain components of the deformation. As we show, computing these projection maps and subsequently the local stiffness matrices, in practice, reduces to the computation of purely geometric quantities. In addition to discussing aspects of implementation of the method, we present several numerical studies in order to verify convergence of the VEM and evaluate its performance for various types of meshes.


T. Zegard and G.H. Paulino. "GRAND — Ground structure based topology optimization for arbitrary 2D domains using MATLAB." Journal of Structural and Multidisciplinary Optimization. Vol. 50, No. 5, pp. 861882. 2014. Link to PDF  Download GRAND v1.0 View/Hide Abstract
The present work describes in detail an imple mentation of the ground structure method for non– orthogonal unstructured and concave domains written in MATLAB, called GRAND — GRound structure ANalysis and Design. The actual computational implementation is provided, and example problems are given for educational and testing purposes. The problem of ground structure generation is translated into a linear algebra approach, which is inspired by the video–game literature. To prevent the ground structure generation algorithm from creating mem bers within geometric entities that no member should inter sect (e.g. holes, passive regions), the concept of “restriction zones” is employed, which is based on collision detection algorithms used in computational geometry and video– games. The aim of the work is to provide an easy–to–use implementation for the optimization of least–weight trusses embedded in any domain geometry.


D. Spring, S. E. Leon, and G. H. Paulino. "Unstructured polygonal meshes with adaptive refinement for the numerical simulation of dynamic cohesive fracture." International Journal of Fracture. Vol. 189, pp. 3357. 2014. Link to PDF View/Hide Abstract
This paper presents a scheme for adaptive mesh refinement on unstructured polygonal meshes to better capture crack patterns in dynamic cohesive fracture simulations. A randomly seeded polygonal mesh leads to an isotropic discretization of the problem domain, which does not bias crack patterns, but restricts the number of paths that a crack may travel at each node. An adaptive refinement scheme is proposed and investigated through a detailed set of geometric studies. The refinement scheme is selectively chosen to optimize the number of paths that a crack may travel, while still maintaining a conforming domain discretization. The details of the refinement scheme are outlined, along with the criterion used to determine the region of refinement and the method of interpolating nodal attributes. Extrinsic cohesive elements are inserted when and where necessary, and follow the constitutive response of the Park–Paulino–Roesler cohesive model. The influence of bulk and cohesive material heterogeneity is investigated through the use of a statistical distribution of material properties. The adaptive mesh modifications are handled through a compact topological data structure. Numerical examples highlight the features of adaptive refinement in capturing physical fracture patterns while addressing computational cost. Thus, the present approach is a step towards obtaining accurate dynamic fracture patterns and fields with polygonal elements.


S. Vatanabe, G. H. Paulino, and E. C. N. Silva. "Maximizing phononic band gaps in piezocomposite materials by means of topology optimization." Journal of the Acoustical Society of America. Vol. 136, pp. 494501. 2014. Link to PDF View/Hide Abstract
Phononic crystals (PCs) can exhibit phononic band gaps within which sound and vibrations at certain frequencies do not propagate. In fact, PCs with large band gaps are of great interest for many applications, such as transducers, elastic/acoustic filters, noise control, and vibration shields. Previous work in the field concentrated on PCs made of elastic isotropic materials; however, band gaps can be enlarged by using nonisotropic materials, such as piezoelectric materials. Because the main property of PCs is the presence of band gaps, one possible way to design microstructures that have a desired band gap is through topology optimization. Thus in this work, the main objective is to maximize the width of absolute elastic wave band gaps in piezocomposite materials designed by means of topology optimization. For band gap calculation, the finite element analysis is implemented with Bloch–Floquet theory to solve the dynamic behavior of twodimensional piezocomposite unit cells. Higher order frequency branches are investigated. The results demonstrate that tunable phononic band gaps in piezocomposite materials can be designed by means of the present methodology.


D. Spring and G. H. Paulino. "A growing library of threedimensional cohesive elements
for use in ABAQUS." Engineering Fracture Mechanics. Vol. 126, pp. 190216. 2014. Link to PDF View/Hide Abstract
In this paper, we present the implementation of a small library of threedimensional cohesive elements. The elements are formatted as userdefined elements, for compatibility with the commercial finite element software ABAQUS. The PPR, potentialbased traction–separation relation is chosen to describe the element’s constitutive model. The intrinsic cohesive formulation is outlined due to its compatibility with the standard, implicit finite element framework present in ABAQUS. The implementation of the cohesive elements is described, along with instructions on how to incorporate the elements into a finite element mesh. Specific areas of the userdefined elements, in which the user may wish to modify the code to meet specific research needs, are highlighted. Numerical examples are provided which display the capabilities of the elements in both small deformation and finite deformation regimes. A sample element source code is provided in an appendix, and the source codes of the elements are supplied through the website of the research group of the authors.


T. Zegard, W. Baker, A. Mazurek, and G. H. Paulino. "Geometrical aspects of lateral bracing systems: Where should the optimal bracing point be?" ASCE Journal of Structural Engineering. Vol. 140, No. 9, 040140631 to 9. 2014. Link to PDF View/Hide Abstract
This paper describes structural optimization formulations for trusses and applies them to size and geometrically optimize lateral bracing systems. The paper provides analytical solutions for simple and limit cases and uses numerical optimization for cases where analytical solutions are not available. The analysis is based on small displacements, constant or no connection costs, static loads, and elastic behavior. This work provides guidance to improve the common engineering practice of locating the bracing point at the middle or at the top of the bay/story panel.


L. L. Aguiar, C. A. Almeida, and G. H. Paulino. "A threedimensional multilayered pipe beam element: Nonlinear
analysis." Computers and Structures. Vol. 138, pp. 142161. 2014. Link to PDF View/Hide Abstract
This paper addresses the behavior of threedimensional multilayered pipe beams with interlayer slip
condition, under general 3D large displacements, in global riser and pipeline analysis applications. A
new finite element model, considering the Timoshenko beam for each element layer, has been formulated
and implemented. It comprises axial, bending and torsional degreesoffreedom, all varying along the
length of the element according to discretization using Hermitian functions: constant axial and torsional
loadings, and linear bending moments. Transverse shear strains due to bending are included in the formulation
by including two generalized constant degreesoffreedom. Nonlinear contact conditions,
together with various friction conditions between the element layers, are also considered. These are
accounted in the model through a proper representation of the constitutive relation for the shear stress
behavior in the binding material. The element formulation and its numerical capabilities are evaluated by
some numerical testing results, which are compared to other numerical or analytical solutions available
in the literature. The tests results show that the new element provides a simple yet robust and reliable
tool for general multilayered piping analyses.


L. L. Beghini, A. Beghini, N. Katz, W. F. Baker, and G. H. Paulino. "Connecting architecture and engineering through structural topology optimization."Engineering Structures. Vol. 59, pp. 716726. 2014. Link to PDF View/Hide Abstract
One of the prevalent issues facing the construction industry in today’s world is the balance between engineering and architecture: traditionally, the goal of the architect has focused more on the aesthetics, or ‘‘form’’ of a structure, while the goal of the engineer has been focused on stability and efﬁciency, or its ‘‘function’’. In this work, we discuss the importance of a close collaboration between these disciplines, and offer an alternative approach to generate new, integrated design ideas by means of a tailored structural topology optimization framework, which can potentially be of beneﬁt to both the architectural and structural engineering communities. Several practical case studies, from actual collaborative design projects, are given to illustrate the successes and limitations of such techniques.


L. L. Beghini, A. Pereira, R. Espinha, I. F. M. Menezes, W. Celes, and G. H. Paulino. "An objectoriented framework for finite element analysis based
on a compact topological data structure." Advances in Engineering Software. Vol. 68, pp. 4048. Feb. 2014. Link to PDF View/Hide Abstract
This paper describes an ongoing work in the development of a finite element analysis system, called
TopFEM, based on the compact topological data structure, TopS
[1,2]
. This new framework was written
to take advantage of the topological data structure together with objectoriented programming concepts
to handle a variety of finite element problems, spanning from fracture mechanics to topology optimiza
tion, in an efficient, but generic fashion. The class organization of the TopFEM system is described and
discussed within the context of other frameworks in the literature that share similar ideas, such as Get
FEM++, deal.II, FEMOOP and OpenSees. Numerical examples are given to illustrate the capabilities of TopS
attached to a finite element framework in the context of fracture mechanics and to establish a benchmark
with other implementations that do not make use of a topological data structure.


A. Cerrone, P. Wawrzynek, G. H. Paulino, and A. Ingraffea. "Implementation and verification of the Park–Paulino–Roesler
cohesive zone model in 3D." Engineering Fracture Mechanics. Vol. 120, pp. 2642. April 2014. Link to PDF View/Hide Abstract
The Park–Paulino–Roesler (PPR) potentialbased model is a cohesive constitutive model
formulated to be consistent under a high degree of modemixity. Herein, the PPR’s generalization
to threedimensions is detailed, its implementation in a finite element framework
is discussed, and its use in singlecore and high performance computing (HPC) applications
is demonstrated. The PPR model is shown to be an effective constitutive model to account
for crack nucleation and propagation in a variety of applications including adhesives,
composites, linepipe steel, and microstructures.


S. E. Leon, E. N. Lages, C. N. de Araujo and G. H. Paulino. "On the effect of constraint parameters on the generalized displacement control method." Mechanics Research Communications. Vol. 56, pp. 153159, 2014. Link to PDF View/Hide Abstract
In this work, we investigate the generalized displacement control method (GDCM) and provide a modification (MGDCM) that results in an equivalent constraint equation as that of the linearized cylindrical arclength control method (LCALCM). Through numerical examples, we illustrate that the MGDCM is more robust than the standard GDCM in capturing equilibrium paths in regions of high curvature. Moreover, we also provide a geometric and physical interpretation of the method, which sheds light on the general class of path following methods in structural mechanics.


C. Talischi, A. Pereira, G. H. Paulino, I. F. M. Menezes and M. S. Carvalho. "Polygonal finite elements for incompressible fluid flow." International Journal for Numerical Methods in Fluids. Vol. 74, No. 2, pp. 134151. 2014. Link to PDF View/Hide Abstract
We discuss the use of polygonal finite elements for analysis of incompressible flow problems. It is wellknown that the stability of mixed finite element discretizations is governed by the socalled infsup condition, which, in this case, depends on the choice of the discrete velocity and pressure spaces. We present a loworder choice of these spaces defined over convex polygonal partitions of the domain that satisfies the infsup condition and, as such, does not admit spurious pressure modes or exhibit locking. Within each element, the pressure field is constant while the velocity is represented by the usual isoparametric transformation of a linearlycomplete basis. Thus, from a practical point of view, the implementation of the method is classical and does not require any special treatment. We present numerical results for both incompressible Stokes and stationary Navier–Stokes problems to verify the theoretical results regarding stability and convergence of the method


C. Talischi and G. H. Paulino. "Addressing integration error for polygonal finite elements through polynomial projections: A patch test connection." Mathematical Models and Methods in Applied Sciences. Vol. 24, No. 8, pp. 17011727. 2014. Link to PDF View/Hide Abstract
Polygonal finite elements generally do not pass the patch test as a result of quadrature
error in the evaluation of weak form integrals. In this work, we examine the consequences of
lack of polynomial consistency and show that it can lead to a deterioration of convergence
of the finite element solutions. We propose a general remedy, inspired by techniques in
the recent literature of mimetic finite differences, for restoring consistency and thereby
ensuring the satisfaction of the patch test and recovering optimal rates of convergence.
The proposed approach, based on polynomial projections of the basis functions, allows for
the use of moderate number of integration points and brings the computational cost of
polygonal finite elements closer to that of the commonly used linear triangles and bilinear
quadrilaterals. Numerical studies of a twodimensional scalar diffusion problem accompany
the theoretical considerations.


Back to Top 
2013  
E. Dave, W. G. Buttlar, S. E. Leon, B. Behnia, and G. H. Paulino. "IlliTC – lowtemperature cracking model for asphalt pavements." Road Materials and Pavement Design. Vol. 14, No. S2, pp. 5778. 2013. Link to PDF View/Hide Abstract
Lowtemperature cracking (LTC) is a major distress and cause of failure for asphalt pavements located in regions with cold climate; however, most pavement design methods do not directly address LTC. The thermal cracking model (TCModel) utilised by American Association of State Highway and Transportation Officials MechanisticEmpirical Pavement Design Guide relies heavily on phenomenological Paris law for crack propagation. The TCModel predictions are primarily based on tensile strength of asphalt mixture and do not account for quasibrittle behaviour of asphalt concrete. Furthermore, TCModel utilises a simplified onedimensional viscoelastic solution for the determination of thermally induced stresses. This article describes a newly developed comprehensive software system for LTC prediction in asphalt pavements. The software system called ‘IlliTC’ utilises a userfriendly graphical interface with a standalone finiteelementbased simulation programme. The system includes a preanalyser and data input generator module that develops a twodimensional finite element (FE) pavement model for the user and which identifies critical events for thermal cracking using an efficient viscoelastic pavement stress simulation algorithm. Cooling events that are identified as critical are rigor ously simulated using a viscoelastic FE analysis engine coupled with a fractureenergybased cohesive zone fracture model. This article presents a comprehensive summary of the compo nents of the IlliTC system. Model verifications, field calibration and preliminary validation results are also presented.


R. Espinha, K. Park, G. H. Paulino and W. Celes. "Scalable parallel dynamic fracture simulation using an extrinsic cohesive
zone model." Computer Methods in Applied Mechanics and Engineering. Vol. 266, pp. 144161. 2013. Link to PDF View/Hide Abstract
In order to achieve realistic cohesive fracture simulation, a parallel computational framework is developed in conjunction with the parallel topology based data structure (ParTopS). Communications with remote partitions are performed by employing proxy nodes, proxy elements and ghost nodes, while synchronizations are identified on the basis of computational patterns (atnode, atelement, nodestoelement, and elementstonode). Several approaches to parallelize a serial code are discussed. An approach combining local computations and replicated computations with stable iterators is proposed, which is shown to be the most efficient one among the approaches discussed in this study. Furthermore, computational experiments demonstrate the scalability of the parallel dynamic fracture simulation framework for both 2D and 3D problems. The total execution time of a test problem remains nearly constant when the number of processors increases at the same rate as the number of elements.


A. L. Gain and G. H. Paulino. "A critical comparative assessment of differential
equationdriven methods for structural
topology optimization." Journal of Structural and Multidisciplinary Optimization. Vol. 48, No. 4, pp. 685710. 2013. Link to PDF View/Hide Abstract
In recent years, differential equationdriven
methods have emerged as an alternate approach for structural topology optimization. In such methods, the design is evolved using special differential equations. Implicit level
set methods are one such set of approaches in which the
design domain is represented in terms of implicit functions
and generally (but not necessarily) using the Hamilton
Jacobi equation as the evolution equation. Another set of
approaches are referred to as phasefield methods; which
generally use a reactiondiffusion equation, such as the
AllenCahn equation, for topology evolution. In this work,
we exhaustively analyze four levelset methods and one
phasefield method, which are representative of the literature. In order to evaluate performance, all the methods
are implemented in MATLAB and studied using two
dimensional compliance minimization problems.


C. Talischi and G. H. Paulino. "An operator splitting algorithm for Tikhonovregularized topology optimization." Computer Methods in Applied Mechanics and Engineering. Vol. 253, No. 1, pp. 599608. 2013. Link to PDF View/Hide Abstract
In this work, we investigate a Tikhonovtype regularization scheme to address the illposedness of the classical compliance minimization problem. We observe that a semiimplicit discretization of the gradient descent ﬂow for minimization of the regularized objective function leads to a convolution of the original gradient descent step with the Green’s
function associated with the modiﬁed Helmholtz equation. The appearance of “ﬁltering” in
this update scheme is different from the current density and sensitivity ﬁltering techniques
in the literature. The next iterate is deﬁned as the projection of this provisional density
onto the space of admissible density functions. For a particular choice of projection mapping, we show that the algorithm is identical to the wellknown forwardbackward splitting
algorithm, an insight that can be further explored for topology optimization. Also of interest is that with an appropriate choice of the projection parameter, nearly all intermediate
densities are eliminated in the optimal solution using the common density material models. We show examples of near binary solutions even for large values of the regularization
parameter.


S. L. Vatanabe, G. H. Paulino, and E. C. N. Silva. "Design of functionally graded piezocomposites using topology
optimization and homogenization – Toward effective energy harvesting
materials." Computer Methods in Applied Mechanics and Engineering. Vol. 266, No. 1, pp. 205218. 2013. Link to PDF View/Hide Abstract
In the optimization of a piezocomposite the objective is to obtain an improvement in its performance characteristics, usually by changing the volume fractions of constituent materials, its properties, shape of inclusions, and the mechanical properties of the polymer matrix in the composite unit cell. Thus, this work proposes a methodology based on topology optimization and homogenization to design functionally graded piezocomposite materials that considers important aspects in the design process aiming at energy harvesting applications, such as the influence of piezoelectric polarization directions and the influence of material gradation between the constituent materials in the unit cell. The influence of the piezoelectric polarization direction is quantitatively verified using the Discrete Material Optimization (DMO) method, which combines gradients with mathematical programming to solve a discrete optimization problem. The homogenization method is implemented using the graded finite element concept which takes into account the continuous gradation inside the finite elements. One of the main questions answered in this work is, quantitatively, how the microscopic stresses can be reduced by combining the functionally graded material (FGM) concept with optimization. In addition, the influence of polygonal elements is investigated, quantitatively, when compared to quadrilateral 4 node finite element meshes which are usually adopted in material design. However, quads exhibit onenode connections and are susceptible to checkerboard patterns in topology optimization applications. To circumvent these problems, Voronoi diagrams are used as an effective means of generating irregular polygonal meshes for piezocomposite design. The present results consist of bidimensional unit cells that illustrate the methodology proposed in this work.


T. Zegard and G. H. Paulino. "Toward GPU accelerated topology optimization on unstructured meshes." Journal of Structural and Multidisciplinary Optimization. Vol. 48, No. 3, pp. 473485. 2013. Link to PDF View/Hide Abstract
The present work investigates the feasibility of finite element methods and topology optimization for unstructured meshes in massively parallel computer architectures, more specifically on Graphics Processing Units or GPUs. Challenges in the parallel implementation, like the parallel assembly race condition, are discussed and solved with simple algorithms, in this case greedy graph coloring. The parallel implementation for every step involved in the topology optimization process is benchmarked and compared against an equivalent sequential implementation. The ultimate goal of this work is to speed up the topology optimization process by means of parallel computing using offtheshelf hardware. Examples are compared with both a standard sequential version of the implementation and a massively parallel version to better illustrate the advantages and disadvantages of this approach.


T. Zegard and G. H. Paulino. "Truss layout optimization within a continuum." Journal of Structural and Multidisciplinary Optimization. Vol. 48, No. 1, pp. 116. 2013. Link to PDF View/Hide Abstract
The present work extends truss layout optimization
by considering the case when it is embedded
in a continuum. Structural models often combine
discrete and continuum members and current requirements
for efficiency and extreme structures push research
in the eld of optimization. Examples of varied
complexity and dimensional space are analyzed and
compared, highlighting the advantages of the proposed
method. The goal of this work is to provide a simple
formulation for the discrete component of the structure,
more specically the truss, to be optimized in presence
of a continuum.


K. Park and G. H. Paulino. "Cohesive zone models: A critical review of tractionseparation relationships across fracture surfaces." Applied Mechanics Reviews. Vol. 64, No. 6, 0608021 to 20. 2011 (published in 2013). Link to PDF View/Hide Abstract
One of the fundamental aspects in cohesive zone modeling is the definition of the tractionseparation relationship across fracture surfaces, which approximates the nonlinear fracture process. Cohesive tractionseparation relationships may be classified as either nonpotentialbased models or potentialbased models. Potentialbased models are of special interest in the present review article. Several potentialbased models display limitations, especially for mixedmode problems, because of the boundary conditions associated with cohesive fracture. In addition, this paper shows that most effective displacementbased models can be formulated under a single framework. These models lead to positive stiffness under certain separation paths, contrary to general cohesive fracture phenomena wherein the increase of separation generally results in the decrease of failure resistance across the fracture surface (i.e., negative stiffness). To this end, the constitutive relationship of mixedmode cohesive fracture should be selected with great caution.


Back to Top 
2012  
M. Alfano, G. Lubineau, F. Furgiuele, and G. H. Paulino. "Study on the role of laser surface irradiation on damage
and decohesion of Al/epoxy joints." International Journal of Adhesion and Adhesives. Vol. 39, pp. 3341. 2012. Link to PDF View/Hide Abstract
In this work we investigate the effect of laser irradiation on the bond toughness of aluminum/epoxy bonded joints. The evolution of substrate surface morphology and wettability, for various sets of laser process parameters (i.e. laser power, line spacing, scan speed), was investigated by means of Scanning Electron Microscopy (SEM) and contact angle measurements. A proper combination of power, line spacing and scan speed was then selected and adhesive bonded Al/epoxy Tpeel joints were prepared and tested. For comparison, similar samples were produced using substrates with classical grit blasting surface treatment. Finally, postfailure SEM analyses of fracture surfaces were performed, and in order to typify the increase in bond toughness of the joints, finite element simulations were carried out using a potential based cohesive zone model of fracture.


E. V. Dave, G. H. Paulino, and W. G. Buttlar. "Viscoelastic functionally graded finite element method with recursive time integration and applications to flexible pavements." International Journal for Numerical and Analytical Methods in Geomechanics. Vol. 36, No. 9, pp. 11941219. 2012. Link to PDF View/Hide Abstract
The finite element method is used for modeling geotechnical and pavement structures exhibiting significant nonhomogeneity. Property gradients generated due to nonhomogeneous distributions of moisture is one such example for geotechnical materials. Aging and temperature induced property gradients are common sources of nonhomogeneity for asphalt pavements. Investigation of time dependent behavior combined with functionally graded property gradation can be accomplished by means of the nonhomogeneous viscoelastic analysis procedure. This paper describes the development of a generalized isoparametric finite element formulation to capture property gradients within elements, and a recursive formulation for solution of hereditary integral equations. The formulation is verified by comparison with analytical and numerical solutions. Two application examples are presented: the first describes stationary crack tip fields for viscoelastic functionally graded materials, and the second example demonstrates the application of the proposed procedures for efficient and accurate simulations of interfaces between layers of flexible pavement.


A. L. Gain and G. H. Paulino. "Phasefield based topology optimization with polygonal elements: a finite volume approach for the evolution equation." Journal of Structural and Multidisciplinary Optimization. Vol. 46, No. 3, pp. 327342. 2012. Link to PDF View/Hide Abstract
Uniform grids have been the common choice of domain discretization in the topology optimization literature.
Overconstraining geometrical features of such spatial discretizations can result in meshdependent, suboptimal
designs. Thus, in the current work, we employ unstructured polygonal meshes constructed using Voronoi tessellations
to conduct structural topology optimization. We utilize the phasefield method, derived from phase transition phenomenon,
which makes use of the AllenCahn differential equation and sensitivity analysis to update the evolving structural topology.
The solution of the AllenCahn evolution equation is accomplished by means of a centroidal Voronoi tessellation (CVT) based
finite volume approach. The unstructured polygonal meshes not only remove mesh bias but also provide greater flexibility in
discretizing complicated (e.g. nonCartesian) domains. The features of the current approach are demonstrated using various
numerical examples for compliance minimization and compliant mechanism problems.


T. H. Nguyen, G. H. Paulino, J. Song, and C. H. Le."Improving multiresolution topology optimization via multiple discretizations." International Journal for Numerical Methods in Engineering. Vol. 92, No. 6, pp. 507530. 2012. Link to PDF View/Hide Abstract
Unlike traditional topology optimization approach that uses the same discretization for finite
element analysis and design optimization, this paper proposes a framework for improving
multiresolution topology optimization (iMTOP) via multiple distinct discretizations for: (1)
finite elements, (2) design variables, and (3) density. This approach leads to high fidelity
resolution with a relatively low computational cost. In addition, an adaptive multiresolution
topology optimization procedure is introduced, which consists of selective adjustment and
refinement of design variable and density fields. Various two and threedimensional
numerical examples demonstrate that the proposed schemes can significantly reduce
computational cost in comparison to the existing elementbased approach.


K. Park and G. H. Paulino. "Computational implementation of the PPR potentialbased cohesive model in ABAQUS: Educational perspective." Engineering Fracture Mechanics. Vol. 93, pp. 239262. 2012. Link to PDF View/Hide Abstract
A potentialbased cohesive zone model, so called the PPR model, is implemented in a commercial software, e.g. ABAQUS, as a userdefined element (UEL) subroutine. The intrinsic cohesive zone modeling approach is employed because it can be formulated within the standard finite element framework. The implementation procedure for a twodimensional linear cohesive element and the algorithm for the PPR potentialbased model are presented indetail. The source code of the UEL subroutine is provided in the Appendix for educational purposes. Three computational examples are investigated to verify the PPR model and its implementation. The computational results of the model agree well with the analytical solutions.


K. Park, G. H. Paulino, W. Celes, and R. Espinha. "Adaptive mesh refinement and coarsening for cohesive dynamic fracture." International Journal for Numerical Methods in Engineering. Vol. 92, No. 1, pp. 135. 2012. Link to PDF View/Hide Abstract
Adaptive mesh refinement and coarsening (AMR&C) schemes are proposed for efficient com putational simulation of dynamic cohesive fracture. The adaptive mesh refinement (AMR) consists of a sequence of edgesplit operators while the adaptive mesh coarsening (AMC) is based on a sequence of vertexremoval (or edgecollapse) operators. Nodal perturbation and edgeswap operators are also employed around a crack tip region, and cohesive sur face elements are adaptively inserted whenever and wherever they are needed (i.e. extrinsic cohesive zone model). Such adaptive mesh modification events are maintained in conjunc tion with a topological data structure (TopS). The socalled PPR (ParkPaulinoRoesler) potentialbased cohesive model is utilized for the constitutive relationship of the cohesive zone model. The examples investigated include mode I fracture, mixedmode fracture and crack branching problems. The computational results using AMR&C are consistent with the results using uniform mesh refinement. The present approach significantly reduces compu tational cost while exhibiting a multiscale effect that captures both global macrocrack and local microcracks.


W. M. Rubio, G. H. Paulino, and E. N. C. Silva. "Analysis, manufacture and characterization of Ni/Cu functionally graded structures." Materials and Design. Vol. 41, pp. 255265. 2012. Link to PDF View/Hide Abstract
In this work, an experimental and numerical analysis and characterization of functionally graded structures
(FGSs) is developed. Nickel (Ni) and copper (Cu) materials are used as basic materials in the numerical
modeling and experimental characterization. For modeling, a MATLAB finite element code is
developed, which allows simulation of harmonic and modal analysis considering the graded finite
element formulation. For experimental characterization, Ni–Cu FGSs are manufactured by using spark
plasma sintering technique. Hardness and Young’s modulus are found by using microindentation and
ultrasonic measurements, respectively. The effective gradation of Ni/Cu FGS is addressed by means of
optical microscopy, energy dispersive spectrometry, scanning electron microscopy and hardness testing.
For the purpose of comparing modeling and experimental results, the hardness curve, along the gradation
direction, is used for identifying the gradation profile; accordingly, the experimental hardness curve is
used for approximating the Young’s modulus variation and the graded finite element modeling is used
for verification. For the first two resonance frequency values, a difference smaller than 1% between
simulated and experimental results is obtained.


L. L. Stromberg, A. Beghini, W. F. Baker, and G. H. Paulino. "Topology optimization for braced frames: Combining continuum and discrete elements." Engineering Structures. Vol. 37, pp. 106124. 2012. Link to PDF View/Hide Abstract
This paper describes an integrated topology optimization technique with concurrent use of both continuum fournode quadrilateral finite elements and discrete twonode beam elements to design structural braced frames that are part of the lateral system of a highrise building. The work explores the analytical aspects of optimal geometry for braced frames to understand the underlying mechanical phenomena and provides a theoretical benchmark to compare numerical results. The influence of the initial assumptions for the interaction between the quadrilaterals and the frame members are discussed. Numerical examples are given to illustrate the present technique on highrise building structures.


C. Talischi, G. H. Paulino, A. Pereira, and I. F. M. Menezes. "PolyMesher: A generalpurpose mesh generator for polygonal elements written in Matlab." Journal of Structural and Multidisciplinary Optimization. Vol. 45, No. 3, pp. 309328. 2012. Link to PDF  Download PolyMesher v1.0  PolyMesher v1.1 View/Hide Abstract
We present a simple and robust Matlab code for polygonal mesh generation that relies on an implicit description of the domain geometry. The mesh generator can provide, among other things, the input needed for ﬁnite element and optimization codes that use linear convex polygons. In topology optimization, polygonal discretizations have been shown not to be susceptible to numerical instabilities such as checkerboard patterns in contrast to lower order triangular and quadrilaterial meshes. Also, the use of polygonal elements makes
possible meshing of complicated geometries with a selfcontained Matlab code. The main ingredients of the present mesh generator are the implicit description of the domain and
the centroidal Voronoi diagrams used for its discretization. The signed distance function provides all the essential information about the domain geometry and oﬀers great ﬂexibility
to construct a large class of domains via algebraic expressions. Examples are provided to illustrate the capabilities of the code, which is compact and has fewer than 135 lines.


C. Talischi, G. H. Paulino, A. Pereira, and I. F. M. Menezes. "PolyTop: a Matlab implementation of a general topology optimization framework using unstructured polygonal finite element meshes." Journal of Structural and Multidisciplinary Optimization. Vol. 45, No. 3, pp. 329357. 2012. Link to PDF  Download PolyTop v1.0  PolyTop v1.1 View/Hide Abstract
We present an efficient Matlab code for structural topology optimization that includes a general finite element routine based on isoparametric polygonal elements which can be viewed as the extension of linear triangles and bilinear quads. The code also features a modular structure in which the analysis routine and the optimization algorithm are sepa rated from the specific choice of topology optimization formulation. Within this framework, the finite element and sensitivity analysis routines contain no information related to the formulation and thus can be extended, developed and modified independently. We address issues pertaining to the use of unstructured meshes and arbitrary design domains in topol ogy optimization that have received little attention in the literature. Also, as part of our examination of the topology optimization problem, we review the various steps taken in casting the optimal shape problem as a sizing optimization problem. This endeavor allows us to isolate the finite element and geometric analysis parameters and how they are related to the design variables of the discrete optimization problem. The Matlab code is explained in detail and numerical examples are presented to illustrate the capabilities of the code.


S. L. Vatanabe, G. H. Paulino, and E. C. N. Silva. "Influence of pattern gradation on the design of piezocomposite energy
harvesting devices using topology optimization." Composites: Part B. Vol. 43, No. 6, pp. 26462654. 2012. Link to PDF View/Hide Abstract
Piezoelectric materials can be used to convert oscillatory mechanical energy into electrical energy.
Energy harvesting devices are designed to capture the ambient energy surrounding the electronics and
convert it into usable electrical energy. The design of energy harvesting devices is not obvious, requiring
optimization procedures. This paper investigates the influence of pattern gradation using topology opti
mization on the design of piezocomposite energy harvesting devices based on bending behavior. The
objective function consists of maximizing the electric power generated in a load resistor. A projection
scheme is employed to compute the element densities from design variables and control the length scale
of the material density. Examples of twodimensional piezocomposite energy harvesting devices are pre
sented and discussed using the proposed method. The numerical results illustrate that pattern gradation
constraints help to increase the electric power generated in a load resistor and guides the problem
toward a more stable solution.


S. E. Leon, G. H. Paulino, A. Pereira, I. F. M. Menezes, and E. N. Lages. "A unified library of nonlinear solution schemes." Applied Mechanics Reviews. Vol. 64, No. 4, Article 0408031 to 26. 2011 (published in 2012). Link to PDF View/Hide Abstract
Nonlinear problems are prevalent in structural and continuum mechanics, and there is high demand for
computational tools to solve these problems. Despite efforts to develop efficient and effective algorithms, one
single algorithm may not be capable of solving any and all nonlinear problems. A brief review of recent
nonlinear solution techniques is first presented. Emphasis, however, is placed on the review of load, displacement,
arclength, work, generalized displacement, and orthogonal residual control algorithms, which are
unified into a single framework. Each of these solution schemes differs in the use of a constraint equation for
the incrementaliterative procedure. The governing finite element equations and constraint equation for each
solution scheme are combined into a single matrix equation, which characterizes the unified approach. This
conceptual model leads naturally to an effective objectoriented implementation. Within the unified framework,
the strengths and weaknesses of the various solution schemes are examined through numerical examples.


Back to Top 
2011  
M. Alfano, F. Furgiuele, L. Pagnotta and G. H. Paulino. "Analysis of fracture in aluminum joints bonded with a bicomponent epoxy adhesive." Journal of Testing and Evaluation. Vol. 39, No. 2. 2011. Link to PDF View/Hide Abstract
Adhesive bonding is a viable alternative to traditional joining systems (e.g., riveting or welding) for a wide class of components
belonging to electronic, automotive, and aerospace industries. However, adhesive joints often contain flaws; therefore, the development of such
technology requires reliable knowledge of the corresponding fracture properties. In the present paper, the candidate mode I fracture toughness of
aluminum/epoxy joints is determined using a double cantilever beam fracture specimen. A proper data reduction scheme for fracture energy calculation
has been selected based on the results of a sensitivity analysis. Furthermore, a scanning electron microscope is used in order to explore the
locus of failure. Finally, the experimental findings are assessed by means of numerical simulations of crack growth carried out using a cohesive zone
model.


M. Alfano, G. Lubineau, F. Furgiuele, and G. H. Paulino. "On the enhancement of bond toughness for Al/epoxy Tpeel
joints with laser treated substrates." International Journal of Fracture. Vol. 171, No. 2, pp. 139150. 2011. Link to PDF View/Hide Abstract
The aim of the present work is to quantify
the enhancement of bond toughness of Al/epoxy joints
associated to substrates laser irradiation. For this reason
a potential based cohesivemodel is employed and cohesive
elements are implemented within the finite element
framework. The influence of the cohesive properties on
the predicted global response of the joints is firstly analyzed.
The coupling between adherents plasticity and
the cohesive properties is then discussed. It is shown
that the global response is mainly affected by cohesive
energy (the bond toughness) and cohesive strength. In
turn, a proper cost function is defined which quantifies
the deviation between numerical and experimental total
dissipated energy. Based on a sensitivity analysis of
the asdefined cost function, it is shown that an accurate
estimation of the bond toughness can be expected
from global data. The situation is different for the cohesive strength, whose estimation could require more
advanced experimental observations or additional tests.
The results reported in the presentwork allowus to conclude,
in a reliable manner, that the laser surface treatment
can lead to a large improvement of bond toughness.


C. A. Almeida, J. C. R. Albino, I. F. M. Menezes, and G. H. Paulino. "Geometric nonlinear analyses of functionally graded beams using a tailored Lagrangian formulation." Mechanics Research Communications. Vol. 38, No. 8, pp. 553559. 2011. Link to PDF View/Hide Abstract
This paper presents a geometric nonlinear analysis formulation for beams of functionally graded crosssections
by means of a Total Lagrangian formulation. The influence of material gradation on the numerical
response is investigated in detail. Two examples are given that illustrate the main features of the formulation,
in which the behavior of beams of graded crosssections is compared with homogeneous material
beams. A motivation for this work is the potential development of functionally graded risers for the
offshore oil exploration industry.


R. C. Carbonari, P. A. MunozRojas, E. Q. Andrade, G. H. Paulino, K. Nishimoto, and E. C. N. Silva. "Design of pressure vessels using shape optimization:
An integrated approach." International Journal of Pressure Vessels and Piping. Vol. 88, No. 57, pp. 198212. 2011. Link to PDF View/Hide Abstract
Previous papers related to the optimization of pressure vessels have considered the optimization of the nozzle independently from the dished end. This approach generates problems such as thickness variation from nozzle to dished end (coupling cylindrical region) and, as a consequence, it reduces the optimality of the final result which may also be influenced by the boundary conditions. Thus, this work discusses shape optimization of axisymmetric pressure vessels considering an integrated approach in which the entire pressure vessel model is used in conjunction with a multiobjective function that aims to minimize the vonMises mechanical stress from nozzle to head. Representative examples are examined and solutions obtained for the entire vessel considering temperature and pressure loading. It is noteworthy that different shapes from the usual ones are obtained. Even though such different shapes may not be profitable considering present manufacturing processes, they may be competitive for future manufacturing technologies, and contribute to a better understanding of the actual influence of shape in the behavior of pressure vessels.


E. V. Dave, A. F. Braham, W. G. Buttlar, and G. H. Paulino. "Development of a flattened indirect tension test for asphalt concrete." Journal of Testing and Evaluation. Vol. 39, No. 3. 2011. Link to PDF View/Hide Abstract
The indirect tension test (IDT) is frequently used in civil engineering because of its benefits over direct tension testing. In the mid1990s, an IDT protocol was developed for evaluating tensile strength and creep properties of asphalt concrete mixtures, as specified by the American Association of State Highway Transportation Officials (AASHTO) in AASHTO T322. However, with the increased use of finer aggregate gradations and polymer modified asphalt binders in asphalt concrete mixtures, the validity of IDT strength results can be questioned in instances where significant crushing occurs under the narrow loading heads. Therefore, a new specimen configuration is proposed for indirect tension testing of asphalt concrete. In place of the standard loading heads, the specimen was trimmed to produce flat planes with parallel faces, creating a “flattened IDT.” A viscoelastic finite element analysis of the flattened configuration was performed to evaluate the optimal trimming width. In addition, the numerically determined geometry was verified by means of laboratory testing of three asphalt concrete mixtures in two flattened configurations. This integrated modeling and testing study showed that when using fine aggregate gradations and compliant asphalt binders, crushing is significantly reduced while maintaining tensile stresses near the center of the specimen. Furthermore, creep compliances were evaluated using the flattened IDT and compared with those obtained following AASHTO T322. Some variation was observed between the creep properties evaluated from the different geometries, particularly for higher compliance values. As a preliminary assessment, the flattened IDT seems to be a suitable geometry for the evaluation of indirect tensile strength of asphalt concrete. Further testing and analysis should be performed on the flattened IDT arrangement for evaluation of the creep compliance. This study provides an initial step towards a possible revision of the current AASHTO standard for IDT testing of asphalt concrete mixtures.


E. V. Dave, G. H. Paulino, and W. G. Buttlar. "Viscoelastic functionally graded finite element method using correspondence principle." ASCE: Journal of Materials in Civil Engineering. Vol. 23, No. 1, pp. 3948. 2011. Link to PDF View/Hide Abstract
The capability to effectively discretize the problem domain makes the finite element method an attractive simulation technique for modeling complicated boundary value problems such as asphalt concrete pavements with material nonhomogeneities. Specialized “graded elements” have been shown to provide an efficient and accurate tool for the simulation of functionally graded materials. Most of the previous research on numerical simulation of functionally graded materials has been limited to elastic material behavior. Thus the current work focuses on finite element analysis of functionally graded viscoelastic materials. The analysis is performed using the elasticviscoelastic correspondence principle, and viscoelastic material gradation is accounted for within the elements by means of the generalized isoparametric formulation. This paper emphasizes viscoelastic behavior of asphalt concrete pavements and several examples, ranging from verification problems to field scale applications, are presented to demonstrate the features of the present approach.


A. L. Gain, J. Carroll, G. H. Paulino, and J. Lambros. "A hybrid experimental/numerical technique to extract cohesive fracture properties for modeI fracture of quasibrittle materials." International Journal of Fracture. Vol. 169, No. 2, pp. 113131. 2011. Link to PDF View/Hide Abstract
We propose a hybrid technique to extract cohesive fracture properties of a quasibrittle (not exhibiting bulk plasticity) material using an inverse numerical analysis and experimentation based on the optical technique of digital image correlation (DIC). Two options for the inverse analysis were used – a shape optimization approach, and a parameter optimization for a potentialbased cohesive constitutive model, the socalled PPR (ParkPaulinoRoesler) model. The unconstrained, derivative free NelderMead algorithm was used for optimization in the inverse analysis. The two proposed schemes were verified for realistic cases of varying initial guesses, and different synthetic and noisy displacement field data. As proof of concept, both schemes were applied to a Polymethylmethacrylate (PMMA) quasistatic crack growth experiment where the near tip displacement field was obtained experimentally by DIC and was used as input to the optimization schemes. The technique was successful in predicting the applied loaddisplacement response of a four point bend edge cracked fracture specimen.


Y. Liu, W. Ye, A. V. Phan, and G. H. Paulino. "Special issue on the advances in mesh reduction methods — In honor of Professor Subrata Mukherjee on the occasion of his 65th birthday." Engineering Analysis with Boundary Elements. Vol. 35, No. 2, p. 158. 2011. (Second of two special issues, Go to first special issue) Link to PDF 

T. H. Nguyen, J. Song, and G. H. Paulino. "Singleloop system reliabilitybased topology optimization considering statistical dependence between limitstates." Structural and Multidisciplinary Optimization. Vol. 44, No. 5, pp. 593611. 2011. Link to PDF View/Hide Abstract
This paper presents a singleloop algorithm for system reliabilitybased topology optimization (SRBTO) that can account for statistical dependence between multiple limitstates, and its applications to computationally demanding topology optimization problems. A singleloop reliabilitybased design optimization (RBDO) algorithm replaces the innerloop iterations to evaluate probabilistic constraints by a noniterative approximation. The proposed singleloop SRBTO algorithm accounts for the statistical dependence between the limitstates by using the matrixbased system reliability (MSR) method to compute the system failure probability and its parameter sensitivities. The SRBTO/MSR approach is applicable to general system events including series, parallel, cutset and linkset systems and provides the gradients of the system failure probability to facilitate gradientbased optimization. In most RBTO applications, probabilistic constraints are evaluated by use of the firstorder reliability method for efficiency. In order to improve the accuracy of the reliability calculations for RBDO or RBTO problems with high nonlinearity, we introduce a new singleloop RBDO scheme utilizing the secondorder reliability method and implement it to the proposed SRBTO algorithm. Moreover, in order to overcome challenges in applying the proposed algorithm to computationally demanding topology optimization problems, we utilize the multiresolution topology optimization (MTOP) method, which achieves computational efficiency in topology optimization by assigning different levels of resolutions to three meshes representing finite element analysis, design variables and material density distribution respectively. The paper provides numerical examples of two and threedimensional topology optimization problems to demonstrate the proposed SRBTO algorithm and its applications. The optimal topologies from deterministic, component and system RBTOs are compared to investigate the impact of optimization schemes on final topologies. Monte Carlo simulations are also performed to verify the accuracy of the failure probabilities computed by the proposed approach.


K. Park and G. H. Paulino. "Parallel computing of wave propagation in three dimensional functionally
graded media."Mechanics Research Communications. Vol. 38, No. 6, pp. 431436. 2011. Link to PDF View/Hide Abstract
Parallel computing techniques are employed to investigate wave propagation in threedimensional
functionally graded media. In order to obtain eﬀective and eﬃcient parallel ﬁnite element mesh representation, a
topologybased data structure (TopS) and a parallel framework for unstructured mesh (ParFUM) are integrated.
The parallel computing framework is veriﬁed by solving a cantilever example, while the Rayleigh wave speed in
functionally graded media is investigated by comparing the results with the homogeneous case. The computational
results illustrate that when the elastic modulus of a graded media increase along the depth direction, the Rayleigh
wave speed of a graded media is higher than the speed of a homogeneous media with the same material properties
on the surface.


W. M. Rubio, G. H. Paulino, and E. C. N. Silva. "Tailoring vibration mode shapes using topology optimization and functionally graded material concepts."Smart Materials and Structures. Vol. 20, No. 2, 025009. 2011. Link to PDF View/Hide Abstract
Tailoring specified vibration modes is a requirement for designing piezoelectric devices aimed
at dynamictype applications. A technique for designing the shape of specified vibration modes
is the topology optimization method (TOM) which finds an optimum material distribution inside
a design domain to obtain a structure that vibrates according to specified eigenfrequencies and
eigenmodes. Nevertheless, when the TOM is applied to dynamic problems, the wellknown
grayscale or intermediate material problem arises which can invalidate the postprocessing of
the optimal result. Thus, a more natural way for solving dynamic problems using TOM is to
allow intermediate material values. This idea leads to the functionally graded material (FGM)
concept. In fact, FGMs are materials whose properties and microstructure continuously change
along a specific direction. Therefore, in this paper, an approach is presented for tailoring
userdefined vibration modes, by applying the TOM and FGM concepts to design functionally
graded piezoelectric transducers (FGPT) and nonpiezoelectric structures (functionally graded
structures—FGS) in order to achieve maximum and/or minimum vibration amplitudes at certain
points of the structure, by simultaneously finding the topology and material gradation function.
The optimization problem is solved by using sequential linear programming. Twodimensional
results are presented to illustrate the method.


B. Shen and G. H. Paulino. "Direct extraction of cohesive fracture properties from digital image correlation: A hybrid inverse technique."Experimental Mechanics. Vol. 51, No. 2, pp. 143163. 2011. Link to PDF View/Hide Abstract
The accuracy of an adopted cohesive zone model (CZM) can affect the simulated fracture response significantly. The CZM has been usually obtained using global experimental response, e.g., load versus either crack opening displacement or loadline displacement. Apparently, deduction of a local material property from a global response does not provide full confidence of the adopted model. The difficulties are: (1) fundamentally, stress cannot be measured directly and the cohesive stress distribution is nonuniform; (2) accurate measurement of the full crack profile (crack opening displacement at every point) is experimentally difficult to obtain. An attractive feature of digital image correlation (DIC) is that it allows relatively accurate measurement of the whole displacement field on a flat surface. It has been utilized to measure the mode I tractionseparation relation. A hybrid inverse method based on combined use of DIC and finite element method is used in this study to compute the cohesive properties of a ductile adhesive, Devcon Plastic Welder II, and a quasibrittle plastic, G10/FR4 Garolite. Fracture tests were conducted on single edgenotched beam specimens (SENB) under fourpoint bending. A fullfield DIC algorithm was employed to compute the smooth and continuous displacement field, which is then used as input to a finite element model for inverse analysis through an optimization procedure. The unknown CZM is constructed using a flexible Bspline without any “a priori” assumption on the shape. The inversely computed CZMs for both materials yield consistent results. Finally, the computed CZMs are verified through fracture simulation, which shows good experimental agreement.


B. Shen and G. H. Paulino. "Identification of cohesive zone model and elastic parameters of fiberreinforced cementitious composites using digital image correlation and a hybrid inverse technique." Cement and Concrete Composites. Vol. 33, No. 5, pp. 572585. 2011. Link to PDF View/Hide Abstract
Traditional methods for the inverse identification of elastic properties and local cohesive zone model (CZM) of solids utilize only global experimental data. In contrast, this paper addresses the inverse identification of elastic properties and CZM of a range of materials, using fullfield displacement through an optimization technique in a finite element (FE) framework. The new experimental–numerical hybrid approach has been applied to fiberreinforced cementitious composites (FRCC). PVA microfibers are used at four volume fractions: 0.5%, 1%, 2% and 3%. Digital image correlation (DIC) technique is used to measure surface displacement fields of the test specimens. Fourpoint bend tests are carried out for the measurement of the modulus of elasticity, E, and the Poisson’s ratio, ν, while single edgenotched beams (SENB) are used for measurement of modeI CZM parameters. A finite element update inverse formulation, which is based on minimization of the difference between measured and computed displacement field, is used for both identification problems. For the identification of E and ν, linearized form of the Hooke’s tensor in plane stress condition has been derived for twodimensional linear elasticity in FE frame, and Newton–Raphson solver is employed for the inverse problem. For the identification of the CZM, generic spline curves have been used for the parameterization of any CZM thus avoiding the need of an assumption of the CZM shape, while derivativefree Nelder–Mead optimization with CZM shape regularization is employed as the solution method, which reduces the complexity of numerical implementation and improves robustness. The computed E and ν are consistent with published results. The computed CZMs of the FRCCs with different fiber volume fractions reveal a strainhardening characteristic. The computed CZM is used in direct problem simulation, the results of which are consistent with the experimental global response.


L. L. Stromberg, A. Beghini, W. F. Baker, and G. H. Paulino. "Application of layout and topology optimization using pattern gradation for the conceptual design of buildings." Structural and Multidisciplinary Optimization. Vol. 43, No. 2, pp. 165180. 2011. Link to PDF View/Hide Abstract
This paper explores the use of manufacturingtype constraints, in particular pattern gradation and repetition, in the context of building layout optimization. By placing constraints on the design domain in terms of number and variable size of repeating patterns along any direction, the conceptual design for buildings is facilitated. To substantiate the potential future applications of this work, examples within the context of highrise building design are presented. Successful development of such ideas can contribute to practical engineering solutions, especially during the building design process. Examples are given to illustrate the ideas developed both in twodimensional and threedimensional building configurations.


Y. Yu, G. H. Paulino, and Y. Luo. "Finite particle method for progressive failure simulation of truss structures."Journal of Structural Engineering. Vol. 137, No. 10, pp. 11681181. 2011. Link to PDF View/Hide Abstract
A structural analysis framework called the finite particle method (FPM) for structure failure simulation is presented in this paper. The traditional finiteelement method is generated from continuum mechanics and the variational principle; vector mechanics form the basis of FPM. It discretizes the domain with finite particles whose motions are described by Newton’s second law. Instead of imposing a global equilibrium of the entire continuous system, FPM enforces equilibrium on each particle. Thus, particles are free to separate from one another, which is advantageous in the simulation of structural failure. One of the features of this approach is that no iterations to follow nonlinear laws are necessary, and no global matrices are formed or solved in this method. A convected material frame is used to evaluate the structure deformation and internal force. The explicit time integration is adopted to solve the equation of motion. To simulate the truss structure failure, a failure criterion on the basis of the ideal plastic constitutive model and a failure modeling algorithm are proposed by using FPM. According to the energy conservation study of a twodimensional (2D) truss, the energy is decomposed and balanced during the failure process. Also, a more complicated threedimensional (3D) structure failure simulation is given. The comparison of the simulation results and the practical failure mode shows the capability of this method.


Back to Top 
2010  
S. R. M. Almeida, G. H. Paulino, and E. C. N. Silva. "Layout and material gradation in topology optimization of functionally graded structures: a globallocal approach." Structural and Multidisciplinary Optimization. Vol. 42, No. 6, pp. 855868. 2010. Link to PDF View/Hide Abstract
By means of continuous topology optimization,
this paper discusses the influence of material gradation
and layout in the overall stiffness behavior of functionally
graded structures. The formulation is associated to symmetry and pattern repetition constraints, including material gradation effects at both global and local levels. For instance,
constraints associated with pattern repetition are applied by
considering material gradation either on the global structure
or locally over the specific pattern. By means of pattern repetition, we recover previous results in the literature which
were obtained using homogenization and optimization of
cellular materials.


R. C. Carbonari, G. H. Paulino, and E. C. N. Silva. "Integral piezoactuator system with optimum placement of functionally graded material  a topology optimization paradigm." Journal of Intelligent Material Systems and Structures. Vol. 21, No. 6, pp. 16531668. 2010. Link to PDF View/Hide Abstract
Piezoactuators consist of compliant mechanisms actuated by two or more piezoceramic devices. During the assembling process, such flexible structures are usually bonded to the piezoceramics. The thin bonding layer(s) between the compliant mechanism and the piezoceramic may induce undesirable behavior, including unusual interfacial nonlinearities. This constitutes a drawback of piezoelectric actuators and, in some applications, such as those associated to vibration control and structural health monitoring (e. g., aircraft industry), their use may become either unfeasible or at least limited. A possible solution to this standing problem can be achieved through the functionally graded material concept and consists of developing 'integral piezoactuators', that is those with no bonding layer(s) and whose performance can be improved by tailoring their structural topology and material gradation. Thus, a topology optimization formulation is developed, which allows simultaneous distribution of void and functionally graded piezoelectric materials (including both piezo and nonpiezoelectric materials) in the design domain in order to achieve certain specified actuation movements. Two concurrent design problems are considered, that is the optimum design of the piezoceramic property gradation, and the design of the functionally graded structural topology. Twodimensional piezoactuator designs are investigated because the applications of interest consist of planar devices. Moreover, material gradation is considered in only one direction in order to account for manufacturability issues. To broaden the range of such devices in the field of smart structures, the design of integral Moonietype functionally graded piezoactuators is provided according to specified performance requirements.


Y. S. Chan, G. H. Paulino, B. F. Feng, and A. Sutradhar. "Dependence of crack tip singularity on loading functions." Mechanics Research Communications. Vol. 37, No. 2, pp. 191197. 2010. Link to PDF View/Hide Abstract
Under the theory of classical linear fracture mechanics, a finite crack sitting in an isotropic and homogeneous medium is considered. We find that the wellknown crack tip singularity, the inverse squareroot singularity 1/sqrt(r), may disappear under certain type of loading traction functions. More specifically, depending on the cracksurface loading function, the behavior of the crack tip field may be shown to be as smooth as possible. The singular integral equation method is used to study the dependence of the crack tip singularity on the loading traction functions. Exact crack opening displacements, stress fields, and their corresponding loading traction functions are provided. Although the method used is somewhat mathematically elementary, the outcome seems to be new and useful.


H. J. Choi and G. H. Paulino. "Interfacial cracking in a graded coating/substrate system loaded by a frictional sliding flat punch." Proceedings of the Royal Society of London. A: Mathematical, Physical and Engineering Sciences. Vol. 466, No. 2115, pp. 853880. 2010. Link to PDF Download Data Supplement View/Hide Abstract
An analysis of a coupled plane elasticity problem of crack/contact mechanics for a coating/substrate system with functionally graded properties is performed, where the rigid flat punch slides over the surface of the coated system that contains a crack. The graded material is treated as a nonhomogeneous interlayer between dissimilar, homogeneous phases of the coated medium and the crack is assumed to exist along the interface between the interlayer and the substrate. Based on the Fourier integral transform method and the transfer matrix approach, formulation of the current coupled mixed boundary value problem lends itself to the derivation of a set of three simultaneous Cauchytype singular integral equations. In the numerical results, the emphasis is placed on the investigation of interactions between the contact stress field and the cracktip behaviour for various combinations of material, geometric and loading parameters of the coated system. Specifically, effects of interfacial cracking on the distributions of the contact pressure and the inplane stress component along the coating surface are examined and the mixedmode stress intensity factors evaluated from the cracktip stress field with the squareroot singularity are provided as a function of punch location. Further addressed is the quantification of the singular character of contact pressure distributions at the trailing and leading edges of the flat punch in terms of the punchedge stress intensity factors. Implicit in this particular analysis of the coupled crack/contact problem presented henceforth is that the crack closure behaviour under the compressive contact stress field is not taken into account, ignoring the influence of crackface contact and friction.


Y. Liu, W. Ye, A. V. Phan, and G. H. Paulino. "Special issue on the advances in mesh reduction methods — In honor of Professor Subrata Mukherjee on the occasion of his 65th birthday." Engineering Analysis with Boundary Elements. Vol. 34, No. 11, pp. 902903. 2010. (First of two special issues, Go to second special issue) Link to PDF 

L. A. M. Mello, E. de Sturler, G. H. Paulino, and E. C. N. Silva. "Recycling Krylov subspaces for efficient largescale electrical
impedance tomography." Computer Methods in Applied Mechanics and Engineering. Vol. 199, No. 4952, pp. 31013110. 2010. Link to PDF View/Hide Abstract
Electrical impedance tomography (EIT) captures images of internal features of a body. Electrodes are attached to the boundary of the body, low intensity alternating currents are applied, and the resulting electric potentials are measured. Then, based on the measurements, an estimation algorithm obtains the threedimensional internal admittivity distribution that corresponds to the image. One of the main goals of medical EIT is to achieve high resolution and an accurate result at low computational cost. However, when the finite element method (FEM) is employed and the corresponding mesh is refined to increase resolution and accuracy, the computational cost increases substantially, especially in the estimation of absolute admittivity distributions. Therefore, we consider in this work a fast iterative solver for the forward problem, which was previously reported in the context of structural optimization. We propose several improvements to this solver to increase its performance in the EIT context. The solver is based on the recycling of approximate invariant subspaces, and it is applied to reduce the EIT computation time for a constant and high resolution finite element mesh. In addition, we consider a powerful preconditioner and provide a detailed pseudocode for the improved iterative solver. The numerical results show the effectiveness of our approach: the proposed algorithm is faster than the preconditioned conjugate gradient (CG) algorithm. The results also show that even on a standard PC without parallelization, a high mesh resolution (more than 150,000 degrees of freedom) can be used for image estimation at a relatively low computational cost.


D. Ngo, K. Park, G. H. Paulino, and Y. Huang. "On the constitutive relation of materials with microstructure using a potentialbased cohesive model for interface interaction." Engineering Fracture Mechanics. Vol. 77, No. 7, pp. 11531174. 2010. Link to PDF View/Hide Abstract
Macroscopic constitutive relationship is estimated by considering the microscopic particle/matrix interfacial debonding. For the interfacial debonding, the PPR potentialbased cohesive model is utilized. The extended Mori–Tanaka model is employed for micromechanics, while a finite elementbased cohesive zone model is used for the computational model. Both models (theoretical and computational) agree well each other in representing the macroscopic constitutive relationship on the basis of the PPR model. The microscopic interfacial cohesive parameters of the PPR model are estimated from macroscopic composite material behavior. In addition, different microscopic debonding processes are observed with respect to different macroscopic constitutive relationships (e.g. hardening, softening, and snapback).


T. H. Nguyen, G. H. Paulino, J. Song, and C. H. Le. "A computational paradigm for multiresolution topology optimization (MTOP)."Structural and Multidisciplinary Optimization. Vol. 41, No. 4, pp. 525539. 2010. Link to PDF View/Hide Abstract
This paper presents a multiresolution topology optimization (MTOP) scheme to obtain high resolution designs with relatively low computational cost. We employ three distinct discretization levels for the topology optimization procedure: the displacement mesh (or finite element mesh) to perform the analysis, the design variable mesh to perform the optimization, and the density mesh (or density element mesh) to represent material distribution and compute the stiffness matrices. We employ a coarser discretization for finite elements and finer discretization for both density elements and design variables. A projection scheme is employed to compute the element densities from design variables and control the length scale of the material density. We demonstrate via various two and threedimensional numerical examples that the resolution of the design can be significantly improved without refining the finite element mesh.


T. H. Nguyen, J. Song, and G. H. Paulino. "Singleloop system reliabilitybased design optimization using matrixbased system reliability method: theory and applications." ASME Journal of Mechanical Design. Vol. 132, No. 1, pp. 011005111. 2010. Link to PDF View/Hide Abstract
This paper proposes a singleloop system reliabilitybased design optimization (SRBDO)
approach using the recently developed matrixbased system reliability (MSR) method. A
singleloop method was employed to eliminate the innerloop of SRBDO that evaluates
probabilistic constraints. The MSR method enables us to compute the system failure
probability and its parameter sensitivities efficiently and accurately through convenient
matrix calculations. The SRBDO/MSR approach proposed in this paper is applicable to
general systems including series, parallel, cutset, and linkset system events. After a
brief overview on SRBDO algorithms and the MSR method, the SRBDO/MSR approach is
introduced and demonstrated by three numerical examples. The first example deals with
the optimal design of a combustion engine, in which the failure is described as a series
system event. In the second example, the crosssectional areas of the members of a
statically indeterminate truss structure are determined for minimum total weight with a
constraint on the probability of collapse. In the third example, the redistribution of the
loads caused by member failures is considered for the truss system in the second example.
The results based on different optimization approaches are compared for further investigation.
Monte Carlo simulation is performed in each example to confirm the accuracy of
the system failure probability computed by the MSR method.


K. Park, G. H. Paulino, and J. Roesler. "Cohesive fracture model for functionally graded fiber reinforced concrete." Cement and Concrete Research. Vol. 40, No. 6, pp. 956965. 2010. Link to PDF View/Hide Abstract
A simple, effective, and practical constitutive model for cohesive fracture of fiber reinforced concrete is proposed by differentiating the aggregate bridging zone and the fiber bridging zone. The aggregate bridging zone is related to the total fracture energy of plain concrete, while the fiber bridging zone is associated with the difference between the total fracture energy of fiber reinforced concrete and the total fracture energy of plain concrete. The cohesive fracture model is defined by experimental fracture parameters, which are obtained through threepoint bending and split tensile tests. As expected, the model describes fracture behavior of plain concrete beams. In addition, it predicts the fracture behavior of either fiber reinforced concrete beams or a combination of plain and fiber reinforced concrete functionally layered in a single beam specimen. The validated model is also applied to investigate continuously, functionally graded fiber reinforced concrete composites.


G. H. Paulino, K. Park, W. Celes, and R. Espinha. "Adaptive dynamic cohesive fracture simulation using nodal
perturbation and edgeswap operators." International Journal for Numerical Methods in Engineering. Vol. 84, No. 11, pp. 13031343. 2010. Link to PDF View/Hide Abstract
Dependence on mesh orientation impacts adversely the quality of computational solutions generated by
cohesive zone models. For instance, when considering crack propagation along interfaces between finite
elements of 4k structured meshes, both extension of crack length and crack angle are biased according
to the mesh configuration. To address mesh orientation dependence in 4k structured meshes and to avoid
undesirable crack patterns, we propose the use of nodal perturbation (NP) and edgeswap (ES) topological
operation. To this effect, the topological data structure TopS (Int. J. Numer. Meth. Engng 2005; 64:
1529–1556), based on topological entities (node, element, vertex, edge and facet), is utilized so that it is
possible to access adjacency information and to manage a consistent data structure in time proportional
to the number of retrieved entities. In particular, the data structure allows the ES operation to be done in
constant time. Three representative dynamic fracture examples using ES and NP operators are provided:
crack propagation in the compact compression specimen, local branching instability, and fragmentation.
These examples illustrate the features of the present computational framework in simulating a range of
physical phenomena associated with cracking.


B. Shen, I. Stanciulescu, and G. H. Paulino. "Inverse computation of cohesive fracture properties from displacement fields."Inverse Problems in Science and Engineering. Vol. 18, No. 8, pp. 11031128. 2010. Link to PDF View/Hide Abstract
Cohesive zone model (CZM) is a key technique for finite element simulation of
fracture of quasibrittle materials, yet it is usually determined empirically from global
measurement. A more convincing way to obtain the cohesive relation is to determine
experimentally the relation between crack separation and crack surface traction. Recent
developments in experimental mechanics such as photoelasticity and digital image
correlation (DIC) enable accurate measurement of whole field surface displacement. The
cohesive stress at the crack surface cannot be measured directly, but may be determined
indirectly through the displacement field near the crack surface. An inverse problem
thereby is formulated in order to extract the cohesive relation by fully utilizing the
measured displacement field. As the focus in this paper is to develop a framework to
solve the inverse problem effectively, synthetic displacement field data obtained from
finite element analysis (FEA) are used. First by assuming the cohesive relation with a few
governing parameters, a direct problem is solved to obtain the complete synthetic
displacement field at certain loading levels. The computed displacement field is then
assumed known, while the cohesive relation is solved in the inverse problem through the
unconstrained, derivativefree NelderMead (NM) optimization method. Linear and
cubic splines with an arbitrary number of control points are used to represent the shape of
the CZM. The unconstrained nature of NM method and the physical validity of the CZM
shape are guaranteed by introducing barrier terms. Comprehensive numerical tests are
carried out to investigate the sensitivity of the inverse procedure to experimental errors.
The results show that even at a high level of experimental error, the computed CZM is
still well estimated, which demonstrates the practical usefulness of the proposed method.
The technique introduced in this work can be generalized to compute other internal or
boundary stresses from whole displacement field using the finite element method.


A. Sutradhar, G. H. Paulino, M. J. Miller and T. H. Nguyen. "Topological optimization for designing patientspecific large craniofacial segmental bone replacements." Proceedings of the National Academy of Sciences. Vol. 107, No. 30, pp. 1322213227. 2010. Link to PDF Link to Supporting Information View/Hide Abstract
Restoring normal function and appearance after massive facial injuries with bone loss is an important unsolved problem in surgery. An important limitation of the current methods is heuristic ad hoc design of bone replacements by the operating surgeon at the time of surgery. This problem might be addressed by incorporating a computational method known as topological optimization into routine surgical planning. We tested the feasibility of using a multiresolution threedimensional topological optimization to design replacements for massive midface injuries with bone loss. The final solution to meet functional requirements may be shaped differently than the natural human bone but be optimized for functional needs sufficient to support full restoration using a combination of soft tissue repair and synthetic prosthetics. Topological optimization for designing facial bone tissue replacements might improve current clinical methods and provide essential enabling technology to translate generic bone tissue engineering methods into specific solutions for individual patients.


C. Talischi, G. H. Paulino, A. Pereira, and I. F. M. Menezes. “Polygonal finite elements for topology optimization: A unifying paradigm.” International Journal for Numerical Methods in Engineering. Vol. 82, No. 6, pp. 671698. 2010. Link to PDF View/Hide Abstract
In topology optimization literature, the parameterization of design is commonly carried out on uniform grids consisting of Lagrangiantype finite elements (e.g. linear quads). These formulations, however, suffer
from numerical anomalies such as checkerboard patterns and onenode connections, which has prompted
extensive research on these topics. A problem less often noted is that the constrained geometry of these
discretizations can cause bias in the orientation of members, leading to meshdependent suboptimal
designs. Thus, to address the geometric features of the spatial discretization, we examine the use of
unstructured meshes in reducing the influence of mesh geometry on topology optimization solutions.
More specifically, we consider polygonal meshes constructed from Voronoi tessellations, which in addition
to possessing higher degree of geometric isotropy, allow for greater flexibility in discretizing complex
domains without suffering from numerical instabilities.


E. E. Theotokoglou, I. H. Stampouloglou, and G. H. Paulino. “An analytical approach for an adhesive layer in a graded elastic wedge.” Mechanics of Advanced Materials and Structures. Vol. 17, No. 6, pp. 393405. 2010. Link to PDF View/Hide Abstract
In this paper the influence of an adhesive layer in a graded elastic wedge consisted of two subwedges radially bonded, is investigated by means of linear elasticity. The adhesive layer in the analytical solution is simulated either by an interface or by an infinitesimal subwedge of very small wedgeangle. The graded character of the wedges is given either by a linearly varying or by an exponentially varying shear modulus. The inhomogeneous anisotropic selfsimilar biwedge and triwedge, loaded by a concentrated force at their apex, are studied analytically under plane strain or generalized plane stress conditions, using the selfsimilarity property. Based on the separation of loading in each subwedge and on the continuity of displacements at the interface, an analytic solution is deduced for the stress and displacements fields. Applications have been made in the case of a graded biwedge and a composite triwedge, in which for specific values of gradation, the stress and displacements fields are determined.


Back to Top 
2009  
M. Alfano, F. Furgiuele, A. Leonardi, C. Maletta, and G. H. Paulino. “Mode I fracture of adhesive joints using tailored cohesive zone models.” International Journal of Fracture. Vol. 157, No. 12, pp. 193204. 2009. Link to PDF View/Hide Abstract
Cohesive zone models are explored in
order to study cleavage fracture in adhesive bonded
joints. A mode I cohesive model is defined which correlates
the tensile traction and the displacement jump
(crack faces opening) along the fracture process zone.
In order to determine the tractionseparation relation,
the main fracture parameters, namely the cohesive
strength and the fracture energy, as well as its shape,
must be specified. However, owing to the difficulties
associated to the direct measurement of the fracture
parameters, very often they are obtained by comparing
a measured fracture property with numerical predictions
based on an idealized traction separation relation.
Likewise in this paper the cohesive strength of an
adhesive layer sandwiched between elastic substrates
is adjusted to achieve a match between simulations
and experiments. To this aim, the fracture energy and
the loaddisplacement curve are adopted as input in
the simulations. In order to assess whether or not the
shape of the cohesive model may have an influence on
the results, three representative cohesive zone models
have been investigated, i.e. exponential, bilinear and
trapezoidal. A good agreement between experiments
and simulations has been generally observed. However,
a slight difference in predicting the loads for damage
onset has been found using the different cohesive
models.


S. R. M. Almeida, G. H. Paulino, and E. C. N. Silva. "A simple and effective inverse projection scheme for void distribution control in topology optimization." Structural and Multidisciplinary Optimization. Vol. 39, No. 4, pp. 359371. 2009. Link to PDF View/Hide Abstract
The ability to control both the minimum size
of holes and the minimum size of structural members
are essential requirements in the topology optimization
design process for manufacturing. This paper addresses
both requirements by means of a unified approach
involving meshindependent projection techniques. An
inverse projection is developed to control the minimum
hole size while a standard direct projection scheme
is used to control the minimum length of structural
members. In addition, a heuristic scheme combining
both contrasting requirements simultaneously is discussed.
Two topology optimization implementations
are contributed: one in which the projection (either
inverse or direct) is used at each iteration; and the
other in which a twophase scheme is explored. In
the first phase, the compliance minimization is carried
out without any projection until convergence. In the
second phase, the chosen projection scheme is applied
iteratively until a solution is obtained while satisfying
either the minimum member size or minimum hole
size. Examples demonstrate the various features of the
projectionbased techniques presented.


R. C. Carbonari, E. C. N. Silva, and G. H. Paulino. “Multiactuated functionally graded piezoelectric microtools design: A multiphysics topology optimization approach.” International Journal for Numerical Methods in Engineering. Vol. 77, No. 2, pp. 301336. 2009.
Link to PDF View/Hide Abstract
Microtools offer significant promise in a wide range of applications such as cell manipulation, microsurgery,
and micro/nanotechnology processes. Such special microtools consist of multiflexible structures
actuated by two or more piezoceramic devices that must generate output displacements and forces at
different specified points of the domain and at different directions. The microtool structure acts as a
mechanical transformer by amplifying and changing the direction of the piezoceramics output displacements.
The design of these microtools involves minimization of the coupling among movements generated
by various piezoceramics. To obtain enhanced microtool performance, the concept of multifunctional and
functionally graded materials is extended by tailoring elastic and piezoelectric properties of the piezoceramics
while simultaneously optimizing the multiflexible structural configuration using multiphysics
topology optimization. The design process considers the influence of piezoceramic property gradation
and also its polarization sign. The method is implemented considering continuum material distribution
with special interpolation of fictitious densities in the design domain. As examples, designs of a single
piezoactuator, an XY nanopositioner actuated by two graded piezoceramics, and a microgripper actuated
by three graded piezoceramics are considered. The results show that material gradation plays an important
role to improve actuator performance, which may also lead to optimal displacements and coupling ratios
with reduced amount of piezoelectric material. The present examples are limited to twodimensional
models because many of the applications for such microtools are planar devices.


R. Espinha, W. Celes, N. Rodriguez, and G. H. Paulino. “ParTopS: compact topological framework for parallel fragmentation simulations.” Engineering
with Computers. Vol. 25, No. 4, pp. 345365. 2009.
Link to PDF View/Hide Abstract
Cohesive models are used for simulation of
fracture, branching and fragmentation phenomena at various
scales. Those models require high levels of mesh
refinement at the crack tip region so that nonlinear behavior
can be captured and physical results obtained. This imposes
the use of large meshes that usually result in computational
and memory costs prohibitively expensive for a single
traditional workstation. If an extrinsic cohesive model is to
be used, support for dynamic insertion of cohesive elements
is also required. This paper proposes a topological
framework for supporting parallel adaptive fragmentation
simulations that provides operations for dynamic insertion
of cohesive elements, in a uniform way, for both two and
threedimensional unstructured meshes. Cohesive elements
are truly represented and are treated like any other regular
element. The framework is built as an extension of a
compact adjacencybased serial topological data structure,
which can natively handle the representation of cohesive
elements. Symmetrical modifications of duplicated entities
are used to reduce the communication of topological
changes among mesh partitions and also to avoid the use of
locks. The correctness and efficiency of the proposed
framework are demonstrated by a series of arbitrary
insertions of cohesive elements into some sample meshes.


F. Evangelista, J. Roesler, and G. H. Paulino. “Numerical simulations of fracture resistance of functionally graded concrete materials.” Transportation Research Record: Journal of the Transportation Research Board No. 2113,
pp. 122131. 2009.
Link to PDF View/Hide Abstract
The concept of grading material composition in a predetermined direction to target multiple objectives and functionality is applicable to the layering and positioning of different materials at specified depths. From a fracture mechanics perspective, this study explores the advantages of using functionally graded concrete materials (FGCMs), that is, plain concrete and fiberreinforced concrete (FRC), in two distinct layers. The fracture energy (G) and residual load capacity (Pdelta) of twolayered concrete beams are investigated by means of numerical simulations with a cohesive zone model (CZM) implemented in a finite element framework. The required fracture parameters for defining the CZM are obtained from individual fracture tests of the plain concrete and FRC materials. The numerical simulations analyzed the effects of FRC thickness and position (whether at the top or bottom of the beam) on the fracture resistance of the twolayered concrete beam. A costbenefit analysis showed that the FRC placed in the bottom lift is more fracture efficient (higher G and Pdeltavalues at lower cost) than when it is placed in the top lift. There is also an optimal FRC thickness in which the benefit in fracture resistance is reduced as the FRC layer is increased. The application of a CZM to predict the fracture behavior of an FGCM beam has demonstrated its potential for also quantifying the effects of FGCMs on the fracture resistance of concrete pavements.


K. Park, J. Pereira, C. A. Duarte, and G. H. Paulino. “Integration of singular enrichment functions in the generalized/extended finite element method for threedimensional problems.” International Journal for Numerical Methods in Engineering. Vol. 78, No. 10, pp. 12201257. 2009.
Link to PDF View/Hide Abstract
A mapping method is developed to integrate weak singularities, which result from enrichment functions
in the generalized/extended finite element method. The integration scheme is applicable to 2D and 3D
problems including arbitrarily shaped triangles and tetrahedra. Implementation of the proposed scheme in
existing codes is straightforward. Numerical examples for 2D and 3D problems demonstrate the accuracy
and convergence properties of the technique.


K. Park, G. H. Paulino, and J. R. Roesler. "A unified potentialbased cohesive model for mixedmode fracture."Journal of the Mechanics and Physics of Solids. Vol. 57, No. 6, pp. 891908. 2009.
Link to PDF View/Hide Abstract
A generalized potentialbased constitutive model for mixedmode cohesive fracture is presented in conjunction with physical parameters such as fracture energy, cohesive strength and shape of cohesive interactions. It characterizes different fracture energies in each fracture mode, and can be applied to various material failure behavior (e.g. quasibrittle). The unified potential leads to both intrinsic (with initial slope indicators to control elastic behavior) and extrinsic cohesive zone models. Path dependence of workofseparation is investigated with respect to proportional and nonproportional paths—this investigation demonstrates consistency of the cohesive constitutive model. The potentialbased model is verified by simulating a mixedmode bending test. The actual potential is named PPR (Park–Paulino–Roesler), after the first initials of the authors’ last names.


G. H. Paulino and C. H. Le. "A modified Q4/Q4 element for topology optimization." Structural and Multidisciplinary Optimization. Vol. 37, No. 3, pp. 255264. 2009.
Link to PDF View/Hide Abstract
Design variable and displacement fields are
distinct. Thus, their respective material and finite element
meshes may also be distinct, as well as their actual
location of nodes for the two fields. The proposed
Q4/Q4M element possesses design variable nodes and
displacement nodes which are not coincident. The
element has been implemented using different approaches,
including continuous approximation of material
distribution and nodal approaches. The results
obtained demonstrate that the element is effective in
generating structural topologies with high resolution.
From a numerical point of view, mesh independent
solutions can be obtained by means of projection.


G. H. Paulino, E. C. N. Silva, and C. H. Le. “Optimal design of periodic functionally graded composites with prescribed properties.” Structural and Multidisciplinary Optimization. Vol. 38, No. 5, pp. 469489. 2009.
Link to PDF View/Hide Abstract
The computational design of a composite
where the properties of its constituents change gradually
within a unit cell can be successfully achieved
by means of a material design method that combines
topology optimization with homogenization. This is an
iterative numerical method, which leads to changes in
the composite material unit cell until desired properties
(or performance) are obtained. Such method has
been applied to several types of materials in the last
few years. In this work, the objective is to extend the
material design method to obtain functionally graded
material architectures, i.e. materials that are graded at
the local level (e.g. microstructural level). Consistent
with this goal, a continuum distribution of the design
variable inside the finite element domain is considered
to represent a fully continuous material variation during
the design process. Thus the topology optimization
naturally leads to a smoothly graded material system.
To illustrate the theoretical and numerical approaches,
numerical examples are provided. The homogenization
method is verified by considering onedimensional material
gradation profiles for which analytical solutions
for the effective elastic properties are available. The
verification of the homogenization method is extended
to two dimensions considering a trigonometric material
gradation, and a material variation with discontinuous
derivatives. These are also used as benchmark examples
to verify the optimization method for functionally
graded material cell design. Finally the influence of
material gradation on extreme materials is investigated,
which includes materials with nearzero shear modulus,
and materials with negative Poisson’s ratio.


W. M. Rubio, E. C. N. Silva, and G. H. Paulino. "Toward optimal design of piezoelectric transducers based on multifunctional and smoothly graded hybrid material systems." Journal of Intelligent Material Systems and Structures. Vol. 20, No. 14, pp. 17251746. 2009.
Link to PDF View/Hide Abstract
This work explores the design of piezoelectric transducers based on functional material gradation, here named functionally graded piezoelectric transducer (FGPT). Depending on the applications, FGPTs must achieve several goals, which are essentially related to the transducer resonance frequency, vibration modes, and excitation strength at specific resonance frequencies. Several approaches can be used to achieve these goals; however, this work focuses on finding the optimal material gradation of FGPTs by means of topology optimization. Three objective functions are proposed: (i) to obtain the FGPT optimal material gradation for maximizing specified resonance frequencies; (ii) to design piezoelectric resonators, thus, the optimal material gradation is found for achieving desirable eigenvalues and eigenmodes; and (iii) to find the optimal material distribution of FGPTs, which maximizes specified excitation strength. To track the desirable vibration mode, a modetracking method utilizing the ‘modal assurance criterion’ is applied. The continuous change of piezoelectric, dielectric, and elastic properties is achieved by using the graded finite element concept. The optimization algorithm is constructed based on sequential linear programming, and the concept of continuum approximation of material distribution. To illustrate the method, 2D FGPTs are designed for each objective function. In addition, the FGPT performance is compared with the nonFGPT one.


C. Talischi, G. H. Paulino, and C. Le. "Honeycomb Wachspress finite elements for structural topology optimization."Structural and Multidisciplinary Optimization. Vol. 37, No. 6, pp. 569583. 2009.
Link to PDF View/Hide Abstract
Traditionally, standard Lagrangiantype finite
elements, such as linear quads and triangles, have
been the elements of choice in the field of topology
optimization. However, finite element meshes with
these conventional elements exhibit the wellknown
“checkerboard” pathology in the iterative solution of
topology optimization problems. A feasible alternative
to eliminate such longstanding problem consists of using
hexagonal (honeycomb) elements with Wachspresstype
shape functions. The features of the hexagonal
mesh include twonode connections (i.e. two elements
are either not connected or connected by two nodes),
and three edgebased symmetry lines per element. In
contrast, quads can display onenode connections,
which can lead to checkerboard; and only have two
edgebased symmetry lines. In addition, Wachspress
rational shape functions satisfy the partition of unity
condition and lead to conforming finite element approximations.
We explore the Wachspresstype hexagonal
elements and present their implementation using
three approaches for topology optimization: elementbased,
continuous approximation of material distribution,
and minimum lengthscale through projection
functions. Examples are presented that demonstrate
the advantages of the proposed element in achieving
checkerboardfree solutions and avoiding spurious finescale
patterns from the design optimization process.


Back to Top 
2008  
Y. S. Chan, G. H. Paulino, and A. C. Fannjiang. "Gradient elasticity theory for mode III fracture in functionally graded materials  Part II: Crack parallel to the material gradation." ASME Journal of Applied Mechanics. Vol. 76, No. 6, 0610151 to 11. 2008. Link to PDF View/Hide Abstract
A ModeIII crack problem in a functionally graded material modeled by anisotropic straingradient elasticity theory is solved by the integral equation method. The gradient elasticity theory has two material characteristic lengths and, which are responsible for volumetric and surface straingradient terms, respectively. The governing differential equation of the problem is derived assuming that the shear modulus G is a function of x, i.e., G = G(x) = G(0)e(beta x), where G(0) and beta are material constants. A hypersingular integrodifferential equation is derived and discretized by means of the collocation method and a Chebyshev polynomial expansion. Numerical results are given in terms of the crack opening displacements, strains, and stresses with various combinations of the parameters l, l', and beta. Formulas for the stress intensity factors, KIII, are derived and numerical results are provided.


Y. Chen, L. J. Struble, and G. H. Paulino. "Using rheology to achieve coextrusion of cementbased materials with graded cellular structures." International Journal of Applied Ceramic Technology. Vol. 5, No. 5, pp. 513521. 2008. Link to PDF View/Hide Abstract
Coextrusion involves simultaneous extrusion of multiple layers and can be used to produce functionally graded materials
whose layers have different properties. Rheological control is vital for successful coextrusion. During extrusion, flow in the
barrel and die land in a ram extruder should be pluglike, while the paste should be sheared and uniformly elongated in the die
entry region. In the barrel of the extruder, the paste flow velocity field was inferred by direct observation of the paste left in the
barrel, and evidence for plug flow in the barrel was seen only at lowextrudate velocities. In the die land, the Benbow nonlinear
model was employed to assess the paste flow behavior, and plug flow was achieved only when the shear stress applied to the
paste by the die land wall was smaller than its yield stress. For coextrusion, a simple method using thinwalled tubes was found
to be effective to prepare layered feedrods. Functionally graded cellular structures of cementbased materials were successfully
coextruded by using a lowextrudate velocity when the paste had decreasing shear viscosity from inner to outer layers.


H. J. Choi and G. H. Paulino. “Thermoelastic contact mechanics for a flat punch sliding over a graded coating/substrate system with frictional heat generation.” Journal of the Mechanics and Physics of Solids. Vol. 56, No. 4, pp. 16731692. 2008. Link to PDF View/Hide Abstract
The problem of thermoelastic contact mechanics for the coating/substrate system with functionally graded properties is
investigated, where the rigid flat punch is assumed to slide over the surface of the coating involving frictional heat
generation. With the coefficient of friction being constant, the inertia effects are neglected and the solution is obtained
within the framework of steadystate plane thermoelasticity. The graded material exists as a nonhomogeneous interlayer
between dissimilar, homogeneous phases of the coating/substrate system or as a nonhomogeneous coating deposited on the
substrate. The material nonhomogeneity is represented by spatially varying thermoelastic moduli expressed in terms of
exponential functions. The Fourier integral transform method is employed and the formulation of the current
thermoelastic contact problem is reduced to a Cauchytype singular integral equation of the second kind for the unknown
contact pressure. Numerical results include the distributions of the contact pressure and the inplane component of the
surface stress under the prescribed thermoelastic environment for various combinations of geometric, loading, and
material parameters of the coated medium. Moreover, in order to quantify and characterize the singular behavior of
contact pressure distributions at the edges of the flat punch, the stress intensity factors are defined and evaluated in terms
of the solution to the governing integral equation.


K. Park, G. H. Paulino and J. R. Roesler. "Determination of the kink point in the bilinear softening model for concrete." Engineering Fracture Mechanics. Vol. 75, No. 13, pp. 38063818. 2008. Link to PDF View/Hide Abstract
The characterization of the softening curve from experimental results is essential for predicting the fracture behavior of
quasibrittle materials like concrete. Among various shapes (e.g. linear, exponential) to describe the softening behavior of
concrete, the bilinear softening relationship has been extensively used and is the model of choice in this work. Currently,
there is no consensus about the location of the kink point in the bilinear softening curve. In this study, the location of the
kink point is proposed to be the stress at the critical crack tip opening displacement. Experimentally, the fracture parameters
required to describe the bilinear softening curve can be determined with the ‘‘twoparameter fracture model” and the
total work of fracture method based on a single concrete fracture test. The proposed location of the kink point compares
well with the range of kink point locations reported in the literature, and is verified by plotting stress profiles along the
expected fracture line obtained from numerical simulations with the cohesive zone model. Finally, prediction of experimental
load versus crack mouth opening displacement curves validate the proposed location of the kink point for different concrete
mixtures and also for geometrically similar specimens with the same concrete mixture. The experiments were
performed on threepoint bending specimens with concrete mixtures containing virgin coarse aggregate, recycled concrete
coarse aggregate (RCA), and a 50–50 blend of RCA and virgin coarse aggregate. The verification and validation studies
support the hypothesis of the kink point occurring at the critical crack tip opening displacement.


K. Park, G. H. Paulino and J. R. Roesler. "Virtual internal pairbond model for quasibrittle materials." ASCE  Journal of Engineering Mechanics. Vol. 134, No. 10, pp. 856866. 2008. Link to PDF View/Hide Abstract
The present multiscale investigation employs the initial and total fracture energy through a virtual internal pairbond (VIPB). model The proposed VIPB model is an extension of the traditional virtual internal bond VIB model. Two different types of potentials,
a steep shortrange potential and a shallow longrange potential, are employed to describe the initial and the total fracture energies,
respectively. The Morse potential function is modified for the virtual bond potential so that it is independent of specific length scales
associated with the lattice geometry. This feature is incorporated in the VIPB model, which uses both fracture energies and cohesive
strength. With respect to the discretization by finite elements, we address the element size dependence in conjunction with the J integral.
Parameters in the VIPB model are evaluated by numerical simulations of a pure tension test in conjunction with measured fracture
parameters. We also validate the VIPB model by predicting load versus crack mouth opening displacement curves for geometrically
similar specimens, and the measured size effect. Finally, we provide an example involving fiberreinforced concrete, which demonstrates
the advantage of the VIPB model over the usual VIB model.


G. H. Paulino, W. Celes, R. Espinha, and Z. (J.) Zhang. “A general topologybased framework for adaptive insertion of cohesive elements in finite element meshes.” Engineering with Computers. Vol. 24, No. 1, pp. 5978. 2008. Link to PDF View/Hide Abstract
Largescale simulation of separation phenomena
in solids such as fracture, branching, and fragmentation
requires a scalable data structure representation of the
evolving model. Modeling of such phenomena can be
successfully accomplished by means of cohesive models of
fracture, which are versatile and effective tools for computational
analysis. A common approach to insert cohesive
elements in finite element meshes consists of adding discrete
special interfaces (cohesive elements) between bulk
elements. The insertion of cohesive elements along bulk
element interfaces for fragmentation simulation imposes
changes in the topology of the mesh. This paper presents a
unified topologybased framework for supporting adaptive
fragmentation simulations, being able to handle two and
threedimensional models, with finite elements of any order.
We represent the finite element model using a compact
and ‘‘complete’’ topological data structure, which is capable
of retrieving all adjacency relationships needed for the
simulation. Moreover, we introduce a new topologybased
algorithm that systematically classifies fractured facets
(i.e., facets along which fracture has occurred). The algorithm
follows a set of procedures that consistently perform
all the topological changes needed to update the model.
The proposed topologybased framework is general and
ensures that the model representation remains always valid
during fragmentation, even when very complex crack
patterns are involved. The framework correctness and
efficiency are illustrated by arbitrary insertion of cohesive
elements in various finite element meshes of selfsimilar
geometries, including both two and threedimensional
models. These computational tests clearly show linear
scaling in time, which is a key feature of the present datastructure
representation. The effectiveness of the proposed
approach is also demonstrated by dynamic fracture analysis
through finite element simulations of actual engineering
problems.


M. J. Pindera and G. H. Paulino. “Honoring professor Erdogan's seminal contributions to mixed boundaryvalue problems of inhomogeneous and functionally graded materials  Foreword.” ASME Journal of Applied Mechanics. Vol. 75, No. 5, 0503011. 2008. Link to PDF Link to News Article View/Hide Abstract
Professor Fazil Erdogan has influenced several generations of
applied and solid mechanicians working in the area of mixed
boundaryvalue problems of inhomogeneous media, most notably
fracture and contact problems. The analytical approaches that he
had developed with his students in the 1960s and 1970s for the
formulation and reduction of fracture mechanics problems involving
layered media to systems of singular integral equations, and
the corresponding solution techniques, have motivated researchers
working in this area throughout the entire world. His subsequent
work on fracture mechanics problems of inhomogeneous media
with smoothly varying elastic moduli has laid the foundation for
applying these techniques to functionally graded materials, which
played key roles in many technologically important applications e.g., spatially tailored structures for the new generation of hypersonic
aircraft, graded cementitious composites for sustainable infrastructure,
highperformance graded components for automobiles,
and graded microtools in mechatronics. Professor
Erdogan’s continuing leadership role and ceaseless contributions
to the fracture and contact mechanics of this new generation of
materials provide guidance and motivation for others to follow.
This special issue honors Professor Erdogan in recognition of his
past and continuing contributions in the area that plays a critical
role in the development of engineered material systems for critical
technological applications, and builds upon a minisymposium under
the above title held at the recent International Conference on
Multiscale and Functionally Graded Materials (MFGM2006) on
Oct. 15–18, 2006, Honolulu, HI.


B. Shen, M. Hubler, G. H. Paulino, and L. J. Struble. "Functionally graded fiberreinforced cement composite: Processing, microstructure, and properties." Cement and Concrete Composites. Vol. 30, No. 8, pp. 663673. 2008. Link to PDF View/Hide Abstract
A functionallygraded material system has been developed with a spatially tailored fiber distribution to produce a fourlayered, functionally
graded fiberreinforced cement composite (FGFRCC). Fiber volume fractions were increased linearly from 0% in the compression
zone to 2% in the tensile zone so that more fibers were available to carry tensile stress in a bending beam experiment. Extrusion was
used to produce single homogeneous layers of constant fiber volume fraction and highly oriented fibers. The FRCC layers with different
fiber volume fractions were stacked according to the desired configuration and then pressed to produce an integrated FGFRCC. Polished
crosssections of the FGFRCC were examined using a scanning electron microscope (SEM) to measure the fiber distribution. Flexural
tests were carried out to characterize the mechanical behavior and to evaluate the effectiveness of the designed fiber distribution. No
delamination between layers was observed in the fractured specimens. Compared to homogeneous FRCC with the same overall fiber
volume fraction, the FGFRCC exhibited about 50% higher strength and comparable work of fracture.


S. H. Song, M. P. Wagoner, G. H. Paulino, and W. G. Buttlar. “delta(25) Crack opening displacement parameter in cohesive zone models: experiments and simulations in asphalt concrete.” Fatigue & Fracture of Engineering Materials and Structures. Vol. 31, No. 10, pp. 850856. 2008. Link to PDF View/Hide Abstract
Recent work with fracture characterization of asphalt concrete has shown that a cohesive
zone model (CZM) provides insight into the fracture process of the materials. However,
a current approach to estimate fracture energy, i.e., in terms of area of force versus crack
mouth opening displacement (CMOD), for asphalt concrete overpredicts its magnitude.
Therefore, the δ25 parameter, which was inspired by the δ5 concept of Schwalbe and
coworkers, is proposed as an operational definition of a crack tip opening displacement
(CTOD). The δ25 measurement is incorporated into an experimental study of validation
of its usefulness with asphalt concrete, and is utilized to estimate fracture energy. The
work presented herein validates the δ25 parameter for asphalt concrete, describes the
experimental techniques for utilizing the δ25 parameter, and presents threedimensional
(3D) CZM simulations with a specially tailored cohesive relation. The integration of the
δ25 parameter and new cohesive model has provided further insight into the fracture
process of asphalt concrete with relatively good agreement between experimental results
and numerical simulations.


H. M. Yin, W. G. Buttlar, G. H. Paulino, and H. Di Benedetto. "Assessment of existing micromechanical models for asphalt mastics considering viscoelastic effects." Road Materials and Pavement Design. Vol. 9, No. 1, pp. 3157. 2008. Link to PDF View/Hide Abstract
Micromechanical models have been directly used to predict the effective complex modulus of asphalt mastics from the mechanical properties of their constituents. Because the micromechanics models traditionally employed have been based on elastic theory, the viscoelastic effects of binders have not been considered. Moreover, due to the unique features of asphalt mastics such as high concentration and irregular shape of filler particles, some micromechanical models may not be suitable. A comprehensive investigation of four existing micromechanical methods is conducted considering viscoelastic effects. It is observed that the selfconsistent model well predicts the experimental results without introducing any calibration; whereas the MoriTanaka model and the generalized selfconsistent model, which have been widely used for asphalt materials, significantly underestimate the complex Young's modulus. Assuming binders to be incompressible and fillers to be rigid, the dilute model and the selfconsistent model provide the same prediction, but they considerably overestimate the complex Young's modulus. The analyses suggest that these conventional assumptions are invalid for asphalt mastics at low temperatures and high frequencies. In addition, contradictory to the assumption of the previous elastic model, it is found that the phase angle of binders produces considerable effects on the absolute value of the complex modulus of mastics.


H. M. Yin, G. H. Paulino, and W. G. Buttlar. "An explicit elastic solution for a brittle film with periodic cracks." International Journal of Fracture. Vol. 153, No. 1, pp. 3952. 2008. Link to PDF View/Hide Abstract
A twodimensional explicit elastic solution is derived for a brittle film bonded to a ductile substrate through either a frictional interface or a fully bonded interface, in which periodically distributed discontinuities are formed within the film due to the applied tensile stress in the substrate and consideration of a "weak form stress boundary condition" at the crack surface. This solution is applied to calculate the energy release rate of threedimensional channeling cracks. Fracture toughness and nominal tensile strength of the film are obtained through the relation between crack spacing and tensile strain in the substrate. Comparisons of this solution with finite element simulations show that the proposed model provides an accurate solution for the film/substrate system with a frictional interface; whereas for a fully bonded interface it produces a good prediction only when the substrate is not overly compliant or when the crack spacing is large compared with the thickness of the film. If the section is idealized as infinitely long, this solution in terms of the energy release rate recovers Beuth's exact solution for a fully cracked film bonded to a semiinfinite substrate. Interfacial shear stress and the edge effect on the energy release rate of an asymmetric crack are analyzed. Fracture toughness and crack spacing are calculated and are in good agreement with available experiments.


H. M. Yin, G. H. Paulino, W. G. Buttlar, and L. Z. Sun. “Effective thermal conductivity of functionally graded particulate nanocomposites with interfacial thermal resistance.” ASME Journal of Applied Mechanics. Vol. 75, No. 5, 0511131 to 6. 2008. Link to PDF View/Hide Abstract
By means of a fundamental solution for a single inhomogeneity embedded in a functionally
graded material matrix, a selfconsistent model is proposed to investigate the effective
thermal conductivity distribution in a functionally graded particulate nanocomposite.
The “Kapitza thermal resistance” along the interface between a particle and the matrix
is simulated with a perfect interface but a lower thermal conductivity of the particle. The
results indicate that the effective thermal conductivity distribution greatly depends on
Kapitza thermal resistance, particle size, and degree of material gradient.


H. M. Yin, G. H. Paulino, W. G. Buttlar, and L. Z. Sun. "Heat flux field for one spherical inhomogeneity embedded in a functionally graded material matrix." International Journal of Heat and Mass Transfer. Vol. 51, Nos. 1112, pp. 30183024. 2008. Link to PDF View/Hide Abstract
The heat flux field for a single particle embedded in a graded material is derived by using the equivalent inclusion method. A linearly
distributed prescribed heat flux field is introduced to represent the material mismatch between the particle and the surrounding graded
materials. By using Green’s function technique, an explicit solution is obtained for the heat flux field in both the particle and the graded
material. Comparison of the present solution with finite element results illustrates the accuracy and limitation of this solution.


Back to Top 
2007  
E. Alfano, F. Furgiuele, A. Leonardi, C. Maletta, and G. H. Paulino. "Cohesive zone modeling of mode I fracture in adhesive bonded joints." Key Engineering Materials. 348349, pp. 1316. 2007. View/Hide Abstract
This paper deals with the application of Cohesive Zone Model (CZM) concepts to study mode I fracture in adhesive bonded joints. In particular, an intrinsic piecewise linear cohesive surface relation is used in order to model fracture in a precracked bonded Double Cantilever Beam (DCB) specimen, Finite element implementation of the CZM is accomplished by means of the user element (UEL) feature available in the FE commercial code ABAQUS. The sensitivity of the cohesive zone parameters (i.e. fracture strength and critical energy release rate) in predicting the overall mechanical response is first examined; subsequently, cohesive parameters are tuned comparing numerical simulations of the loaddisplacement curve with experimental results retrieved from literature.


R. C. Carbonari, E. C. N. Silva, and G. H. Paulino. “Topology optimization design of functionally graded bimorphtype piezoelectric actuators.” Smart Materials & Structures. Vol. 16, No. 6, pp. 26052620. 2007. Link to PDF View/Hide Abstract
The concept of a functionally graded material (FGM) is useful for
engineering advanced piezoelectric actuators. For instance, it can lead
to locally improved properties, and to increased lifetime of bimorph
piezoelectric actuators. By selectively grading the elastic,
piezoelectric, and/or dielectric properties along the thickness of a
piezoceramic, the resulting gradation of electromechanical properties
influences the behavior and performance of piezoactuators. In this
work, topology optimization is applied to find the optimum gradation
and polarization sign variation in piezoceramic domains in order to
improve actuator performance measured in terms of selected output
displacements. A bimorphtype actuator is emphasized, which is designed
by maximizing the output displacement or output force at selected
location(s) (e.g. the tip of the actuator). The numerical
discretization is based on the graded finite element concept such that
a continuum approximation of material distribution, which is
appropriate to model FGMs, is achieved. The present results consider
twodimensional models with a planestrain assumption. The material
gradation plays an important role in improving the actuator performance
when distributing piezoelectric (PZT5A) and nonpiezoelectric (gold)
materials in the design domain; however, the performance is not
improved when distributing two types of similar piezoelectric material.
In both cases, the polarization sign change did not play a significant
role in the results. However, the optimizer always finds a solution
with opposite polarization (as expected).


J. H. Kim and G. H. Paulino. “On fracture criteria for mixedmode crack propagation in functionally graded materials.” Mechanics of Advanced Materials and Structures. Vol. 14, No. 4, pp. 227244. 2007. Link to PDF View/Hide Abstract
Recently, the functionally graded material (FGM) concept has been
explored in piezoelectric materials to improve properties and to
increase the lifetime of bimorph piezoelectric actuators. For instance,
elastic, piezoelectric, and/or dielectric properties may be graded
along the thickness of a piezoceramic. Thus, the gradation of
piezoceramic properties influences the performance of piezoactuators.
The usual FGM modelling using traditional finite element formulation
and discrefization into layers gives a highly discontinuous stress
distribution, which is undesirable. In this work, we focus on
nonhomogeneous piezoelectric materials using a generalized
isoparametric formulation based on the graded finite element concept,
in which the properties change smoothly inside the element. This
approach provides a continuum material distribution, which is
appropriate to model FGMs. Both fournode quadrilaterals and eightnode
quadrilaterals for piezoelectric FGMs were implemented using the graded
finite element concept. A closed form twodimensional analytical model
of piezoelectric FGMs is also developed to check the accuracy of these
finite elements and to assess the influence of material property
gradation on the behavior of piezoelectric FGMs. The paper discusses
and compares the behavior of piezoelectric graded elements under four
loading conditions with respect to the analytical solutions (derived in
this work) considering exponential variation of elastic, piezoelectric,
and dielectric properties separately. The analytical solutions provide
benchmark problems to verify numerical procedures (such as the finite
element method and the boundary element method).


J. R. Roesler, G. H. Paulino, C. Gaedicke, A. Bordelon, and K. Park. "Fracture behavior of functionally graded concrete materials for rigid pavements." Transportation Research Record 3047, pp. 4049. 2007. Link to PDF View/Hide Abstract
Currently, in concrete pavements, a single concrete mixture design and structural surface layer are selected to resist mechanical loading without an attempt to affect concrete pavement shrinkage, ride quality, or noise attenuation adversely. An alternative approach is to design sublayers within the concrete pavement surface that have specific functions and thus to achieve higher performance at a lower cost. The objective of this research was to address the structural benefits of functionally graded concrete materials (FGCMs) for rigid pavements by testing and modeling the fracture behavior of different combinations of layered plain concrete materials and concrete materials reinforced with synthetic fibers. The threepoint bendingbeam test was used to obtain the softening behavior and fracture parameters of each FGCM. The peak loads and initial fracture energy between the plain, fiberreinforced, and FGCMs were similar; this signified similar crack initiation. The total fracture energy clearly indicated the improvements in fracture behavior of FGCM relative to fulldepth plain concrete. The fracture behavior of FGCM depended on the position of the fiberreinforced layer relative to the starter notch. The fracture parameters of both the fiberreinforced and plain concrete were embedded into a finite elementbased cohesive zone model. The model successfully captured the experimental behavior of the FGCMs and now can be implemented to predict the fracture behavior of proposed FGCM configurations and structures such as rigid pavements. This integrated approach (testing and modeling) is promising and demonstrates the viability of FGCM for designing layered concrete pavement systems.


G. H. Paulino and J. H. Kim. "The weak patch test for nonhomogeneous materials modeled with graded finite elements." Journal of the Brazilian Society of Mechanical Sciences and Engineering. Vol. XXIX, No. 1, pp. 6381. 2007. Link to PDF View/Hide Abstract
Functionally graded materials have an additional length scale
associated to the spatial variation of the material property field
which competes with the usual geometrical length scale of the boundary
value problem. By considering the length scale of nonhomogeneity, this
paper presents the weak patch test (rather than the standard one) of
the graded element for nonhomogeneous materials to assess convergence
of the finite element method (FEM). Both consistency (as the size of
elements approach zero, the FEM approximation represents the exact
solution) and stability (spurious mechanisms are avoided) conditions
are addressed. The specific graded elements considered here are
isoparametric quadrilaterals (e.g. 4, 8 and 9node) considering two
dimensional plane and axisymmetric problems. The finite element
approximate solutions are compared with exact solutions for
nonhomogeneous materials.


J. Roesler, G. H. Paulino, K. Park and C. Gaedicke. "Concrete fracture prediction using bilinear softening." Cement & Concrete Composites. Vol. 29, No. 4, pp. 300312. 2007. Link to PDF View/Hide Abstract
A finite elementbased cohesive zone model was developed using bilinear
softening to predict the monotonic load versus crack mouth opening
displacement curve of geometrically similar notched concrete specimens.
The softening parameters for concrete material are based on concrete
fracture tests, total fracture energy (G(F)), initial fracture energy
(G(f)), and tensile strength (f(t)'), which are obtained from a
threepoint bending configuration. The features of the finite element
model are that bulk material elements are used for the uncracked
regions of the concrete, and an intrinsicbased tractionopening
constitutive relationship for the cracked region. Size effect
estimations were made based on the material dependent properties (G(f)
and f(t)') and the size dependent property (G(F)). Experiments using
the threepoint bending configuration were completed to verify that the
model predicts the peak load and softening behavior of concrete for
multiple specimen depths. The fracture parameters, based oil the size
effect method or the twoparameter fracture model, were found to
adequately characterize the bilinear softening model.


E. C. N. Silva, R. C. Carbonari, and G. H. Paulino. “On graded elements for multiphysics applications.” Smart Materials & Structures. Vol. 16, No. 6, pp. 24082428. 2007. Link to PDF View/Hide Abstract
Recently, the functionally graded material (FGM) concept has been
explored in piezoelectric materials to improve properties and to
increase the lifetime of bimorph piezoelectric actuators. For instance,
elastic, piezoelectric, and/or dielectric properties may be graded
along the thickness of a piezoceramic. Thus, the gradation of
piezoceramic properties influences the performance of piezoactuators.
The usual FGM modelling using traditional finite element formulation
and discrefization into layers gives a highly discontinuous stress
distribution, which is undesirable. In this work, we focus on
nonhomogeneous piezoelectric materials using a generalized
isoparametric formulation based on the graded finite element concept,
in which the properties change smoothly inside the element. This
approach provides a continuum material distribution, which is
appropriate to model FGMs. Both fournode quadrilaterals and eightnode
quadrilaterals for piezoelectric FGMs were implemented using the graded
finite element concept. A closed form twodimensional analytical model
of piezoelectric FGMs is also developed to check the accuracy of these
finite elements and to assess the influence of material property
gradation on the behavior of piezoelectric FGMs. The paper discusses
and compares the behavior of piezoelectric graded elements under four
loading conditions with respect to the analytical solutions (derived in
this work) considering exponential variation of elastic, piezoelectric,
and dielectric properties separately. The analytical solutions provide
benchmark problems to verify numerical procedures (such as the finite
element method and the boundary element method).


S. H. Song, G. H. Paulino and W. G. Buttlar. Erratum: "Simulation of crack propagation in asphalt concrete using an intrinsic cohesive zone model." ASCE Journal of Engineering Mechanics. Vol. 133, No.7, p. 85312151223. 2007. Link to PDF 

F. V. Stump, E. C. N. Silva, and G. H. Paulino. "Optimization of material distribution in functionally graded structures with stress constraints." Communications in Numerical Methods in Engineering. Vol. 23, No.6, pp. 535551. 2007. Link to PDF View/Hide Abstract
This work describes a topology optimization framework to design the
material distribution of functionally graded structures considering
mechanical stress constraints. The problem of interest consists in
minimizing the volumetric density of a material phase subjected to a
global stress constraint. Due to the existence of microstructure, the
microlevel stress is considered, which is computed by means of a
mechanical concentration factor using a pnorm of the Von Mises stress
criterium (applied to the microlevel stress). Because a 01
(voidsolid) material distribution is not being sought, the singularity
phenomenon of stress constraint does not occur as long as the material
at any point of the medium does not vanish and it varies smoothly
between material 1 and material 2. To design a smoothly graded material
distribution, a material model based on a nonlinear interpolation of
the HashinStrikhman upper and lower bounds is considered. Consistently
with the framework adopted here in, the socalled 'continuous
approximation of material distribution' approach is employed, which
considers a continuous distribution of the design variable inside the
finite element. As examples, the designs of functionally graded disks
subjected to centrifugal body force are presented. The method generates
smooth material distributions, which are able to satisfy the stress
constraint.


H. Tan, Y. Huang, C. Liu, G. Ravichandran, and G. H. Paulino. “Constitutive behavior of composites with interface debonding: the extended MoriTanaka method for uniaxial tension.” International Journal of Fracture. Vol. 146, No. 3, pp. 139148. 2007. Link to PDF View/Hide Abstract
Debonding of particle/matrix interfaces can significantly affect the
macroscopic behavior of composite materials. We have used a nonlinear
cohesive law for particle/matrix interfaces to study the effect of
interface debonding on the macroscopic behavior of particlereinforced
composite materials subject to uniaxial tension. The MoriTanaka
method, which is suitable for composites with high particle volume
fraction, is extended to account for interface debonding. At a fixed
particle volume fraction, small particles lead to the hardening
behavior of the composite while large particles yield softening. The
interface sliding may contribute significantly to the macroscopic
behavior of the composite.


S. Wang, E. de Sturler and G. H. Paulino. "Largescale topology optimization using preconditioned Krylov subspace methods with recycling." International Journal for Numerical Methods in Engineering. Vol. 69, No. 12 , pp. 24412468. 2007. Link to PDF View/Hide Abstract
The computational bottleneck of topology optimization is the solution
of a large number of linear systems arising in the finite element
analysis. We propose fast iterative solvers for large threedimensional
topology optimization problems to address this problem. Since the
linear systems in the sequence of optimization steps change slowly from
one step to the next, we can significantly reduce the number of
iterations and the runtime of the linear solver by recycling selected
search spaces from previous linear systems. In addition, we introduce a
MINRES (minimum residual method) version with recycling (and a
shortterm recurrence) to make recycling more efficient for symmetric
problems. Furthermore, we discuss preconditioning to ensure fast
convergence. We show that a proper rescaling of the linear systems
reduces the huge condition numbers that typically occur in topology
optimization to roughly those arising for a problem with constant
density. We demonstrate the effectiveness of our solvers by solving a
topology optimization problem with more than a million unknowns on a
fast PC.


H. M. Yin, G. H. Paulino, W. G. Buttlar and L. Z. Sun. "Micromechanicsbased thermoelastic model for functionally graded particulate materials with particle interactions." Journal of the Mechanics and Physics of Solids. Vol. 55, No. 1, pp. 132160. 2007. Link to PDF View/Hide Abstract
Thermoelastic behavior of functionally graded particulate materials is
investigated with a micromechanical approach. Based on a special
representative volume element constructed to represent the graded
microstructure of a macroscopic material point, the relation between
the averaged strains of the particle and matrix phases is derived with
pairwise particle interactions, and a set of governing equations for
the thermoelastic behavior of functionally graded materials is
presented. The effective coefficient of thermal expansion at a material
point is solved through the overall averaged strain of two phases
induced by temperature change under the stressfree condition, and is
shown to exhibit a weak anisotropy due to the particle interactions
within the graded microstructure. When the material gradient is
eliminated, the proposed model predicts the effective coefficient of
thermal expansion for uniform composites as expected. If the particle
interactions are disregarded, the proposed model recovers the Kerner
model. The proposed semianalytical scheme is consistent and general,
and can handle any thermal loading variation. As examples, the thermal
stress distributions of graded thermal barrier coatings are solved for
two types of thermal loading: uniform temperature change and
steadystate heat conduction in the gradation direction.


H. M. Yin, W. G. Buttlar and G. H. Paulino. "Simplified solution for periodic thermal discontinuities in an asphalt overlay bonded to rigid pavements." Journal of Transportation Engineering. Vol. 133, No. 1, pp. 3946. 2007. Link to PDF View/Hide Abstract
This work investigates the elastic fields which develop in an overlay
bonded to a rigid substrate when the system is subjected to thermally
induced stress. A twodimensional solution of the displacement field is
derived for periodic discontinuities distributed in a hot mix asphalt
overlay bonded to a rigid pavement, where the length of the pavement
before cracking develops is much larger than its layer thickness. A
series form solution is obtained, requiring calibration due to the
limitation of the basis functions used. The formulation allows thermal
cracks of variable depth to be considered, and its accuracy is verified
through comparisons with numerical solutions obtained with ABAQUS.
Energy release rates are calculated from the model for topdown plane
strain cracking and threedimensional channeling. By comparing the
energy release rates with the fracture toughness of the overlay,
conditions for crack initiation and an estimation of crack depth for a
given temperature change can be estimated. Although several simplifying
assumptions are made in the current approach, it is shown to be more
general and therefore more widely applicable as compared to existing
closedform solutions. The solutions are valuable to the pavement
analyst who seeks to understand the general mechanisms of thermally
induced pavement deterioration and for the researcher wishing to
perform early stage verification of more complex pavement models.


Z. Zhang, G. H. Paulino and W. Celes. "Extrinsic cohesive zone modeling of dynamic fracture and microbranching instability in brittle materials." International Journal for Numerical Methods in Engineering. Vol. 72, No. 8, pp. 893923. 2007. Link to PDF View/Hide Abstract
Dynamic crack microbranching processes in brittle materials are
investigated by means of a computational fracture mechanics approach
using the finite element method with special interface elements and a
topological data structure representation. Experiments indicate
presence of a limiting crack speed for dynamic crack in brittle
materials as well as increasing fracture resistance with crack speed.
These phenomena are numerically investigated by means of a cohesive
zone model (CZM) to characterize the fracture process. A critical
evaluation of intrinsic versus extrinsic CZMs is briefly presented,
which highlights the necessity of adopting an extrinsic approach in the
current analysis. A novel topologybased data structure is employed to
enable fast and robust manipulation of evolving mesh information when
extrinsic cohesive elements are inserted adaptively. Compared to
intrinsic CZMs, which include an initial hardening segment in the
tractionseparation curve, extrinsic CZMs involve additional issues
both in implementing the procedure and in interpreting simulation
results. These include time discontinuity in stress history, fracture
pattern dependence on time step control, and numerical energy balance.
These issues are investigated in detail through a 'quasisteadystate'
crack propagation problem in poly m ethyl methacrylate. The simulation
results compare reasonably well with experimental observations both
globally and locally, and demonstrate certain advantageous features of
the extrinsic CZM with respect to the intrinsic CZM.


Z. Zhang and G. H. Paulino. “Wave propagation and dynamic analysis of smoothly graded heterogeneous continua using graded finite elements.” International Journal of Solids and Structures. Vol. 44, No. 1112, pp. 36013626. 2007. Link to PDF View/Hide Abstract
The dynamic behavior of smoothly graded heterogeneous materials is
investigated using the finite element method. The global variation of
material properties (e.g., Young's modulus, Poisson's ratio and mass
density) is treated at the element level using a generalized
isoparametric formulation. Three classes of examples are presented to
illustrate this approach and to investigate the influence of material
inhomogeneity on the characteristics of wave propagation pattern and
stress redistribution. First, a cantilever beam example is presented
for verification purposes. Emphasis is placed on the comparison of
numerical results with analytical ones, as well as modal analysis for
beams with different material gradation profiles. Second, wave
propagation patterns are explored for a fixedfree slender bar
considering homogeneous, bimaterial, trilayered and smoothly graded
materials (steel/alumina), which also provide further verification of
the numerical procedures. Comparison of stress histories in these
samples indicates that the smooth transition of material gradation
considerably alleviates the stress discontinuity in the bimaterial
system (with sharp interface). Third, a threepointbending epoxy/glass
graded beam specimen is investigated for validation purposes. The beam
is graded along the height direction. Stress evolution history at a
location of interest is analyzed in detail, which not only reveals the
dependence of stress evolution on material gradation direction, but
also provides information predictive of potential material failure time
for graded beams with different material gradation profiles. Jointly,
these three classes of examples provide proper verification and
validation for the present numerical techniques.


Back to Top 
2006  
W. G. Buttlar, G. H. Paulino and S. H. Song. "Application of graded finite elements for asphalt pavement." ASCE Journal of Engineering Mechanics. Vol. 132, No. 3, pp. 240249. 2006. Link to PDF View/Hide Abstract
Asphalt paving layers, particularly the surface course, exhibit vertically graded material properties. This grading is caused
primarily by temperature gradients and aging related stiffness gradients. Most conventional existing analysis models do not directly
account for the continuous grading of properties in flexible pavement layers. As a result, conventional analysis methods may lead to
inaccurate prediction of pavement responses and distress under traffic and environmental loading. In this paper, a theoretical formulation
for the graded finite element method is provided followed by an implementation using the user material subroutine (UMAT) capability of
the finite element software ABAQUS. Numerical examples using the UMAT are provided to illustrate the benefits of using graded elements
in pavement analysis.


Y. S. Chan, G. H. Paulino and A. C. Fannjiang. "Change of constitutive relations due to interaction between straingradient effect and material gradation." ASME Journal of Applied Mechanics. Vol. 73, No. 5, pp. 871875. 2006. Link to PDF View/Hide Abstract
For classical elasticity, the constitutive equations (Hooke’s law) have the same functional form for both homogeneous and nonhomogeneous materials. However, for straingradient elasticity, such is not the case. This paper shows that for straingradient elasticity with volumetric and surface energy (Casal’s continuum), extra terms appear in the constitutive equations which are associated with the interaction between the material gradation and the nonlocal effect of strain gradient. The corresponding governing partial differential equations are derived and their solutions are discussed.


Y. Chen, L. J. Struble, and G. H. Paulino. "Extrudability of cementbased materials." American Ceramic Society Bulletin 85 (6), pp. X1X5. 2006. View/Hide Abstract
The extrudability of cementitious materials with the use of penetration resistance (PR) test is investigated. The extrusion of cementitious pastes requires controlling the rheology to obtain goodquality products. PR is used to measure the setting time of concrete in the ASTM Standard C 403 and PR has been related to the rheological parameters of a cementitious paste. A Type I Portland cement was used and methyl hydroxyethyl cellulose was used to produce extrudable pastes. The test was performed at room temperature and was consolidated into a rectangular container. A series of needles was used to penetrate the paste and a smaller needle was used as the PR increased. The vertical force was applied and the force was recorded at each penetration against hydration time. It was observed that the upper limit for PR corresponding to current extrusion setup was above 2000 kPa and the paste was to dry to be extruded. After mixing with water paste was extrudable and in the range of 102000 kPa.


G. H. Paulino, H. M. Yin and L. Z. Sun. "Micromechanicsbased interfacial debonding model for damage of functionally graded materials with particle interactions." International Journal of Damage Mechanics, special issue on "Microstructural Mechanics and Damage Mechanics of Materials." Vol. 15, No. 3, pp. 267288. 2006. Link to PDF View/Hide Abstract
A micromechanical damage model is developed for twophase functionally graded materials (FGMs) considering the interfacial debonding of particles and pairwise interactions between particles. Given an applied mechanical loading on the upper and lower boundaries of an FGM, in the particle–matrix zones, interactions from all other particles over the representative volume element (RVE) are integrated to calculate the homogenized elastic fields. A transition function is constructed to solve the elastic field in the transition zone. The progressive damage process is dependent on the applied loading and is represented by the debonding angles which are obtained from the relation between particle stress and interfacial strength. In terms of the elastic equivalency, debonded, isotropic particles are replaced by perfectly bonded, orthotropic particles. Correspondingly, the effective elasticity distribution in the gradation direction is solved. The computational implementation is discussed and numerical simulations are provided to illustrate the capability of the proposed model.


G. H. Paulino and A. Sutradhar. "The simple boundary element method for multiple cracks in functionally graded media governed by potential theory: a threedimensional Galerkin approach." International Journal for Numerical Methods in Engineering. Vol. 65, No. 12, pp. 20072034. 2006. Link to PDF View/Hide Abstract
The simple boundary element method consists of recycling existing codes for homogeneous media to solve problems in nonhomogeneous media while maintaining a purely boundaryonly formulation. Within this scope, this paper presents a ‘simple’ Galerkin boundary element method for multiple cracks in problems governed by potential theory in functionally graded media. Steadystate heat conduction is investigated for thermal conductivity varying either parabolically, exponentially, or trigonometrically in one or more coordinates. A threedimensional implementation which merges the dual boundary integral equation technique with the Galerkin approach is presented. Special emphasis is given to the treatment of crack surfaces and boundary conditions. The test examples simulated with the present method are verified with finite element results using graded finite elements. The numerical examples demonstrate the accuracy and efficiency of the present method especially when multiple interacting cracks are involved.


D. J. Shim, G. H. Paulino and R. H. Dodds Jr. "A boundary layer framework considering material gradation effects." Engineering Fracture Mechanics. Vol. 73, No. 5, pp. 593615. 2006. Link to PDF View/Hide Abstract
This paper describes the development and application of a novel modified boundary layer (MBL) model for graded nonhomogeneous materials, e.g. functionally graded materials (FGMs). The proposed model is based on a middlecrack tension, M(T), specimen with traction boundary conditions applied to the top and lateral edges of the model. Finite element analyses are performed using WARP3D, a fracture mechanics research finite element code, which incorporates elements with graded elastic and plastic properties. Elastic cracktip fields obtained from the proposed MBL model show excellent agreement with those obtained from full models of the cracked component for homogeneous and graded nonhomogeneous materials. The K–T dominance of FGMs is investigated by comparing the actual stress fields with the asymptotic stress fields (the Williams solution). The examples investigated in the present study consider a crack parallel to the material gradient. Results of the present study provide insight into the K–T dominance of FGMs and also show the range of applicability of the proposed MBL model. The MBL model is applied to analyze the elastic–plastic cracktip response of a Ti/TiB FGM SE(T) specimen. The numerical results demonstrate that the proposed MBL model captures the effect of Tstress on plastic zone size and shape, constraint effects, etc. for such configurations.


D. J. Shim, G. H. Paulino and R. H. Dodds Jr. "Effect of material gradation on Kdominance of fracture specimens." Engineering Fracture Mechanics. Vol. 73, No. 5, pp. 643648. 2006. Link to PDF View/Hide Abstract
This work describes the effect of material gradation (parallel to the crack plane) on stress intensity factors and Kdominance, i.e. the dominance of the singular region, of fracture specimens; SE(T), SE(B) and C(T). The extent of Kdominance is investigated by comparing the actual stress field with the Williams asymptotic stress field. Linearelastic finite element analyses are performed using graded elements which incorporate graded material properties at the element level. Material gradation and crack geometry are systematically varied to perform the parametric study. Results reveal that the effect of material gradation on KI is most pronounced when a short crack is located on the stiffer side of the fracture specimen. For a given specimen and crack geometry, the extent of Kdominance yields a curve with a peak point at a certain material gradation. Results of the present study provide valuable insight into the Kdominance of FGMs.


D. J. Shim, G. H. Paulino and R. H. Dodds Jr. "J resistance behavior in functionally graded materials using cohesive zone and modified boundary layer models." International Journal of Fracture. Vol. 139, No. 1, pp. 91117. 2006. Link to PDF View/Hide Abstract
This paper describes elastic–plastic crack growth resistance simulation in a ceramic/metal functionally graded material (FGM) under mode I loading conditions using cohesive zone and modified boundary layer (MBL) models. For this purpose, we first explore the applicability of two existing, phenomenological cohesive zone models for FGMs. Based on these investigations, we propose a new cohesive zone model. Then, we perform crack growth simulations for TiB/Ti FGM SE(B) and SE(T) specimens using the three cohesive zone models mentioned above. The crack growth resistance of the FGM is characterized by the J integral. These results show that the two existing cohesive zone models overestimate the actual J value, whereas the model proposed in the present study closely captures the actual fracture and crack growth behaviors of the FGM. Finally, the cohesive zone models are employed in conjunction with the MBL model. The two existing cohesive zone models fail to produce the desired K–T stress field for the MBL model. On the other hand, the proposed cohesive zone model yields the desired K–T stress field for the MBL model, and thus yields JR curves that match the ones obtained from the SE(B) and SE(T) specimens. These results verify the application of the MBL model to simulate crack growth resistance in FGMs.


E. C. N. Silva, M. C. Walters and G. H. Paulino. "Modeling bamboo as a functionally graded material: lessons for the analysis of affordable materials." Journal of Materials Science, special issue on "Materials for Affordable Housing and Infrastructure." Vol. 41, No. 21, pp. 69917004. 2006. Link to PDF View/Hide Abstract
Natural fibers are promising for engineering applications due to their low cost. They are abundantly available in tropical and subtropical regions of the world, and they can be employed as construction materials. Among natural fibers, bamboo has been widely used for housing construction around the world. Bamboo is an optimized composite that exploits the concept of Functionally Graded Material (FGM). Biological structures such as bamboo have complicated microstructural shapes and material distribution, and thus the use of numerical methods such as the finite element method, and multiscale methods such as homogenization, can help to further understanding of the mechanical behavior of these materials. The objective of this work is to explore techniques such as the finite element method and homogenization to investigate the structural behavior of bamboo. The finite element formulation uses graded finite elements to capture the varying material distribution through the bamboo wall. To observe bamboo behavior under applied loads, simulations are conducted under multiple considerations such as a spatially varying Young’s modulus, an averaged Young’s modulus, and orthotropic constitutive properties obtained from homogenization theory. The homogenization procedure uses effective, axisymmetric properties estimated from the spatially varying bamboo composite. Threedimensional models of bamboo cells were built and simulated under tension, torsion, and bending load cases.


S. H. Song, G. H. Paulino, and W. G. Buttlar. "A bilinear cohesive zone model tailored for fracture of asphalt concrete considering viscoelastic bulk material." Engineering Fracture Mechanics. Vol. 73, No.18, pp. 28292848. 2006. Link to PDF View/Hide Abstract
A bilinear cohesive zone model (CZM) is employed in conjunction with a viscoelastic bulk (background) material to investigate fracture behavior of asphalt concrete. An attractive feature of the bilinear CZM is a potential reduction of artificial compliance inherent in the intrinsic CZM. In this study, finite material strength and cohesive fracture energy, which are cohesive parameters, are obtained from laboratory experiments. Finite element implementation of the CZM is accomplished by means of a usersubroutine which is employed in a commercial finite element code (e.g., UEL in ABAQUS). The cohesive parameters are calibrated by simulation of mode I diskshaped compact tension results. The ability to simulate mixedmode fracture is demonstrated. The singleedge notched beam test is simulated where cohesive elements are inserted over an area to allow cracks to propagate in any general direction. The predicted mixedmode crack trajectory is found to be in close agreement with experimental results. Furthermore, various aspects of CZMs and fracture behavior in asphalt concrete are discussed including: compliance, convergence, and energy balance.


S. H. Song and G. H. Paulino. "Dynamic stress intensity factors for homogeneous and smoothly heterogeneous materials using the interaction integral method." International Journal of Solids and Structures. Vol. 43, No. 16, pp. 48304866. 2006. Link to PDF View/Hide Abstract
Dynamic stress intensity factors (DSIFs) are important fracture parameters in understanding and predicting dynamic fracture behavior of a cracked body. To evaluate DSIFs for both homogeneous and nonhomogeneous materials, the interaction integral (conservation integral) originally proposed to evaluate SIFs for a static homogeneous medium is extended to incorporate dynamic effects and material nonhomogeneity, and is implemented in conjunction with the finite element method (FEM). The technique is implemented and verified using benchmark problems. Then, various homogeneous and nonhomogeneous cracked bodies under dynamic loading are employed to investigate dynamic fracture behavior such as the variation of DSIFs for different material property profiles, the relation between initiation time and the domain size (for integral evaluation), and the contribution of each distinct term in the interaction integral.


S. H. Song, G. H. Paulino and W. G. Buttlar. "Simulation of crack propagation in asphalt concrete using a intrinsic cohesive zone model." ASCE Journal of Engineering Mechanics. Vol. 132, No. 11, pp. 12151223. 2006. Link to PDF View/Hide Abstract
This is a practical paper which consists of investigating fracture behavior in asphalt concrete using an intrinsic cohesive zone model CZM. The separation and traction response along the cohesive zone ahead of a crack tip is governed by an exponential cohesive law specifically tailored to describe cracking in asphalt pavement materials by means of softening associated with the cohesive law. Finiteelement implementation of the CZM is accomplished by means of a user subroutine using the user element capability of the ABAQUS software, which is verified by simulation of the double cantilever beam test and by comparison to closedform solutions. The cohesive parameters of finite material strength and cohesive fracture energy are calibrated in conjunction with the singleedge notched beam [SE(B)] test. The CZM is then extended to simulate mixedmode crack propagation in the SE(B) test. Cohesive elements are inserted over an area to allow cracks to propagate in any direction. It is shown that the simulated crack trajectory compares favorably with that of experimental results.


M. C. Walters, G. H. Paulino, and R. H. Dodds, Jr. "Computation of mixedmode stress intensity factors for cracks in threedimensional functionallygraded solids." ASCE Journal of Engineering Mechanics. Vol. 132, No. 1, pp. 115. 2006. Link to PDF View/Hide Abstract
This work applies a twostate interaction integral to obtain stress intensity factors along cracks in threedimensional functionally graded materials. The procedures are applicable to planar cracks with curved fronts under mechanical loading, including crackface tractions. Interactionintegral terms necessary to capture the effects of material nonhomogeneity are identical in form to terms that arise due to crackfront curvature. A discussion reviews the origin and effects of these terms, and an approximate interactionintegral expression that omits terms arising due to curvature is used in this work to compute stress intensity factors. The selection of terms is driven by requirements imposed by material nonhomogeneity in conjunction with appropriate mesh discretization along the crack front. Aspects of the numerical implementation with (isoparametric) graded finite elements are addressed, and examples demonstrate the accuracy of the proposed method.


Back to Top 
2005  
W. Celes, G. H. Paulino and R. Espinha. "A compact adjacencybased topological data structure for finite element mesh representation." International Journal for Numerical Methods in Engineering. Vol. 64, No. 11, pp. 15291556. 2005. Link to PDF View/Hide Abstract
This paper presents a novel compact adjacencybased topological data structure for finite element mesh representation. The proposed data structure is designed to support, under the same framework, both two and threedimensional meshes, with any type of elements defined by templates of ordered nodes. When compared to other proposals, our data structure reduces the required storage space while being ‘complete’, in the sense that it preserves the ability to retrieve all topological adjacency relationships in constant time or in time proportional to the number of retrieved entities. Element and node are the only entities explicitly represented. Other topological entities, which include facet, edge, and vertex, are implicitly represented. In order to simplify accessing topological adjacency relationships, we also define and implicitly represent oriented entities, associated to the use of facets, edges, and vertices by an element. All implicit entities are represented by concrete types, being handled as values, which avoid usual problems encountered in other reduced data structures when performing operations such as entity enumeration and attribute attachment. We also extend the data structure with the use of ‘reverse indices’, which improves performance for extracting adjacency relationships while maintaining storage space within reasonable limits. The data structure effectiveness is demonstrated by two different applications: for supporting fragmentation simulation and for supporting volume rendering algorithms.


W. Celes, G. H. Paulino and R. Espinha. "Efficient handling of implicit entities in reduced mesh representations." Journal of Computing and Information Science in Engineering, special issue on "MeshBased Geometric Data Processing." Vol. 5, No. 4, pp. 348359. 2005. Link to PDF View/Hide Abstract
Stateoftheart numerical analyses require mesh representation with a data structure that provides topological information. Due to the increasing size of the meshes currently used for simulating complex behaviors with finite elements or boundary elements (e.g., adaptive and/or coupled analyses), several researchers have proposed the use of reduced mesh representations. In a reduced representation, only a few types of the defined topological entities are explicitly represented; all the others are implicit and retrieved “onthefly,” as required. Despite being very effective in reducing the memory space needed to represent large models, reduced representations face the challenge of ensuring the consistency of all implicit entities when the mesh undergoes modifications. As implicit entities are usually described by references to explicit ones, modifying the mesh may change the way implicit entities (which are not directly modified) are represented, e.g., the referenced explicit entities may no longer exist. We propose a new and effective strategy to treat implicit entities in reduced representations, which is capable of handling transient nonmanifold configurations. Our strategy allows, from the application point of view, explicit and implicit entities to be interchangeably handled in a uniform and transparent way. As a result, the application can list, access, attach properties to, and hold references to implicit entities, and the underlying data structure ensures that all such information remains valid even if the mesh is modified. The validity of the proposed approach is demonstrated by running a set of computational experiments on different models subjected to dynamic remeshing operations.


J. H. Kim and G. H. Paulino. "Consistent formulations of the interaction integral method for fracture of functionally graded materials." ASME Journal of Applied Mechanics. Vol. 72, No. 3, pp. 351364. 2005. Link to PDF View/Hide Abstract
The interaction integral method provides a unified framework for evaluating fracture parameters (e.g., stress intensity factors and T stress) in functionally graded materials. The method is based on a conservation integral involving auxiliary fields. In fracture of nonhomogeneous materials, the use of auxiliary fields developed for homogeneous materials results in violation of one of the basic relations of mechanics, i.e., equilibrium, compatibility or constitutive, which naturally leads to three independent formulations: “nonequilibrium,” “incompatibility,” and “constantconstitutivetensor.” Each formulation leads to a consistent form of the interaction integral in the sense that extra terms are added to compensate for the difference in response between homogeneous and nonhomogeneous materials. The extra terms play a key role in ensuring path independence of the interaction integral. This paper presents a critical comparison of the three consistent formulations and addresses their advantages and drawbacks. Such comparison is made both from a theoretical point of view and also by means of numerical examples. The numerical implementation is based on finite elements which account for the spatial gradation of material properties at the element level (graded elements).


J. H. Kim and G. H. Paulino. "On accurate numerical evaluation of stress intensity factors and Tstress in functionally graded materials." Materials Science Forum 492493, pp. 403408. 2005. Link to PDF View/Hide Abstract
This paper revisits the interaction integral method to evaluate both the mixedmode stress intensity factors and the Tstress in functionally graded materials under mechanical loading. A nonequilibrium formulation is developed in an equivalent domain integral form, which is naturally suitable to the finite element method. Graded material properties are integrated into the element stiffness matrix using the generalized isoparametric formulation. The type of material gradation considered includes continuum functions, such as an exponential function, but the present formulation can be readily extended to micromechanical models. This paper presents a fracture problem with an inclined center crack in a plate and assesses the accuracy of the present method compared with available semianalytical solutions.


J. H. Kim and G. H. Paulino. "Mixedmode crack propagation in functionally graded materials." Materials Science Forum. Vols. 492493, pp. 409414. 2005. Link to PDF View/Hide Abstract
This paper presents numerical simulation of mixedmode crack propagation in functionally graded materials by means of a remeshing algorithm in conjunction with the finite element method. Each step of crack growth simulation consists of the calculation of the mixedmode stress intensity factors by means of a nonequilibrium formulation of the interaction integral method, determination of the crack growth direction based on a specific fracture criterion, and local automatic remeshing along the crack path. A specific fracture criterion is tailored for FGMs based on the assumption of local homogenization of asymptotic cracktip fields in FGMs. The present approach uses a userdefined crack increment at the beginning of the simulation. Crack trajectories obtained by the present numerical simulation are compared with available experimental results.


G. H. Paulino and Z. Zhang. "Dynamic fracture of functionally graded composites using an intrinsic cohesive zone model." Materials Science Forum. Vols. 492493, pp. 447452. 2005. Link to PDF View/Hide Abstract
This paper presents a Cohesive Zone Model (CZM) approach for investigating dynamic failure processes in homogeneous and Functionally Graded Materials (FGMs). The failure criterion is incorporated in the CZM using both a ßnite cohesive strength and work to fracture in the material description. A novel CZM for FGMs is explored and incorporated into a finite element framework. The material gradation is approximated at the element level using a graded element formulation. A numerical example is provided to demonstrate the eácacy of the CZM approach, in which the inàuence of the material gradation on the crack branching pattern is studied.


G. H. Paulino and E. C. N. Silva. "Design of functionally graded structures using topology optimization." Materials Science Forum. Vols. 492493, pp. 435440. 2005. Link to PDF View/Hide Abstract
The concept of functionally graded materials (FGMs) is closely related to the concept of topology optimization, which consists in a design method that seeks a continuum optimum material distribution in a design domain. Thus, in this work, topology optimization is applied to design FGM structures considering a minimum compliance criterion. The present approach applies the socalled “continuous topology optimization” formulation where a continuous change of material properties is considered inside the design domain by using the graded finite element concept. A new design is obtained where distribution of the graded material itself is considered in the design domain, and the material properties change in a certain direction according to a specified variation, leading to a structure with asymmetric stiffness properties.


G. H. Paulino and A. Sutradhar. "A simple Galerkin boundary element method for threedimensional crack problems in functionally graded materials." Materials Science Forum. Vols. 492493, pp. 367372. 2005. Link to PDF View/Hide Abstract
This paper presents a Galerkin boundary element method for solving crack problems governed by potential theory in nonhomogeneous media. In the simple boundary element method, the nonhomogeneous problem is reduced to a homogeneous problem using variable transformation. Cracks in heat conduction problem in functionally graded materials are investigated. The thermal conductivity varies parabolically in one or more coordinates. A three dimensional boundary element implementation using the Galerkin approach is presented. A numerical example demonstrates the eáciency of the method. The result of the test example is in agreement with ßnite element simulation results.


S. H. Song, G. H. Paulino, and W. G. Buttlar. "Cohesive zone simulation of mode I and mixedmode crack propagation in asphalt concrete." Geotechnical Special Publication (130142), pp. 189198. 2005. View/Hide Abstract
A cohesive zone model (CZM) is employed to investigate fracture behavior of asphalt concrete. The separation and traction response along the cohesive zone ahead of a crack tip is modeled by an exponential cohesive law specifically tailored to describe the cracking in asphalt pavement materials by means of a softening cohesive law. This exponential cohesive model is implemented into a userdefined element (UEL) of the ABAQUS software. Using the CZM, first, crack propagation in a mode I singleedge notched beam (SE(B)) test is simulated such that the cohesive parameters of finite material strength and fracture energy are calibrated based upon the experimental results. Then, the mixedmode SE(B) test is simulated using the calibrated cohesive parameters. The crack trajectory of the numerical simulation is found to compare favorably with experimental results.


A. Sutradhar, G. H. Paulino and L. J. Gray. "On hypersingular surface integrals in the symmetric Galerkin boundary element method: Application to heat conduction in exponentially graded materials." International Journal for Numerical Methods in Engineering. Vol. 62, No. 1, pp.122157. 2005. Link to PDF View/Hide Abstract
A symmetric Galerkin formulation and implementation for heat conduction in a threedimensional functionally graded material is presented. The Green’s function of the graded problem, in which the thermal conductivity varies exponentially in one coordinate, is used to develop a boundaryonly formulation without any domain discretization. The main task is the evaluation of hypersingular and singular integrals, which is carried out using a direct ‘limit to the boundary’ approach. However, due to complexity of the Green’s function for graded materials, the usual direct limit procedures have to be modified, incorporating Taylor expansions to obtain expressions that can be integrated analytically. Several test examples are provided to verify the numerical implementation. The results of test calculations are in good agreement with exact solutions and corresponding finite element method simulations.


E. E. Theotokoglou and G. H. Paulino. "A crack in the homogeneous half plane interacting with a crack at the interface between the nonhomogeneous coating and the homogeneous halfplane." International Journal of Fracture. Vol. 134, No. 1, pp. L11L18. 2005. Link to PDF View/Hide Abstract
A fundamental solution is established for a crack in a homogeneous halfplane interacting with a crack at the interface between the homogeneous elastic halfplane and the nonhomogeneous elastic coating in which the shear modulus varies exponentially with one coordinate. The problem is solved under plane strain or generalized plane stress conditions using the Fourier integral transform method. The stress field in the homogeneous half plane is evaluated by the superposition of two states of stresses, one of which is associated with a local coordinate system in the infinite fractured plate, while the other one in the infinite half plane defined in a structural coordinate system.


E. E. Theotokoglou and G. H. Paulino. "Interaction between an embedded crack and an interface crack in nonhomogeneous coating system." Materials Science Forum. Vols. 492493, pp. 397402. 2005. Link to PDF View/Hide Abstract
A general methodology is constructed for the fundamental solution of a crack in the homogeneous halfplane interacting with a crack at the interface between the homogeneous elastic halfplane and the nonhomogeneous elastic coating in which the shear modulus varies exponentially with one coordinate. The problem is solved under plane strain or generalized plane stress condition using the Fourier integral transform method. The stress field in the homogeneous half plane is evaluated by the superposition of two states of stresses, one is associated with a local coordinate system in the infinite fractured plate, while the other in the infinite half plane defined in a structural coordinate system.


M. P. Wagoner, W. G. Buttlar, and G. H. Paulino. "Development of a singleedge notched beam test for asphalt concrete mixtures." Journal of Testing and Evaluation. Vol. 33, No. 6, pp. 452460. 2005. Experimental Mechanics. Vol. 45, No. 3, pp. 270277. 2005. Link to PDF View/Hide Abstract
This paper describes the development of a fracture test for determining the fracture energy of asphalt concrete. The test will be used in combination with numerical analysis and field studies to obtain a better understanding of the mechanisms of reflective cracking in asphalt concrete overlays. A review of the literature revealed that a singleedge notched beam (SE(B)) test specimen was the most promising fracture test for the objectives of the reflective cracking study. Existing servohydraulic testing equipment was modified to perform the SE(B) test along with new loading fixtures, sensors, data collection, and analysis procedures. Preliminary tests were conducted to develop test procedures, to obtain a better understanding of crackfront characteristics, to investigate test repeatability, to examine variations of fracture energy with temperature, and to investigate mixedmode fracture. The results from the tests follow expected trends and test variability appears to be within a range typical for asphalt concrete fracture testing.


M. P. Wagoner, W. G. Buttlar, and G. H. Paulino. "Development of a singleedge notched beam test for the study of asphalt concrete fracture." Geotechnical Special Publication (130142), pp. 137149. 2005.
View/Hide Abstract
The fundamental mechanisms behind the initiation and propagation of reflective cracks in asphalt concrete overlays are not well understood. A novel approach is currently underway to investigate reflective cracking by integrating laboratory tests and numerical analysis to develop and validate an approach to predict reflective cracking in pavement structures. Asphalt overlays are subjected to thermal and mechanical loadings, which induce a combination of tension and shear in the overlay. The selection of a laboratory fracture test that has the ability to induce a range of modemixity at the crack tip is paramount for the integrated approach to be successful. A promising test geometry that has successfully been used for obtaining Mode I (tensile opening) fracture properties of asphalt concrete is the singleedge notched beam (SE(B)). With a slight modification to the notch location, mixedmode (combination of opening and shearing) fracture can be investigated with the notched beam geometry using the same testing fixtures, instrumentation, etc. Test procedures have been developed for the SE(B) test that allow investigation of the crack front profile, discrete crack growth using crack detection gages, and variation of fracture energy with temperature. A single mixture was tested using the proposed mixedmode test geometry with various notch offsets. Through the initial stages of the SE(B) test development, the test provided satisfactory results in terms of the expected trends in fracture energy.


M. P. Wagoner, W. G. Buttlar and G. H. Paulino. "Diskshaped compact tension test for asphalt concrete fracture." Experimental Mechanics. Vol. 45, No. 3, pp. 270277. 2005. Link to PDF View/Hide Abstract
A diskshaped compact tension (DC(T)) test has been developed as a practical method for obtaining the fracture energy of asphalt concrete. The main purpose of the development of this specimen geometry is the ability to test cylindrical cores obtained from inplace asphalt concrete pavements or gyratorycompacted specimens fabricated during the mixture design process. A suitable specimen geometry was developed using the ASTM E399 standard for compact tension testing of metals as a starting point. After finalizing the specimen geometry, a typical asphalt concrete surface mixture was tested at various temperatures and loading rates to evaluate the proposed DC(T) configuration. The variability of the fracture energy obtained from the DC(T) geometry was found to be comparable with the variability associated with other fracture tests for asphalt concrete. The ability of the test to detect changes in the fracture energy with the various testing conditions (temperature and loading rate) was the benchmark for determining the potential of using the DC(T) geometry. The test has the capability to capture the transition of asphalt concrete from a brittle material at low temperatures to a more ductile material at higher temperatures. Because testing was conducted on ungrooved specimens, special care was taken to quantify deviations of the crack path frorn the pure mode I crack path. An analysis of variance of test data revealed that the prototype DC(T) can detect statistical differences Jn fracture energy resulting for tests conducted across a useful range of test temperatures and loading rates. This specific "analysis also indicated that fracture energy is not correlated to crack deviation angle. This paper also provides an overvwew of ongoing work integrating experimental results and observations with numerical analysis by means of a cohesive zone model tailored for asphalt concrete fracture behavior.


M. P. Wagoner, W. G. Buttlar, G. H. Paulino, and P. Blankenship. "Investigation of the fracture resistance of hotmix asphalt concrete using a diskshaped compact tension test." Transportation Research Record (1929), pp. 183192. 2005. Link to PDF View/Hide Abstract
In recent years the transportation materials research community has focused a great deal of attention on the development of testing and analysis methods to shed light on fracture development in asphalt pavements. Recently it has been shown that crack initiation and propagation in asphalt materials can be realistically modeled with cuttingedge computational fracture mechanics tools. However, much more progress is needed toward the development of practical laboratory fracture tests to support these new modeling approaches. The goal of this paper is twofold: (a) to present a diskshaped compact tension [DC(T)] test, which appears to be a practical method for determining lowtemperature fracture properties of cylindrically shaped asphalt concrete test specimens, and (b) to illustrate how the DC(T) test can be used to obtain fracture properties of asphalt concrete specimens obtained from Held cores following dynamic modulus and creep compliance tests performed on the same specimens. Testing four mixtures with varied composition demonstrated that the DC(T) test could detect the transition from quasibrittle to brittle fracture by testing at several low temperatures selected to span across the glass transition temperatures of the asphalt binder used. The tendency toward brittle fracture with increasing loading rate was also detected. Finally, the DC(T) test was used in a forensic study to investigate premature reflective cracking of an isolated portion of pavement in Rochester, New York. One benefit of the DC(T) test demonstrated during testing of field samples was the ability to obtain mixture fracture properties as part of an efficient suite of tests performed on cylindrical specimens.


M. C. Walters, G. H. Paulino and R. H. Dodds, Jr. "Interactionintegral procedures for 3D curved cracks including surface tractions." Engineering Fracture Mechanics. Vol. 72, No. 11, pp. 16351775. 2005. Link to PDF View/Hide Abstract
This study examines a twostate interaction integral for the direct computation of mixedmode stress intensity factors along curved cracks under remote mechanical loads and applied crackface tractions. We investigate the accuracy of stress intensity factors computed along planar, curved cracks in homogeneous materials using a simplified interaction integral that omits terms to reflect specifically the effects of local crackfront curvature. We examine the significance of the crackface traction term in the interaction integral, and demonstrate the benefit of a simple, exact numerical integration procedure to evaluate the integral for one class of threedimensional elements. The work also discusses two approaches to compute auxiliary, interaction integral quantities along cracks discretized by linear and curved elements. Comparisons of numerical results with analytical solutions for stress intensity factors verify the accuracy of the proposed interaction integral procedures.


H. M. Yin, G. H. Paulino, W. G. Buttlar and L. Z. Sun. "Effective thermal conductivity of twophase functionally graded particulate composites." Journal of Applied Physics. Vol. 98, No. 6, Article 063704 (9 pages). 2005. Link to PDF View/Hide Abstract
A multiscale modeling method is proposed to derive effective thermal conductivity in twophase graded particulate composites. In the particlematrix zone, a graded representative volume element is constructed to represent the random microstructure at the neighborhood of a material point. At the steady state, the particle’s averaged heat flux is solved by integrating the pairwise thermal interactions from all other particles. The homogenized heat flux and temperature gradient are further derived, through which the effective thermal conductivity of the graded medium is calculated. In the transition zone, a transition function is introduced to make the homogenized thermal fields continuous and differentiable. By means of temperature boundary conditions, the temperature profile in the gradation direction is solved. When the material gradient is zero, the proposed model can also predict the effective thermal conductivity of uniform composites with the particle interactions. Parametric analyses and comparisons with other models and available experimental data are presented to demonstrate the capability of the proposed method.


H. M. Yin, L. Z. Sun, and G. H. Paulino. "A multiscale framework for elastic deformation of functionally graded composites." Materials Science Forum. Vols. 492493, pp. 391396. 2005. Link to PDF View/Hide Abstract
A micromechanicsbased elastic model is developed for twophase functionally graded composites with locally pairwise particle interactions. In the gradation direction, there exist two microstructurally distinct zones: particlematrix zone and transition zone. In the particlematrix zone, the homogenized elastic fields are obtained by integrating the pairwise interactions from all other particles over the representative volume element. In the transition zone, a transition function is constructed to make the homogenized elastic fields continuous and differentiable in the gradation direction. The averaged elastic fields are solved for transverse shear loading and uniaxial loading in the gradation direction.


Z. Zhang and G. H. Paulino. "Cohesive zone modeling of dynamic failure in homogeneous and functionally graded materials." International Journal of Plasticity, special issue on "Inelastic Response of Multiphase Materials." Vol. 21, No. 6, pp. 11951254. 2005. Link to PDF View/Hide Abstract
This work investigates dynamic failure processes in homogeneous and functionally graded materials (FGMs). The failure criterion is incorporated in the cohesive zone model (CZM) using both a finite cohesive strength and work to fracture in the material description. A novel CZM for FGMs is explored and incorporated into a finite element framework. The material gradation is approximated at the element level using a graded element formulation. Examples are provided to verify the numerical approach, and to investigate the influence of material gradation on crack initiation and propagation in ModeI as well as in mixedmode fracture problems. The examples include spontaneous rapid crack growth in homogeneous and FGM strips, dynamic crack propagation in actual monolithic and epoxy/glass FGM beams (threepoint bending) under impact loading, and mixedmode crack propagation in precracked steel and graded plates.


Back to Top 
2004  
Y. S. Chan, L. J. Gray, T. Kaplan and G. H. Paulino. "Green's function for a twodimensional exponentially graded elastic medium." Proceedings of the Royal Society of London. A: Mathematical, Physical and Engineering Sciences. Vol. 460, No. 2046, pp. 16891706. 2004. Link to PDF View/Hide Abstract
The freespace Green function for a twodimensional exponentially graded elastic medium is derived. The shear modulus p is assumed to be an exponential function of the Cartesian coordinates (x, y), i.e. mu
equivalent to mu(x, y) = mu(0)e(2(beta1x+beta2y)), where mu(0), beta(1)
and beta(2) are material constants, and the Poisson ratio is assumed
constant. The Green function is shown to consist of a singular part, involving modified Bessel func tions, and a nonsingular term. The nonsingular component is expressed in terms of onedimensional Fouriertype integrals that can be computed by the fast Fourier transform.


N. A. Dumont, R. A. P Chaves, and G. H. Paulino. "The hybrid boundary element method applied to problems of potential theory in nonhomogeneous materials." International Journal of Computational Engineering Science. Vol. 5, No. 4, pp. 863891. 2004. Link to PDF View/Hide Abstract
Since the introduction of the hybrid boundary element method in 1987, it has been applied to various problems of elasticity and potential theory, including timedependent problems. This paper focuses on establishing the conceptual framework for applying both the variational formulation and a simplified version of the hybrid boundary element method to nonhomogeneous materials. Several classes of fundamental solutions for problems of potential are derived. Thus, the boundaryonly feature of the method is preserved even with a spatially varying material property. Several numerical examples are given in terms of an efficient patch test including irregularly bounded, unbounded, and multiply connected regions submitted to high gradients.


J. H. Kim and G. H. Paulino. "On the Poisson's ratio effect on mixedmode stress intensity factors and Tstress in functionally graded materials." International Journal of Computational Engineering Science, special issue on “Modelling of Functionally Graded Materials." Vol. 5, No. 4, pp. 833861. 2004. Link to PDF View/Hide Abstract
Poisson's ratio is an important factor for fracture of functionally graded materials (FGMs). It may have significant influence on fracture parameters (e.g. stress intensity factors and Tstress) for a crack in FGMs under mixedmode loading conditions, while its effect on such parameters is negligible in homogeneous materials. For instance, when tension load is applied in the direction parallel to material gradation, the fracture parameters may show significant influence on the Poisson's ratio. This paper uses a new formulation, socalled nonequilibrium formulation, of the interaction integral method. It also presents a few numerical examples where Poisson's ratio is assumed either constant or linearly varying function, and Young's modulus is assumed to be exponential or hyperbolictangent function.


J. H. Kim and G. H. Paulino. “Simulation of crack propagation in functionally graded materials under mixedmode and nonproportional loading.” International Journal of Mechanics and Materials in Design. Vol. 1, No. 1, pp. 6394. 2004. Link to PDF View/Hide Abstract
Automatic simulation of crack propagation in homogeneous and functionally graded materials is performed by means of a remeshing algorithm in conjunction with the finite element method. The crack propagation is performed under mixedmode and nonproportional loading. Each step of crack growth simulation consists of calculation of mixedmode stress intensity factors by means of a novel formulation of the interaction integral method, determination of crack growth direction based on a specific fracture criterion, and local automatic remeshing along the crack path. The present approach requires a userdefined crack increment at the beginning of the simulation. Crack trajectories obtained by the present numerical simulation are compared with available experimental results.


J. H. Kim and G. H. Paulino. "Tstress in orthotropic functionally graded materials: Lekhnitskii and Stroh formalisms." International Journal of Fracture. Vol. 126, No. 4, pp. 345384. 2004. Link to PDF View/Hide Abstract
A new interaction integral formulation is developed for evaluating the elastic Tstress for mixedmode crack problems with arbitrarily oriented straight or curved cracks in orthotropic nonhomogeneous materials. The development includes both the Lekhnitskii and Stroh formalisms. The former is physical and relatively simple, and the latter is mathematically elegant. The gradation of orthotropic material properties is integrated into the element stiffness matrix using a “generalized isoparametric formulation” and (special) graded elements. The specific types of material gradation considered include exponential and hyperbolictangent functions, but micromechanics models can also be considered within the scope of the present formulation. This paper investigates several fracture problems to validate the proposed method and also provides numerical solutions, which can be used as benchmark results (e.g. investigation of fracture specimens). The accuracy of results is verified by comparison with analytical solutions.


G. H. Paulino. "Guest editorial: Modeling of functionally graded materials." International Journal of Computational Engineering Science 5 (4), pp. iii. 2004. 

G. H. Paulino and J. H. Kim. "A new approach to compute Tstress in functionally graded materials by means of the interaction integral method." Engineering Fracture Mechanics. Vol. 71, Nos. 1314, pp.19071950. 2004. Link to PDF View/Hide Abstract
A ‘‘nonequilibrium’’ formulation is developed for evaluating T stress in functionally graded materials with mixedmode cracks. The T stress is evaluated by means of the interaction integral (conservation integral) method in conjunction with the finite element method. The gradation of material properties is integrated into the element stiffness matrix using the socalled ‘‘generalized isoparametric formulation’’. The types of material gradation considered include exponential, linear, and radially graded exponential functions; however, the present formulation is not limited to specific functions and can be readily extended to micromechanics models. This paper investigates several fracture problems (including both homogeneous and functionally graded materials) to verify the proposed formulation, and also provides numerical solutions to various benchmark problems. The accuracy of numerical results is discussed by comparison with available analytical, semianalytical, or numerical solutions. The revisited interaction integral method is shown to be an accurate and robust scheme for evaluating T stress in functionally graded materials.


A. Sutradhar and G. H. Paulino. "A simple boundary element method for problems of potential problems in nonhomogeneous media." International Journal for Numerical Methods in Engineering. Vol. 60, No. 13, pp. 22032230. 2004. Link to PDF View/Hide Abstract
A simple boundary element method for solving potential problems in nonhomogeneous media is presented. A physical parameter (e.g. heat conductivity, permeability, permittivity, resistivity, magnetic permeability) has a spatial distribution that varies with one or more coordinates. For certain classes of material variations the nonhomogeneous problem can be transformed to known homogeneous problems such as those governed by the Laplace, Helmholtz and modified Helmholtz equations. A threedimensional Galerkin boundary element method implementation is presented for these cases. However, the present development is not restricted to Galerkin schemes and can be readily extended to other boundary integral methods such as standard collocation. A few test examples are given to verify the proposed formulation. The paper is supplemented by an Appendix, which presents an ABAQUS usersubroutine for graded finite elements. The results from the finite element simulations are used for comparison with the present boundary element solutions.


A. Sutradhar and G. H. Paulino. "The simple boundary element method for transient heat conduction in functionally graded materials." Computer Methods in Applied Mechanics and Engineering. Vol. 193, No. 4244, pp. 45114539. 2004. Link to PDF View/Hide Abstract
This paper presents a ‘‘simple’’ boundary element method for transient heat conduction in functionally graded materials, which leads to a boundaryonly formulation without any domain discretization. For a broad range of functional material variation (quadratic, exponential and trigonometric) of thermal conductivity and specific heat, the nonhomogeneous problem can be transformed into the standard homogeneous diffusion problem. A threedimensional boundary element implementation, using the Laplace transform approach and the Galerkin approximation, is presented. The time dependence is restored by numerically inverting the Laplace transform by means of the Stehfest algorithm. A number of numerical examples demonstrate the efficiency of the method. The results of the test examples are in excellent agreement with analytical solutions and finite element simulation results.


A. Sutradhar and G. H. Paulino. "Symmetric Galerkin boundary element computation of Tstress and stress intensity factors for mixedmode cracks by the interaction integral method." Engineering Analysis with Boundary Elements. Vol. 28, No. 11, pp. 13351350. 2004. Link to PDF View/Hide Abstract
An interaction integral method for evaluating Tstress and mixedmode stress intensity factors (SIFs) for twodimensional crack problems using the symmetric Galerkin boundary element method is presented. By introducing a known auxiliary field solution, the SIFs for both mode I and mode II are obtained from a pathindependent interaction integral entirely based on the Jð; J1Þ integral. The Tstress can be calculated directly from the interaction integral by simply changing the auxiliary field. The numerical implementation of the interaction integral method is relatively simple compared to other approaches and the method yields accurate results. A number of numerical examples with available analytical and numerical solutions are examined to demonstrate the accuracy and efficiency of the method.


M. C. Walters, G. H. Paulino, and R. H. Dodds. "Stressintensity factors for surface cracks in functionally graded materials under modeI thermomechanical loading." International Journal of Solids and Structures 41 (34), pp. 10811118. 2004. Link to PDF View/Hide Abstract
This paper describes the development and application of a general domain integral method to obtain Jvalues along crack fronts in threedimensional configurations of isotropic, functionally graded materials (FGMs). The present work considers modeI, linearelastic response of cracked specimens subjected to thermomechanical loading, although the domain integral formulation accommodates elasticplastic behavior in FGMs. Finite element solutions and domain integral Jvalues for a twodimensional edge crack show good agreement with available analytical solutions for both tension loading and temperature gradients. A displacement correlation technique provides pointwise stressintensity values along semielliptical surface cracks in FGMs for comparison with values derived from the proposed domain integral. Numerical implementation and mesh refinement issues to maintain path independent Jvalues are explored. The paper concludes with a parametric study that provides a set of stressintensity factors for semielliptical surface cracks covering a practical range of crack sizes, aspect ratios and material property gradations under tension, bending and spatiallyvarying temperature loads.


H. M. Yin, L. Z. Sun, and G. H Paulino. "Micromechanical modeling of functionally graded composites." American Society of Mechanical Engineers, Applied Mechanics Division, AMD 255, pp. 18. 2004. View/Hide Abstract
A micromechanicsbased elastic model is developed for twophase functionally graded materials with locally pairwise interactions between particles. While the effective material properties change gradually along the gradation direction, there exist two microstructurally distinct zones: particlematrix zone and transition zone. In the particlematrix zone, pairwise interactions between particles are employed using a modified Green's function method. By integrating the interactions from all other particles over the representative volume element, the homogenized elastic fields are obtained. The effective stiffness distribution over the gradation direction is further derived. In the transition zone, a transition function is constructed to make the homogenized elastic fields continuous and differentiable in the gradation direction. The model prediction is compared with other models and experimental data to demonstrate the capability of the proposed method.


H. M. Yin, L. Z. Sun and G. H. Paulino. "Micromechanicsbased elastic model for functionally graded materials with particle interactions." Acta Materialia. Vol. 52, No. 12, pp. 3535–3543. 2004. Link to PDF View/Hide Abstract
A micromechanicsbased elastic model is developed for twophase functionally graded materials with locally pairwise interactions between particles. While the effective material properties change gradually along the gradation direction, there exist two microstructurally distinct zones: particle–matrix zone and transition zone. In the particle–matrix zone, pairwise interactions between particles are employed using a modified Green’s function method. By integrating the interactions from all other particles over the representative volume element, the homogenized elastic fields are obtained. The effective stiffness distribution over the gradation direction is further derived. In the transition zone, a transition function is constructed to make the homogenized elastic fields continuous and differentiable in the gradation direction. The model prediction is compared with other models and experimental data to demonstrate the capability of the proposed method.


Back to Top 
2003  
Y. S. Chan, A. C. Fannjiang, and G. H. Paulino. "Integral equations with hypersingular kernels  theory and applications to fracture mechanics." International Journal of Engineering Science. Vol. 41, No. 7, pp. 683720. 2003. Link to PDF View/Hide Abstract
Hypersingular integrals were investigated for general integers a (positive) and m (nonnegative). The integrals were evaluated to calculate the stress intensity factors. Examples involving crack problems were given and discussed with emphasis on the linkage between mathematics and mechanics of fracture. Analysis shows that as material property variation in space and higher order graded continuum theories are considered, the formulation of the crack problem and the associated kernels becomes quite involved.


L. J. Gray, T. Kaplan, G. H. Paulino and J. D Richardson. "Green's functions and boundary integral analysis for exponentially graded materials: Heat conduction." ASME Journal of Applied Mechanics. Vol. 70, No. 4, pp. 543549. 2003. Link to PDF View/Hide Abstract
Free space Green’s functions are derived for graded materials in which the thermal conductivity varies exponentially in one coordinate. Closedform expressions are obtained for the steadystate diffusion equation, in two and three dimensions. The corresponding boundary integral equation formulations for these problems are derived, and the threedimensional case is solved numerically using a Galerkin approximation. The results of test calculations are in excellent agreement with exact solutions and finite element simulations.


L. J. Gray, A. V. Phan, G. H. Paulino and T. Kaplan. "Improved quarterpoint crack tip element." Engineering Fracture Mechanics. Vol. 70, No. 2, pp. 269283. 2003. Link to PDF View/Hide Abstract
We present a modification to the quarterpoint crack tip element and employ this element in twodimensional boundary integral fracture analysis. The standard singular element is adjusted so that the neartip crack opening displacement satisfies a known constraint: the coefficient of the term which is linear in the distance to the tip must vanish. Stress intensity factors calculated with the displacement correlation technique are shown to be highly accurate, and significantly more accurate than with the standard element. The improvements are especially dramatic for mixedmode problems involving curved and interacting cracks.


Z. H. Jin, G. H. Paulino and R. H. Dodds. "Cohesive fracture modeling of elasticplastic crack growth in functionally graded materials." Engineering Fracture Mechanics. Vol. 70, No. 14, pp. 18851912. 2003. Link to PDF View/Hide Abstract
This work investigates elastic–plastic crack growth in ceramic/metal functionally graded materials (FGMs). The study employs a phenomenological, cohesive zone model proposed by the authors and simulates crack growth by the gradual degradation of cohesive surfaces ahead of the crack front. The cohesive zone model uses six materialdependent parameters (the cohesive energy densities and the peak cohesive tractions of the ceramic and metal phases, respectively, and two cohesive gradation parameters) to describe the constitutive response of the material in the cohesive zone. A volume fraction based, elastic–plastic model (extension of the original Tamura–Tomota–Ozawa model) describes the elastic–plastic response of the bulk background material. The numerical analyses are performed using WARP3D, a fracture mechanics research finite element code, which incorporates solid elements with graded elastic and plastic properties and interfacecohesive elements coupled with the functionally graded cohesive zone model. Numerical values of volume fractions for the constituents specified at nodes of the finite element model set the spatial gradation of material properties with isoparametric interpolations inside interface elements and background solid elements to define pointwise material property values. The paper describes applications of the cohesive zone model and the computational scheme to analyze crack growth in a singleedge notch bend, SE(B), specimen made of a TiB/Ti FGM. Cohesive parameters are calibrated using the experimentally measured load versus average crack extension (across the thickness) responses of both Ti metal and TiB/Ti FGM SE(B) specimens. The numerical results show that with the calibrated cohesive gradation parameters for the TiB/Ti system, the load to cause crack extension in the FGM is much smaller than that for the metal. However, the crack initiation load for the TiB/Ti FGM with reduced cohesive gradation parameters (which may be achieved under different manufacturing conditions) could compare to that for the metal. Crack growth responses vary strongly with values of the exponent describing the volume fraction profile for the metal. The investigation also shows significant crack tunneling in the Ti metal SE(B) specimen. For the TiB/Ti FGM system, however, crack tunneling is pronounced only for a metalrich specimen with relatively smaller cohesive gradation parameter for the metal.


J. H. Kim and G. H. Paulino. "An accurate scheme for mixedmode fracture analysis of functionally graded materials using the interaction integral and micromechanics models." International Journal for Numerical Methods in Engineering. Vol. 58, No. 10, pp. 14571497. 2003. Link to PDF View/Hide Abstract
The interaction integral is a conservation integral that relies on two admissible mechanical states for evaluating mixedmode stress intensity factors (SIFs). The present paper extends this integral to functionally graded materials in which the material properties are determined by means of either continuum functions (e.g. exponentially graded materials) or micromechanics models (e.g. selfconsistent, Mori–Tanaka, or threephase model). In the latter case, there is no closedform expression for the materialproperty variation, and thus several quantities, such as the explicit derivative of the strain energy density, need to be evaluated numerically (this leads to several implications in the numerical implementation). The SIFs are determined using conservation integrals involving known auxiliary solutions. The choice of such auxiliary elds and their implications on the solution procedure are discussed in detail. The computational implementation is done using the nite element method and thus the interaction energy contour integral is converted to an equivalent domain integral over a nite region surrounding the crack tip. Several examples are given which show that the proposed method is convenient, accurate, and computationally efficient.


J. H. Kim and G. H. Paulino. "The interaction integral for fracture of orthotropic functionally graded materials: Evaluation of stress intensity factors." International Journal of Solids and Structures. Vol.40, No. 15, pp. 39674001. 2003. Link to PDF View/Hide Abstract
The interaction integral is an accurate and robust scheme for evaluating mixedmode stress intensity factors. This paper extends the concept to orthotropic functionally graded materials and addresses fracture mechanics problems with arbitrarily oriented straight and/or curved cracks. The gradation of orthotropic material properties are smooth functions of spatial coordinates, which are integrated into the element stiffness matrix using the socalled ‘‘generalized isoparametric formulation’’. The types of orthotropic material gradation considered include exponential, radial, and hyperbolictangent functions. Stress intensity factors for mode I and mixedmode twodimensional problems are evaluated by means of the interaction integral and the finite element method. Extensive computational experiments have been performed to validate the proposed formulation. The accuracy of numerical results is discussed by comparison with available analytical, semianalytical, or numerical solutions.


J. H. Kim and G. H. Paulino. "Mixedmode Jintegral formulation and implementation using graded elements for fracture analysis of nonhomogeneous orthotropic materials." Mechanics of Materials. Vol. 35, No. 12, pp. 107128. 2003. Link to PDF View/Hide Abstract
The pathindependent J*k integral, in conjunction with the finite element method (FEM), is presented for mode I and mixedmode crack problems in orthotropic functionally graded materials (FGMs) considering plane elasticity. A general procedure is presented where the crack is arbitrarily oriented, i.e. it does not need to be aligned with the principal orthotropy directions. Smooth spatial variations of the independent engineering material properties are incorporated into the element stiffness matrix using a ‘‘generalized isoparametric formulation’’, which is natural to the FEM. Both exponential and linear variations of the material properties are considered. Stress intensity factors and energy release rates for pure mode I and mixedmode boundary value problems are numerically evaluated by means of the equivalent domain integral especially tailored for orthotropic FGMs. Numerical results are discussed and validated against available theoretical and numerical solutions.


J. H. Kim and G. H. Paulino. "Erratum: Mixedmode Jintegral formulation and implementation using graded elements for fracture analysis of nonhomogeneous orthotropic materials." Mechanics of Materials. Vol. 35, No. 12, pp. 107128, 2003. Link to PDF 

J. H. Kim and G. H. Paulino. "TStress, mixedmode stress intensity factors, and crack initiation angles in functionally graded materials: a unified approach using the interaction integral method." Computer Methods in Applied Mechanics and Engineering. Vol. 192, No. 1112, pp. 14631494. 2003. Link to PDF View/Hide Abstract
For linear elastic functionally graded materials (FGMs), the fracture parameters describing the crack tip fields include not only stress intensity factors (SIFs) but also Tstress (nonsingular stress). These two fracture parameters are important for determining the crack initiation angle under mixedmode loading conditions in brittle FGMs (e.g. ceramic/ceramic such as TiC/SiC). In this paper, the mixedmode SIFs and Tstress are evaluated by means of the interaction integral, in the form of an equivalent domain integral, in combination with the finite element method. In order to predict the crack initiation angle in brittle FGMs, this paper makes use of a fracture criterion which incorporates the Tstress effect. This type of criterion involves the mixedmode SIFs, the Tstress, and a physical length scale rc (representative of the fracture process zone size). Various types of material gradations are considered such as continuum models (e.g. exponentially graded material) and micromechanics models (e.g. selfconsistent model). Several examples are given to show the accuracy and efficiency of the interaction integral scheme for evaluating mixedmode SIFs, Tstress, and crack initiation angle. The techniques developed provide a basic framework for quasistatic crack propagation in FGMs.


G. Li, G. H. Paulino, and N. R. Aluru. "Coupling of the meshfree finite cloud method with the boundary element method: a collocation approach." Computer Methods in Applied Mechanics and Engineering. Vol. 192, No. 2021 , pp. 23552375. 2003. Link to PDF View/Hide Abstract
Meshless and meshbased methods are among the tools frequently applied in the numerical treatment of partial differential equations (PDEs). This paper presents a coupling of the meshless finite cloud method (FCM) and the standard (meshbased) boundary element method (BEM), which is motivated by the complementary properties of both methods. While the BEM is appropriate for solving linear PDEs with constant coefficients in bounded or unbounded domains, the FCM is appropriate for either homogeneous, inhomogeneous or even nonlinear problems in bounded domains. Both techniques (FCM and BEM) use a collocation procedure in the numerical approximation. No mesh is required in the FCM subdomain and its nodal point distribution is completely independent of the BEM subdomain. The coupling approach is demonstrated for linear elasticity problems. Because both FCM and BEM employ traction displacement relations, no transformations between forces and tractions (typical of BEM and finite element coupling) are needed. Several numerical examples are given to demonstrate the proposed approach.


S. Mukherjee and G. H. Paulino. "The elasticviscoelastic correspondence principle for functionally graded materials, revisited." Journal of Applied Mechanics, Transactions ASME. Vol. 70, No. 3, pp. 359363. 2003. Link to PDF View/Hide Abstract
Paulino and Jin [Paulino, G. H., and Jin, Z.H., 2001, "Correspondence Principle in Viscoelastic Functionally Graded Materials," ASME J. Appl. Mech., 68, pp. 129132], have recently shown that the viscoelastic correspondence principle remains valid for a linearly isotropic viscoelastic functionally graded material with separable relaxation (or creep) functions in space and time. This paper revisits this issue by addressing some subtle points regarding this result and examines the reasons behind the success or failure of the correspondence principle for viscoelastic functionally graded materials. For the inseparable class of nonhomogeneous materials, the correspondence principle fails because of an inconsistency between the replacements of the moduli and of their derivatives. A simple but informative onedimensional example, involving an exponentially graded material, is used to further clarify these reasons.


G. H. Paulino, A. Sutradhar, and L. J. Gray. "Boundary element methods for functionally graded materials." Boundary Elements 4, pp. 137146. 2003. View/Hide Abstract
Functionally graded materials (FGMs) possess a smooth variation of material properties due to gradual change in microstructural details. For example, the material gradation may change gradually from a pure ceramic to a pure metal. This work focuses on threedimensional potential (both steady state and transient) and elasticity problems for inhomogeneous materials. The Green's function (GF) for these materials (exponentially graded) are expressed as the GF for the homogeneous material plus additional terms due to material gradation. The numerical implementations are performed using a Galerkin (rather than collocation) approximation. The results of test calculations are in good agreement with available analytical solution and finite element simulations.


G. H. Paulino, A. C. Fannjiang and Y. S. Chan. "Gradient elasticity theory for mode III fracture in functionally graded materialsPart I: crack perpendicular to the material gradation." ASME Journal of Applied Mechanics. Vol. 70, No. 4, pp. 531542. 2003. Link to PDF View/Hide Abstract
Anisotropic strain gradient elasticity theory is applied to the solution of a mode III crack in a functionally graded material.The
theory possesses two material characteristic lengths, l and l', which
describe the size scale effect resulting from the underlining
microstructure, and are associated to volumetric and surface strain
energy, respectively. The governing differential equation of the
problem is derived assuming that the shear modulus is a function of the
Cartesian coordinate y, i.e., G G (y) = G(0)e(gammay), where G(0) and y
are material constants. The crack boundary value problem is solved by means of Fourier transforms and the hypersingular integrodifferential equation method. The integral equation is discretized using the collocation method and a Chebyshev polynomial expansion. Formulas for stress intensity factors, KIII , are derived, and numerical results of KIII for various combinations of C, l', and gamma are provided. Finally, conclusions are inferred and potential extensions of this work are discussed.


A. Sutradhar, G. H. Paulino, and L. J. Gray. "Erratum: Transient heat conduction in homogeneous and nonhomogeneous materials by the Laplace transform Galerkin boundary element method (Engineering Analysis with Boundary Elements (2002) 26 (119132)." Engineering Analysis with Boundary Elements. Vol. 27, No. 6, pp. 639. 2003. 

Back to Top 
2002  
Y. S. Chan, A. C. Fannjiang, and G. H. Paulino. "Integral equations with hypersingular kernels  Theory and applications to fracture mechanics." International Journal of Engineering Science. Vol. 41, No. 7, pp. 683720. 2002. Link to PDF View/Hide Abstract
Hypersingular integrals of the type
Ialpha(Tn,m,r) = integral(1)(1) Tn(s)(1s(2))(m1/2)/(sr)(alpha) ds, r < 1 and Ialpha(Un,m,r) = integral(1)(1) Un(s)(1s(2))(m1/2)/(sr)(alpha) ds, r < 1 are investigated for general integers alpha (positive) and m (nonnegative), where Tn(s) and Un(s) are the Chebyshev polynomials of the first and second kinds, respectively. Exact formulas are derived for the cases alpha = 1, 2, 3, 4 and m = 0, 1, 2, 3; most of them corresponding to new solutions derived in this paper. Moreover, a systematic approach for evaluating these integrals when a > 4 and m > 3 is provided. The integrals are also evaluated as r > 1 in order to calculate the stress intensity factors. Examples involving crack problems are given and discussed with emphasis on the linkage between mathematics and mechanics of fracture. The examples include classical linear elastic fracture mechanics, functionally graded materials, and gradient elasticity theory. An appendix, with an alternative derivation of the formulae, supplements the paper. 

N. A. Dumont, R. A. P. Chaves, and G. H. Paulino. "The hybrid boundary element method applied to functionally graded materials." International Series on Advances in Boundary Elements. Vol. 5, No. 4, pp. 267276. 2002. View/Hide Abstract
The hybrid boundary element method has been successfully applied to various problems of elasticity and potential theory. This paper focuses on establishing the conceptual framework for applying the hybrid boundary element method to functionally graded materials. Several classes of fundamental solutions for problems of potential are derived. In particular, a recently developed fundamental solution is employed to model heat conduction in materials with exponentially varying thermal conductivity. Thus, the boundaryonly feature of the method is preserved even with such spatially varying material property. Two numerical examples are given in terms of a patch test for irregular regions submitted to high gradients.


A. C. Fannjiang, Y. S. Chan and G. H. Paulino. "Strain gradient elasticity for mode III cracks: A hypersingular integrodifferential equation approach." SIAM Journal on Applied Mathematics. Vol. 62, No. 3, pp. 10661091. 2002. Link to PDF View/Hide Abstract
We consider Casal’s strain gradient elasticity with two material lengths L, L associated with volumetric and surface energies, respectively. For a Mode III finite crack we formulate a hypersingular integrodifferential equation for the crack slope supplemented with the natural cracktip conditions. The fullfield solution is then expressed in terms of the crack profile and the Green function, which is obtained explicitly. For L = 0, we obtain a closed form solution for the crack profile. The case of small L is shown to be a regular perturbation. The question of convergence, as L, L' > 0, is studied in detail both analytically and numerically.


M. R. Hill, R. D. Carpenter, G. H. Paulino, Z. A. Munir, and J. C. Gibeling. "Fracture testing of a layered functionally graded material." ASTM Special Technical Publication (1409), pp. 169184. 2002.
View/Hide Abstract
This paper describes measurements of the K1 Rcurve for a layered ceramicmetallic functionally graded material (FGM) composed of Ti and TiB phases. Single edge notch bend specimens were fabricated for crack propagation perpendicular to the graded layers from the brittle to the ductile side of the FGM. The precracking method and residual stresses affected the measured toughness. A new reverse bending method produced a sharp precrack without damaging the material. A representative sample indicates Rcurve behavior rising from Kmeas, = 7 MPa·m1/2 in the initial cracktip location (38% Ti and 62% TiB) to 31 MPa·m1/2 in pure Ti. Residual stresses acting normal to the crackplane were measured in the parent plate using the crack compliance technique. The residual stress intensity factor in the SE(B) specimen determined from the weight function was 2.4 MPa·m1/2 at the initial precrack length. Thus, the actual material toughness at the start of the test was 34% less than the measured value.


Z. H. Jin, G. H. Paulino, and R. H. Dodds, Jr. "Finite element investigation of quasistatic crack growth in functionally graded materials using a novel cohesive zone fracture model." ASME Journal of Applied Mechanics. Vol. 69, No. 3, pp. 370379. 2002. Link to PDF View/Hide Abstract
This work studies mode I crack growth in ceramic/metal functionally graded materials (FGMs) using threedimensional interfacecohesive elements based upon a new phenomenological cohesive fracture model. The local separation energies and peak tractions for the metal and ceramic constituents govern the cohesive fracture process. The model formulation introduces two cohesive gradation parameters to control the transition of fracture behavior between the constituents. Numerical values of volume fractions for the constituents specified at nodes of the finite element model set the spatial gradation of material properties with standard isoparametric interpolations inside interface elements and background solid elements to define pointwise material property values. The paper describes applications of the cohesive fracture model and computational scheme to analyze crack growth in compact tension, C(T), and singleedge notch bend, SE(B), specimens with material properties characteristic of a TiB/Ti FGM. Young’s modulus and Poisson’s ratio of the background solid material are determined using a selfconsistent method (the background material remains linear elastic). The numerical studies demonstrate that the load to cause crack extension in the FGM compares to that for the metal and that crack growth response varies strongly with values of the cohesive gradation parameter for the metal. These results suggest the potential to calibrate the value of this parameter by matching the predicted and measured crack growth response in standard fracture mechanics specimens.


Z. H. Jin and G. H. Paulino. "A viscoelastic functionally graded strip containing a crack subjected to inplane loading." Engineering Fracture Mechanics, special issue on "Fracture of Functionally Graded Materials." Vol. 69, No. 1416, pp. 17691790. 2002. Link to PDF View/Hide Abstract
In this paper, a crack in a viscoelastic strip of a functionally graded
material (FGM) is studied under tensile loading conditions. The
extensional relaxation modulus is assumed as E = E(0) exp(beta*y/h)f(t), where h is a scale length and f(t) is a nondimensional function of
time t either having the form f(t) = E(infinity)/E(0) + (1 
E(infinity)/E(0)) exp(t/t(0)) for a linear standard solid or f(t) =
(t(0)/t)(q) for a power law material model, where E(0), E(infinity),
beta, t(0) and q are material constants. An extensional relaxation
function in the form E = E(0) exp(beta*y/h)[t(0) exp(delta*y/h)/t](q) is
also considered, in which the relaxation time depends on the Cartesian
coordinate y exponentially with delta being a material constant
describing the gradation of the relaxation time. The Poisson's ratio is
assumed to have the form nu = nu(0)(1 + gamma*y/h) exp(beta*y/h)g(t),
where nu(0) and gamma are material constants, and g(t) is a
nondimensional function of time t. An elastic FGM crack problem is first solved and
the ‘‘correspondence principle’’ is used to obtain both mode I and mode II stress intensity factors, and the crack
opening/sliding displacements for the viscoelastic FGM considering various material models.


J. H. Kim and G. H. Paulino. "Finite element evaluation of mixedmode stress intensity factors in functionally graded materials." International Journal for Numerical Methods in Engineering. Vol. 53, No. 8, pp. 19031935. 2002. Link to PDF View/Hide Abstract
This paper is directed towards finite element computation of fracture parameters in functionally graded material (FGM) assemblages of arbitrary geometry with stationary cracks. Graded finite elements are developed where the elastic moduli are smooth functions of spatial coordinates which are integrated into the element stiffness matrix. In particular, stress intensity factors for mode I and mixedmode twodimensional problems are evaluated and compared through three di=erent approaches tailored for FGMs: pathindependent Jk*integral, modi:ed crackclosure integral method, and displacement correlation technique. The accuracy of these methods is discussed based on comparison with available theoretical, experimental or numerical solutions.


J. H. Kim and G. H. Paulino. "Isoparametric graded finite elements for nonhomogeneous isotropic and orthotropic materials." ASME Journal of Applied Mechanics. Vol. 69, No. 4, pp. 502514. 2002. Link to PDF View/Hide Abstract
Graded finite elements are presented within the framework of a generalized isoparametric formulation. Such elements possess a spatially varying material property field, e.g. Young’s modulus (E) and Poisson’s ratio (v) for isotropic materials; and principal Young’s moduli (E11 ,E22), inplane shear modulus (G_12) , and Poisson’s ratio (v_12) for orthotropic materials. To investigate the influence of material property variation, both exponentially and linearly graded materials are considered and compared. Several boundary value problems involving continuously nonhomogeneous isotropic and orthotropic materials are solved, and the performance of graded elements is compared to that of conventional homogeneous elements with reference to analytical solutions. Such solutions are obtained for an orthotropic plate of infinite length and finite width subjected to various loading conditions. The corresponding solutions for an isotropic plate are obtained from those for the orthotropic plate. In general, graded finite elements provide more accurate local stress than conventional homogeneous elements, however, such may not be the case for fournode quadrilateral (Q4) elements. The framework described here can serve as the basis for further investigations such as thermal and dynamic problems in functionally graded materials.


J. H. Kim and G. H. Paulino. "Mixedmode fracture of orthotropic functionally graded materials using finite elements and the modified crack closure method." Engineering Fracture Mechanics, special issue on "Fracture of Functionally Graded Materials." Vol. 69, No. 1416, pp. 15571586. 2002. Link to PDF View/Hide Abstract
A finite element methodology is developed for fracture analysis of orthotropic functionally graded materials (FGMs) where cracks are arbitrarily oriented with respect to the principal axes of material orthotropy. The graded and orthotropic material properties are smooth functions of spatial coordinates, which are integrated into the element stiffness matrix using the isoparametric concept and special graded finite elements. Stress intensity factors (SIFs) for mode I and mixedmode twodimensional problems are evaluated and compared by means of the modified crack closure (MCC) and the displacement correlation technique (DCT) especially tailored for orthotropic FGMs. An accurate technique to evaluate SIFs by means of the MCC is presented using a simple twostep (predictor–corrector) process in which the SIFs are first predicted (e.g. by the DCT) and then corrected by Newton iterations. The effects of boundary conditions, crack tip mesh discretization and material properties on fracture behavior are investigated in detail. Many numerical examples are given to validate the proposed methodology. The accuracy of results is discussed by comparison with available (semi) analytical or numerical solutions.


S. Mukherjee and G. H. Paulino. "The elasticviscoelastic correspondence principle for functionally graded materials, revisited." ASME Journal of Applied Mechanics. Vol. 70, No. 3, pp. 359363. 2002. Link to PDF View/Hide Abstract
Paulino and Jin [Paulino, G. H., and Jin, Z.H., 2001, ‘‘Correspondence Principle in Viscoelastic Functionally Graded Materials,’’ ASME J. Appl. Mech., 68, pp. 129–132], have recently shown that the viscoelastic correspondence principle remains valid for a linearly isotropic viscoelastic functionally graded material with separable relaxation (or creep) functions in space and time. This paper revisits this issue by addressing some subtle points regarding this result and examines the reasons behind the success or failure of the correspondence principle for viscoelastic functionally graded materials. For the inseparable class of nonhomogeneous materials, the correspondence principle fails because of an inconsistency between the replacements of the moduli and of their derivatives. A simple but informative onedimensional example, involving an exponentially graded material, is used to further clarify these reasons.


G. H. Paulino. "Fracture of functionally graded materials." Engineering Fracture Mechanics. Editorial, Vol. 69, Nos. 1416, pp. 15191520. 2002. Link to PDF View/Hide Abstract
A finite element methodology is developed for fracture analysis of orthotropic functionally graded materials (FGMs) where cracks are arbitrarily oriented with respect to the principal axes of material orthotropy. The graded and orthotropic material properties are smooth functions of spatial coordinates, which are integrated into the element stiffness matrix using the isoparametric concept and special graded finite elements. Stress intensity factors (SIFs) for mode I and mixedmode twodimensional problems are evaluated and compared by means of the modified crack closure (MCC) and the displacement correlation technique (DCT) especially tailored for orthotropic FGMs. An accurate technique to evaluate SIFs by means of the MCC is presented using a simple twostep (predictor–corrector) process in which the SIFs are first predicted (e.g. by the DCT) and then corrected by Newton iterations. The effects of boundary conditions, crack tip mesh discretization and material properties on fracture behavior are investigated in detail. Many numerical examples are given to validate the proposed methodology. The accuracy of results is discussed by comparison with available (semi) analytical or numerical solutions.


A. Sutradhar, G. H. Paulino and L. J. Gray. "Transient heat conduction in homogeneous and nonhomogeneous materials by the Laplace Transform Galerkin boundary element method." Engineering Analysis with Boundary Elements. Vol. 26, No. 2, pp. 119132. 2002. Link to PDF View/Hide Abstract
The Green's function for threedimensional transient heat conduction
(diffusion equation) for functionally graded materials (FGMs) is
derived. The thermal conductivity and heat capacitance both vary
exponentially in one coordinate. In the process of solving this
diffusion problem numerically, a Laplace transform (LT) approach is
used to eliminate the dependence on time. The fundamental solution in
Laplace space is derived and the boundary integral equation formulation
for the Laplace Transform boundary element method (LTBEM) is obtained.
The numerical implementation is performed using a Galerkin
approximation, and the timedependence is restored by numerical
inversion of the LT. Two numerical inversion techniques have been
investigated: a Fourier series method and Stehfest's algorithm, the
latter being preferred. A number of test problems have been examined,
and the results are in excellent agreement with available analytical
solutions.


Back to Top 
2001  
E. M. CarrilloHeian, R. D. Carpenter, G. H. Paulino, J. C. Gibeling, and Z. A. Munir. "Dense layered MoSi2/SiC functionally graded composites formed by field activated synthesis." Journal of the American Ceramic Society. Vol. 84, No. 5, pp. 962968. 2001. Link to PDF View/Hide Abstract
Dense, layered, single and gradedcomposition composites of MoSi2 and SiC were formed from elemental powders in one step, using the fieldactivated pressureassisted combustion method. Compositions ranging from 100% MoSi2 to 100% SiC were synthesized, with relative densities ranging from 99% to 76%, respectively. Xray diffractometry results indicated the formation of pure phases when the concentration of MoSi2 was high and the appearance of a ternary phase, Mo4.8Si3C0.6, when the concentration of SiC was high. Electron microprobe analysis results showed the formation of stoichiometric and uniformly distributed phases. A layertolayer variation in composition of 10 mol% was sufficient to prevent thermal cracking during formation of the layered functionally graded materials.


E. M. CarrilloHeian, C. Unuvar, J. C. Gibeling, G. H. Paulino, and Z. A. Munir. "Simultaneous synthesis and densification of Niobium Silicide/Niobium composites." Scripta Materialia. Vol. 45, No. 4, pp. 405412. 2001. Link to PDF 

Y. S. Chan, G. H. Paulino, and A. C. Fannjiang. "The crack problem for nonhomogeneous materials under antiplane shear loading  A displacement based formulation." International Journal of Solids and Structures. Vol. 38, No. 17, pp. 29893005. 2001. Link to PDF View/Hide Abstract
This paper presents a displacement based integral equation formulation
for the mode III crack problem in a nonhomogeneous medium with a
continuously differentiable shear modulus, which is assumed to be an
exponential function, i.e., G(x) = G(0) exp(beta*x). This formulation
leads naturally to a hypersingular integral equation. The problem is
solved for a finite crack and results are given for crack profiles and
stress intensity factors. The results are affected by the parameter
beta describing the material nonhomogeneity. This study is motivated by
crack problems in straingradient elasticity theories where higher
order singular integral equations naturally arise even in the
slopebased formulation.


M. K. Chati, S. Mukherjee, and G. H. Paulino. "The meshless hypersingular boundary node method for threedimensional potential theory and linear elasticity problems." Engineering Analysis with Boundary Elements. Vol. 25, No. 8, pp. 639653. 2001. Link to PDF View/Hide Abstract
The Boundary Node Method (BNM) represents a coupling between Boundary
Integral Equations (BIEs) and Moving Least Squares (MLS) approximants.
The main idea here is to retain the dimensionality advantage of the
former and the meshless attribute of the latter. The result is a
'meshfree' method that decouples the mesh and the interpolation
procedures. The BNM has been applied to solve 2D and 3D problems in
potential theory and linear elasticity. The Hypersingular Boundary
Element Method (HBEM) has diverse important applications in areas such
as fracture mechanics, wave scattering, error analysis and adaptivity,
and to obtain a symmetric Galerkin boundary element formulation. The
present work presents a coupling of Hypersingular Boundary Integral
Equations (HBIEs) with MLS approximants, to produce a new meshfree
methodthe Hypersingular Boundary Node Method (HBNM). Numerical results
from this new method, for selected 3D problems in potential theory and
in linear elasticity, are presented and discussed in this paper.


M. K. Chati, G. H. Paulino, and S. Mukherjee. "The meshless standard and hypersingular boundary node methods  Applications to error estimation and adaptivity in threedimensional problems." International Journal for Numerical Methods in Engineering. Vol. 50, No. 9, pp. 22332269. 2001. Link to PDF View/Hide Abstract
The standard (singular) boundary node method (BNM) and the novel hypersingular boundary node method (HBNM) are employed for the usual and adaptive solutions of threedimensional potential and elasticity problems. These methods couple boundary integral equations with moving leastsquares interpolants while retaining the dimensionality advantage of the former and the meshless attribute of the latter. The `hypersingular residuals', developed for error estimation in the meshbased collocation boundary element method (BEM) and symmetric Galerkin BEM by Paulino et al., are extended to the meshless BNM setting. A simple `a posteriori' error estimation and an e
ective adaptive renement procedure are presented. The implementation of all the techniques involved in this work are discussed, which includes aspects regarding parallel implementation of the BNM and HBNM codes. Several numerical examples are given and discussed in detail. Conclusions are inferred and relevant extensions of the methodology introduced in this work are provided.


Z. H. Jin and G. H. Paulino. "Transient thermal stress analysis of an edge crack in a functionally graded material." International Journal of Fracture. Vol. 107, No. 1, pp. 7398. 2001. Link to PDF View/Hide Abstract
An edge crack in a strip of a functionally graded material (FGM) is studied under transient thermal loading conditions. The FGM is assumed having constant Young’s modulus and Poisson’s ratio, but the thermal properties of the material vary along the thickness direction of the strip. Thus the material is elastically homogeneous but thermally nonhomogeneous. This kind of FGMs include some ceramic/ceramic FGMs such as TiC/SiC, MoSi2/Al2O3 and MoSi2/SiC, and also some ceramic/metal FGMs such as zirconia/nickel and zirconia/steel. A multilayered material model is used to solve the temperature field. By using the Laplace transform and an asymptotic analysis, an analytical first order temperature solution for short times is obtained. Thermal stress intensity factors (TSIFs) are calculated for a TiC/SiC FGM with various volume fraction profiles of the constituent materials. It is found that the TSIF could be reduced if the thermally shocked cracked edge of the FGM strip is pure TiC, whereas the TSIF is increased if the thermally shocked edge is pure SiC.


G. H. Paulino, G. Menon and S. Mukherjee. "Error estimation using hypersingular integrals in boundary element methods for linear elasticity." Engineering Analysis with Boundary Elements. Vol. 25, No. 7, pp. 523534. 2001. Link to PDF View/Hide Abstract
A natural measure of the error in the boundary element method rests on
the use of both the standard boundary integral equation (BIE) and the
hypersingular BIE (HBIE). An approximate (numerical) solution can be
obtained using either one of the BIEs. One expects that the residual,
obtained when such an approximate solution is substituted to the other
BIE is related to the error in the solution. The present work is
developed for vector field problems of linear elasticity. In this
context, suitable 'hypersingular residuals' are shown, under certain
special circumstances, to be globally related to the error. Further,
heuristic arguments are given for general mixed boundary value
problems. The calculated residuals are used to compute element error
indicators, and these error indicators are shown to compare well with
actual errors in several numerical examples, for which exact errors are
known. Conclusions are drawn and potential extensions of the present
error estimation method are discussed.


G. H. Paulino and Z. H. Jin. "Correspondence principle in viscoelastic functionally graded materials." ASME Journal of Applied Mechanics. Vol. 68, No. 1, pp. 129132. 2001. Link to PDF View/Hide Abstract
This paper presents an extension of the correspondence principle (as applied to homogeneous viscoelastic solids) to nonhomogeneous viscoelastic solids under the assumption that the relaxation (or creep) moduli be separable functions in space and time. A few models for graded viscoelastic materials are presented and discussed. The revisited correspondence principle extends to specific instances of thermoviscoelasticity and fracture of functionally graded materials.


G. H. Paulino and Z. H. Jin. "A crack in a viscoelastic functionally graded material layer embedded between two dissimilar homogeneous viscoelastic layers Antiplane shear analysis." International Journal of Fracture. Vol. 111, No. 3, pp. 283303. 2001. Link to PDF View/Hide Abstract
A crack in a viscoelastic functionally graded material (FGM) layer
sandwiched between two dissimilar homogeneous viscoelastic layers is
studied under antiplane shear conditions. The shear relaxation modulus
of the FGM layer follows the power law of viscoelasticity, i.e., mu (0)
= mu (0) exp(betay/ h)[t(0) exp(delta (y)/h)/t](q), where h is a scale
length, and mu (0), mu (0), beta, delta and q are material constants.
Note that the FGM layer has positiondependent modulus and relaxation
time. The shear relaxation functions of the two homogeneous
viscoelastic layers are ft = mu (1) (t(1)/t)(q) for the bottom layer
and mu = mu (2)(t(2)/t)(q) for the top layer, where mu (1) and mu (2)
are material constants, and t(1) and t(2) are relaxation times. An
elastic crack problem of the composite structure is first solved and
the 'correspondence principle' is used to obtain stress intensity
factors (SIFs) for the viscoelastic system. Formulae for SIFs and crack
displacement profiles are derived. Several examples are given which
include interface cracking between a viscoelastic functionally graded
interlayer and a viscoelastic homogeneous material coating. Moreover, a
parametric study is conducted considering various material and
geometric parameters and loading conditions.


G. H. Paulino, J. C. Gibeling, R. D. Carpenter, W. W. Liang and Z. A. Munir. "Fracture testing and finite element modeling of pure Titanium." Engineering Fracture Mechanics. Vol. 68, No. 12, pp. 14171432. 2001. Link to PDF View/Hide Abstract
This paper presents the results of fracture experiments and
corresponding finite element analyses (FEA) of pure titanium. This
investigation was motivated by the desire to develop a JR testing
protocol and numerical procedures that are applicable to a
titanium/titanium boride layered functionally graded material. Tensile
tests and a twodimensional axisymmetric finite element model were used
to determine the plasticity data for the titanium. Crack growth
experiments were conducted in threepoint bending using single edge
notched bend specimens. Threedimensional FEA of crack growth
initiation and twodimensional FEA with automatic crack propagation
were performed. Two crack propagation conditions based on experimental
data were used: (a) crack length versus loadline displacement and (b)
crack length versus crack mouth opening displacement. The subsequent
predictions of the nonlinear finite element models are in reasonable
agreement with the measured value of J at initiation and with the
rising JR data during crack propagation.


G. H. Paulino and Y. Liu. "Implicit consistent and continuum tangent operators in elastoplastic boundary element formulations." Computer Methods in Applied Mechanics and Engineering. Vol. 190, Nos. 1517, pp. 21572179. 2001. Link to PDF View/Hide Abstract
This paper presents an assessment and comparison of boundary element method (BEM) formulations for elastoplasticity using both the consistent tangent operator (CTO) and the continuum tangent operator (CON). These operators are integrated into a single computational implementation using linear or quadratic elements for both boundary and domain discretizations. This computational setting is also used in the development of a method for calculating the J integral, which is an important parameter in (nonlinear) fracture mechanics. Various twodimensional examples are given and relevant response parameters such as the residual norm, computational processing time, and results obtained at various load and iteration steps, are provided. The examples include fracture problems and J integral evaluation. Finally, conclusions are inferred and extensions of this work are discussed.


G. H. Paulino and Z.H Jin. "Viscoelastic functionally graded materials subjected to antiplane shear fracture." ASME Journal of Applied Mechanics. Vol. 68, No. 2, pp. 284293. 2001. Link to PDF View/Hide Abstract
In this paper, a crack in a strip of a viscoelastic functionally graded
material is studied under antiplane shear conditions. The shear
relaxation function of the material is assumed as mu = mu (0)exp
(betay/h)f(t), where h is a length scale and f(t) is a nondimensional
function of time t having either the form f(t) = mu (infinity)/mu (0) +
(1  mu (0))exp (t/t(0)) for a linear standard solid, or f(t) =
(t(0)/t)q for a powerlaw material model. We also consider the shear
relaxation function mu = mu (0)exp(betay/h)[t(0)exp(deltay/h)/t](q) in
which the relaxation time depends on the Cartesian coordinat y
exponentially. Thus this latter model represents a powerlaw material
with positiondependent relaxation time. In this above expression, the
parameters beta, mu (0), mu (infinity), t(0); delta, q are material
constants. An elastic crack problem is first solved and the
correspondence principle (revisited) is used to obtain stress intensity
factors for the viscoelastic functionally graded material. Formulas for
stress intensity factors and crack displacement profiles are derived.
Results for these quantities are discussed considering various material
models are loading conditions.


Back to Top 
2000 and Prior  
R. D. Carpenter, G. H. Paulino, Z. A. Munir and J. C. Gibeling. "A novel technique to generate sharp cracks in metalic/ceramic functionally graded materials by reverse 4point bending." Scripta Materialia. Vol. 43, No. 6, pp. 547552. 2000. Link to PDF 

Z. H. Jin and G. H. Paulino. "Transient thermal stress analysis of an interior crack in functionally graded materials." American Society of Mechanical Engineers, Aerospace Division (Publication) AD 60, pp. 121125. 2000. View/Hide Abstract
An internal crack in a strip of a functionally graded material (FGM) is studied under transient thermal loading conditions. The FGM is assumed having constant Young's modulus and Poisson's ratio, but the thermal properties of the material vary along the thickness direction of the strip. Thus the material is elastically homogeneous but thermally nonhomogeneous. This kind of FGMs include some ceramic/ceramic FGMs such as TiC/SiC and MoSi2/Al2O3, and also some ceramic/metal FGMs such as zirconia/nickel and zirconia/steel. Thermal stress intensity factors (TSIFs) are calculated for a TiC/SiC FGM with various volume fraction profiles of the constituent materials. It is found that the TSIF could be reduced if the thermally shocked cracked edge of the FGM strip is pure TiC.


L. Yong, G. H. Paulino, and L. Lihua. "Elastoviscoplastic consistent tangent operator conceptbased implicit boundary element methods." Science in China, Series E: Technological Sciences 43 (2), pp. 154164. 2000.
Link to PDF View/Hide Abstract
An elastoviscoplastic consistent tangent operator (CTO) conceptbased implicit algorithm for nonlinear boundary element methods is presented. Both kinematic and isotropic strain hardening are considered. The elastoviscoplastic radial return algorithm (RRA) and the elastoviscoplastic CTO and its related scheme are developed. In addition, the limit cases (e.g. elastoplastic problem) of viscoplastic RRA and CTO are discussed. Finally, numerical examples, which are compared with the latest FEM results of Ibrahimbegovic et al. and ABAQUS results, are provided.


R. D. Carpenter, W. W. Liang, G. H. Paulino, J. C. Gibeling, and Z. A. Munir. "Fracture testing and analysis of a layered functionally graded Ti/TiB beam in 3point bending." Materials Science Forum 308311, pp. 837842. 1999. Link to PDF View/Hide Abstract
This paper presents a methodology for testing and analyzing a layered functionally graded beam consisting of Ti and TiB phases. The test specimen is a single edge notch bending (SENB), with composition varying from a TiBrich layer at the bottom to commercially pure (CP) Ti at the top. The crack front is aligned parallel to the material gradation and it grows along this direction (mode I). For this crack orientation, 'average' Rcurve behavior is investigated in order to understand the mechanics of crack growth. The beam is modeled using twodimensional finite elements considering elastoplastic behavior. Properties of the mixed phase interlayers are approximated using a ruleofmixtures approach, and the validity of this approach is discussed. Characteristic response quantities are investigated such as the mechanical fields, crack mouth opening displacement, and J integral.


E. N. Lages, G. H. Paulino, I. F. M. Menezes, and R. R. Silva. "Nonlinear finite element analysis using an objectoriented philosophy  Application to beam elements and to the Cosserat continuum." Engineering with Computers. Vol. 15, No. 1, pp. 7389. 1999. Link to PDF View/Hide Abstract
An ObjectOriented Programming (OOP) framework is presented for solving nonlinear structural mechanics problems by means of the Finite Element Method (FEM). Emphasis is placed on engineering applications (geometrically nonlinear beam model, and elastoplastic Cosserat continuum), and OOP is employed as an effective tool, which plays an important role in the FEM treatment of such applications. The implementation is based on computational abstractions of both mathematical and physical concepts associated to structural mechanics problems involving geometrical and material nonlinearities. The overall class organization for nonlinear mechanics modeling is discussed in detail. All the analyses rely on a generic control class where several classical and modern nonlinear solution schemes are available. Examples which explore, demonstrate and validate the main features of the overall computational system are presented and discussed.


G. Menon, G. H. Paulino, and S. Mukherjee. "Analysis of hypersingular residual error estimates in boundary element methods for potential problems." Computer Methods in Applied Mechanics and Engineering. Vol. 173, Nos. 34, pp. 449473. 1999. Link to PDF View/Hide Abstract
A novel iteration scheme, using boundary integral equations, is
developed for error estimation in the boundary element method. The
iteration scheme consists of using the boundary integral equation for
solving the boundary value problem and iterating this solution with the
hypersingular boundary integral equation to obtain a new solution. The
hypersingular residual r is consistently defined as the difference in
the derivative quantities on the boundary, i.e.
r = partial derivative phi((1))/partial derivative n  partial derivative phi((2))/partial derivative n where phi is the potential and (partial derivative phi/partial derivative n)((i)), i = i, 2, is the flux obtained by solution (i). Here, i = 1 refers to the boundary integral equation, and i = 2 refers to the hypersingular boundary integral equation. The hypersingular residual is interpreted in the sense of the iteration scheme defined above and it is shown to provide an error estimate for the boundary value problem. Errorhypersingular residual relations are developed for Dirichlet and Neumann problems, which are shown to be limiting cases of the more general relation for the mixed boundary value problem. These relations lead to global bounds on the error. Four numerical examples, involving Galerkin boundary elements, are given, and one of them involves a physical singularity on the boundary and preliminary adaptive calculations. These examples illustrate important features of the hypersingular residual error estimate proposed in this paper. 

G. H. Paulino, A. C. Fannjiang, and Y. S. Chan. "Gradient elasticity theory for a mode III crack in a functionally graded material." Materials Science Forum 308311, pp. 971976. 1999. Link to PDF View/Hide Abstract
Anisotropic strain gradient elasticity theory is applied to the solution of a mode III crack in a functionally graded material (FGM). The theory includes both volumetric and surface energy terms, and a particular form of the moduli variation. The crack boundary value problem is solved by means of Fourier transform and a hypersingular integral equation of the Fredholm type. The solution naturally leads to a cusping crack, which is consistent with Barenblatt's 'cohesive zone' theory, but without the assumption regarding existence of interatomic forces. The numerical implementation is discussed and examples are given, which provide insight into the cracking phenomenon in FGMs governed by strain gradient elasticity with characteristic lengths associated to volumetric and surface strain energy terms.


G. H. Paulino and L. J. Gray. "Galerkin Residuals for adaptive symmetricGalerkin boundary element methods." ASCE Journal of Engineering Mechanics. Vol. 125, No. 5, pp. 575585. 1999. Link to PDF View/Hide Abstract
This paper presents a simple a posteriori error estimator and an effective adaptive mesh refinement procedure for the symmetric Galerkin boundary element method. The ‘‘hypersingular residuals,’’ developed for error estimation in a standard collocation BEM, are extended to the symmetric Galerkin setting. This leads to the formulation of ‘‘Galerkin residuals,’’ which are intrinsic to the symmetric Galerkin boundary integral approach and form the basis of the present error estimation scheme. Several computational experiments are conducted to test both the accuracy and the reliability of the proposed technique. These experiments involve potential theory and various problem configurations including mixed boundary conditions, corners, and nonconvex domains. The numerical results indicate that reliable solutions to practical engineering problems can be obtained with this method.


G. H. Paulino, I. F. M. Menezes, J. B. C. Neto, and L. F. R. C. Martha. "A methodology for adaptive finite element analysis: towards an integrated computational environment." Computational Mechanics. Vol. 23, Nos. 56, pp. 361388. 1999. Link to PDF View/Hide Abstract
This work introduces a methodology for selfadaptive numerical procedures, which relies on the various components of an integrated, objectoriented, computational environment involving pre, analysis, and postprocessing modules. A basic platform for numerical experiments and further development is provided, which allows implementation of new elements/error estimators and sensitivity analysis. A general implementation of the Superconvergent Patch Recovery (SPR) and the recently proposed Recovery by Equilibrium in Patches (REP) is presented. Both SPR and REP are compared and used for error estimation and for guiding the adaptive remeshing process. Moreover, the SPR is extended for calculating sensitivity quantities of ®rst and higher orders. The mesh (re)generation process is accomplished by means of modern methods combining quadtree and Delaunay triangulation techniques. Surface mesh generation in arbitrary domains is performed automatically (i.e. with no user intervention) during the selfadaptive analysis using either quadrilateral or triangular elements. These ideas are implemented in the Finite Element System Technology in Adaptivity (FESTA) software. The effectiveness and versatility of FESTA are demonstrated by representative numerical examples illustrating the interconnections among finite element analysis, recovery procedures, error estimation/adaptivity and automatic mesh generation.


L. J. Gray and G. H. Paulino. "Crack tip interpolation, revisited." SIAM Journal on Applied Mathematics. Vol. 58, No. 2, pp. 428455. 1998. Link to PDF View/Hide Abstract
It is well known that the near tip displacement weld on a crack surface can be represented in a power series in the variable sqrt(r), where r is the distance to the tip. It is shown herein that the coefficients of the linear terms on the two sides of the crack are equal. Equivalently, the linear term in the crack opening displacement vanishes. The proof is a completely general argument, valid for an arbitrary (e.g., multiple, nonplanar) crack conguration and applied boundary conditions. Moreover, the argument holds for other equations, such as Laplace. A limit procedure for calculating the surface stress in the form of a hypersingular boundary integral equation is employed to enforce the boundary conditions along the crack faces. Evaluation of the nite surface stress and examination of potentially singular terms lead to the result. Inclusion of this constraint in numerical calculations should result in a more accurate approximation of the displacement and stress elds in the tip region, and thus a more accurate evaluation of stress intensity factors.


K. Sivathasan, G. H. Paulino, X. S. Li, and K. Arulanandan. "Validation of site characterization method for the study of dynamic pore pressure response." Geotechnical Special Publication (75 I), pp. 469481. 1998. View/Hide Abstract
An attempt is being made to develop procedures to obtain initial state and constitutive model parameters by insitu testing. A centrifuge model test is used in the present study for comparative purposes. A nondestructive electrical method is used to characterize the soil. Representative electrical parameters to the centrifuge model are hypothesized. A fully coupled nonlinear effective stress analysis procedure SUMDES is used in the analysis with soil behavior simulated by a bounding surface hypoplasticity model for sand. The development of liquefaction and the pore water pressure generation is investigated in the fully coupled nonlinear analysis. The excess pore water pressures, calculated from the analysis, are in close agreement with the measured excess pore water pressure during the centrifuge test.


L. J. Gray and G. H. Paulino. "Symmetric Galerkin boundary integral fracture analysis for plane orthotropic elasticity." Computational Mechanics. Vol. 20, Nos..12, pp. 2633. 1997. Link to PDF View/Hide Abstract
This paper discusses the formulation and implementation of the symmetric Galerkin boundary integral method for two dimensional linear elastic orthotropic fracture analysis. For the usual case of a tractionfree crack, the symmetry of the coefficient matrix can be effectively exploited to significantly reduce the computational work required to construct the linear system. Inaddition, computation time is reduced by employing efficient analytic integration formulas for the analysis of the orthotropic singular and hypersingular integrals. Preliminary test calculations indicate that the method is both accurate and efficient.


L. J. Gray and G. H. Paulino. "Symmetric Galerkin boundary integral formulation for interface and multizone problems."International Journal for Numerical Methods in Engineering. Vol. 40, No. 16, pp. 30853101. 1997. Link to PDF View/Hide Abstract
Domains containing an `internal boundary', such as a bimaterial interface, arise in many applications, e.g. composite materials and geophysical simulations. This paper presents a symmetric Galerkin boundary integral method for this important class of problems. In this situation, the physical quantities are known to satisfy continuity conditions across the interface, but no boundary conditions are specied. The algorithm described herein achieves a symmetric matrix of reduced size. Moreover, the symmetry can also be invoked to lessen the numerical work involved in constructing the system of equations, and thus the method is computationally very ecient. A prototype numerical example, with several variations in the boundary conditions and material properties, is employed to validate the formulation and corresponding numerical procedure. The boundary element results are compared with analytical solutions and with numerical results obtained with the finite element method.


S. H. Hsieh, G. H. Paulino and J. F. Abel. "Evaluation of automatic domain partitioning algorithms for parallel finite element analysis." International Journal for Numerical Methods in Engineering. Vol. 40, No. 6, pp. 10251051. 1997. Link to PDF View/Hide Abstract
This paper studies and compares the domain partitioning algorithms presented by Farhat,1 AlNasra and Nguyen,2 Malone,3 and Simon4/Hsieh et al.5, 6 for load balancing in parallel Þnite element analysis. Both the strengths and weaknesses of these algorithms are discussed. Some possible improvements to the partitioning algorithms are also suggested and studied. A new approach for evaluating domain partitioning algorithms is described. Direct numerical comparisons among the considered partitioning algorithms are then conducted using this suggested approach with both regular and irregular Þnite element meshes of di¤erent order and dimensionality. The test problems used in the comparative studies along with the results obtained provide a set of benchmark examples for other researchers to evaluate both new and existing partitioning algorithms. In addition, interactive graphics tools used in this work to facilitate the evaluation and comparative studies are presented.


K. M. Mosalam and G. H. Paulino. "Evolutionary characteristic length method for smeared cracking finite element methods." Finite Elements in Analysis and Design. Vol. 27, No. 1, pp. 99108. 1997. Link to PDF View/Hide Abstract
The fixed smeared crack concept with strain decomposition is reformulated utilizing a selfadaptive strategy at the constitutive level. This formulation focuses on continuous adaptation of the crack band width based on the incremental finite element solution and the idea of nonlocal continuum. The required nonlocal forms are obtained by means of the superconvergent patch recovery procedure. Comparison with experimental results indicates the superiority of the present formulation over the standard smeared cracking.


G. H. Paulino, F. Shi, S. Mukherjee and P. S. Ramesh. "Nodal sensitivities as error estimates in computational mechanics." Acta Mechanica. Vol. 121, Nos. 14, pp. 191213. 1997. Link to PDF View/Hide Abstract
This paper proposes the use of special sensitivities, called nodal sensitivities, as error indicators
and estimators for numerical analysis in mechanics. Nodal sensitivities are defined as rates of change of
response quantities with respect to nodal positions. Direct analytical differentiation is used to obtain the
sensitivities, and the infinitesimal perturbations of the nodes are forced to lie along the elements. The idea
proposed here can be used in conjunction with general purpose computational methods such as the Finite
Element Method (FEM), the Boundary Element Method (BEM) or the Finite Difference Method (FDM);
however, the BEM is the method of choice in this paper. The performance of the error indicators is evaluated
through two numerical examples in linear elasticity.


G. H. Paulino, L. J. Gray and V. Zarikian. "Hypersingular residuals  A new approach for error estimation in the boundary element method." International Journal for Numerical Methods in Engineering. Vol. 39, No. 12, pp. 20052029. 1996. Link to PDF View/Hide Abstract
This paper presents a new approach for a posteriori 'pointwise' error estimation in the boundary element method. The estimator relies upon evaluation of the residual of hypersingular integral equations, and is therefore intrinsic to the boundary integral equation approach. A methodology is developed for approximating the error on the boundary as well as in the interior of the domain. Extensive computational experiments have been performed for the twodimensional Laplace equation and the numerical results indicate that the error estimates successfully track the form of the exact error curve. Moreover, a reasonable estimate of the magnitude of the actual error is also predicted.


S. H. Hsieh, G. H. Paulino, and J. F. Abel. "Recursive Spectral algorithms for automatic domain partitioning in parallel finite element analysis." Computer Methods in Applied Mechanics and Engineering. Vol. 121, Nos. 14, pp. 137162. 1995. Link to PDF View/Hide Abstract
Recently, several domain partitioning algorithms have been proposed to effect loadbalancing among processors in parallel finite element analysis. The recursive spectral bisection (RSB) algorithm [l] has been shown to be effective. However, the bisection nature of the RSB results in partitions of an integer power of two, which is too restrictive for computing environments consisting of an arbitrary number of processors. This paper presents two recursive spectral partitioning algorithms, both of which generalize the RSB algorithm for an arbitrary number of partitions. These algorithms are based on a graph partitioning approach which includes spectral techniques and graph representation of finite element meshes. The ‘algebraic connectivity vector’ is introduced as a parameter to assess the quality of the partitioning results. Both nodebased and elementbased partitioning strategies are discussed. The spectral algorithms are also evaluated and compared for coarsegrained partitioning using different types of structures modelled by lD, 2D and 3D finite elements.


G. H. Paulino, I. F. M. Menezes, M. Gattass and S. Mukherjee. "A new algorithm for finding a pseudoperipheral vertex or the endpoints of a pseudodiameter in a graph." Communications in Numerical Methods in Engineering. Vol. 10, No. 11, pp. 913926. 1994. Link to PDF View/Hide Abstract
Based on the concept of the Laplacian matrix of a graph, this paper presents the SGPD (spectral graph pseudoperipheral and pseudodiameter) algorithm for finding a pseudoperipheral vertex or the endpoints of a pseudodiameter in a graph. This algorithm is compared with the ones by Grimes et al. (1990), George and Liu (1979), and Gibbs et al. (1976). Numerical results from a collection of benchmark test problems show the effectiveness of the proposed algorithm. Moreover, it is shown that this algorithm can be efficiently used in conjunction with heuristic algorithms for ordering sparse matrix equations. Such heuristic algorithms, of course, must be the ones which use the pseudoperipheral vertex or pseudodiameter concepts.


G. H. Paulino, I. F. M. Menezes, M. Gattass and S. Mukherjee. "Node and element resequencing using the Laplacian of a finite element graph. Part I: General concepts and algorithm." International Journal for Numerical Methods in Engineering. Vol. 37, No. 9, pp. 15111530. 1994. Link to PDF View/Hide Abstract
A Finite Element Graph (FEG) is defined here as a nodal graph (G), a dual graph (G*), or a communication graph (G') associated with a generic finite element mesh. The Laplacian matrix ((L(G),L(G*) or L(G')), used for the study of spectral properties of an FEG, is constructed from usual vertex and edge connectivities of a graph. An automatic algorithm, based on spectral properties of an FEG (G, G* or G O ) , is proposed to reorder the nodes and/or elements of the associated finite element mesh. The new algorithm is called Spectral FEG Resequencing (SFR). This algorithm uses global information in the graph, it does not depend on a pseudoperipheral vertex in the resequencing process, and it does not use any kind of level structure of the graph. Moreover, the SFR algorithm is of special advantage in computing environments with vector and parallel processing capabilities.
Nodes or elements in the mesh can be reordered depending on the use of an adequate graph representation associated with the mesh. If G is used, then the nodes in the mesh are properly reordered for achieving profile and wavefront reduction of the finite element stiffness matrix. If either G* or G' is used, then the elements in the mesh are suitably reordered for a finite element frontal solver. A unified approach involving FEGs and finite element concepts is presented. Conclusions are inferred and possible extensions of this research are pointed out.
In Part II of this work,' the computational implementation of the SFR algorithm is described and several numerical examples are presented. The examples emphasize important theoretical, numerical and practical aspects of the new resequencing method.


G. H. Paulino, I. F. M. Menezes, M. Gattass, and S. Mukherjee. "Node and element resequencing using the Laplacian of a finite element graph. Part II: Implementation and numerical results." International Journal for Numerical Methods in Engineering. Vol. 37, No. 9, pp. 15311555. 1994. Link to PDF View/Hide Abstract
In Part I of this work, Paulino et a/.' have presented an algorithm for profile and wavefront reduction of large sparse matrices of symmetric configuration. This algorithm is based on spectral properties of a Finite Element Graph (FEG). An FEG has been defined as a nodal graph G, a dual graph G* or a communication graph C' associated with a generic finite element mesh. The novel algorithm has been called Spectral FEG Resequencing (SFR). This algorithm has specific features that distinguish it from previous algorithms. These features include (1) use of global information in the graph, (2) no need of a pseudoperipheral vertex or the endpoints of a pseudodiameter, and (3) no need of any type of level structure of the FEG. To validate this algorithm in a numerical sense, extensive computational testing on a variety of problems is presented here. This includes algorithmic performance evaluation using a library of benchmark test problems which contains both connected and nonconnected graphs, study of the algebraic connectivity (lambda_2) of an FEG, eigensolver convergence verification, running time performance evaluation and assessment of the algorithm on a set of practical finite element examples. It is shown that the SFR algorithm is effective in reordering nodes and/or elements of generic finite element meshes. Moreover, it computes orderings which compare favourably with the ones obtained by some previous algorithms that have been published in the technical literature.


M. Gattass, G. H. Paulino, and J. C. Gortaire. "Geometrical and topological consistency in interactive graphical preprocessors of threedimensional framed structures." Computers and Structures. Vol. 46, No. 1, pp. 99124. 1993. Link to PDF View/Hide Abstract
This paper shows how interactive graphical preprocessors can treat inconsistent data to generate consistent finite element and structural models. Algorithms are described for checking the geometrical and topological consistency of threedimensional meshes made up of onedimensional finite elements. With the consistency guaranteed, the user of a graphical preprocessor can freely create the geometry and topology of a structure with his or her attention concentrated only on the image shown on the screen. Inconsistencies like repeated nodes or partially overlapped bars are left to the system to solve.


G. H. Paulino, M. T. A. Saif, and S. Mukherjee. "A finite elastic body with a curved crack loaded in antiplane shear." International Journal of Solids and Structures. Vol. 30, No. 8, pp. 10151037. 1993. Link to PDF View/Hide Abstract
This paper presents a Boundary Integral Equation Method (BIEM) for an arbitrarily shaped, linearly elastic, homogeneous and isotropic body with a curved crack loaded in antiplane shear. The crack must be modeled as an arc of a circle and wholly inside the solidotherwise its position and orientation with respect to the boundary of the body is arbitrary.
The effect of the crack on the stress field is incorporated in an augmented kernel developed for the mode III crack problem such that discretization of the cutout boundary is no longer necessary. This modification of the kernel of the integral equation leads to solutions on and near the cutout with great accuracy. An asymptotic analysis is conducted in order to derive the Stress Intensity Factor (SIF) Kin, at each crack tip, in closed form. In this formulation, a straight crack can be viewed as a particular case of the more general curved crack. In particular, attention is paid,to the influence of crack curvature and edge effect on the stress intensity factors at the right and left crack tips. A rigorous mathematical formulation is developed, the main aspects of the numerical implementation are discussed and several representative numerical examples are presented in this paper.


Back to Top 