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 three-dimensional 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 density-based 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.

A robust ground structure topology optimization framework is presented to handle the uncertainty of load direction and design for the worst case compliance scenario. The deterministic optimization framework is formulated by a min-max compliance objective to first determine the critical load angle corresponding to the worst case compliance and then to design the topology for compliance minimization. The first optimization problem, based on our load definition, is shown to be equivalent to a maximum eigenvalue function, thus causing significant drawbacks in gradient-based optimization approaches in the case of eigenvalue coalescence. Here, we propose a method to treat the non-differentiability of the maximum eigenvalue optimization problem by a smooth maximum regularization function; hence, presenting a framework for optimizing ground structure networks considering an infinite number of load directions. The results achieved demonstrate that the proposed framework provides solutions with low compliance in all possible loading directions leading to robust structural designs.

The concepts of origami and kirigami have often been presented separately. Here, we put forth a synergistic approach—the folded kirigami—in which kirigami assemblies are complemented by means of folding, typical of origami patterns. Besides the emerging patterns themselves, the synergistic approach also leads to topological mechanical metamaterials. While kirigami metamaterials have been fabricated by various methods, such as 3D printing, cutting, casting, and assemblage of building blocks, the “folded kirigami” claim their distinctive properties from the universal folding protocols. For a target kirigami pattern, we design an extended high-genus pattern with appropriate sets of creases and cuts, and proceed to fold it sequentially to yield the cellular structure of a 2D lattice endowed with finite out-of-plane thickness. The strategy combines two features that are generally mutually exclusive in canonical methods: fabrication involving a single piece of material and realization of nearly ideal intercell hinges. We test the approach against a diverse portfolio of triangular and quadrilateral kirigami configurations. We demonstrate a plethora of emerging metamaterial functionalities, including topological phase-switching reconfigurability between polarized and nonpolarized states in kagome kirigami, and availability of nonreciprocal mechanical response in square-rhombus kirigami.

Origami-based designs refer to the application of the ancient art of origami to solve engineering problems of different nature. Despite being implemented at dimensions that range from the nano to the meter scale, origami-based designs are always defined by the laws that govern their geometrical properties at any scale. It is thus not surprising to notice that the study of their applications has become of cross-disciplinary interest. This article aims to review recent origami-based applications in engineering, design methods and tools, with a focus on research outcomes from 2015 to 2020. First, an introduction to origami history, mathematical background and terminology is given. Origami-based applications in engineering are reviewed largely in the following fields: biomedical engineering, architecture, robotics, space structures, biomimetic engineering, fold-cores, and metamaterials. Second, design methods, design tools, and related manufacturing constraints are discussed. Finally, the article concludes with open questions and future challenges.

Origami foldcores, especially the blockfold cores, have emerged as promising components of high-performance sandwich composites. Inspired by the blockfold origami, we propose the axisymmetric blockfold origami (ABO), which is composed of both rectangular and trapezoidal panels. The ABO inherits the non-flat-foldability of the blockfold origami, and furthermore, displays self-locking mechanisms and enhanced stiffness. The geometry and folding kinematics of the ABO are formulated with respect to the geometric parameters and the folding angle of the assembly. The mathematical conditions are derived for the existence of self-locking mechanisms. We perform compression test simulations to demonstrate enhanced stiffness and increased critical loads. We find that the existence of rectangular panels not only dominates the non-flat-foldability of the ABO, but also contributes to the enhancement of the stiffness. Our results suggest the potential applications of the ABO for building load-bearing structures with rotational symmetry. Moreover, we discuss the prospects of designing tightly assembled multi-layered origami structures with prestress induced by the mismatch of successive layers to enlighten future research.

Soft robots often draw inspiration from nature to navigate different environments. Although the inching motion and crawling motion of caterpillars have been widely studied in the design of soft robots, the steering motion with local bending control remains challenging. To address this challenge, we explore modular origami units which constitute building blocks for mimicking the segmented caterpillar body. Based on this concept, we report a modular soft Kresling origami crawling robot enabled by electrothermal actuation. A compact and lightweight Kresling structure is designed, fabricated, and characterized with integrated thermal bimorph actuators consisting of liquid crystal elastomer and polyimide layers. With the modular design and reprogrammable actuation, a multiunit caterpillar-inspired soft robot composed of both active units and passive units is developed for bidirectional locomotion and steering locomotion with precise curvature control. We demonstrate the modular design of the Kresling origami robot with an active robotic module picking up cargo and assembling with another robotic module to achieve a steering function. The concept of modular soft robots can provide insight into future soft robots that can grow, repair, and enhance functionality.

From a kinematics perspective, a Kresling origami cell couples axial displacement (contraction/expansion) with twist, leading to non-rigid origami behavior. From an energy landscape perspective, the selection of the Kresling origami geometry, together with its fabrication process and material, lead to energy envelopes allowing single or multiple stable states. In this context, this paper explores the Kresling origami mechanics through mathematical modeling integrated with experimental testing. On the theoretical mechanics front, we present a comprehensive model incorporating the representative geometrical parameters of the Kresling origami cell into the corresponding energy function in order to capture its essential mechanical behavior. On the experimental mechanics front, we create two fixtures that demonstrate the ability to control axial displacement (contraction/expansion) and twist independently, without imposing any constraints on the chiral arrangement of individual cells within the Kresling origami array (composed of multiple cells). Finally, we show the coexistence of multiple mechanical and morphological configurations within the same Kresling array by programming its loading modes, i.e., compression or twist. The fundamental nature of this work makes it applicable to several fields of engineering, including soft robotics and mechanical computing.

Origami provides a method to transform a flat surface into complex three dimensional geometries, which has applications in deployable structures, meta-materials, robotics, and beyond. The Miura-ori and the Eggbox are two fundamental planar origami patterns. Both patterns have been studied closely, and have become the basis for many engineering applications and derivative origami patterns. Here, we study the Hybrid structure formed by combining unit cells of the Miura-ori and Eggbox patterns. We find the compatibility constraints required to form the Hybrid structure and derive properties of its kinematics such as self-locking and Poisson's ratio. We then compare the aforementioned properties of the Miura-Eggbox Hybrid with those of the Morph pattern, another generalization of the Miura-ori and Eggbox patterns. In addition, we study the structure formed by combining all three unit cells of the Miura-ori, Eggbox, and Morph. Our results show that such patterns have tunable self-locking states and Poisson's ratio beyond their constituent components. Hybrid patterns formed by combining different origami patterns are an avenue to derive more functionality from simple constituents for engineering applications.

The development of strut-and-tie models (STM) for the design of reinforced concrete deep beams considering a general multi-material and multi-volume topology optimization framework is presented. The general framework provides flexibility to control location/inclination/length scale of the ties according to practical design requirements. Well-understood optimality conditions are applied to evaluate the performance of the optimized STM layouts. Specifically, the Michell number Z (or load path) is used as a simple and effective criterion to quantify the STM. The experimental results verify that the STM layout with the lowest load path Z achieves the highest ultimate load. Moreover, significantly reduced cracking is observed in the optimized STM layouts compared to the traditional STM layout. This observation implies that the optimized layouts may require less crack control reinforcement, which would lower the total volume of steel required for the deep beams.

Topology optimization problems typically consider a single load case or a small, discrete number of load cases; however, practical structures are often subjected to infinitely many load cases that may vary in intensity, location and/or direction (e.g. moving/rotating loads or uncertain fixed loads). The variability of these loads significantly influences the stress distribution in a structure and should be considered during the design. We propose a locally stress-constrained topology optimization formulation that considers loads with continuously varying direction to ensure structural integrity under more realistic loading conditions. The problem is solved using an Augmented Lagrangian method, and the continuous range of load directions is incorporated through a series of analytic expressions that enables the computation of the worst-case maximum stress over all possible load directions. Variable load intensity is also handled by controlling the magnitude of load basis vectors used to derive the worst-case load. Several two- and three-dimensional examples demonstrate that topology-optimized designs are extremely sensitive to loads that vary in direction. The designs generated by this formulation are safer, more reliable, and more suitable for real applications, because they consider realistic loading conditions.

Numerous topology optimization formulations have been proposed in order to enhance structural resistance to material failure. A clear line can often be drawn between those methods which attempt to constrain local failure criteria and those that explicitly model the failure physics during the optimization process. In this work the former method is extended in a manner inspired by the mathematical form of typical gradient-enhanced damage models. Importantly, the proposed formulation relies on linear physics during the optimization procedure, which greatly increases its speed and robustness, both of which are essential in industrial applications where large numerical models are typically used. The size effect introduced by using such a numerical model is further investigated and select observations are provided, such as spurious ‘‘fin-like’’ patterns that emerge depending on the type of structure and loading conditions. Finally, the load capacity is verified for each optimized design through a post-optimization verification procedure which is unaffected by the density-based design parameterization and associated material interpolation schemes.

The multi-degree of freedom Resch pattern forms most of its surfaces by tuning the folding angles of its creases. In this research, we program the triangular Resch pattern to naturally achieve surfaces with various curvatures by predefining the neutral angle and stiffness of the creases for the whole pattern. We simulate the free unfolding of tessellations by combining a “bar-and-hinge” model with an explicit meshless method, namely, the finite particle method. The effects of the damping factor, crease stiffness, and neutral angle on the unfolding process are investigated. The neutral angle of the creases plays a critical role in determining the final stable shape of the pattern. Then, we break the natural symmetry of the tessellations by changing the neutral angle and applying specific constraints to active creases to create stable unfolding surfaces with various curvatures. This study provides a foundation for the development of programmed curvatures of metamaterials that can be folded into origami patterns with multiple degrees of freedom.

The high computational cost of topology optimization has prevented its widespread use as a generative design tool. To reduce this computational cost, we propose an artificial intelligence approach to drastically accelerate topology optimization without sacrificing its accuracy. The resulting AI-driven topology optimization can fully capture the underlying physics of the problem. As a result, the machine learning model, which consists of a convolutional neural network with residual links, is able to generalize what it learned from the training set to solve a wide variety of problems with different geometries, boundary conditions, mesh sizes, volume fractions and filter radius. We train the machine learning model separately from the topology optimization, which allows us to achieve a considerable speedup (up to 30 times faster than traditional topology optimization). Through several design examples, we demonstrate that the proposed AI-driven topology optimization framework is effective, scalable and efficient. The speedup enabled by the framework makes topology optimization a more attractive tool for engineers in search of lighter and stronger structures, with the potential to revolutionize the engineering design process. Although this work focuses on compliance minimization problems, the proposed framework can be generalized to other objective functions, constraints and physics.

Geometrical-frustration-induced anisotropy and inhomogeneity are explored to achieve unique properties of metamaterials that set them apart from conventional materials. According to Neumann's principle, to achieve anisotropic responses, the material unit cell should possess less symmetry. Based on such guidelines, a triclinic metamaterial system of minimal symmetry is presented, which originates from a Trimorph origami pattern with a simple and insightful geometry: a basic unit cell with four tilted panels and four corresponding creases. The intrinsic geometry of the Trimorph origami, with its changing tilting angles, dictates a folding motion that varies the primitive vectors of the unit cell, couples the shear and normal strains of its extrinsic bulk, and leads to an unusual Poisson effect. Such an effect, associated with reversible auxeticity in the changing triclinic frame, is observed experimentally, and predicted theoretically by elegant mathematical formulae. The nonlinearities of the folding motions allow the unit cell to display three robust stable states, connected through snapping instabilities. When the tristable unit cells are tessellated, phenomena that resemble linear and point defects emerge as a result of geometric frustration. The frustration is reprogrammable into distinct stable and inhomogeneous states by arbitrarily selecting the location of a single or multiple point defects. The Trimorph origami demonstrates the possibility of creating origami metamaterials with symmetries that are hitherto nonexistent, leading to triclinic metamaterials with tunable anisotropy for potential applications such as wave propagation control and compliant microrobots.

Origami metamaterials are known to display highly tunable Poisson’s ratio depending on their folded state. Most studies on the Poisson effects in deployable origami tessellations are restricted to theory and simulation. Experimental realization of the desired Poisson effects in origami metamaterials requires special attention to the boundary conditions to enable nonlinear deformations that give rise to tunability. In this work, we present a novel experimental setup suitable to study the Poisson effects in 2D origami tessellations that undergo simultaneous deformations in both the applied and transverse directions. The setup comprises a gripping mechanism, which we call a Saint-Venant fixture, to eliminate Saint-Venant end effects during uniaxial testing. Using this setup, we conduct Poisson’s ratio measurements of the Morph origami pattern whose configuration space combines features of the Miura-ori and Eggbox parent patterns. We experimentally observe the Poisson’s ratio sign switching capability of the Morph pattern, along with its ability to display either completely positive or negative values of Poisson’s ratio by virtue of topological transformations. To demonstrate the versatility of the novel setup we also perform experiments on the standard Miura-ori and the standard Eggbox patterns. Our results demonstrate the agreement between the theory, the simulations, and the experiments on the Poisson’s ratio measurement and its tunability in origami metamaterials. The proposed experimental technique can be adopted for investigating other tunable properties of origami metamaterials in static and in dynamic regimes, such as elastic thermal expansion, and wave propagation control.

A system and method for accelerating topology optimization of a design includes a topology optimization module con figured to determine state variables of the topology using a two - scale topology optimization using design variables for a coarse - scale mesh and a fine - scale mesh for a number of optimization steps . A machine learning module includes a fully connected deep neural network having a tunable num ber of hidden layers configured to execute an initial training of a machine learning - based model using the history data , determine a predicted sensitivity value related to the design variables using the trained machine learning model , execute an online update of the machine learning - based model using updated history data , and update the design variables based on the predicted sensitivity value . The model predictions reduce the number of two - scale optimizations for each optimization step to occur only for initial training and for online model updates.

As a series of tessellation origami patterns consisting of more than one type of polygons, Resch patterns are generally rigid foldable but with a large number of degrees of freedom (DOFs). The patterns are not flat foldable, instead the fully folded configuration is a sandwich style in flat profile. In order to achieve one-DOF forms of triangular Resch pattern units, the thick-panel technique is employed at vertices to replace spherical linkages with spatial linkages. The compatibility among all the vertices is studied by a systematic kinematic analysis of the thickpanel Resch pattern. Finally, two design schemes are obtained to form a one-DOF origami structure between the initial planar structure and the fully folded 3D sandwich structure with dome-curved intermediate folded configurations. The one-DOF Resch pattern structure can be used to generate reliable shape transformation.

Geometric compatibility constraints dictate the mechanical response of soft systems that can be utilized for the design of mechanical metamaterials such as the negative Poisson’s ratio Miura-ori origami crease pattern. Here, we develop a formalism for linear compatibility that enables explicit investigation of the interplay between geometric symmetries and functionality in origami crease patterns. We apply this formalism to a particular class of periodic crease patterns with unit cells composed of four arbitrary parallelogram faces and establish that their mechanical response is characterized by an anticommuting symmetry. In particular, we show that the modes are eigenstates of this symmetry operator and that these modes are simultaneously diagonalizable with the symmetric strain operator and the antisymmetric curvature operator. This feature reveals that the anticommuting symmetry defines an equivalence class of crease pattern geometries that possess equal and opposite in-plane and out-of-plane Poisson’s ratios. Finally, we show that such Poisson’s ratios generically change sign as the crease pattern rigidly folds between degenerate ground states and we determine subfamilies that possess strictly negative in-plane or out-of-plane Poisson’s ratios throughout all configurations.

The present study introduces a formulation for topology optimization of structures with constraints on the frst principal stress. We solve the problem considering local stress constraints via the augmented Lagrangian method, which enables the solution of large-scale problems without the need for ad hoc aggregation schemes and clustering methods. Numerical examples are provided which demonstrate the efectiveness of the framework for practical problems with numerous (e.g., in the range of million(s)) local constraints imposed on the maximum principal stress. One of the examples is a three-dimensional antenna support bracket, which represents a realistic engineering design problem. This example, which has more than one million constraints, is proposed as a benchmark problem for stress-constrained topology optimization.

Biomimetic soft robotic crawlers have attracted extensive attention in various engineering fields, owing to their adaptivity to different terrains. Earthworm-like crawlers realize locomotion through in-plane contraction, while inchworm-like crawlers exhibit out-of-plane bending-based motions. Although in-plane contraction crawlers demonstrate effective motion in confined spaces, miniaturization is challenging because of limited actuation methods and complex structures. Here, we report a magnetically actuated small-scale origami crawler with in-plane contraction. The contraction mechanism is achieved through a four-unit Kresling origami assembly consisting of two Kresling dipoles with two-level symmetry. Magnetic actuation is used to provide appropriate torque distribution, enabling a small-scale and untethered robot with both crawling and steering capabilities. The crawler can overcome large resistances from severely confined spaces by its anisotropic and magnetically tunable structural stiffness. The multifunctionality of the crawler is explored by using the internal cavity of the crawler for drug storage and release. The magnetic origami crawler can potentially serve as a minimally invasive device for biomedical applications.

Spinodal architected materials with tunable anisotropy unify optimal design and manufacturing of multiscale structures. By locally varying the spinodal class, orientation, and porosity during topology optimization, a large portion of the anisotropic material space is exploited such that material is efficiently placed along principal stress trajectories at the microscale. Additionally, the bicontinuous, nonperiodic, unstructured, and stochastic nature of spinodal architected materials promotes mechanical and biological functions not explicitly considered during optimization (e.g., insensitivity to imperfections, fluid transport conduits). Furthermore, in contrast to laminated composites or periodic, structured architected materials (e.g., lattices), the functional representation of spinodal architected materials leads to multiscale, optimized designs with clear physical interpretation that can be manufactured directly, without special treatment at spinodal transitions. Physical models of the optimized, spinodal-embedded parts are manufactured using a scalable, voxel-based strategy to communicate with a masked stereolithography (m-SLA) 3D printer.

As a series of tessellation origami patterns consisting of more than one type of polygons, Resch patterns are generally rigid foldable but with a large number of degrees of freedom (DOFs). The patterns are not flat foldable, instead the fully folded configuration is a sandwich style in flat profile. In order to achieve one-DOF forms of triangular Resch pattern units, the thick-panel technique is employed at vertices to replace spherical linkages with spatial linkages. The compatibility among all the vertices is studied by a systematic kinematic analysis of the thickpanel Resch pattern. Finally, two design schemes are obtained to form a one-DOF origami structure between the initial planar structure and the fully folded 3D sandwich structure with dome-curved intermediate folded configurations. The one-DOF Resch pattern structure can be used to generate reliable shape transformation.

A wide range of modern science and engineering applications are formulated as optimization problems with a system of partial differential equations (PDEs) as constraints. These PDE-constrained optimization problems are typically solved in a standard discretize-then-optimize approach. In many industry applications that require high-resolution solutions, the discretized constraints can easily have millions or even billions of variables, making it very slow for the standard iterative optimizer to solve the exact gradients. In this work, we propose a general framework to speed up PDE-constrained optimization using online neural synthetic gradients (ONSG) with a novel twoscale optimization scheme. We successfully apply our ONSG framework to computational morphogenesis, a representative and challenging class of PDE-constrained optimization problems. Extensive experiments have demonstrated that our method can significantly speed up computational morphogenesis (also known as topology optimization), and meanwhile maintain the quality of final solution compared to the standard optimizer. On a large-scale 3D optimal design problem with around 1,400,000 design variables, our method achieves up to 7.5x speedup while producing optimized designs with comparable objectives.

We analyze the folding kinematics of a recently proposed origami-based tessellated structure called the Morph pattern, using thin, rigid panel assumptions. We discuss the geometry of the Morph unit cell that can exist in two characteristic modes differing in the mountain/valley assignment of a degree-four vertex and explain how a single tessellation of the Morph structure can undergo morphing through rigid origami kinematics resulting in multiple hybrid states. We describe the kinematics of the tessellated Morph pattern through multiple branches, each path leading to different sets of hybrid states. We study the kinematics of the tessellated structure through local and global Poisson’s ratios and derive an analytical condition for which the global ratio switches between negative and positive values. We show that the interplay between the local and global kinematics results in folding deformations in which the hybrid states are either locked in their current modes or are transformable to other modes of the kinematic branches, leading to a reprogrammable morphing behavior of the system. Finally, using a bar-and-hinge model-based numerical framework, we simulate the nonlinear folding behavior of the hybrid systems and verify the deformation characteristics that are predicted analytically.

We present an augmented Lagrangian-based approach for stress-constrained topology optimization of structures subjected to general dynamic loading. The approach renders structures that satisfy the stress constraints locally at every time step. To solve problems with a large number of stress constraints, we normalize the penalty term of the augmented Lagrangian function with respect to the total number of constraints (i.e., the number of elements in the mesh times the number of time steps). Moreover, we solve the stress-constrained problem effectively by penalizing constraints associated with high stress values more severely than those associated with low stress values. We integrate the equations of motion using the HHTmethod and conduct the sensitivity analysis consistently with this method via the “discretize-then-differentiate” approach. We present several numerical examples that elucidate the effectiveness of the approach to solve dynamic, stress-constrained problems under several loading scenarios including loads that change in magnitude and/or direction and loads that change in position as a function of time.

An optimization-based approach is proposed to design elastostatic cloaking devices in twodimensional (2D) lattices. Given an elastic lattice with a defect, i.e. a circular or elliptical hole, a small region (cloak) around the hole is designed to hide the effect of the hole on the elastostatic response of the lattice. Inspired by the direct lattice transformation approach to elastostatic cloaking in 2D lattices, the lattice nodal positions in the design region are obtained using a coordinate transformation of the reference (undisturbed) lattice nodes. Subsequently, additional connectivity (i.e. a ground structure) is defined in the design region and the stiffness properties of these elements are optimized to mimic the global stiffness characteristics of the reference lattice. A weighted least-squares objective function is proposed, where the weights have a physical interpretation—they are the design-dependent coefficients of the design lattice stiffness matrix. The formulation leads to a convex objective function that does not require a solution to an additional adjoint system. Optimization-based cloaks are designed considering uniaxial tension in multiple directions and are shown to exhibit approximate elastostatic cloaking, not only when subjected to the boundary conditions they were designed for but also for uniaxial tension in directions not used in design and for shear loading.

Inspired by the embodied intelligence observed in octopus arms, we introduce magnetically controlled origami robotic arms based on Kresling patterns for multimodal deformations, including stretching, folding, omnidirectional bending, and twisting. The highly integrated motion of the robotic arms is attributed to inherent features of the reconfigurable Kresling unit, whose controllable bistable deploying/folding and omnidirectional bending are achieved through precise magnetic actuation. We investigate single- and multiple-unit robotic systems, the latter exhibiting higher biomimetic resemblance to octopus’ arms. We start from the single Kresling unit to delineate the working mechanism of the magnetic actuation for deploying/folding and bending. The two-unit Kresling assembly demonstrates the basic integrated motion that combines omnidirectional bending with deploying. The four-unit Kresling assembly constitutes a robotic arm with a larger omnidirectional bending angle and stretchability. With the foundation of the basic integrated motion, scalability of Kresling assemblies is demonstrated through distributed magnetic actuation of double-digit number of units, which enables robotic arms with sophisticated motions, such as continuous stretching and contracting, reconfigurable bending, and multiaxis twisting. Such complex motions allow for functions mimicking octopus arms that grasp and manipulate objects. The Kresling robotic arm with noncontact actuation provides a distinctive mechanism for applications that require synergistic robotic motions for navigation, sensing, and interaction with objects in environments with limited or constrained access. Based on small-scale Kresling robotic arms, miniaturized medical devices, such as tubes and catheters, can be developed in conjunction with endoscopy, intubation, and catheterization procedures using functionalities of object manipulation and motion under remote control.

Structural instability, once a catastrophic phenomenon to be avoided in engineering applications, is being harnessed to improve functionality of structures and materials, and has been a catalyst of substantial research in the field. One important application is to create functional metamaterials that deform their internal structure to adjust performance, resembling phase transformations in natural materials. In this paper, we propose a novel origami pattern, named the Shrimp pattern, with application to multi-phase architected metamaterials whose phase transition is achieved mechanically by snap-through. The Shrimp pattern consists of units that can be easily tessellated in two dimensions, either periodically with homogeneous local geometry or non-periodically with heterogeneous local geometries. We can use a few design parameters to program the unit cell to become either monostable or tune the energy barrier between the bistable states. By tessellating these unit cells into an architected metamaterial, we can create complex yet navigable energy landscapes, leading to multiple metastable phases of the material. As each phase has different geometries, the metamaterial can switch between different mechanical properties and shapes. The geometric origin of the multi-stable behavior implies that, conceptually, our designs are scale-independent, making them candidates for a variety of innovative applications, including reprogramable materials, reconfigurable acoustic waveguides, and microelectronic mechanical systems and energy storage systems.

We present a Matlab implementation for topology optimization of structures subjected to dynamic loads. The code, which we name PolyDyna, is built on top of PolyTop—a Matlab code for static compliance minimization based on polygonal finite elements. To solve the structural dynamics problem, we use the HHT-α method, which is a generalization of the classical Newmark-β method. In order to handle multiple regional volume constraints efficiently, PolyDyna uses a variation of the ZPR design variable update scheme enhanced by a sensitivity separation technique, which enables it to solve non-self-adjoint topology optimization problems. We conduct the sensitivity analysis using the adjoint method with the “discretize-then-differentiate” approach, such that the sensitivity analysis is consistently evaluated on the discretized system (both in space and time). We present several numerical examples, which are explained in detail and summarized in a library of benchmark problems. PolyDyna is intended for educational purposes and the complete Matlab code is provided as electronic supplementary material.

As a ubiquitous paper folding art, origami has promising applications in science and engineering. Many software and parameterized methods have been proposed to draw, analyze and design origami patterns. Here we focus on the shape grammar formalism and the Shape Machine, a shape grammar interpreter that has managed to automate the seamless shape calculations that the shape grammar formalism advocates. Different from other origami pattern generation methods, shape grammars generate origami patterns through recursive applications of shape rewriting rules using lines and curves. Based on this concept, the transformations between some common origami patterns are reorganized following visual cues and reasoning. Four examples of generating origami pattern are presented to show the capability of Shape Machine in origami design, including construction, modification and programming of an origami pattern. The new origami designs inspired by this work prove that shape grammars and Shape Machine provide a perspective and modeling technique for creating origami tessellation patterns.

We present a B-bar formulation of the virtual element method (VEM) for the analysis of both nearly incompressible and compressible materials. The material stiffness is decomposed into dilatational and deviatoric parts, and only the deviatoric part of the material stiffness is utilized for stabilization of the element stiffness matrix. A feature of the formulation is that locking behavior for nearly incompressible materials is successfully removed by the spectral decomposition of the material stiffness. The eigenvalue analysis demonstrates that the method eliminates higher energy modes associated with locking behavior for nearly incompressible materials, while capturing constant strain energy modes for both compressible and nearly incompressible materials. The convergence and accuracy of the B-bar VEM are discussed using 2D and 3D examples with various element shapes (convex and non-convex).

Because of increased geometric freedom at a widening range of length scales and access to a growing material space, additive manufacturing has spurred renewed interest in topology optimization of parts with spatially varying material properties and structural hierarchy. Simultaneously, a surge of micro/nanoarchitected materials have been demonstrated. Nevertheless, multiscale design and micro/nanoscale additive manufacturing have yet to be sufficiently integrated to achieve free-form, multiscale, biomimetic structures. We unify design and manufacturing of spatially varying, hierarchical structures through a multimicrostructure topology optimization formulation with continuous multimicrostructure embedding. The approach leads to an optimized layout of multiple microstructural materials within an optimized macrostructure geometry, manufactured with continuously graded interfaces. To make the process modular and controllable and to avoid prohibitively expensive surface representations, we embed the microstructures directly into the 3D printer slices. The ideas provide a critical, interdisciplinary link at the convergence of material and structure in optimal design and manufacturing.

We present PolyStress, a Matlab implementation for topology optimization with local stress constraints considering linear and material nonlinear problems. The implementation of PolyStress is built upon PolyTop, an educational code for compliance minimization on unstructured polygonal finite elements. To solve the nonlinear elasticity problem, we implement a Newton-Raphson scheme, which can handle nonlinear material models with a given strain energy density function. To solve the stress-constrained problem, we adopt a scheme based on the augmented Lagrangian method, which treats the problem consistently with the local definition of stress without employing traditional constraint aggregation techniques. The paper discusses several theoretical aspects of the stress-constrained problem, including details of the augmented Lagrangian-based approach implemented herein. In addition, the paper presents details of the Matlab implementation of PolyStress, which is provided as electronic supplementary material. We present several numerical examples to demonstrate the capabilities of PolyStress to solve stress-constrained topology optimization problems and to illustrate its modularity to accommodate any nonlinear material model. Six appendices supplement the paper. In particular, the first appendix presents a library of benchmark examples, which are described in detail and can be explored beyond the scope of the present work.

We put forward a general machine learning-based topology optimization framework, which greatly accelerates the design process of large-scale problems, without sacrifice in accuracy. The proposed framework has three distinguishing features. First, a novel online training concept is established using data from earlier iterations of the topology optimization process. Thus, the training is done during, rather than before, the topology optimization. Second, a tailored two-scale topology optimization formulation is adopted, which introduces a localized online training strategy. This training strategy can improve both the scalability and accuracy of the proposed framework. Third, an online updating scheme is synergistically incorporated, which continuously improves the prediction accuracy of the machine learning models by providing new data generated from actual physical simulations. Through numerical investigations and design examples, we demonstrate that the aforementioned framework is highly scalable and can efficiently handle design problems with a wide range of discretization levels, different load and boundary conditions, and various design considerations (e.g., the presence of non-designable regions).

This paper examines the structural behavior of a steel catenary riser (SCR) with a segment composed of a functionally graded material (FGM) at the touchdown zone (TDZ). A co-rotational beam formulation is employed to evaluate the influence of material gradation on the static and dynamic nonlinear behavior of the riser. A finite element model is developed to predict the riser response considering environmental and operational conditions such as self-weight, buoyancy, current, soil contact, waves, buoys, and ship motion. Moreover, a comparison of SCR configurations is made between those constructed from FGM and those constructed from standard homogeneous materials. The study indicates that an SCR configuration with an incorporated FGM segment improves the curvatures near the TDZ, where designing to mitigate fatigue damage is critical.

To consistently coarsen arbitrary unstructured meshes, a computational morphogenesis process is built in conjunction with a numerical method of choice, such as the virtual element method with adaptive meshing. The morphogenesis procedure is performed by clustering elements based on a posteriori error estimation. Additionally, an edge straightening scheme is introduced to reduce the number of nodes and improve accuracy of solutions. The adaptive morphogenesis can be recursively conducted regardless of element type and mesh generation counting. To handle mesh modification events during the morphogenesis, a topology‐based data structure is employed, which provides adjacent information on unstructured meshes. Numerical results demonstrate that the adaptive mesh morphogenesis effectively handles mesh coarsening for arbitrarily shaped elements while capturing problematic regions such as those with sharp gradients or singularity.

Origami structures demonstrate great theoretical potential for creating metamaterials with exotic properties. However, there is a lack of understanding of how imperfections influence the mechanical behaviour of origami-based metamaterials, which, in practice, are inevitable. For conventional materials, imperfection plays a profound role in shaping their behaviour. Thus, this paper investigates the influence of small random geometric imperfections on the nonlinear compressive response of the representative Miura-ori, which serves as the basic pattern for many metamaterial designs. Experiments and numerical simulations are used to demonstrate quantitatively how geometric imperfections hinder the foldability of the Miura-ori, but on the other hand, increase its compressive stiffness. This leads to the discovery that the residual of an origami foldability constraint, given by the Kawasaki theorem, correlates with the increase of stiffness of imperfect origami-based metamaterials. This observation might be generalizable to other flat-foldable patterns, in which we address deviations from the zero residual of the perfect pattern; and to non-flat-foldable patterns, in which we would address deviations from a finite residual.

We address material nonlinear topology optimization problems considering the Drucker–Prager strength criterion by means of a surrogate nonlinear elastic model. The nonlinear material model is based on a generalized J2 deformation theory of plasticity. From an algorithmic viewpoint, we consider the topology optimization problem subjected to prescribed energy, which leads to robust convergence in nonlinear problems. The objective function of the optimization problem consists of maximizing the strain energy of the system in equilibrium subjected to a volume constraint. 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 through direct minimization of the structural strain energy using Newton’s method with an inexact line search strategy. Four numerical examples demonstrate features of the proposed nonlinear topology optimization framework considering the Drucker–Prager strength criterion.

In order to capture nonlinear crack tip behavior and account for mixed-mode fatigue crack growth, a cohesive zone based fatigue crack growth model is proposed in conjunction with the Park-Paulino-Roesler (PPR) traction-separation relationship. The model clearly defines five stages during arbitrary fatigue loading: softening, unloading, reloading, contact, and complete failure. The cohesive traction-separation relationship is based on the PPR fracture potential, while the fatigue damage is accumulated by introducing two conjugate damage measures. One damage measure is associated with the rate of separation, while the other is related to the rate of traction (or local stress). Additionally, two model constants are introduced to control the normal contact condition, which may be associated with physical conditions such as crack closure, crack face roughness, and oxidation of fracture surface. Furthermore, computational simulations of the proposed fatigue crack growth model are performed for a simple mode-I test, double cantilever beam test, modified mixed-mode bending test, and three-point bending test. The observed computational results lead to stable and consistent fatigue crack growth.

We introduce a general multi-material topology optimization framework for large deformation problems that effectively handles an arbitrary number of candidate hyperelastic materials and addresses three major associated challenges: material interpolation, excessive distortion of low-density elements, and computational efficiency. To account for many nonlinear elastic materials, we propose a material interpolation scheme that, instead of interpolating multiple material parameters (such as Young’s modulus), interpolates multiple nonlinear stored-energy functions. To circumvent convergence difficulties caused by excessive distortions of low-density elements under large deformations, an energy interpolation scheme is revisited to account for multiple candidate hyperelastic materials. Computational efficiency is addressed from both structural analysis and optimization perspectives. To solve the nonlinear state equations efficiently, we employ the lower-order Virtual Element Method in conjunction with tailored adaptive mesh refinement and coarsening strategies. To efficiently update the design variables of the multi-material system, we exploit the separable nature and improve the ZPR (Zhang–Paulino–Ramos) update scheme to account for positive sensitivities and update the design variables associated with each volume constraint in parallel. Four design examples with three types of nonlinear material models demonstrate the efficiency and effectiveness of the proposed framework.

Deployability, multifunctionality, and tunability are features that can be explored in the design space of origami engineering solutions. These features arise from the shape-changing capabilities of origami assemblies, which require effective actuation for full functionality. Current actuation strategies rely on either slow or tethered or bulky actuators (or a combination). To broaden applications of origami designs, we introduce an origami system with magnetic control. We couple the geometrical and mechanical properties of the bistable Kresling pattern with a magnetically responsive material to achieve untethered and local/distributed actuation with controllable speed, which can be as fast as a tenth of a second with instantaneous shape locking. We show how this strategy facilitates multimodal actuation of the multicell assemblies, in which any unit cell can be independently folded and deployed, allowing for on-the-fly programmability. In addition, we demonstrate how the Kresling assembly can serve as a basis for tunable physical properties and for digital computing. The magnetic origami systems are applicable to origami-inspired robots, morphing structures and devices, metamaterials, and multifunctional devices with multiphysics responses.

This paper presents a consistent topology optimization formulation for mass minimization with local stress constraints by means of the augmented Lagrangian method. To solve problems with a large number of constraints in an effective way, we modify both the penalty and objective function terms of the augmented Lagrangian function. The modification of the penalty term leads to consistent solutions under mesh refinement and that of the objective function term drives the mass minimization towards black and white solutions. In addition, we introduce a piecewise vanishing constraint, which leads to results that outperform those obtained using relaxed stress constraints. Although maintaining the local nature of stress requires a large number of stress constraints, the formulation presented here requires only one adjoint vector, which results in an efficient sensitivity evaluation. Several 2D and 3D topology optimization problems, each with a large number of local stress constraints, are provided.

The microstructural configuration of a material affects its macroscopic viscoelastic behavior, which suggests that materials can be engineered to achieve a desired viscoelastic behavior over a range of frequencies. To this end, we leverage topology optimization to find the optimized topology of a multi-phase viscoelastic composite to tailor its energy dissipation behavior as a function of frequency. To characterize the behavior of each material phase, we use a fractional viscoelastic constitutive model. This type of material model uses differential operators of non-integer order, which are appropriate to represent hereditary phenomena with long- and short-term memory. The topology optimization formulation aims to find the lightest microstructure that minimizes the sum of squared loss modulus residuals for a given set of target frequencies. This leads to the design of materials with either maximized loss modulus for a given target frequency or tailored loss modulus for a predefined set of frequencies. We present several numerical examples, both in 2D and 3D, which demonstrate that the microstructural configuration of multi-phase materials affects its macroscopic viscoelastic behavior. Thus, if properly designed, the material behavior can be tailored to dissipate energy for a desired frequency (maximized loss modulus) or for a range of frequencies (tailored energy dissipation behavior).

We present a virtual element method (VEM)-based topology optimization framework using polyhedral elements, which allows for convenient handling of non-Cartesian 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. In the structural problem, the VEM is adopted to solve the three-dimensional elasticity equation. Compared to the finite element method, the VEM does not require numerical integration (when linear elements are used) 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 the VEM basis functions. When compared to the standard element-wise constant approximation, the continuous approximation enriches the geometrical representation 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 non-Cartesian domains, demonstrating the main features of the proposed VEM-based topology optimization framework. The source code for a MATLAB implementation of the proposed work, named PolyTop3D, is available in the (electronic) Supplementary Material accompanying this publication.

Mechanical metamaterials inspired by the Japanese art of paper folding have gained considerable attention because of their potential to yield deployable and highly tunable assemblies. The inherent foldability of origami structures enlarges the material design space with remarkable properties such as auxeticity and high deformation recoverability and deployability, the latter being key in applications where spatial constraints are pivotal. This work integrates the results of the design, 3D direct laser writing fabrication, and in situ scanning electron microscopic mechanical characterization of microscale origami metamaterials, based on the multimodal assembly of Miura-Ori tubes. The origami-architected metamaterials, achieved by means of microfabrication, display remarkable mechanical properties: stiffness and Poisson’s ratio tunable anisotropy, large degree of shape recoverability, multistability, and even reversible auxeticity whereby the metamaterial switches Poisson’s ratio sign during deformation. The findings here reported underscore the scalable and multifunctional nature of origami designs, and pave the way toward harnessing the power of origami engineering at small scales.

Structures containing tension-only members, i.e., cables, are widely used in engineered structures (e.g., suspension and cable-stayed bridges, tents, and bicycle wheels) and are also found in nature (e.g., spider webs). We seek to use the ground structure method to obtain optimal cable network configurations. The structures are modeled using principles of nonlinear elasticity that allow for large displacements, i.e., global configuration changes, and large deformations. The material is characterized by a hyperelastic constitutive relation in which the strain energy is nonzero only when the axial stretch of a member is greater than or equal to one (i.e., tension-only behavior). We maximize the stationary potential energy of the equilibrated system, which avoids the need for an additional adjoint equation in computing the derivatives needed for the solution of the optimization problem. Several examples demonstrate the capabilities of the proposed formulation for topology optimization of cable networks. Motivated by nature, a spider web–inspired cable net is designed.

An interesting, yet challenging problem in topology optimization consists of finding the lightest structure that is able to withstand a given set of applied loads without experiencing local material failure. Most studies consider material failure via the von Mises criterion, which is designed for ductile materials. To extend the range of applications to structures made of a variety of different materials, we introduce a unified yield function that is able to represent several classical failure criteria including von Mises, Drucker–Prager, Tresca, Mohr–Coulomb, Bresler–Pister and Willam–Warnke, and use it to solve topology optimization problems with local stress constraints. The unified yield function not only represents the classical criteria, but also provides a smooth representation of the Tresca and the Mohr–Coulomb criteria—an attribute that is desired when using gradient-based optimization algorithms. The present framework has been built so that it can be extended to failure criteria other than the ones addressed in this investigation. We present numerical examples to illustrate how the unified yield function can be used to obtain different designs, under prescribed loading or design-dependent loading (e.g. self-weight), depending on the chosen failure criterion.

We present a simple, effective, and scalable approach for significantly accelerating the convergence in Topology Optimization simulations. Specifically, treating the design process as a fixed-point 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.

This paper presents a density-based topology optimization formulation for the design of multi-material thermoelastic structures. The formulation is written in the form of a multi-objective 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 non-self-adjoint topology optimization problems with an arbitrary number of volume constraints, which can be imposed either to a subset of the candidate materials, or to sub-regions 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.

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.

When conventional filtering schemes are used in reliability-based 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 potential-energy-based 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 reliability-based truss sizing optimization for the purpose of topology optimization under uncertainties. Furthermore, a single-loop approach is adopted to enhance the computational efficiency of the proposed RBTO method. Numerical examples of two- and three-dimensional 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 first-order reliability method and Monte Carlo simulations.

Origami offers an avenue to program three-dimensional shapes via scale-independent and non-destructive 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 non-uniform (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 multi-stable metasurfaces with programmable non-Euclidean geometries.

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 non-convex discretizations of elastodynamic problems in 2D and 3D using an explicit time integration scheme. In the formulation, a diagonal matrix-based 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 non-simply connected elements, which can lead to convenient modeling of arbitrarily-shaped inclusions in composites.

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 Class-1 (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.

Rate-dependent 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 mixed-mode rate-dependent 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 non-integer order, leading to power-law-type 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 self-healing. 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 rate-dependent fracture process of rubber-like materials for a wide range of loading rates and that it can capture rate-dependent mixed-mode fracture processes accurately. Results from the last example indicate that the order of the fractional differential operators acts as a memory-like parameter that allows for the fracture modeling of long- and short-term 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.

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 line-search strategy. Four numerical examples demonstrate the features of the proposed material nonlinear topology optimization framework for approximating standard von Mises plasticity.

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 well-known tensegrities and new tensegrities from the ground structure network, some of which are prototyped with the aid of additive manufacturing.

Exploring the configurational space of specific origami patterns [e.g., Miura-ori (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 four-vertex 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 in-plane stretching and out-of-plane 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.

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.

We investigate and experimentally observe the existence of topologically protected interface modes in a one-dimensional mechanical lattice and we report on the effect of nonlinearities on topological protection. The lattice consists of a one-dimensional array of spinners with nearest-neighbor 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 edge-to-bulk 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.

This paper introduces a general recovery-based 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.

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 first-passage probability, often needs to be evaluated. In this paper, a new method is proposed to incorporate constraints on the first-passage probability into reliability-based optimization of structural design or topology. For efficient evaluations of first-passage 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 first-passage probability to facilitate the use of efficient gradient-based optimization algorithms. The proposed method is successfully demonstrated by numerical examples of a space truss and building structures subjected to stochastic earthquake ground motions.

This paper presents a set of deployable origami tube structures that can create smooth functional surfaces while simultaneously maintaining a high out-of-plane 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. Three-point 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 outof-plane 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.

The tremendous increase in the number of components in typical electrical and communication modules requires low-cost, 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 continuous-state ranges from folded (compact) to unfolded (large surface) configurations. Special “bridge-like” 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 Miura-Ori pattern is used to fabricate spatial filters—frequency-selective surfaces (FSSs) with dipole resonant elements placed along the fold lines. The electrical length of the dipole elements in these structures changes when the Miura-Ori is folded, which facilitates tunable frequency response for the proposed shape-reconfigurable FSS structure. Higherorder spatial filters are realized by creating multilayer Miura-FSS configurations, which further increase the overall bandwidth of the structure. Such multilayer Miura-FSS structures feature the unprecedented capability of on-the-fly 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.

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 multi-material 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 non-convex 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 Zhang-Paulino-Ramos 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.

This paper contributes a scheme for two-dimensional (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 Park-Paulino-Roesler potential-based 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 I-II, 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 FEM-LB coupling demonstrate the potential of this numerical procedure to model hydraulic fracturing processes in geomaterials of complex geometries.

Origami engineering principles have recently been applied to a wide range of applications, including soft robots, stretchable electronics, and mechanical metamaterials. In order to achieve the 3D nature of engineered structures (e.g. load-bearing capacity) and capture the desired kinematics (e.g., foldability), many origami-inspired engineering designs are assembled from smaller parts and often require binding agents or additional elements for connection. Attempts at direct fabrication of 3D origami structures have been limited by available fabrication technologies and materials. Here, we propose a new method to directly 3D print origami assemblages (that mimic the behavior of their paper counterparts) with acceptable strength and load-bearing capacity for engineering applications. Our approach introduces hinge-panel elements, where the hinge regions are designed with finite thickness and length. The geometrical design of these hinge-panels, informed by both experimental and theoretical analysis, provides the desired mechanical behavior. In order to ensure foldability and repeatability, a novel photocurable elastomer system is developed and the designs are fabricated using digital light processing-based 3D printing technology. Various origami assemblages are produced to demonstrate the design flexibility and fabrication efficiency offered by our 3D printing method for origami structures with enhanced load bearing capacity and selective deformation modes.

This paper considers a co-rotational 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 two-node Hermitian beam elements which are referred to a co-rotation coordinate system attached to the element local frame of coordinates. The geometric non-linear effects of the beam are considered under large displacement and rotations, but under small-strain conditions. The deflections of the riser result from forces caused by self-weight, buoyancy, sea currents, waves, the action of floaters, seabed-structure 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.

A framework is presented for multi-material compliance minimization in the context of continuum based topology optimization. We adopt the common approach of finding 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 sufficiently general volume/mass constraints can be specified, 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 sub-region 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.

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 bar-and-hinge 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 Miura-ori 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.

Multi-material 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 cross-sectional 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 multi-material formulation that handles an arbitrary number of volume constraints and candidate materials. To efficiently handle such volume constraints, we employ the ZPR (Zhang-Paulino-Ramos) design update scheme for multi-material 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 Karush-Kuhn-Tucker (KKT) conditions. Through numerical examples, we demonstrate that the proposed algorithm, which can be readily implemented in multi-material optimization, along with the ZPR update scheme, is robust and effective for selecting a single preferred material among multiple candidate materials.

The structural performance of a grid-shell depends directly on the geometry of the design. Form-finding methods, which are typically based on the search for bending-free configurations, aid in achieving structurally efficient geometries. This manuscript proposes two form-finding methods for grid-shells: 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 well-known force density method, are evaluated and compared using numerical examples. The advantages and drawbacks of the methods are reviewed, compared and highlighted.

Rate-dependent 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 rate-dependent fracture process. The latter approach better represents the physics driving the rate-dependent fracture process, yet it is computationally more expensive. This work presents a methodology for studying mixed-mode rate-dependent fracture that is both efficient and approximates the viscoelastic material behavior in front of the crack tip. In this mixed-mode approach, we approximate the viscoelastic behavior in front of the crack tip using two rate-dependent 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 rate-dependent springs using two damage parameters extracted from the Park-Paulino-Roesler (PPR) cohesive fracture model. The rate-dependent 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.

There are few methods capable of capturing the full spectrum of pervasive fracture behavior in three-dimensions. 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 three-dimensional 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.

Large facial bone loss usually requires patient-specific bone implants to restore the structural integrity and functionality that also affects the appearance of each patient. Titanium alloys (e.g., Ti-6Al-4V) 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 stress-shielding. 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 well-known 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..

Multi-material 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 multi-material 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 Ogden-based, bilinear, and linear materials, we demonstrate that the proposed multi-material 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 (Zhang-Paulino-Ramos Jr.).

Origami-inspired designs possess attractive applications to science and engineering (e.g. deployable, self-assembling, 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 non-rigid origami structures. In this paper, we propose a fully nonlinear, displacement-based implicit formulation for performing static/quasi-static analyses of nonrigid origami structures based on ‘bar-and-hinge’ 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.

The ground structure method seeks to approximate Michell optimal solutions for real-world 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 post-processing 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 end-filter in the traditional plastic and nested elastic approaches of the ground structure method.

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.

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.

Tensegrity structures with detached struts are naturally suitable for deployable applications, both in terrestrial and outer-space structures, as well as morphing devices. Composed of discontinuous struts and continuous cables, such systems are only structurally stable when self-stress 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 shape-change 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 component-wise 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 light-weight structures, shape-changing soft robots, morphing antenna and RF devices, and biomedical devices.

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 reduced-order modeling, and demonstrate its application to two- and three-dimensional 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 Newton-Raphson approach with line search, which is essential for convergence. In addition, the use of reduced-order 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.

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 trial-and-error to obtain local fracture properties such that the DEMmodels approximate the experimental load-crack 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 single-edge notched beam test, SE(B), under three-point 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 fracture-related 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 Newton-like 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.

We present a general virtual element method (VEM) framework for finite elasticity, which emphasizes two issues: element-level 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 two-field mixed and an equivalent displacement-based, which is free of volumetric locking. Convergence and accuracy of the VEM are verified by means on numerical examples, and engineering applications are demonstrated.

Although the Park-Paulino-Roesler (PPR) potenital-based 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.

The definition of a traction-separation 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 traction-separation 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 traction-separation relationship in Abaqus can lead to non-physical 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 traction-separation relationship should be developed and selected with great caution, especially under mixed-mode conditions.

We investigate a recently proposed variational principle with rigid-body 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 rigid-body 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.

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 Editor-in-Chief with highest dedication until his death. He had a highly productive and long research life.

Despite being an effective and a general method to obtain optimal solutions, topology optimization generates solutions with complex geometries, which are neither cost-effective nor practical from a manufacturing (industrial) perspective. Manufacturing constraint techniques based on a unified projection-based 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 pseudo-density 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 projection-based approach.

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 well-established 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 low-order 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.

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, Voronoi-based and structured quadrilateral discretizations, and proposes two simple and effective ground structure generation approaches: the macro-element and macro-patch 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.

This paper proposes an efficient gradient-based optimization approach for reliability-based 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 first-order 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 multi-point 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 FORM-based approaches. The present approach typically does not lead to a fully stressed design, and thus this feature will be verified by numerical examples.

We use the graphical processing unit (GPU) to perform dynamic fracture simulation using adaptively refined and coarsened finite elements and the inter-element 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 reduced-scale 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.

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, patient-specific bone replacement shapes can be designed for different clinical cases that ensure revival of efficient load transfer mechanisms in the mid-face. In this work, four clinical cases are introduced and their respective patient-specific 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 full-field 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 mid-face reconstruction.

This paper proposes an efficient method named segmental multi-point 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 gradient-based routines for solving reliability-based design optimization problems where accuracy of the parameter sensitivity is essential for convergence to an optimal solution.

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 well-defined 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.

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 Voronoi-type 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, non-convergent 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 physically-based examples including a study of filled elastomers with interphasial effect and a qualitative comparison with cavitation experiments for fiber reinforced elastomers.

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 non-matching 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.

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 density-based topology optimization. The results obtained from these topology optimization methods require some degree of post-processing 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 three-dimensional density-based topology optimization, a tool suitable for result inspection and generation of additive manufacturing output is also provided.

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 single-field and mixed finite elements of arbitrary order. For demonstration purposes, applications of the proposed rigid-body-constraint formulation are illustrated within the context of elastomers, reinforced with periodic and random distributions of rigid filler particles, undergoing finite deformations.

Thin sheets can be assembled into origami tubes to create a variety of deployable, reconfigurable and mechanistically unique three-dimensional structures. We introduce and explore origami tubes with polygonal, translational symmetric cross-sections 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 cross-section 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 two-scale 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 new concrete fracture geometry is presented, which can quantify multiple fracture properties from a single specimen test. The disk-shaped 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 crack-tip opening displacement, CTOD_c, and initial fracture energy, G_f, as well as the specimen-dependent total fracture energy, G_F. The DCT-based fracture properties have the same experimental variation as the single-edge notched beam test. The experimentally derived fracture parameters were implemented into a cohesive zone model, which enabled estimation of concrete tensile strength from field-extracted cores.

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.

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 large-scale problems, and uses its object-oriented 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 ‘matrix-free’ 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.

This paper presents an integrated theoretical and computational investigation into the macroscopic behavior of composite materials containing multi-phase reinforcing particles with simultaneous nonlinear debonding along the microconstituent interfaces. The interfacial debonding is characterized by the nonlinear Park-Paulino-Roesler potential-based cohesive zone model. The extended Mori-Tanaka 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 three-dimensional finite element-based 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, snap-back and softening stages.

This paper describes an integrated topology optimization framework using discrete and continuum elements for buckling and stiffness optimization of high-rise buildings. The discrete (beam/truss) elements are optimized based on their cross-sectional 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 high-rise building design.

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 zipper-coupled 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.

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/herringbone-base folded sheet metamaterials such as Miura-ori. 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 Miura-ori, 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 in-plane Poisson’s ratios. By introducing a class of zigzag-base materials in the current study, we unify the concept of in-plane Poisson’s ratio for similar materials in the literature and extend it to the class of zigzagbase folded sheet materials.

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 pre-processing 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 ‘on-the-fly’, when needed and where needed. The main challenge for implementing a GPU-based 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 pre-processing 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.

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 carbon-black reinforced rubber, an interphase or ``bound rubber'' phase often occurs at the particle-matrix 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 Park-Paulino-Roesler traction-separation relation, are inserted between each particle and its corresponding interphase to account for debonding. To account for friction, we present a new, coupled cohesive-friction 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.

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 two-dimensions.

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 three-dimensional 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 first-order 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.

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 fine-scale 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 non-Gaussian 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 finite-element simulations of the large deformation response of Gaussian and non-Gaussian 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.

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 large-size, 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 super-component representing the set are computed respectively using the CSP method and combined by the chain-rule. 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 first-passage probability of a building structure under stochastic excitations, modeled by use of finite elements.

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 non-convergent. 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 pseudo-derivative method developed by Belytschko and co-workers and, for linear elliptic problems, with our previous remedy that involves splitting of polynomial and non-polynomial 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. 

Failure of adhesive bonded structures often occurs concurrent with the formation of a non-negligible 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 traction-separation relation. Yet it is usually determined empirically, by using calibration procedures combining experimental data, such as load-displacement or crack length data, with finite element simulation of fracture. Thanks to the recent progress in image processing, and the availability of low-cost 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 so-called 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.

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, Saint-Venant, 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.

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 ill-posedness 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 two-metric 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.

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 non-convex stored-energy 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 displacement-based and mixed polygonal finite element approximations, the latter of which consist of a piecewise constant pressure field and a linearly-complete 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 physically-based examples illustrate the potential of polygonal finite elements in studying and modeling nonlinear elastic materials with complex microstructures under finite deformations.

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.

We explore the recently-proposed 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 three-dimensions and elaborate upon the key concepts underlying the first-order 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.

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.

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.

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 non-isotropic 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 two-dimensional 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.

In this paper, we present the implementation of a small library of three-dimensional cohesive elements. The elements are formatted as user-defined elements, for compatibility with the commercial finite element software ABAQUS. The PPR, potential-based 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 user-defined 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.

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.

This paper addresses the behavior of three-dimensional multilayered pipe beams with interlayer slip condition, under general 3-D 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 degrees-of-freedom, 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 degrees-of-freedom. 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.

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 efficiency, 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 benefit 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.

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 object-oriented 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.

The Park–Paulino–Roesler (PPR) potential-based model is a cohesive constitutive model formulated to be consistent under a high degree of mode-mixity. Herein, the PPR’s generalization to three-dimensions is detailed, its implementation in a finite element framework is discussed, and its use in single-core 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.

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 arc-length 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.

We discuss the use of polygonal finite elements for analysis of incompressible flow problems. It is well-known that the stability of mixed finite element discretizations is governed by the so-called inf-sup condition, which, in this case, depends on the choice of the discrete velocity and pressure spaces. We present a low-order choice of these spaces defined over convex polygonal partitions of the domain that satisfies the inf-sup 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 linearly-complete 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

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 two-dimensional scalar diffusion problem accompany the theoretical considerations.

Low-temperature 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 Mechanistic-Empirical 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 quasi-brittle behaviour of asphalt concrete. Furthermore, TCModel utilises a simplified one-dimensional 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 user-friendly graphical interface with a stand-alone finite-element-based simulation programme. The system includes a preanalyser and data input generator module that develops a two-dimensional 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 fracture-energy-based 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.

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 (at-node, at-element, nodes-to-element, and elements-to-node). 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.

In recent years, differential equation-driven 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 phase-field methods; which generally use a reaction-diffusion equation, such as the Allen-Cahn equation, for topology evolution. In this work, we exhaustively analyze four level-set methods and one phase-field 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.

In this work, we investigate a Tikhonov-type regularization scheme to address the illposedness of the classical compliance minimization problem. We observe that a semiimplicit discretization of the gradient descent flow 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 modified Helmholtz equation. The appearance of “filtering” in this update scheme is different from the current density and sensitivity filtering techniques in the literature. The next iterate is defined 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 well-known forward-backward 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.

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 one-node 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 bi-dimensional unit cells that illustrate the methodology proposed in this work.

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 off-the-shelf 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.

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 specically the truss, to be optimized in presence of a continuum.

One of the fundamental aspects in cohesive zone modeling is the definition of the traction-separation relationship across fracture surfaces, which approximates the nonlinear fracture process. Cohesive traction-separation relationships may be classified as either nonpotential-based models or potential-based models. Potential-based models are of special interest in the present review article. Several potential-based models display limitations, especially for mixed-mode problems, because of the boundary conditions associated with cohesive fracture. In addition, this paper shows that most effective displacement-based 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 mixed-mode cohesive fracture should be selected with great caution.

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 T-peel joints were prepared and tested. For comparison, similar samples were produced using substrates with classical grit blasting surface treatment. Finally, post-failure 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.

The finite element method is used for modeling geotechnical and pavement structures exhibiting significant non-homogeneity. Property gradients generated due to non-homogeneous distributions of moisture is one such example for geotechnical materials. Aging and temperature induced property gradients are common sources of non-homogeneity for asphalt pavements. Investigation of time dependent behavior combined with functionally graded property gradation can be accomplished by means of the non-homogeneous 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.

Uniform grids have been the common choice of domain discretization in the topology optimization literature. Over-constraining geometrical features of such spatial discretizations can result in mesh-dependent, sub-optimal designs. Thus, in the current work, we employ unstructured polygonal meshes constructed using Voronoi tessellations to conduct structural topology optimization. We utilize the phase-field method, derived from phase transition phenomenon, which makes use of the Allen-Cahn differential equation and sensitivity analysis to update the evolving structural topology. The solution of the Allen-Cahn 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. non-Cartesian) domains. The features of the current approach are demonstrated using various numerical examples for compliance minimization and compliant mechanism problems.

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 three-dimensional numerical examples demonstrate that the proposed schemes can significantly reduce computational cost in comparison to the existing element-based approach.

A potential-based cohesive zone model, so called the PPR model, is implemented in a commercial software, e.g. ABAQUS, as a user-defined 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 two-dimensional linear cohesive element and the algorithm for the PPR potential-based model are presented in-detail. 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.

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 edge-split operators while the adaptive mesh coarsening (AMC) is based on a sequence of vertex-removal (or edge-collapse) operators. Nodal perturbation and edge-swap 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 so-called PPR (Park-Paulino-Roesler) potential-based cohesive model is utilized for the constitutive relationship of the cohesive zone model. The examples investigated include mode I fracture, mixed-mode 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 macro-crack and local micro-cracks.

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.

This paper describes an integrated topology optimization technique with concurrent use of both continuum four-node quadrilateral finite elements and discrete two-node beam elements to design structural braced frames that are part of the lateral system of a high-rise 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 high-rise building structures.

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 finite 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 self-contained 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 offers great flexibility 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.

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.

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 two-dimensional 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.

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, arc-length, 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 incremental-iterative 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 object-oriented implementation. Within the unified framework, the strengths and weaknesses of the various solution schemes are examined through numerical examples.

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.

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 as-defined 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.

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 cross-sections 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.

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 multi-objective function that aims to minimize the von-Mises 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.

The indirect tension test (IDT) is frequently used in civil engineering because of its benefits over direct tension testing. In the mid-1990s, 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.

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 non-homogeneities. 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 elastic-viscoelastic correspondence principle, and viscoelastic material gradation is accounted for within the elements by means of the generalized iso-parametric 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.

We propose a hybrid technique to extract cohesive fracture properties of a quasi-brittle (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 potential-based cohesive constitutive model, the so-called PPR (Park-Paulino-Roesler) model. The unconstrained, derivative free Nelder-Mead 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 Polymethyl-methacrylate (PMMA) quasi-static 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 load-displacement response of a four point bend edge cracked fracture specimen.

This paper presents a single-loop algorithm for system reliability-based topology optimization (SRBTO) that can account for statistical dependence between multiple limit-states, and its applications to computationally demanding topology optimization problems. A single-loop reliability-based design optimization (RBDO) algorithm replaces the inner-loop iterations to evaluate probabilistic constraints by a non-iterative approximation. The proposed single-loop SRBTO algorithm accounts for the statistical dependence between the limit-states by using the matrix-based 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, cut-set and link-set systems and provides the gradients of the system failure probability to facilitate gradient-based optimization. In most RBTO applications, probabilistic constraints are evaluated by use of the first-order 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 single-loop RBDO scheme utilizing the second-order 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 three-dimensional 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.

Parallel computing techniques are employed to investigate wave propagation in three-dimensional functionally graded media. In order to obtain effective and efficient parallel finite element mesh representation, a topology-based data structure (TopS) and a parallel framework for unstructured mesh (ParFUM) are integrated. The parallel computing framework is verified 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.

Tailoring specified vibration modes is a requirement for designing piezoelectric devices aimed at dynamic-type 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 well-known grayscale or intermediate material problem arises which can invalidate the post-processing 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 user-defined vibration modes, by applying the TOM and FGM concepts to design functionally graded piezoelectric transducers (FGPT) and non-piezoelectric 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. Two-dimensional results are presented to illustrate the method.

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 load-line 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 non-uniform; (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 traction-separation 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 quasi-brittle plastic, G-10/FR4 Garolite. Fracture tests were conducted on single edge-notched beam specimens (SENB) under four-point bending. A full-field 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 B-spline 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.

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 full-field displacement through an optimization technique in a finite element (FE) framework. The new experimental–numerical hybrid approach has been applied to fiber-reinforced 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. Four-point bend tests are carried out for the measurement of the modulus of elasticity, E, and the Poisson’s ratio, ν, while single edge-notched beams (SENB) are used for measurement of mode-I 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 two-dimensional 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 derivative-free 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 strain-hardening characteristic. The computed CZM is used in direct problem simulation, the results of which are consistent with the experimental global response.

This paper explores the use of manufacturing-type 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 high-rise 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 two-dimensional and three-dimensional building configurations.

A structural analysis framework called the finite particle method (FPM) for structure failure simulation is presented in this paper. The traditional finite-element 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 two-dimensional (2D) truss, the energy is decomposed and balanced during the failure process. Also, a more complicated three-dimensional (3D) structure failure simulation is given. The comparison of the simulation results and the practical failure mode shows the capability of this method.

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.

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 non-piezoelectric 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. Two-dimensional 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 Moonie-type functionally graded piezoactuators is provided according to specified performance requirements.

Under the theory of classical linear fracture mechanics, a finite crack sitting in an isotropic and homogeneous medium is considered. We find that the well-known crack tip singularity, the inverse square-root singularity 1/sqrt(r), may disappear under certain type of loading traction functions. More specifically, depending on the crack-surface 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.

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 non-homogeneous 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 Cauchy-type singular integral equations. In the numerical results, the emphasis is placed on the investigation of interactions between the contact stress field and the crack-tip 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 in-plane stress component along the coating surface are examined and the mixed-mode stress intensity factors evaluated from the crack-tip stress field with the square-root 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 punch-edge 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 crack-face contact and friction.

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 three-dimensional 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.

Macroscopic constitutive relationship is estimated by considering the microscopic particle/matrix interfacial debonding. For the interfacial debonding, the PPR potential-based cohesive model is utilized. The extended Mori–Tanaka model is employed for micromechanics, while a finite element-based 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 snap-back).

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 three-dimensional numerical examples that the resolution of the design can be significantly improved without refining the finite element mesh.

This paper proposes a single-loop system reliability-based design optimization (SRBDO) approach using the recently developed matrix-based system reliability (MSR) method. A single-loop method was employed to eliminate the inner-loop 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, cut-set, and link-set 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 cross-sectional 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.

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 three-point 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.

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 edge-swap (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.

Cohesive zone model (CZM) is a key technique for finite element simulation of fracture of quasi-brittle 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, derivative-free Nelder-Mead (N-M) 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 N-M 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.

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 three-dimensional 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.

In topology optimization literature, the parameterization of design is commonly carried out on uniform grids consisting of Lagrangian-type finite elements (e.g. linear quads). These formulations, however, suffer from numerical anomalies such as checkerboard patterns and one-node 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 mesh-dependent sub-optimal 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.

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 wedge-angle. The graded character of the wedges is given either by a linearly varying or by an exponentially varying shear modulus. The inhomogeneous anisotropic self-similar bi-wedge and tri-wedge, loaded by a concentrated force at their apex, are studied analytically under plane strain or generalized plane stress conditions, using the self-similarity 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 bi-wedge and a composite tri-wedge, in which for specific values of gradation, the stress and displacements fields are determined.

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 traction-separation 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 load-displacement 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.

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 mesh-independent 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 two-phase 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 projection-based techniques presented.

Micro-tools offer significant promise in a wide range of applications such as cell manipulation, microsurgery, and micro/nanotechnology processes. Such special micro-tools consist of multi-flexible 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 micro-tool structure acts as a mechanical transformer by amplifying and changing the direction of the piezoceramics output displacements. The design of these micro-tools involves minimization of the coupling among movements generated by various piezoceramics. To obtain enhanced micro-tool performance, the concept of multifunctional and functionally graded materials is extended by tailoring elastic and piezoelectric properties of the piezoceramics while simultaneously optimizing the multi-flexible 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 nano-positioner actuated by two graded piezoceramics, and a micro-gripper 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 two-dimensional models because many of the applications for such micro-tools are planar devices.

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 three-dimensional 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 adjacency-based 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.

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 fiber-reinforced concrete (FRC), in two distinct layers. The fracture energy (G) and residual load capacity (P-delta) of two-layered 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 two-layered concrete beam. A cost-benefit analysis showed that the FRC placed in the bottom lift is more fracture efficient (higher G- and P-delta-values 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.

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.

A generalized potential-based constitutive model for mixed-mode 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. quasi-brittle). The unified potential leads to both intrinsic (with initial slope indicators to control elastic behavior) and extrinsic cohesive zone models. Path dependence of work-of-separation is investigated with respect to proportional and non-proportional paths—this investigation demonstrates consistency of the cohesive constitutive model. The potential-based model is verified by simulating a mixed-mode bending test. The actual potential is named PPR (Park–Paulino–Roesler), after the first initials of the authors’ last names.

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.

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 one-dimensional 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 near-zero shear modulus, and materials with negative Poisson’s ratio.

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 mode-tracking 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 non-FGPT one.

Traditionally, standard Lagrangian-type 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 well-known “checkerboard” pathology in the iterative solution of topology optimization problems. A feasible alternative to eliminate such long-standing problem consists of using hexagonal (honeycomb) elements with Wachspresstype shape functions. The features of the hexagonal mesh include two-node connections (i.e. two elements are either not connected or connected by two nodes), and three edge-based symmetry lines per element. In contrast, quads can display one-node connections, which can lead to checkerboard; and only have two edge-based symmetry lines. In addition, Wachspress rational shape functions satisfy the partition of unity condition and lead to conforming finite element approximations. We explore the Wachspress-type hexagonal elements and present their implementation using three approaches for topology optimization: elementbased, continuous approximation of material distribution, and minimum length-scale through projection functions. Examples are presented that demonstrate the advantages of the proposed element in achieving checkerboard-free solutions and avoiding spurious finescale patterns from the design optimization process.

A Mode-III crack problem in a functionally graded material modeled by anisotropic strain-gradient 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 strain-gradient 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.

Co-extrusion 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 co-extrusion. During extrusion, flow in the barrel and die land in a ram extruder should be plug-like, 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 low-extrudate 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 co-extrusion, a simple method using thin-walled tubes was found to be effective to prepare layered feedrods. Functionally graded cellular structures of cement-based materials were successfully co-extruded by using a low-extrudate velocity when the paste had decreasing shear viscosity from inner to outer layers.

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 steady-state 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 Cauchy-type singular integral equation of the second kind for the unknown contact pressure. Numerical results include the distributions of the contact pressure and the in-plane 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.

The characterization of the softening curve from experimental results is essential for predicting the fracture behavior of quasi-brittle 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 ‘‘two-parameter 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 three-point 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.

The present multiscale investigation employs the initial and total fracture energy through a virtual internal pair-bond
(VIPB). model The proposed VIPB model is an extension of the traditional virtual internal bond
VIB model. Two different types of potentials, a steep short-range potential and a shallow long-range 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 fiber-reinforced concrete, which demonstrates the advantage of the VIPB model over the usual VIB model.

Large-scale 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 topology-based framework for supporting adaptive fragmentation simulations, being able to handle two- and three-dimensional 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 topology-based 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 topology-based 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 self-similar geometries, including both two- and three-dimensional 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.

Professor Fazil Erdogan has influenced several generations of applied and solid mechanicians working in the area of mixed boundary-value 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, high-performance 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.

A functionally-graded material system has been developed with a spatially tailored fiber distribution to produce a four-layered, functionally- graded fiber-reinforced 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 cross-sections 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.

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 co-workers, 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 three-dimensional (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.

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 self-consistent model well predicts the experimental results without introducing any calibration; whereas the Mori-Tanaka model and the generalized self-consistent 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 self-consistent 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.

A two-dimensional 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 three-dimensional 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 semi-infinite 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.

By means of a fundamental solution for a single inhomogeneity embedded in a functionally graded material matrix, a self-consistent 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.

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.

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 piece-wise linear cohesive surface relation is used in order to model fracture in a pre-cracked 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 load-displacement curve with experimental results retrieved from literature.

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 bimorph-type 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 two-dimensional models with a plane-strain assumption. The material gradation plays an important role in improving the actuator performance when distributing piezoelectric (PZT5A) and non-piezoelectric (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).

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 four-node quadrilaterals and eight-node quadrilaterals for piezoelectric FGMs were implemented using the graded finite element concept. A closed form two-dimensional 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).

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 three-point bending-beam test was used to obtain the softening behavior and fracture parameters of each FGCM. The peak loads and initial fracture energy between the plain, fiber-reinforced, and FGCMs were similar; this signified similar crack initiation. The total fracture energy clearly indicated the improvements in fracture behavior of FGCM relative to full-depth plain concrete. The fracture behavior of FGCM depended on the position of the fiber-reinforced layer relative to the starter notch. The fracture parameters of both the fiber-reinforced and plain concrete were embedded into a finite element-based 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.

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 9-node) considering two dimensional plane and axisymmetric problems. The finite element approximate solutions are compared with exact solutions for nonhomogeneous materials.

A finite element-based 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 three-point 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 intrinsic-based traction-opening 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 three-point 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 two-parameter fracture model, were found to adequately characterize the bilinear softening model.

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 four-node quadrilaterals and eight-node quadrilaterals for piezoelectric FGMs were implemented using the graded finite element concept. A closed form two-dimensional 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).

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 micro-level stress is considered, which is computed by means of a mechanical concentration factor using a p-norm of the Von Mises stress criterium (applied to the micro-level stress). Because a 0-1 (void-solid) 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 non-linear interpolation of the Hashin-Strikhman upper and lower bounds is considered. Consistently with the framework adopted here in, the so-called '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.

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 particle-reinforced composite materials subject to uniaxial tension. The Mori-Tanaka 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.

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 three-dimensional 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 short-term 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.

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 pair-wise 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 stress-free 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 semi-analytical 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 steady-state heat conduction in the gradation direction.

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 two-dimensional 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 top-down plane strain cracking and three-dimensional 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 closed-form 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.

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 topology-based 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 traction-separation 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 'quasi-steady-state' 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.

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 fixed-free slender bar considering homogeneous, bi-material, tri-layered 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 bi-material system (with sharp interface). Third, a three-point-bending 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.

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.

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 strain-gradient 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.

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 good-quality 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 10-2000 kPa.

A micromechanical damage model is developed for two-phase functionally graded materials (FGMs) considering the interfacial debonding of particles and pair-wise 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.

The simple boundary element method consists of recycling existing codes for homogeneous media to solve problems in non-homogeneous media while maintaining a purely boundary-only 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. Steady-state heat conduction is investigated for thermal conductivity varying either parabolically, exponentially, or trigonometrically in one or more co-ordinates. A three-dimensional 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.

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 middle-crack 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 crack-tip 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 crack-tip response of a Ti/TiB FGM SE(T) specimen. The numerical results demonstrate that the proposed MBL model captures the effect of T-stress on plastic zone size and shape, constraint effects, etc. for such configurations.

This work describes the effect of material gradation (parallel to the crack plane) on stress intensity factors and K-dominance, i.e. the dominance of the singular region, of fracture specimens; SE(T), SE(B) and C(T). The extent of K-dominance is investigated by comparing the actual stress field with the Williams
asymptotic stress field. Linear-elastic 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 K-dominance yields a curve with a peak point at a certain material gradation. Results of the present study provide valuable insight into the K-dominance of FGMs.

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.

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. Three-dimensional models of bamboo cells were built and simulated under tension, torsion, and bending load cases.

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 user-subroutine which is employed in a commercial finite element code (e.g., UEL in ABAQUS). The cohesive parameters are calibrated by simulation of mode I disk-shaped compact tension results. The ability to simulate mixed-mode fracture is demonstrated. The single-edge notched beam test is simulated where cohesive elements are inserted over an area to allow cracks to propagate in any general direction. The predicted mixed-mode 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.

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 non-homogeneous materials, the interaction integral (conservation integral) originally proposed to evaluate SIFs for a static homogeneous medium is extended to incorporate dynamic effects and material non-homogeneity, and is implemented in conjunction with the finite element method (FEM). The technique is implemented and verified using benchmark problems. Then, various homogeneous and non-homogeneous 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.

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. Finite-element 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 closed-form solutions. The cohesive parameters of finite material strength and cohesive fracture energy are calibrated in conjunction with the single-edge notched beam [SE(B)] test. The CZM is then extended to simulate mixed-mode 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.

This work applies a two-state interaction integral to obtain stress intensity factors along cracks in three-dimensional functionally graded materials. The procedures are applicable to planar cracks with curved fronts under mechanical loading, including crack-face tractions. Interaction-integral terms necessary to capture the effects of material nonhomogeneity are identical in form to terms that arise due to crack-front curvature. A discussion reviews the origin and effects of these terms, and an approximate interaction-integral 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.

This paper presents a novel compact adjacency-based topological data structure for finite element mesh representation. The proposed data structure is designed to support, under the same framework, both two- and three-dimensional 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.

State-of-the-art 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 “on-the-fly,” 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.

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 “constant-constitutive-tensor.” 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).

This paper revisits the interaction integral method to evaluate both the mixed-mode stress intensity factors and the T-stress 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 semi-analytical solutions.

This paper presents numerical simulation of mixed-mode 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 non-equilibrium 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 crack-tip fields in FGMs. The present approach uses a user-defined crack increment at the beginning of the simulation. Crack trajectories obtained by the present numerical simulation are compared with available experimental results.

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.

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 so-called “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.

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.

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 user-defined element (UEL) of the ABAQUS software. Using the CZM, first, crack propagation in a mode I single-edge 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 mixed-mode 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 symmetric Galerkin formulation and implementation for heat conduction in a three-dimensional functionally graded material is presented. The Green’s function of the graded problem, in which the thermal conductivity varies exponentially in one co-ordinate, is used to develop a boundary-only 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.

A fundamental solution is established for a crack in a homogeneous half-plane interacting with a crack at the interface between the homogeneous elastic half-plane 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.

A general methodology is constructed for the fundamental solution of a crack in the homogeneous half-plane interacting with a crack at the interface between the homogeneous elastic half-plane 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.

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 single-edge 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 crack-front characteristics, to investigate test repeatability, to examine variations of fracture energy with temperature, and to investigate mixed-mode fracture. The results from the tests follow expected trends and test variability appears to be within a range typical for asphalt concrete fracture testing.

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 mode-mixity 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 single-edge notched beam (SE(B)). With a slight modification to the notch location, mixed-mode (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 mixed-mode 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.

A disk-shaped 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 in-place asphalt concrete pavements or gyratory-compacted 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.

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 cutting-edge 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 disk-shaped compact tension [DC(T)] test, which appears to be a practical method for determining low-temperature 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 quasi-brittle 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.

This study examines a two-state interaction integral for the direct computation of mixed-mode stress intensity factors along curved cracks under remote mechanical loads and applied crack-face 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 crack-front curvature. We examine the significance of the crack-face 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 three-dimensional 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.

A multiscale modeling method is proposed to derive effective thermal conductivity in two-phase graded particulate composites. In the particle-matrix 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.

A micromechanics-based elastic model is developed for two-phase functionally graded composites with locally pair-wise particle interactions. In the gradation direction, there exist two microstructurally distinct zones: particle-matrix zone and transition zone. In the particle-matrix zone, the homogenized elastic fields are obtained by integrating the pair-wise 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.

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 Mode-I as well as in mixed-mode 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 (three-point bending) under impact loading, and mixed-mode crack propagation in pre-cracked steel and graded plates.

The free-space Green function for a two-dimensional 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 non-singular term. The non-singular component is expressed in terms of one-dimensional Fourier-type integrals that can be computed by the fast Fourier transform.

Since the introduction of the hybrid boundary element method in 1987, it has been applied to various problems of elasticity and potential theory, including time-dependent 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 boundary-only 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.

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 T-stress) for a crack in FGMs under mixed-mode 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, so-called non-equilibrium 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.

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 mixed-mode and non-proportional loading. Each step of crack growth simulation consists of calculation of mixed-mode 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 user-defined crack increment at the beginning of the simulation. Crack trajectories obtained by the present numerical simulation are compared with available experimental results.

A new interaction integral formulation is developed for evaluating the elastic T-stress for mixed-mode 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 hyperbolic-tangent 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.

A ‘‘non-equilibrium’’ 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 so-called ‘‘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, semi-analytical, 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 simple boundary element method for solving potential problems in non-homogeneous 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 co-ordinates. For certain classes of material variations the non-homogeneous problem can be transformed to known homogeneous problems such as those governed by the Laplace, Helmholtz and modified Helmholtz equations. A three-dimensional 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 user-subroutine for graded finite elements. The results from the finite element simulations are used for comparison with the present boundary element solutions.

This paper presents a ‘‘simple’’ boundary element method for transient heat conduction in functionally graded materials, which leads to a boundary-only formulation without any domain discretization. For a broad range of functional material variation (quadratic, exponential and trigonometric) of thermal conductivity and specific heat, the non-homogeneous problem can be transformed into the standard homogeneous diffusion problem. A three-dimensional 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.

An interaction integral method for evaluating T-stress and mixed-mode stress intensity factors (SIFs) for two-dimensional 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 path-independent interaction integral entirely based on the Jð; J1Þ integral. The T-stress 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.

This paper describes the development and application of a general domain integral method to obtain J-values along crack fronts in three-dimensional configurations of isotropic, functionally graded materials (FGMs). The present work considers mode-I, linear-elastic response of cracked specimens subjected to thermomechanical loading, although the domain integral formulation accommodates elastic-plastic behavior in FGMs. Finite element solutions and domain integral J-values for a two-dimensional edge crack show good agreement with available analytical solutions for both tension loading and temperature gradients. A displacement correlation technique provides pointwise stress-intensity values along semi-elliptical surface cracks in FGMs for comparison with values derived from the proposed domain integral. Numerical implementation and mesh refinement issues to maintain path independent J-values are explored. The paper concludes with a parametric study that provides a set of stress-intensity factors for semi-elliptical surface cracks covering a practical range of crack sizes, aspect ratios and material property gradations under tension, bending and spatially-varying temperature loads.

A micromechanics-based elastic model is developed for two-phase functionally graded materials with locally pair-wise 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, pair-wise 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.

A micromechanics-based elastic model is developed for two-phase functionally graded materials with locally pair-wise 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, pair-wise 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.

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.

Free space Green’s functions are derived for graded materials in which the thermal conductivity varies exponentially in one coordinate. Closed-form expressions are obtained for the steady-state 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.

We present a modification to the quarter-point crack tip element and employ this element in two-dimensional boundary integral fracture analysis. The standard singular element is adjusted so that the near-tip 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.

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 material-dependent 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 interface-cohesive 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 single-edge 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 metal-rich specimen with relatively smaller cohesive gradation parameter for the metal.

The interaction integral is a conservation integral that relies on two admissible mechanical states for evaluating mixed-mode 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. self-consistent, Mori–Tanaka, or three-phase model). In the latter case, there is no closed-form expression for the material-property 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.

The interaction integral is an accurate and robust scheme for evaluating mixed-mode 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 so-called ‘‘generalized isoparametric formulation’’. The types of orthotropic material gradation considered include exponential, radial, and hyperbolic-tangent functions. Stress intensity factors for mode I and mixed-mode two-dimensional 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, semi-analytical, or numerical solutions.

The path-independent J*k -integral, in conjunction with the finite element method (FEM), is presented for mode I and mixed-mode 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 mixed-mode 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.

For linear elastic functionally graded materials (FGMs), the fracture parameters describing the crack tip fields include not only stress intensity factors (SIFs) but also T-stress (nonsingular stress). These two fracture parameters are important for determining the crack initiation angle under mixed-mode loading conditions in brittle FGMs (e.g. ceramic/ceramic such as TiC/SiC). In this paper, the mixed-mode SIFs and T-stress 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 T-stress effect. This type of criterion involves the mixed-mode SIFs, the T-stress, 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. self-consistent model). Several examples are given to show the accuracy and efficiency of the interaction integral scheme for evaluating mixed-mode SIFs, T-stress, and crack initiation angle. The techniques developed provide a basic framework for quasi-static crack propagation in FGMs.

Meshless and mesh-based 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 (mesh-based) 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.

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 one-dimensional example, involving an exponentially graded material, is used to further clarify these reasons.

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 three-dimensional 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.

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.

Hypersingular integrals of the type

I-alpha(T-n,m,r) = integral(-1)(1) T-n(s)(1-s(2))(m-1/2)/(s-r)(alpha) ds, |r| < 1

and I-alpha(U-n,m,r) = integral(-1)(1) U-n(s)(1-s(2))(m-1/2)/(s-r)(alpha) ds, |r| < 1

are investigated for general integers alpha (positive) and m (nonnegative), where T-n(s) and U-n(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.

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 boundary-only 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.

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 full-field 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.

This paper describes measurements of the K1 R-curve for a layered ceramic-metallic 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 R-curve behavior rising from Kmeas, = 7 MPa·m1/2 in the initial crack-tip location (38% Ti and 62% TiB) to 31 MPa·m1/2 in pure Ti. Residual stresses acting normal to the crack-plane 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.

This work studies mode I crack growth in ceramic/metal functionally graded materials (FGMs) using three-dimensional interface-cohesive 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 single-edge 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 self-consistent 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.

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.

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 co-ordinates which are integrated into the element stiffness matrix. In particular, stress intensity factors for mode I and mixed-mode two-dimensional problems are evaluated and compared through three di=erent approaches tailored for FGMs: path-independent Jk*-integral, modi:ed crack-closure integral method, and displacement correlation technique. The accuracy of these methods is discussed based on comparison with available theoretical, experimental or numerical solutions.

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), in-plane 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 four-node 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.

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 mixed-mode two-dimensional 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 two-step (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.

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 one-dimensional example, involving an exponentially graded material, is used to further clarify these reasons.

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 mixed-mode two-dimensional 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 two-step (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.

The Green's function for three-dimensional 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 time-dependence 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.

Dense, layered, single- and graded-composition composites of MoSi2 and SiC were formed from elemental powders in one step, using the field-activated pressure-assisted combustion method. Compositions ranging from 100% MoSi2 to 100% SiC were synthesized, with relative densities ranging from 99% to 76%, respectively. X-ray 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 layer-to-layer variation in composition of 10 mol% was sufficient to prevent thermal cracking during formation of the layered functionally graded materials.

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 strain-gradient elasticity theories where higher order singular integral equations naturally arise even in the slope-based formulation.

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 2-D and 3-D 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 method-the Hypersingular Boundary Node Method (HBNM). Numerical results from this new method, for selected 3-D problems in potential theory and in linear elasticity, are presented and discussed in this paper.

The standard (singular) boundary node method (BNM) and the novel hypersingular boundary node method (HBNM) are employed for the usual and adaptive solutions of three-dimensional potential and elasticity problems. These methods couple boundary integral equations with moving least-squares 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 mesh-based 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 renement 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.

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 multi-layered 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.

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.

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.

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 position-dependent 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.

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 J-R testing protocol and numerical procedures that are applicable to a titanium/titanium boride layered functionally graded material. Tensile tests and a two-dimensional axisymmetric finite element model were used to determine the plasticity data for the titanium. Crack growth experiments were conducted in three-point bending using single edge notched bend specimens. Three-dimensional FEA of crack growth initiation and two-dimensional FEA with automatic crack propagation were performed. Two crack propagation conditions based on experimental data were used: (a) crack length versus load-line displacement and (b) crack length versus crack mouth opening displacement. The subsequent predictions of the non-linear finite element models are in reasonable agreement with the measured value of J at initiation and with the rising J-R data during crack propagation.

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 two-dimensional 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.

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 power-law 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 power-law material with position-dependent 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.

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.

An elasto-viscoplastic consistent tangent operator (CTO) concept-based implicit algorithm for nonlinear boundary element methods is presented. Both kinematic and isotropic strain hardening are considered. The elasto-viscoplastic radial return algorithm (RRA) and the elasto-viscoplastic 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.

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 TiB-rich 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' R-curve behavior is investigated in order to understand the mechanics of crack growth. The beam is modeled using two-dimensional finite elements considering elastoplastic behavior. Properties of the mixed phase interlayers are approximated using a rule-of-mixtures 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.

An Object-Oriented 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.

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. Error-hypersingular 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.

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.

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.

This work introduces a methodology for selfadaptive numerical procedures, which relies on the various components of an integrated, object-oriented, computational environment involving pre-, analysis, and post-processing 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 self-adaptive 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.

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 conguration 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.

An attempt is being made to develop procedures to obtain initial state and constitutive model parameters by in-situ testing. A centrifuge model test is used in the present study for comparative purposes. A non-destructive 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.

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 traction-free 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.

Domains containing an `internal boundary', such as a bi-material 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 specied. 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.

This paper studies and compares the domain partitioning algorithms presented by Farhat,1 Al-Nasra 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.

The fixed smeared crack concept with strain decomposition is reformulated utilizing a self-adaptive 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.

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.

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 two-dimensional 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.

Recently, several domain partitioning algorithms have been proposed to effect load-balancing 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 node-based and element-based partitioning strategies are discussed. The spectral algorithms are also evaluated and compared for coarse-grained partitioning using different types of structures modelled by l-D, 2-D and 3-D finite elements.

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 end-points 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.

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.

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 non-connected 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.

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 three-dimensional meshes made up of one-dimensional 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.

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 anti-plane shear. The crack must be modeled as an arc of a circle and wholly inside the solid-otherwise 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.