STRUCTURAL INSTABILITIES FOR THE DESIGN AND CONTROL OF COMPUTATIONAL MATERIALS AND STRUCTURES By Talal Husain Ibrahem Salem A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Civil Engineering – Doctor of Philosophy 2021 ABSTRACT STRUCTURAL INSTABILITIES FOR THE DESIGN AND CONTROL OF COMPUTATIONAL MATERIALS AND STRUCTURES By Talal Husain Ibrahem Salem The main hypothesis in this work is that structural instabilities, in single microbeams and assembled structures, can be tailored to build in computational capabilities into structural systems, thus mimicking advanced circuitry functionalities using only mechanical loading in structural elements. In this thesis, we demonstrate the possibilities of solving complicated fourth order partial differential equations using potential energy optimization in axially loaded beams. A similar concept to the minimum energy path in analog circuitry. The proposed approach relies on controlling the post-buckling instability of elastic Euler- Bernoulli beams subjected to gradually increasing loading. The focus in the first few chapters in this thesis, is to study novel methods to tailor and control the mechanical response of anisotropic bilaterally constrained beams with arbitrary cross-section geometries. In particular, theoretical models are developed, using the small deformation theory assumptions, and the static post- buckling behavior of the elastic beams is analyzed using the energy method. Non- prismatic and Functionally-Graded Materials (FGM) beams are studied with respect to several parameters and effects (i.e., beam shape configurations, materials, geometric properties). Effective control levels over the post-buckling behavior, of the bi-walled beams, were achieved. In addition, the manufacturing of the structural system was studied using 3D printing techniques. It was demonstrated that by tuning the beams’ parameters (material and geometry), the snap-through transition events between post-buckling configurations (i.e., sudden energy release), can be configured to happen at specific times (during the loading stage), and that the generated energy output can be controlled. This principle is the main concept used to implement multifunctional capabilities in structural systems. ACKNOWLEDGEMENTS I would like to express my greatest gratitude to my advisor, Dr. Nizar Lajnef, for his consistent support and encouragement. I appreciate the support of the members of my PhD committee: Prof. Neeraj Buch, Dr. Sara Roccabianca and Dr. Weiyi Lu. Throughout this research, I benefited from discussions with members of my research team, including Dr. Imen Zaabar, and other Postdoc, PhD and MS students working with Dr. Lajnef. Also, I would like to thank my parents for their continued support to chase my dreams and reach the stars. My siblings Dr. Farah, Ms. Nada and my love Ms. Sarah deserve my wholehearted thanks as well. iv TABLE OF CONTENTS LIST OF TABLES ...................................................................................................................... vii LIST OF FIGURES ................................................................................................................... viii Chapter 1: Dissertation Overview – Motivation and Vision ..................................................... 1 1.1 Outline .............................................................................................................................. 1 1.2 Summary of Research ...................................................................................................... 2 1.2.1 Development of Bi-walled Nonuniform Beams ....................................................... 2 1.2.2 Bi-Walled Mechanical Metamaterials with Patternable Anisotropy ........................ 2 1.2.3 Carbon Nanotube and FG Plate-like 2D Mechanical Metamaterials........................ 3 1.2.4 Postbuckling-based Partial Differential Equations Solver ........................................ 4 Chapter 2: Programmable Assembly of Bi-walled Nonuniform Beams: Concept, Modeling and Performance ........................................................................................................................... 5 2.1 Introduction ...................................................................................................................... 5 2.2 Design Principle and Problem Statement ......................................................................... 7 2.2.1 Geometry Variation – Single Beam Analysis ........................................................... 9 2.2.2 Material Variation– Single Beam Analysis ............................................................ 12 2.2.3 Beam Assembly Analysis ....................................................................................... 14 2.3 Experimental and Numerical Validations ...................................................................... 17 2.4 Control of Assembly System Response using the Number of Beams, Stiffness, Length and Thickness ............................................................................................................................ 20 2.4.1 Influence on the Postbuckling Response ................................................................ 20 2.4.2 Influence on the Postbuckling Response ................................................................ 23 2.5 Conclusions .................................................................................................................... 23 Chapter 3: Functionally Graded Materials Beams Subjected to Bilateral Constraints: Structural Instability and Material Topology .......................................................................... 25 3.1 Introduction .................................................................................................................... 25 3.2 Fabrication and Experiments of the FGM beams........................................................... 28 3.2.1 Fabrication of FGM Beams Using 3D Printing Technique .................................... 28 3.2.2 Experimental Setup and Testing ............................................................................. 30 3.3 Theoretical and Numerical Modeling of the Bi-walled FGM Beams ............................ 33 3.3.1 Theoretical Modeling for Postbuckling Analysis ................................................... 33 3.3.2 Numerical Modeling and Results Comparison ....................................................... 39 3.4 Parametric study ............................................................................................................. 40 3.5 Material Topology to Design the FGM Beams with Maneuverable Postbuckling Response.................................................................................................................................... 42 3.6 Additional Studies .......................................................................................................... 46 3.7 Conclusions .................................................................................................................... 49 v Chapter 4: Maneuverable Postbuckling of Extensible Mechanical Metamaterials Using Functionally Graded Materials and Carbon Nanotubes......................................................... 50 4.1 Introduction .................................................................................................................... 50 4.2 Problem Formulation of the Microscale FGM-MM and CNT-MM .............................. 52 4.3 Effective Materials Properties of the FGM-MM and CNT-MM ................................... 54 4.3.1 Structural Characteristics of MM............................................................................ 54 4.3.2 Material Characteristics of FGMs and CNTs ......................................................... 55 4.4 Postbuckling Analysis of the FGM/CNT-MM Subjected to Bilateral Constraints........ 58 4.4.1 General Governing Equations ................................................................................. 58 4.4.2 Postbuckling Analysis of the FGM-MM and CNT-MM ........................................ 62 4.4.3 Energy Method to Solve the Postbuckling Response ............................................. 68 4.5 Numerical Modelling and Validation of the FGM-MM and CNT-MM ........................ 72 4.5.1 Numerical Setup and Modelling ............................................................................. 72 4.5.2 Comparison of the Theoretical and Numerical Models .......................................... 74 4.5.3 Comparison of the Theoretical and Numerical Results with the Existing Study ... 77 4.6 Maneuverability of the Postbuckling Response for the FGM-MM and CNT-MM ....... 79 4.6.1 Material Influence of FGMs and CNTs .................................................................. 79 4.6.2 Geometric Influence of the Cylindrical Corrugation in MM .................................. 80 4.7 Additional Studies .......................................................................................................... 85 4.7.1 User-Defined Subroutine USDFLD for the FGM-MM .......................................... 86 4.8 Conclusions .................................................................................................................... 86 Chapter 5: Structural Instability for the Design and Control of Computational Materials and Structures ............................................................................................................................. 88 5.1 Introduction .................................................................................................................... 88 5.2 Conceptual Design ......................................................................................................... 89 5.3 Theoretical formulation for PDE solver ......................................................................... 91 5.4 Trigonometric Functions ................................................................................................ 91 5.5 Exponential and Polynomial Functions.......................................................................... 95 5.6 PDE solver Fabrication and Experimental Validation ................................................... 96 5.7 Experimental Validation ................................................................................................ 97 5.8 Conclusions .................................................................................................................. 107 Chapter 6: Conclusions and Future Work ............................................................................. 109 6.1 Conclusions .................................................................................................................. 109 6.1.1 Post-Buckling Analysis of Bi-walled Non-Uniform Beams Using Small Deformation Assumptions ................................................................................................... 109 6.1.2 Post-Buckling Analysis of Bi-walled Functionally Graded Material Beams Using Small Deformation Assumptions ........................................................................................ 109 6.1.3 Design and Control of Partial Differential Equations-Solver ............................... 110 6.2 Future Work ................................................................................................................. 110 6.2.1 Optimization of the PDE-solver Parameters to Account for Different PDEs ....... 110 6.2.2 Optimization of the Algorithm to Include Different Loading Conditions ............ 110 REFERENCES ........................................................................................................................ 112 vi LIST OF TABLES Table 1. Geometry functions considered in the theoretical model. ................................................ 9 Table 2. Materials properties of the FGM beams. ........................................................................ 30 Table 3. Geometric properties of the FGM specimens. ................................................................ 32 Table 4. Material distribution functions in the theoretical modelling and the axial force for the postbuckling snap-through in the 3rd, 5th and 7th modes (all in mm and N).................................. 47 Table 5. Buckling mode vs. computing time. ............................................................................... 48 Table 6. Comparison of the material functions between the FGM-MM and CNT-MM. ............. 53 Table 7. Geometric and material properties, element size and type, and loading condition of the FE models for the CNT-MM and FGM-MM. .............................................................................. 76 Table 8. Geometric and material properties used in the comparison between the presented and existing studies. ............................................................................................................................. 78 Table 9. PDE-solver (trigonometric functions), parameter definition. ......................................... 95 Table 10. PDE-solver (exponential and polynomial), parameter definition. ................................ 96 vii LIST OF FIGURES Figure 1. (a) Studied cases of the nonuniform beams with the linear width, sinusoidal width, and radical width. (b) Illustration of the assembly of nonuniform beams under different postbuckling modes. (c) Side view of the postbuckling deformation configurations of the nonuniform beams subjected to bilateral constraints. (d) Symbols defined in the theoretical modelling of the bi-walled beam. ............................................................................................................................................... 8 Figure 2. Comparison of the force-displacement relationship between the current and existing models. .......................................................................................................................................... 18 Figure 3. (a) Experimental setup. (b) Aluminum frames for the bilateral constraints and the fixtures for the clamped-clamped boundary conditions. Nonuniform beam samples tested in the (c) first (i.e., PLA) and (d) second (i.e., polycarbonate) beam assemblies. (e) Deformation configurations of the constant, radical and sinusoidal PLA beams. (f) Comparison between the theoretical, experimental and numerical results for the polycarbonate and PLA beam assemblies. ............... 18 Figure 4. Influence of the (a) beam number, (b) bending stiffness, (c) length and (d) thickness on the postbuckling behavior of the beam assembly. ........................................................................ 21 Figure 5. Effect of the beam width and thickness on travelling distance between mode transitions for (a) linear width, (b) sinusoidal width and (c) radical width. ................................................... 22 Figure 6. (a) 3D printing device used to manufacture the FGM beams (Ultimaker s5 3D printer). (b) Illustrative demonstration of the 3D printing technique used to fabricate the FGM beams. The material distributions of the FGM beams were designed with six segments using PLA, ABS and nylon. (A: PLA, B: ABS/Nylon)................................................................................................... 31 Figure 7. (a) Experimental setup of the FGM beams subjected to the bilateral constraints. (b) Postbuckling mode shape configurations of the PLA/ABS beams in buckling modes 1, 3 and 5 (i.e., Φ1, Φ3 and Φ5) (Green demonstrates the ABS material and blue demonstrates the PLA material). ....................................................................................................................................... 31 Figure 8. Schematic illustration of the bi-walled FGM beam postbuckled to the 7th buckling mode under the axial compression.......................................................................................................... 33 Figure 9. (a) Numerical simulations of the FGM beams defined as bi-layered lamination and the postbuckling shape configurations in mode 1, 3 and 5 (i.e., Φ1, Φ3 and Φ5). (b) Comparison of the normalized force-displacement relationship between the theoretical, numerical and experimental results for the PLA/ABS and PLA/nylon beams. ................................................... 38 Figure 10. (a) Effects of the material functions on the Young’s modulus of the FGM beams, and the (b) postbuckling resistance with respect to the material functions. (c) Effect of the volume fraction r on the Young’s modulus, and the (d) postbuckling resistance with respect to the volume fraction r. (e) Effect of the length L on the Young’s modulus, and the (f) postbuckling resistance viii with respect to different length L (gap, width and thickness of the FGM beams are fixed as 4 mm, 25 mm and 1.32 mm, respectively, for all cases).......................................................................... 41 Figure 11. Flowchart of the optimization of the FGM beams subjected to bilateral constraints, including the postbuckling analysis for structural instability and the material topology for material function optimization. ................................................................................................................... 43 Figure 12. (a) Drop-forces prior and posterior to buckling mode transitions for the FGM beams with the power material function. Distributions of the drop-forces in buckling modes Φ3, Φ5, Φ7 and the summation with respect to E2 and r for the FGM beams with the material functions defined in (b) power, (c) cosine square, (d) decimal logarithm, (e) natural exponential, and (f) sinusoidal functions (L = 200 mm, b = 25 mm, t = 1.82 mm and EA = 3.85 GPa for all cases). ............ 44 Figure 13. Total drop-forces posterior to buckling mode transitions for the FGM beams with power, cosine square, decimal logarithm, natural exponential, and sinusoidal functions. ........... 46 Figure 14. Comparison of the postbuckling responses of the bilaterally constrained beam obtained using the 10-, 20- and 30-mode theoretical model and experiments. ........................................... 48 Figure 15. Principles of the maneuverable, extensible FGM-MM and CNT-MM designed using the structural method of MM and material method of FGMs and CNTs. .................................... 54 Figure 16. (a) Demonstration of the postbuckling deformation for the extensible CNT-MM and FGM-MM subjected to bilateral constraints. (b) Deformation analysis of an arbitrary segment in the postbuckled CNT-MM and FGM-MM. .................................................................................. 56 Figure 17. Flowchart of numerically solving the postbuckling response of FGM-MM and CNT- MM using the energy minimization. ............................................................................................. 72 Figure 18. (a) Meshed CNT-MM and FGM-MM with cylindrical corrugation under the axial force and bilateral constraints. (b) Postbuckling shape configurations of the CNT-MM, FGM1-MM and FGM2-MM from the first buckling mode (Φ1), flatten first buckling mode (Φ1 flatten), third buckling mode (Φ3), flatten third buckling mode (Φ3 flatten), and fifth buckling mode (In this study, FGM1-MM denote the MM designed with the first type of FGMs, and FGM2-MM are designed by the second FGMs). .................................................................................................... 74 Figure 19. Comparison of the force-displacement relations and postbuckling shape configurations at buckling mode transitions between the theoretical and numerical results for the (a) CNT-MM and (b) FGM1-MM and FGM2-MM. ........................................................................................... 77 Figure 20. Comparison of the force-displacement relations between the reduced CNT-MM and FGM-MM models in this study, and the theoretical results in the existing study [63]. ............... 78 Figure 21. Influence of the material functions on the postbuckling response of the CNT-MM, FGM1-MM and FGM2-MM. Force-displacement relations of the (a) CNT-MM affected by the volume fraction V, (b) FGM1-MM affected by the index b, and (c) FGM2-MM affected by the index b. Postbuckling shape configurations of the (d) CNT-MM in Φ1 for V = 2.8%, Φ3 for V = 12%, and Φ5 for V = 28%, and the (e) FGM1-MM and (f) FGM2-MM in Φ1 for b = 0.3, Φ3 for ix b = 0.6 , and Φ5 for b = 0.9 ( L = 2 mm , Wb = 0.4 mm , tMM = 10 μm , h0 = 30 μm , and EM = 130 GPa). ........................................................................................................................... 80 Figure 22. Influence of the cylindrical corrugation on the postbuckling response of the CNT-MM, FGM1-MM and FGM2-MM. Force-displacement relations of the (a) CNT-MM, (b) FGM1-MM and(c) FGM2-MM affected by the diameter-to-height ratio DH. Postbuckling shape configurations of the (d) CNT-MM, (e) FGM1-MM and (f) FGM2-MM in Φ1 for DH = 0.2, Φ3 for DH = 1, and Φ5 for DH = 5. (DH are obtained from the cases of 20 μm-to-100 μm, 20 μm-to-20 μm and 100 μm-to-20 μm). .............................................................................................................................. 82 Figure 23. Comparison of the material and geometric influences on the force-displacement relations for the (a) CNT-MM with DH = 0.2 and V = 2.8%, DH = 1 and V = 12%, and DH = 5 and V = 28%, and (b) FGM1-MM and (c) FGM2-MM with DH = 0.2 and b = 0.3, DH = 1 and b = 0.6, and DH = 5 and b = 0.9. Distributions of A1, A3, A5 and A7 for the (d) CNT-MM with DH = 0.2 and b = 0.3, and the (e) FGM1-MM and (f) FGM2-MM with DH = 0.2 and b = 0.3. ................................................................................................................................................ 84 Figure 24. Variation of the maximum axial force in buckling mode transitions for the (a) CNT- MM, (b) FGM1-MM and (c) FGM2-MM. ................................................................................... 85 Figure 25. Demonstration of the (a) FGM-MM and CNT-MM designed with cylindrical corrugation, and the (b) hollow beam that has the same length, width, height and thickness as the FGM-MM and CNT-MM. ............................................................................................................ 86 Figure 26. Postbuckling shape configurations in mode 1st, 3rd, 5th and 7th. .................................. 90 Figure 27. PDE solver arbitrary cross-section. ............................................................................. 91 Figure 28. 3D printer for the fabrication of the PDE solver. ........................................................ 97 Figure 29. Comparison of the force-displacement relationship between the theoretical and experimental results for the material and geometric variation cases. ......................................... 107 x Chapter 1: Dissertation Overview – Motivation and Vision This work describes the efforts undertaken towards implementing, studying and controlling the mechanism of post-buckling in simple and assembled structural systems. The tuning and control of snap-through events in deformed constrained strips, was demonstrated as a proof-of- concept approach to embed computational capabilities in structural systems, in particular for this work, the ability to solve complex fourth order differential equations under real loading in real time. An example of possible practical use can be the design of a mechanical structural apparatus that implements the partial differential equations of a real physical problem (for example a vibrating component on a bridge). The apparatus can solve and simulate the component response in real time under real loading, alleviating the need for any complex finite element or other heavy modeling approaches. This concept can be easily extended to simulate in real time any physical phenomena described by partial differential equations, or can be directly used to endow materials with the capability to compute and analyze their response under service loading, with the purpose of real time structural analysis. Practical manufacturing and experimental considerations were also studied in this work. 1.1 Outline The dissertation can be drawn as: Chapter 1 summarizes the dissertation chapters and sections and provide the main motivation of this work. Chapter 2 present the successful attempts to control the postbuckling behavior of bi-walled beams by tuning its geometric properties. Chapters 3 and 4 reports the tailoring of the stored potential energy for the constrained system by varying the materials’ properties. Chapter 5 shows the implementation of the mechanism of structural instability for solving high order PDEs. Chapter 6 presents the conclusions and the suggested future work. 1 1.2 Summary of Research 1.2.1 Development of Bi-walled Nonuniform Beams Postbuckling structural instabilities has been shown to have useful mechanical characteristics such as well deformation resistance and recovery. However, the difficulties in control and programmability, due to the complex interconnected sensitivities to geometric and material properties, severely hinder the phenomenon use in multifunctional structural applications. In this chapter, we propose a concept of nonuniform beam assembly to increase the controllability of the postbuckling response. Bilaterally constrained, nonuniform beams are theoretically investigated to obtain the buckling instability, and the predictions are compared with the experimental and numerical results with satisfactory agreements. Parametric studies are carried out to demonstrate the tunability of the reported beam assembly with respect to the geometric properties and material parameters (i.e., Young’s modulus) of the nonuniform beams. Finally, the use of the proposed beam assembly method is investigated for novel applications as mechanical triggers and deformation detectors. This study demonstrates an exciting approach to tune the mechanical characteristics of engineered assembly structures for novel applications, such as material embedded mechanical sensing. 1.2.2 Bi-Walled Mechanical Metamaterials with Patternable Anisotropy In recent years, the study of Functionally graded materials (FGM) opened exciting new venues for the control and manipulation of engineered materials and structures. In this chapter, we investigate bilaterally constrained FGM beams with programmable material functions. The FGM beams are fabricated using 3D printing techniques, and tested to understand the behavior of structural instability (i.e., postbuckling) under the bilateral confinements. Theoretical and numerical models are developed to investigate the postbuckling response, and the results are 2 compared to experimental observations with satisfactory agreements. Material topology optimization is then carried out to investigate the influences of the material functions on the release of the stored energy during the bucking mode transitions in the FGM beams. It is found that stored energy variations can be used to optimize the material functions, which allows for the guided design of bi-walled FGM beams with well-defined controllability over the structural instability. The reported bilaterally constrained FGM beams with optimized material functions can be used in a multitude of different applications. 1.2.3 Carbon Nanotube and FG Plate-like 2D Mechanical Metamaterials Architected, structural materials have been reported with promising enhancement of mechanical performance using the structural method (e.g., mechanical metamaterials MM) and material method (e.g., composite materials). Here, we develop the extensible, plate-like mechanical metamaterials at the microscale using functionally graded materials (FGM-MM) and carbon nanotubes (CNT-MM) to obtain the advanced structural materials with well maneuverability over the postbuckling response. Theoretical models are developed to investigate the postbuckling response of the FGM-MM and CNT-MM subjected to bilateral constraints, and numerical simulations are carried out to validate the theoretical results. The theoretical models are then used to investigate the maneuverability of the postbuckling behaviors with respect to the material properties (i.e., FGMs and CNTs) and geometric properties (i.e., corrugated microstructures). The findings show that more significant corrugation in MM and more obvious composition in FGMs and CNTs provide higher controllability on the buckling mode transitions. The reported CNT-MM and FGM-MM provide a novel direction of programming mechanical response of artificial materials. 3 1.2.4 Postbuckling-based Partial Differential Equations Solver Solving nonlinear partial differential equations (PDEs) is complex, time consuming and may be error-prone. Thus, they typically require high-performance computers to be solved. Generally, solving PDEs analytically is based on finding a change of variable to transform the equation into soluble form. Here, a novel design concept of PDEs analog solver is reported based on structural instabilities (i.e., postbuckling) of a bilaterally confined beam. 3D-printed beams are fabricated in the centimeter scale to consider both design maneuverability and manufacturing cost. In this approach, the potential energy balance after the snap-through energy release events, represent the solution for a fourth order PDE equation. Analytical models are developed to obtain the response generated by the PDE-solver under gradually increasing static axial compressive load. 4 Chapter 2: Programmable Assembly of Bi-walled Nonuniform Beams: Concept, Modeling and Performance 2.1 Introduction Buckling has been usually considered as a limit state to avoid in typical structural design [1- 3]. However, in recent years significant interest has been raised toward the potential positive usages of the buckling and postbuckling states in smart materials and structures [4]. For instance, instability of slender elements (i.e., buckling) has been designed in advanced structures to obtain predominant mechanical response such as negative Poisson’s ratio [5], enhancement of deformation resistance [6], or deformation recovery caused by geometric nonlinearity [7]. Buckling-enabled monostable, bistable, and multistable systems have been reported in multifunctional devices [8]. Moreover, the instability behavior of a stiff thin film device integrated on a cylindrical substrate was investigated for potential uses in smart wearable devices (e.g., smartwatch, wristband, etc.) [9]. The buckling behavior of mechanical metamaterials was also studied to develop energy absorption devices by maximizing its snap-through response [10]. In oreder to maximize the generated electrical energy, the control and optimization of the postbuckled beams (generators) mechanical characteristics, is required beams. Four control approaches have been reported for the tuning of the postbuckling response: (i) Material altering, specifically the switching of isotropic and orthotropic properties [11-13]; (ii) Geometric approach absed on investigating non-uniform and non-prismatic geometries [14, 15]; (iii) Constraint changing strategies, from regular confinements to irregularly distributed constraints [16, 17]; and (iv) our investigated assembly strategies based combining multiple easy to build and design uniform beams to form a system with multi-variable control possibilities 5 Theoretical, numerical and experimental studies were conducted to investigate the influence of material and geometric properties on the beams. It was shown that the material variables are less influential when compared with the effects of geometry changes. More importantly, modifying a single beam element is insuffient for effective optimization of the overall postbuckling mode transitions. Therefore, the constraint strategy was used to change the rigid and fixed walls into movable and spring-constrained walls [18, 19]. In addition, Borchani at al. used the assembly strategy to investigate a multi-beam system that stacks the elasticas in a parallel configuration, thus tailoring the postbuckling response of the system using the number and geometric properties of the uniform beams [20]. However, using only uniform beams limits the controllability of the postbuckling behavior (i.e., unable to achieve certain postbuckling mode transitions). In this chapter, postbuckled non-uniform beam assembly systems are studied theoretically, numerically, and experimentally. This work combines all of the discussed control strategies to achieve a precise control over the postbuckling performance (i.e., snap-through mode transitions) of bi-walled elasticas. In principle, placing multiple nonuniform beams between parallel confinements (i.e., varying the geometric properties and number of non-prismatic beams) provides effective tunability over the postbuckling mode transitions. Tuning the stacked beams in the assembly leads to more accurate controllability that is beneficial for the considered applications. To investigate the mechanical characteristics of the nonuniform beam assembly, we first develop a discretized theoretical model that accounts for the buckling-induced instability of multiple bi- constrained beams under quasi-static loading. The total potential energy of the beam system is used to determine the postbuckling response of the beams with linear, sinusoidal and radical width. Experiments and numerical simulations are then conducted to validate the theoretical results, and satisfactory agreements are obtained. Parametric studies are carried out using the theoretical model 6 to investigate the sensitivity of the postbuckling response tunability with respect to the geometric strategy (i.e., geometric variation of the beams) and assembly strategy (i.e., number and arrangement of the beams). It is shown that nonuniform beams assembly offers an effective approach to control the postbuckling response, which is required for the envisioned multiscale applications, such as mechanical triggers and damage detectors. 2.2 Design Principle and Problem Statement In this chapter, we combine the beam shape and assembly methods. We develop a system assembly approach of multiple nonuniformly shaped beams, in varying configurations, in order to effectively tune the postbuckling mode transitions (Φ defines the buckling mode) under a range of different excitations, e.g., external force, environmental temperature fluctuations, etc. Effectively controlling the mode transitions is critical, which can be verified by the characteristics of the postbuckling response. The axial stiffness (i.e., slope of the force-displacement relation) is maintained between different postbuckling modes. The axial displacement can be converted into transverse motion, which can be used to a activate piezoelectric patch to generate electrical signals for different purposes such as energy harvesting [14, 16]. Tuning the mode transitions significantly influences the generated electrical signals. For example, it is desirable to maximize the harvested electrical power for energy harvesting purposes, while the objective in sensing applications is to match the response amplitude to classified damage events. Given the difficulties, reported in previous studies, of precisely controlling the postbuckling mode transitions, we introduce the nonuniform beam assembly approach to demonstrate a programmable postbuckling system response. Three beam shapes are investigated, linear, sinusoidal and radical, as shown in Figure 1(a). Each individual beam ‘k’ is fully defined by its length 𝐿𝑘 , elastic modulus 𝐸𝑘 , net gap between the 7 bilateral constraints ℎ𝑘 (𝑥) thickness 𝑡𝑘 (𝑥), width 𝑏𝑘 (𝑥), nonuniform cross section area 𝐴𝑘 (𝑥) and moment of inertia 𝐼𝑘 (𝑥). Figure 1(b) illustrates the nonuniform beam assembly system that will exhibit different buckling mode transitions for the different beams under an external applied load on the system. Figure 1(c) provides the side view of the postbuckling shape configurations for the nonuniform beams subjected to bilateral constraints. The change in the buckling modes with increasing displacement can be demonstrated in four stages: (1) original shape, (2) point contact, (3) line contact, and (4) snap-through. The beams switch from (1) to (2) when the buckling limit is exceeded under the axial displacement. Increasing the displacement, the deformation is more severe, and the point contact is enlarged to a line contact. Further, the deformed beam snaps to a higher buckling mode. Figure 1. (a) Studied cases of the nonuniform beams with the linear width, sinusoidal width, and radical width. (b) Illustration of the assembly of nonuniform beams under different postbuckling modes. (c) Side view of the postbuckling deformation configurations of the nonuniform beams subjected to bilateral constraints. (d) Symbols defined in the theoretical modelling of the bi-walled beam. 8 2.2.1 Geometry Variation – Single Beam Analysis The postbuckling response of a bi-walled single elastic beam with nonuniform cross section under axial force control was studied and investigated [14]. For a single beam, the system consists of clamped-clamped nonuniform beam positioned between two bilateral fixed frictionless walls as presented in Figure 1(d) The beam will be subjected to axial displacement ∆, and the beam has a length 𝐿, elastic modulus 𝐸, net gap between the bilateral constraints ℎ(𝑥) thickness 𝑡(𝑥), width 𝑏(𝑥)𝑡 3 (𝑥) 𝑏(𝑥), nonuniform cross section area 𝐴(𝑥) = 𝑏(𝑥)𝑡(𝑥) and area moment of inertia 𝐼(𝑥) = 12 that vary continuously along the axial direction. Small deformation assumptions are adopted in Euler–Bernoulli Table 1. Geometry functions considered in the theoretical model. Case Geometry function* Linear width 𝑏(𝑋) = 𝑏𝑡𝑜𝑝 + ( 𝑏𝑏𝑜𝑡 − 𝑏𝑡𝑜𝑝 )𝑋 Sinusoidal width 𝑏(𝑋) = 𝑏𝑏𝑜𝑡 + 2𝐴̅𝑠𝑖𝑛(2𝜋𝑚𝑋) Radical width 𝑏(𝑋) = 𝑏𝑏𝑜𝑡 − 2𝐵̅ √𝑋𝐿 * : 𝐴̅ is the sine amplitude and 𝐵̅ is the radical width variation parameter. beam theory to model the beam under consideration. According to the conducted studies [21, 22], the simple local balance of moments equation that governs the buckling behavior of beam with nonuniform cross section, under an axial compressive force 𝑝̂ , may be written as 𝑑2 𝑀(𝑋) + 𝑁 𝑀(𝑋) = 0; (2-1) 𝑑𝑋 2 𝑝̂ where, 𝑁 = 𝐸𝐼(𝑋) and the normalized transverse deflection 𝑊(𝑋) is used to define the curvature of the beam 𝑀(𝑋), as 𝑑2 𝑊(𝑋) 𝑀(𝑋) = ; (2-2) 𝑑𝑋 2 9 𝑥 ̂ (𝑋𝐿) 𝑤 where 𝑋 = 𝐿 and the normalized transverse displacement 𝑊(𝑋) = . The boundary ℎ conditions of the governing equation are 𝑊(0) = 𝑊(1) = 0; {𝑑𝑊(𝑋) 𝑑𝑊(𝑋) (2-3) |𝑋=0 = 𝑑𝑋 |𝑋=1 = 0; 𝑑𝑋 and the general solution of Eq. (2-2) can be expressed as follows, 𝑀(𝑋) = 𝐶1 Ω1 (𝑋) + 𝐶2 Ω2 (𝑋); (2-4) where 𝐶1 and 𝐶2 are the unknown integral constants, and Ω1 (𝑋) and Ω2 (𝑋) represent the linearly independent special solutions for different cross section area configurations. Integrating the curvature of the beam 𝑀(𝑋) leads to the general solution of the diverse beam as, 𝑑𝑊(𝑋) = 𝐶1 ∫ Ω1 (𝑋)𝑑𝑋 + 𝐶2 ∫ Ω2 (𝑋)𝑑𝑋 + 𝐶3 ; (2-5) 𝑑𝑋 𝑊(𝑋) = 𝐶1 ∬ Ω1 (𝑋)𝑑𝑋𝑑𝑋 + 𝐶2 ∬ Ω2 (𝑋)𝑑𝑋𝑑𝑋 + 𝐶3 𝑋 + 𝐶4 . (2-6) Higher order expansion improves the solution. Therefore, direct numerical integration was used and since different independent specific solutions will result from different beam shape configurations, three possible cases have been studied, (i) linear width; (ii) sinusoidal width; and (iii) radical width. The beam geometry varies with respect to the normalized longitudinal coordinate X . Table 1 summarized the functions which the beams follow in varying their geometry. Since the linearly independent special solutions Ω1 (𝑋) and Ω2 (𝑋) depend on the distribution of the beam flexural stiffness in Eq. (2-7), therefore the first step is to evaluate the variable quantities 𝛼, 𝛽 and ζ for each beam shape configurations with respect to its geometry function. 𝐸𝐼(𝑋) = 𝛼(1 + 𝛽𝑓(𝑋))ζ . (2-7) According to previous conducted study, the linearly independent special solutions for each of the prementioned cases can be expressed as follows, 10 2.2.1.1 Linear width Ω1 (𝑋) = √1 + 𝛽𝑋 𝐽1 [𝜑]𝐿𝑛𝑟𝑤 ; Ω2 (𝑋) = √1 + 𝛽𝑋 𝑌1 [𝜑]𝐿𝑛𝑟𝑤 ; Ω(𝑋)𝐿𝑛𝑟𝑤 = 1 (2-8) 𝜑 𝐿𝑛𝑟𝑤 = 2𝑛(1 + 𝛽𝑋)2 ; 𝑝̂𝐿2 𝑛 = √𝛼𝛽2 ; { 𝐸𝑏𝑡𝑜𝑝 𝑡 3 𝑏𝑏𝑜𝑡 where, 𝛼 = ;𝛽= − 1; ζ = 1 and 𝐽; 𝑌 are the first and second type of Bessel function. 12 𝑏𝑡𝑜𝑝 Substituting Eq. (2-8) into Eqs. (2-5) and (2-6), to get the general deflection functions, then the integration constants C3 and C4 can be evaluated and expressed in terms of C1 and C2 by using the given boundary conditions in Eq. (2-3) as follows, 0 3 −1 3 5 𝑑𝑊(0)𝐿𝑛𝑟𝑤 2𝐾2 𝑝 𝐹𝑞 [ ;(2, ); −𝐾] 2 2 −0.5 0.5 𝐾 0= = 𝐶1 5 + 𝐶2 𝐺𝑚 𝑛 . 5 + 𝐶3 ; (2-9) 𝑑𝑋 3𝛽𝑛2 √𝜋 𝑝 𝑞 −1 −1 𝛽𝑛2 √𝜋 √𝐾 ( 0.5) 0 5 −1 3 7 𝐿𝑛𝑟𝑤 4𝐾2 𝑝 𝐹𝑞 [ ;(2, ); −𝐾] 2 2 −0.5 0.5 𝐾 0 = 𝑊(0) = 𝐶1 5 + 𝐶2 𝐺𝑚 𝑛 . 9 + 𝐶4 ; (2-10) 15𝛽 2 𝑛2 √𝜋 𝑝 𝑞 −2 −1 𝛽 2 𝑛 2 √𝜋 √𝐾 ( 0.5) 0 5 −1 3 7 𝐿𝑛𝑟𝑤 4𝐾2 𝑝 𝐹𝑞 [ ;(2, ); −𝐾] 2 2 −0.5 0.5 𝐾 0 = 𝑊(1) = 𝐶1 + 𝐶2 𝐺𝑚 𝑛 . +𝐶 + 𝐶4 ; (2-11) 5 15𝛽 2 𝑛2 √𝜋 𝑝 𝑞 −2 −1 𝛽2 𝑛92 √𝜋 3 √𝐾 ( 0.5) 11 where, 𝑝 𝐹𝑞 and 𝐺𝑚 𝑛 are the generalized hypergeometric and Meijer G functions. With all 𝑝 𝑞 integration constants defined in terms of C4, eigenvalues and buckling mode shapes can be obtained numerically. 2.2.1.2 Sinusoidal width Ω1 (𝑋) = √1 + 𝛽sin (2𝜋𝑚𝑋) 𝐽1 [𝜑 𝑆𝑖𝑛𝑢𝑤 ]; Ω2 (𝑋) = √1 + 𝛽sin (2𝜋𝑚𝑋) 𝑌1 [𝜑 𝑆𝑖𝑛𝑢𝑤 ]; Ω(𝑋)𝑆𝑖𝑛𝑢𝑤 = 1 (2-12) 𝜑 𝑆𝑖𝑛𝑢𝑤 = 2𝑛(1 + 𝛽𝑋)2 ; 𝑝̂𝐿2 𝑛 = √𝛼𝛽2 ; { 𝐸𝑏𝑡𝑜𝑝 𝑡 3 2𝐴̅ where, 𝛼 = ;𝛽= ; and ζ = 1. 12 𝑏𝑏𝑜𝑡 2.2.1.3 Radical width Ω1 (𝑋) = √1 + 𝛽𝑋 𝐽2 [𝜑 𝑅𝑎𝑑𝑤 ]; 3 Ω2 (𝑋) = √1 + 𝛽X 𝑌2 [𝜑 𝑅𝑎𝑑𝑤 ]; 3 Ω(𝑋)𝑅𝑎𝑑𝑤 = 𝑅𝑎𝑑𝑤 4 3 (2-13) 𝜑 = 3 𝑛(1 + 𝛽𝑋)4 ; 𝑝̂𝐿2 𝑛 = √𝛼𝛽2 ; { 𝐸√ 𝑏𝑏𝑜𝑡 𝑡 3 2𝐵̅ 1 where, 𝛼 = ; 𝛽=− ; and ζ = 2 . 12 𝑏𝑏𝑜𝑡 By adopting the same procedure followed in the linear width case for other cases, the shape functions of the cases under consideration can be obtained. 2.2.2 Material Variation– Single Beam Analysis For an anisotropic beam defined in Cartesian coordinates, the beam has a length 𝐿 in the 𝑥- direction, width 𝑏 in the 𝑦-direction, a thickness 𝑡 in the 𝑧-direction, a cross section area 𝐴 = 𝑏 𝑡, bt3 an area moment of inertia 𝐼 = and a Young’s modulus 𝐸(𝑥) that varies continuously along the 12 longitudinal direction 𝑥 in accordance with the rule of mixtures. Different distribution functions 12 (i.e., simple power law, natural exponential and decimal logarithm) can be used to control the variation of the Young’s modulus in the longitudinal direction, from pure material 1 at (𝑥 = 0) to pure material 2 at (𝑥 = 𝐿), 𝑥 𝑟 𝐸1 − [𝐸1 − 𝐸2 ] (𝐿 ) ; Simple power law [𝐸1 −𝐸2 ] 𝐸(𝑥) = 𝐸1 + 𝑥𝑟 ( ) ; Natural exponential (2-14) 𝑒 𝐿 𝑥 2 {𝐸1 + [𝐸1 − 𝐸2 ] log (𝐿 ) 𝑟 ; Decimal logarithm 𝑉 where, 𝑟 = 𝑉1 describes the variation in the volume fraction, 0 < 𝑥 < 𝐿 , 𝐸1 is defined as the 2 Young’s modulus of the stiffer material, [𝐸1 − 𝐸2 ] is the difference between Young’s moduli of the two materials being used, and 𝐸2 is the softer material’s Young’s modulus. The governing equation of the buckling behavior of the presented anisotropic beam, may be re-written as 𝑑2 𝑀(𝑋) + Ň 𝑀(𝑋) = 0; (2-15) 𝑑𝑋 2 12𝑝̂ where, Ň = 𝑏𝑡 3 𝐸(𝑋) and integration constants, eigenvalues, and buckling shape functions can be obtained numerically by adopting the same prescribed boundary conditions. According to the previous study [16], the obtained buckling shape functions form an orthogonal basis and thus, the superposition method is used to express the transverse deflection as a linear combination of buckling modes as 𝑊(𝑋) = ∑∞ 𝑗=1,3,5… 𝐴𝑗 𝜓𝑗 (𝑋) ; (2-16) where 𝐴𝑗 are the coefficients that express each buckling mode influence on the transverse deflection. 13 2.2.3 Beam Assembly Analysis In the case of an assembly system, the previously mentioned analytical formulation is extended to account for multiple different beams subjected to axial displacement simultaneously. The system consists of N anisotropic/nonuniform beams stacked in parallel. All beams will be subjected to the specified axial displacement ∆𝑘 simultaneously. We use displacement control since it helps to increase the controllability of loading on the beams, thus leading to more precise tuning over the post-buckling process. The post buckling responses for the beams with anisotropic/nonuniform sections are analyzed by applying the energy method. The weight coefficients are obtained through minimizing the total energy between the bilateral constraints. Since the beams are deflected quasi-statically between frictionless bilateral constraints, the kinetic energy can be ignored, and the potential energy will signify the total energy. The summation of the bending strain energy 𝑢𝑏𝑘 , compressive strain energy 𝑢𝑐𝑘 , and energy of external work 𝑢𝑝𝑘 represent the total potential energy. The beams are subjected to a gradually increasing axial displacement. The axial force 𝑝̂𝑘 due to the applied axial displacement can be expressed as 𝐿 𝑑𝑤̂ (𝑥) 2 𝑘 ∆𝑘 ∫0 𝑘 ( 𝑑𝑥 ) 𝑑𝑥 GV: 𝑝̂𝑘 = 𝐸𝑘 [ 𝐿 1 − 𝐿𝑘 1 ]; ∫0 𝑘 𝐴 (𝑥) 𝑑𝑥 2 ∫0 𝐴 (𝑥) 𝑑𝑥 𝑘 𝑘 ̂ 𝑘 (𝑥) 2 (2-17) 𝐿𝑘 𝑑𝑤 ∆𝑘 ∫0 ( 𝑑𝑥 ) 𝑑𝑥 MV: 𝑝̂ 𝑘 = 𝐴𝑘 [ 𝐿𝑘 1 − 𝐿 1 ]; ∫0 𝐸 (𝑥) 𝑑𝑥 2 ∫0 𝑘 (𝑥) 𝑑𝑥 { 𝑘 𝐸𝑘 where GV and MV denote the geometry variation for nonuniform beams and the material variation for anisotropic beams, respectively, The prescribed energy quantities can be expressed in terms of axial displacement ∆𝑘 as follows, 14 2 𝐸𝑘 𝐿 𝑑2 𝑤̂ 𝑘 (𝑥) GV: 𝑢𝑏𝑘 = 𝑘 ∫0 𝐼𝑘 (𝑥) ( ) 𝑑𝑥 ; 2 𝑑𝑥 2 { 2 (2-18) 𝑏𝑘 𝑡𝑘3 𝐿𝑘 𝑑2 𝑤̂ 𝑘 (𝑥) MV: 𝑢𝑏𝑘 = ∫ 𝐸𝑘 (𝑥) ( ) 𝑑𝑥 ; 24 0 𝑑𝑥 2 𝑆𝑘 ∆𝑘 𝑢𝑐𝑘 = 𝑐 ; (2-19) 2 𝑝̂𝑘 ∆𝑘 𝑢𝑝𝑘 = ; (2-20) 2 where 𝑆𝑘 and ∆𝑘𝑐 , respectively, refer to the axial compressive force and deformation under the geometry and material considerations as, 𝐿 𝑑𝑤̂ (𝑥) 2 𝑘 ∆𝑘 ∫0 𝑘 ( 𝑑𝑥 ) 𝑑𝑥 GV: S𝑘 = 𝐸𝑘 [ 𝐿 1 − 𝐿𝑘 1 ]; ∫0 𝑘 𝐴 (𝑥) 𝑑𝑥 2 ∫0 𝐴 (𝑥) 𝑑𝑥 𝑘 𝑘 ̂ 𝑘 (𝑥) 2 (2-21) 𝐿 𝑑𝑤 ∆𝑘 ∫0 𝑘( 𝑑𝑥 ) 𝑑𝑥 MV: S𝑘 = 𝐴𝑘 [ 𝐿𝑘 1 − 𝐿𝑘 1 ]; { ∫0 𝐸 (𝑥) 𝑑𝑥 2 ∫0 𝐸𝑘 (𝑥) 𝑑𝑥 𝑘 and 𝑝̂ 𝐿 1 GV: ∆𝑘𝑐 = 𝐸𝑘 ∫0 𝑘 𝐴 𝑑𝑥 ; 𝑘 𝑘 (𝑥) { 𝑝̂𝑘 𝐿𝑘 1 (2-22) MV: ∆𝑘𝑐 = ∫ 0 𝑑𝑥 . 𝑏𝑘 𝑡𝑘 𝐸𝑘 (𝑥) Note that ∆𝑘 is the variation of the beam length given as 1 𝐿 𝑑𝑤̂ 𝑘 (𝑥) 2 ∆𝑘𝑐 = ∆𝑘 − 2 ∫0 𝑘 ( ) 𝑑𝑥. (2-23) 𝑑𝑥 Substituting Eqs. (2-21) to (2-23) into Eqs. (2-18) to (2-20), the potential energy is the summation of the partial potential energies as, 15 𝛱(𝐴𝑘𝑗 ) = ∑𝑁 𝑘 𝑘 𝑘 𝑘=1 𝑢𝑏 (𝐴𝑗 ) + 𝑢𝑐 (𝐴𝑗 ) − 𝑢𝑝 (𝐴𝑗 ) ; 𝑘 𝑘 𝑘 (2-24) which can be written as, 2 𝐿 𝑑𝑤̂ 𝑘(𝑥) 2 𝐿 𝑑𝑤 ̂ (𝑥) 2 2 𝐸𝑘 (2∆𝑘 ∫0 𝑘 ( ) 𝑑𝑥−[∫0 𝑘( 𝑘 ) 𝑑𝑥] ) 𝐸𝑘 𝐿𝑘 𝑑2 𝑤 ̂ 𝑘 (𝑥) 𝑑𝑥 𝑑𝑥 𝑘 𝛱 = ∫ 𝐼 (𝑥) ( ) 𝑑𝑥 − ; (2-25a) 2 0 𝑘 𝑑𝑥 2 𝐿 8 ∫0 𝑘 1 𝑑𝑥 𝐴𝑘 (𝑥) 2 2 𝐿 𝑑𝑤 ̂ 𝑘 (𝑥) 2 𝐿 𝑑𝑤 ̂ (𝑥) 2 𝐿 𝑑2 𝑤 ̂ 𝑘 (𝑥) 𝑏𝑘 𝑡𝑘 (2∆𝑘 ∫0 𝑘 ( ) 𝑑𝑥−[∫0 𝑘 ( 𝑘 ) 𝑑𝑥] ) 𝑏𝑘 𝑡𝑘3 ∫0 𝑘 ( 2 ) 𝑑𝑥 𝑑𝑥 𝑑𝑥 𝑑𝑥 𝛱𝑘 = [ 𝐿 1 ]− 𝐿 1 ; (2-25b) 24 ∫0 𝑘 𝐸 (𝑥)𝑑𝑥 8 ∫0 𝑘 𝐸𝑘 (𝑥) 𝑑𝑥 𝑘 Substituting the boundary conditions and Eqs. (2-16) and (2-17) into Eq. (2-25), the total potential energy can be expressed as, 2 1 𝑑𝜓𝑗 (𝑋) 2 1 𝑑𝜓𝑗 (𝑋) 2 2 2∆𝑘 ∫0 ℎ𝑘2 (𝑋)[∑∞ 𝑗=1 𝐴𝑗 𝑘 ] 𝑑𝑋−[∫0 ℎ𝑘2 (𝑋)[∑∞ 𝑘 𝑗=1 𝐴𝑗 𝑑𝑋 ] 𝑑𝑋] 𝐿𝑘 1 2 𝑑2𝜓 𝑗 (𝑋) 𝑑𝑋 GV: 𝛱 𝑘 = 𝐸𝑘 ( ∫ ℎ (𝑋)𝑏𝑘 (𝑋)𝑡𝑘3 (𝑋) [∑∞ 24 0 𝑘 𝑗=1 𝐴𝑗 𝑘 𝑑𝑋 2 ] 𝑑𝑋 − 1 1 ); 8𝐿𝑘 ∫0 𝑑𝑋 𝑏𝑘 (𝑋)𝑡𝑘 (𝑋) 2 𝑑2 𝜓𝑗 (𝑋) 2 1 𝑑𝜓𝑗 (𝑋) 2 1 𝑑𝜓𝑗 (𝑋) 2 1 ∫0 [∑∞ 𝑘 𝑏𝑘 𝑡𝑘 (2∆𝑘 ∫0 [∑∞ 𝑘 𝑗=1 𝐴𝑗 ] 𝑑𝑋−[∫0 [∑∞ 𝑗=1 𝐴𝑗 𝑘 ] 𝑑𝑋] ) 𝑗=1 𝐴𝑗 ] 𝑑𝑋 𝑏𝑘 𝑡𝑘3 𝐿𝑘 𝑑𝑋2 𝑑𝑋 𝑑𝑋 MV: 𝛱 𝑘 = ℎ𝑘2 ( [ 1 1 ]− 1 1 ); 24 ∫0 𝐸 (𝑋)𝑑𝑋 8𝐿𝑘 ∫0 𝐸 (𝑋)𝑑𝑋 𝑘 𝑘 { (2-26) where, 𝑗 = 1, … … . , ∞, 𝑘 = 1, … … . , 𝑁. The unknown weight coefficients 𝐴𝑗𝑘 can be determined by solving the following constrained minimization problem, which represent the physical confinement of the beam between the bilateral boundaries. Min[𝛱(𝐴𝑘𝑗 )]; 𝑗 = 1, … , ∞, 𝑘 = 1, … , 𝑁 { (2-27) 0 ≤ 𝑊𝑘 (𝑋) ≤ 1. The Nelder–Mead algorithm is adopted to solve the energy minimization problem in Eq. (2-27). It is worthwhile noting that we consider 30 buckling modes to solve the numerical minimization problem in this study. This is mainly due to the considerations of accuracy and computational cost. The 30-mode model offers the most accurate theoretical results when compared with the experimental results, especially for the postbuckling mode transitions (i.e., Φ3, 16 Φ5 and Φ7). On the other hand, the computational cost of solving the 30-mode theoretical model is approximately two times of that for the 20-mode model. 2.3 Experimental and Numerical Validations The theoretical model was simplified to account for a single uniform beam under displacement control to validate the results with the existing model [20]. Figure 2 presents the theoretical force- displacement response for a constant width beam subjected to the quasi-static axial force. The beam has the length of 250 mm, thickness of 2.34 mm, total separated gap of 4 mm, Young’s modulus of 2.31 GPa and constant width of 30 mm. The third, fifth and seventh buckling mode shapes are labeled as Φ3, Φ5 and Φ7, respectively, indicating the buckling shapes after the snap- through events. Both elastic curves match with the differences of 1.22, 1.77 and 1.73% at Φ3, Φ5 and Φ7, respectively. The difference between both theoretical models, could be due to the numerical errors while minimizing the total potential energy. However, since the range of error is very small (<2%), the theoretical model can be adopted to investigate the effect of other parameters on the postbuckling response of a nonuniform cross-section beams assembly subjected to an axial displacement loading. To validate the theoretical postbuckling response of the beams, the beam assemblies were experimentally investigated. The experimental setup was similar to the adopted setup presented in previous conducted study, the procedure was adjusted to account for multiple non-uniform beams (i.e., different lengths and thicknesses) simultaneously, as shown in Figure 3(a). 17 Figure 2. Comparison of the force-displacement relationship between the current and existing models. Figure 3. (a) Experimental setup. (b) Aluminum frames for the bilateral constraints and the fixtures for the clamped-clamped boundary conditions. Nonuniform beam samples tested in the (c) first (i.e., PLA) and (d) second (i.e., polycarbonate) beam assemblies. (e) Deformation configurations of the constant, radical and sinusoidal PLA beams. (f) Comparison between the theoretical, experimental and numerical results for the polycarbonate and PLA beam assemblies. 18 The experimental setup structure consists of aluminum frames with 13 mm thickness, in order to hold the beams firmly in place and emphasize the independent response of each beam, Softer restraints may influence the responses due to the contact interactions between each the buckled beams and the adjacent walls. Transverse reactions (i.e., forces, moments) will accelerate the development of the higher buckling modes, resulting in earlier than expected mode transitions. The beam ends were fixed and positioned between the two rigid constraining walls, while the gaps between the walls were free to adjust which allows the use of varied geometries. The aluminum frames provided the bilateral constraints, and the fixed-fixed boundary conditions were induced by the top and bottom fixtures, as shown in Figure 3(b). In this study, the bilateral constraints are flat, rigid and fixed with limited gap distance, such that the postbuckling mode transitions of the reported beam assembly can be more tunable. Note that the aluminum constraints were lubricated to omit the influence of friction in the experiments. We only considered the clamped-clamped boundary conditions in the beam assembly because of their higher controllability in an experimental setup, compared to other types of boundary conditions such as pinned-pinned or pinned-clamped. An MTS Flextest 40 and a loading frame unit model370 were used for displacement-control testing. Gradually increasing axial force was applied to the top of the rigid beam placed on the top of the assembly frame. The loading stage was limited to a maximum total shorting of 6 mm, and the loading period was set to be 40 seconds. Before loading the nonuniform beam assembly, cyclic load test was carried out on a constant polycarbonate beam, and repeatability was verified. We tested two sets of samples. For the first set, three polycarbonate beams (𝐸 = 2.25 GPa) with different geometries were tested. It includes two linear and one constant width beams (Figure 3(c)). The beams length and Young’s modulus are 250 mm and 2.25 GPa, respectively. In the 19 second set of experiments, three beams (i.e., constant, radical and sinusoidal beams with different length and thickness) were 3D printed, using an UltiMaker 3D printer, with 0.1 mm layer height polylactic acid (PLA) material (Figure 3(d)). The beams width and Young’s modulus are set to 25 mm and 3.5 GPa, respectively. Figure 3(e) displays the postbuckling shape configurations of the 3D printed constant, radical and sinusoidal PLA beams. Finite Element models were developed in Abaqus v16.1. The exact same beam assemblies, discussed above, were simulated. The contact interaction was defined to consider the effect of the bilateral constraints. The walls and microbeams were simulated using shell elements (S4R). Figure 3(f) compares the theoretical, experimental and numerical results for the first (i.e., polycarbonate) and second (i.e., PLA) beam assemblies. Satisfactory agreements are obtained between the theoretical and experimental results. 2.4 Control of Assembly System Response using the Number of Beams, Stiffness, Length and Thickness 2.4.1 Influence on the Postbuckling Response The developed theoretical model is used to analyze the postbuckling behavior of the multi- nonuniform beams assembly to investigate the effects on the assembly’s elastic response of the length L, thickness t, Young’s modulus E and the number of the beams. We consider the second set from the experiments (i.e., assembly of 3D printed beams). The gap between the bilateral 20 constraints was fixed at 4 mm. In the case of sinusoidal and radical width, the amplitude A ̅ and the width variation parameter B ̅ are set to be 0.5 and 0.315, respectively. Figure 4. Influence of the (a) beam number, (b) bending stiffness, (c) length and (d) thickness on the postbuckling behavior of the beam assembly. 21 Figure 5. Effect of the beam width and thickness on travelling distance between mode transitions for (a) linear width, (b) sinusoidal width and (c) radical width. Figure 4(a) shows the influence mechanism of changing the number of beams as a controlling parameter on the assembly elastic behavior. Three, six (2 linear, 2 sinusoidal, 2 radical), and nine (3 linear, 3 sinusoidal, 3 radical) beams assembly configurations were investigated. The length, thickness and Young’s modulus were 250 mm, 2.14 mm and 2.25 GPa, respectively for all beams in this comparison. The results show that increasing the number of beams in the assembly will lead to an axial force increment with respect to the number of beams being used, i.e., the axial forces were twice the values in the case of a using six beams and three times the values in the case of nine beams, however the axial displacement at the transition were the same. Figure 4(b) shows the effect of changing the elastic modulus, in the range 1.65 ≤ 𝐸 ≤ 2.85 GPa, on the assembly’s force-displacement response. The three beams have a length of 250 mm and a thickness of 2.14 mm. The axial forces at each mode transition, for the three different configurations, are considered. 22 The results show that varying material properties (elastic modulus) has no effects on the transitions displacements or stiffness. This implies that it will result in significant increase in the forces at each mode transitions, but fails in controlling the gap between these mode transitions. Figure 4(c) investigates the influence of the beam length, and Figure 4(d) shows the effect beam thickness variation. It is observed that decreasing the beam length and increasing the beam thickness will decrease the gap between the mode transitions, which leads to a better resolution. 2.4.2 Influence on the Postbuckling Response In this section, we define the gap between the mode transitions as the ‘travelling distance’. The aim is to investigate ways to pre-define it and to design an assembly system with linear, sinusoidal and radical beams that will exhibit the desired travel distance response. Although the travelling distance can be varied using the assembly of multiple uniform beams, we aim here to exactly tune it to a set value using the complex nonuniform beams assembly configuration. Figure 5 shows the variation of the travelling distance, between the buckling mode transitions, with respect to the beam cross-sectional area (combined effect of beam thickness and width). Figure 5(a) presents the results for a nonuniform beam with linear width. Figure 5(b) shows the sinusoidal beam results and Figure 5(c) displays the radical beam output. The results show that beam thickness and width can be used to tune the traveling distance between postbuckling mode transitions. Note that the sinusoidal beam shows a significant influence in Figure 5(b), which indicates that sinusoidal shape affects the travelling distance more significant than the other two beam patterns. 2.5 Conclusions This study developed an approach for nonuniform beam assembly to tune and control the postbuckling response of a loaded system. A theoretical model was developed to investigate the postbuckling performance of the nonuniform beams subjected to bilateral constraints. 23 Experimental and numerical simulations were conducted to validate the theoretical predictions. Satisfactory agreements were obtained. Parametric studies were carried out to investigate the tunability of the reported assembly system response in terms of the geometric and material attributes of the nonuniform beams. The reported assembly of nonuniform elasticas offers a novel approach to tune the mechanical characteristics of architected structures for use in novel applications such as energy harvesting and extremely low frequency force/displacement sensing. 24 Chapter 3: Functionally Graded Materials Beams Subjected to Bilateral Constraints: Structural Instability and Material Topology 3.1 Introduction Functionally graded materials (FGM), a type of functional composites focused on exploring mechanical responses through material design and optimization, have marked their debut decades ago and have been experiencing rapid developments ever since then, which rekindles the popularity of research on composite materials [23]. Compared to their traditional counterparts, FGM have several advantages and unusual properties such as enhanced deformation resistance, enhanced toughness, ultra-light, well recoverability, etc. [24-26]. Surpassing the mechanical characteristics of natural materials, FGM are reported with superior response behavior due to their effectively tuned and engineered material properties [27]. A research field has therefore emerged aiming to design, characterize and harness the functional patterns of FGMs in order to obtain desirable performances for specific applications. This is a highly multidisciplinary research community that comprises structural analysis, material science, mechanics and engineering. As a consequence, research efforts have been dedicated to designing composites’ functionally graded material properties, such that to obtain advanced physical and mechanical responses that are otherwise not accessible in nature materials [28-31]. A significant number of studies have been conducted in recent years to characterize and design FGMs under different conditions. Static, free vibration and wave propagation of bi-material beams fused with an FGM layer were investigated analytically using the first-order shear deformation theory [32]. Structural behavior of functionally graded tapered clamped free Euler–Bernoulli beams was studied to determine the critical buckling loads and free frequencies in the longitudinal and transverse directions [33]. General formulas were proposed to obtain the effective stiffness 25 coefficients of the elastic beams made of FGM and piecewise homogeneous materials [34]. Modified couple stress theory and von Kármán geometric nonlinearity were used to investigate the nonlinear vibration of FGM microbeams, where the material properties were assumed to be graded in the thickness direction [35, 36]. The authors discussed the influences of the length scale parameter and material properties on the nonlinear mechanical response of the FGM microbeams. Higher-order shear deformation beam theories were developed for bending and free vibration of FGM beams [37]. The authors reported that varying the power law index in the material functions significantly affected the stiffness of the FGM beams. Effects of material composition on the thermal buckling and vibration of FGM beams were also investigated [38-41]. Higher order shear deformation theory was adopted to develop finite element models for refined mixed beam element [42]. The mixed finite beam model was proposed to explore the vibrational behavior of FGM beams, in which material properties were described using a power law distribution [43]. The mechanical behavior of porous beams, made of FGMs, was inspected by different researchers to explore their principles and features (i.e., abundance of delamination) compared with the conventional materials [44]. The dynamic behavior of the porous FGM beams was studied while the Young's moduli and mass density were nonlinearly graded along the thickness direction [45- 47]. The authors reported the effects of the porosity distributions and porosity coefficients on the dynamic behavior of the FGM porous beams. In addition to investigating the influence of material functions on the mechanical response of FGM structures, recent research studies have also been dedicated to optimizing the material functions in order to tune and control the response. Two optimization strategies have been proposed to tune FGMs for enhanced mechanical or thermal behavior, including varying the material functions with respect to the 1) volume fraction distribution of FGM, and 2) relationships between material gradation and property variations, 26 which leads to various potential applications (e.g., dentistry, bone implants, energy harvesting, etc.) [48-52]. Tailoring material properties of FGM to achieve specific desired goals has received considerable attention in recent years [53-58]. To date, the majority of the studies has been focusing on the material properties of FGM. However, only a limited number of reported studies have investigated FGM systems optimization for complex performance with structural instabilities (e.g., postbuckling behavior). Structures with controlled instabilities (i.e., postbuckling) in bilaterally confined FGM beams, have been investigated for deployment for smart applications, in several recent studies, including innovative devices and techniques for energy harvesting, structural health monitoring (SHM), artificial muscles, and soft robotics [8, 59, 60]. It is critical to effectively tune the postbuckling response of the bi-walled FGM beams in those devices. Therefore, studies have been particularly focused on investigating the postbuckling mechanisms. Buckling mode transitions of bi-walled beams were investigated using the geometric equilibrium conditions [61, 62]. More recently, an energy method was developed to determine the postbuckling response by minimizing the total energy of the deformed beams between the bilateral confinements [63, 64]. To improve the overall performance of the devices in those applications, studies were conducted on the control of postbuckled beams using structure-related methods. The postbuckling response was tuned by varying the geometric properties of bi-walled beams (i.e., nonuniform beams) [65], assembling the beams into controllable systems (i.e., beam assembly) [20], and changing the lateral confinements from rigid to flexible discontinuous (i.e., irregular constraints) [66]. However, studies have not been proposed on the bi-walled FGM beams for controllarble postbuckling. Behavior that could be harnessed for multiple applications. To address this research gap, we investigate in this chapter, bi-walled FGM beams with a focus on optimizing the functional patterns. Here, we propose bilaterally confined FGM beams with 27 programmable postbuckling behavior by optimizing the material functions. The FGM beams are fabricated using 3D printing techniques, and the experiments are conducted to test the postbuckling response under axial loading. Theoretical and numerical models are developed to validate the experimental results and satisfactory agreements are obtained. The material topology model is then developed to optimize the material pattern such that to accurately control the postbuckling response. This study develops a promising approach to design and optimize FGM structures for multipurpose applications. For example, taking advantage of programmable postbuckling responses which convert quasi-static excitations into snap-through high accelerations, the reported bi-walled FGM beams can be applied to trigger piezoelectric materials for energy harvesting under low-frequency excitations [67, 68]. Identifying the electric signals generated by piezoelectric transducer as limit states, the FGM beams are also expected to be used in SHM methods to detect thermal damage under low-frequency ambient temperature variation [69]. The chapter is outlined as follow: Section 2 introduces the fabrication of the FGM beams using 3D printing technique, and discusses the experimental setups for postbuckling under bilateral constraints. Section 3 presents the analytical and numerical models of the bi-walled FGM beams and compare the results with the experiments. Section 4 discusses the material topology used to optimize the functions of the reported bi-walled FGM beams for maneuverable structural instability. Section 5 summarizes the main findings in this study. 3.2 Fabrication and Experiments of the FGM beams 3.2.1 Fabrication of FGM Beams Using 3D Printing Technique In this study, the FGM beams were fabricated using a 3D printing technique, and the bilateral constraints were made of alumina. The Ultimaker-S5 dual nozzle printer implementing the fused filament fabrication technology was used in the experiments. This provided flexibility and 28 maintained structural integrity over the FGM samples. Fig. 1 illustrates the principle and printing process of the dual nozzle 3D printing. The polymer filaments were fed and heated by the 3D printer to convert the materials from the glassy to rubbery state, then extruded through the nozzle to the platform such that to construct the desired FGM beam distributions. The extrusion head consisted of two 0.4 mm diameter brass nozzles, which can be moved in 360○. In particular, the nozzles were fed with different materials, which were moved in the x − y plane to deposit the FGM beam designed in AutoCAD 2019 and Siemens NX. The direction of action of the compressive load was in the x-axis, the layer build direction was in the z-axis, and the raster was moving in the x-y plane to deposit the material. The platform was moved in the z direction layer by layer until the designed thickness was reached [70]. Note that the thickness tolerance was approximately 0.06 mm. FGM specimens were fabricated using the polylactic acid (PLA), acrylonitrile butadiene styrene (ABS) and Nylon mainly since the difference of their failure modes (i.e., Nylon and ABS typically experienced elastic failure while PLA was brittle) [71]. In 3D printing, the dual extrusion head rate was set to be 25 mm/s, and the printing bed temperature was varied from 0-85 oC depending on the melt temperatures of the feeding materials. The temperatures of the nozzles for the material pattern of PLA/ABS was 200/230 oC and PLA/nylon pattern was 200/245 oC. The standby temperatures (which the inactive nozzle cools down to, while switching nozzles in dual extrusion machines) for PLA and nylon were fixed as 175 oC and the temperature for ABS was 85 oC. The fan speeds used to print PLA, ABS and nylon were set as 100%, 2% and 40 %, respectively. The FGM beams were printed with the layer tolerance of 0.06 mm. The material properties of the PLA, nylon and ABS used in the experiments are listed in Table 2. 29 The FGM beams were divided into six segments with the length of 28.3 mm, and the material fraction was varied in each segment in the 𝑥 −direction, as shown in Figure 6(b). The volume distribution of the FGM beams can be described as the power, cosine squared, decimal logarithm, natural exponential, and sinusoidal functions in the longitudinal direction, 𝑥 𝑟 𝑉1 − (𝑉1 − 𝑉2 ) (𝐿 ) ; 𝟐𝛑𝑥 𝑟 𝑉1 + (𝑉1 − 𝑉2 ) cos2 ( ) ; 𝐿 𝑥 𝑉(𝑥) = 𝑉1 + (𝑉1 − 𝑉2 ) log (𝐿 ) 𝑟 2 ; (0 < 𝑥 < 𝐿) (3-1) (𝑉1 −𝑉2 ) 𝑉1 + 𝑥𝑟 ( ) ; 𝑒 𝐿 𝛑𝑥 𝑟 { 𝑉1 + (𝑉1 − 𝑉2 ) sin ( 2𝐿 ) ; where the total beam length and volume fraction r are 170 mm and 0.76, respectively. 3.2.2 Experimental Setup and Testing The postbuckling response of the bilaterally constrained FGM beams was experimentally investigated. Figure 7 presents the experimental setup and testing results of the deformed Table 2. Materials properties of the FGM beams. PLA Nylon ABS Density (g/cm3) 1.24 1.14 1.1 Young’s modulus (GPa) 3.47 0.889 2.07 Elongation at break (%) 5.2 210 4.8 Flexural modulus/tensile modulus 1.34 0.8 1.28 Filament length/weight (mm/g) 126.67 137.33 142.67 30 Figure 6. (a) 3D printing device used to manufacture the FGM beams (Ultimaker s5 3D printer). (b) Illustrative demonstration of the 3D printing technique used to fabricate the FGM beams. The material distributions of the FGM beams were designed with six segments using PLA, ABS and nylon. (A: PLA, B: ABS/Nylon). Figure 7. (a) Experimental setup of the FGM beams subjected to the bilateral constraints. (b) Postbuckling mode shape configurations of the PLA/ABS beams in buckling modes 1, 3 and 5 (i.e., Φ1, Φ3 and Φ5) (Green demonstrates the ABS material and blue demonstrates the PLA material). 31 Table 3. Geometric properties of the FGM specimens. FGM beam Width w (mm) 25 Length L (mm) 170 Thickness t (mm) 2.89 Volume fraction r 0.76 shape configurations. Figure 7(a) shows the bilateral constraints made in alumina for the postbuckling response of the FGM beams. The testing protocol consisted of applying the quasi- static axial displacement in three continuous steps (i.e., initial calibration, loading, and unloading) using the MTS Flextest 40 mechanical testing machine [14]. The FGM beams were placed between the adaptable, frictionless aluminum frames under the clamped-clamped boundary conditions using the testing fixtures. The gap between the aluminum frames was 5 mm and the gradually increasing loading was applied to the top of the samples within the loading time period of 40 s. The aluminum frames were designed to test three samples parallelly, and therefore, the testing in the middle frames can effectively prevent probable experimental imperfections (e.g., vibration or rotation). In addition, the frames were lubricated to ensure the negligibility of the friction between the FGM beam and the constraints. Two types of the FGM beams (i.e., PLA/ABS and PLA/nylon) were tested in the experiments Figure 7(b) shows the postbuckling shape configurations of the 3D printed FGM beams under bilateral constraints in buckling modes 1, 3 and 5 (i.e., Φ1, Φ3 and Φ5). 32 3.3 Theoretical and Numerical Modeling of the Bi-walled FGM Beams 3.3.1 Theoretical Modeling for Postbuckling Analysis In the theoretical modelling, the FGM beams are considered with the length 𝐿 in the 𝑥 - direction, width 𝑏 in the 𝑦-direction, thickness 𝑡 in the 𝑧-direction, cross-sectional area 𝐴 = 𝑏𝑡, bt3 and moment of inertia 𝐼 = , as shown in Figure 8. The FGM beams are under the clamped- 12 clamped boundary conditions. The net gap ℎ is defined as the distance from the FGM beams to the constraint at the far end, i.e., ℎ = ℎ0 − 𝑡, where ℎ0 is the gap between the constraints. The beams are subjected to axial displacement ∆, and the small deformation assumptions are adopted into the Euler–Bernoulli beam theory to develop the theoretical model since the ratio of the Figure 8. Schematic illustration of the bi-walled FGM beam postbuckled to the 7th buckling mode under the axial compression. ℎ net gap-to-length is relatively small (𝐿 ≪ 1). The material functions of the FGM beams are defined by the Young’s modulus, which is varied from complete material A at one end (𝑥 = 0) to complete material B at the other end (𝑥 = 𝐿). The effective Young’s moduli are defined as the power, cosine squared, decimal logarithm, natural exponential, and sinusoidal functions in the longitudinal direction, which are given as 33 𝑥 𝑟 𝐸𝐴 − |𝐸𝐴 − 𝐸𝐵 | (𝐿 ) ; 𝟐𝛑𝑥 𝑟 𝐸𝐴 + |𝐸𝐴 − 𝐸𝐵 | cos2 ( ) ; 𝐿 𝑥 𝐸(𝑥) = 𝐸𝐴 + |𝐸𝐴 − 𝐸𝐵 | log (𝐿 ) 𝑟 2 ; (0 < 𝑥 < 𝐿) (3-2) |𝐸𝐴 −𝐸𝐵 | 𝐸𝐴 + 𝑥𝑟 ( ) ; 𝑒 𝐿 𝛑𝑥 𝑟 { 𝐸𝐴 + |𝐸𝐴 − 𝐸𝐵 | sin ( 2𝐿 ) ; 𝑉𝐴 where, 𝑟 = describes the variation in the volume fraction, 𝐸𝐴 is defined as the Young’s modulus 𝑉𝐵 of the stiffer material, |𝐸𝐴 − 𝐸𝐵 | is the absolute difference between the Young’s moduli of the two materials, and 𝐸𝐵 is the Young’s modulus of the softer material. The second order differential equation that governs the buckling behavior of the FGM beams, under an axial compressive force 𝑝̂ , can be written as [21, 22] d2 𝑀(𝑋) + 𝑁(𝑋) 𝑀(𝑋) = 0; (3-3) d𝑋 2 12𝑝̂ where 𝑁(𝑋) = 𝑏𝑡 3 𝐸(𝑋). Note that the normalized governing equation accounts for the postbuckling response of the FGM beams with variable geometries (i.e., the differential equation is valid for variable coefficients), since 𝑁(𝑋) is defined with respect to the variable modulus 𝐸(𝑋). The normalized transverse deflection 𝑊(𝑋) is used to define the curvature 𝑀(𝑋) as d2 𝑊(𝑋) 𝑀(𝑋) = ; (3-4) d𝑋 2 𝑥 ̂ (𝑋𝐿) 𝑤 where 𝑋 = 𝐿 and 𝑊(𝑋) = . The boundary conditions are ℎ 𝑊(0) = 𝑊(1) = 0 {d𝑊(𝑋) d𝑊(𝑋) (3-5) |𝑋=0 = d𝑋 |𝑋=1 = 0. d𝑋 34 𝑀(𝑋) = 𝐶1 sin(𝛽𝑋) + 𝐶2 cos(𝛽𝑋); (3-6) 12𝑝̂ where 𝐶1 and 𝐶2 are the unknown integral constants, and 𝛽 = √𝑏𝑡 3 𝐸(𝑋). Integrating 𝑀(𝑋) leads to obtain the general solution for the FGM beams as d𝑊(𝑋) cos(𝛽𝑋) sin(𝛽𝑋) = −𝐶1 + 𝐶2 + 𝐶3 ; (3-7) d𝑋 𝛽 𝛽 sin(𝛽𝑋) cos(𝛽𝑋) 𝑊(𝑋) = −𝐶1 − 𝐶2 + 𝐶3 𝑋 + 𝐶4 ; (3-8) 𝛽 𝛽 where 𝐶1 ,𝐶2 ,𝐶3 and 𝐶4 are the integration constants that can be determined using the boundary conditions. Substituting Eq. (3-5) into Eqs. (3-7) and (3-8), the four algebraic equations are obtained as 𝐶 − 𝛽1 + 𝐶3 = 0 cos(𝛽) sin(𝛽) −𝐶1 + 𝐶2 + 𝐶3 = 0 𝛽 𝛽 𝐶 ; (3-9) − 𝛽2 + 𝐶4 = 0 𝑠𝑖𝑛(𝛽) 𝑐𝑜𝑠(𝛽) {−𝐶1 𝛽 − 𝐶2 𝛽 + 𝐶3 + 𝐶4 = 0 where 𝐶1 ,𝐶2 and 𝐶3 are defined in terms of 𝐶4 . The eigenvalues and buckling shape functions 𝜓𝑗 can be obtained numerically. The general form of the buckling shape function can be expressed as 𝛽 sin(0.5𝛽)+1.4142𝛽 cos(0.7071𝛽) sin(𝛽𝑋) 𝜓(𝑋) = 1 − cos(𝛽𝑋) + cos(0.7071𝛽)+1.4142𝛽 sin(0.7071𝛽)−cos(0.5𝛽) [ − 𝑋]. (3-10) 𝛽 35 According to [20], the obtained buckling shape functions form an orthogonal basis, and thus, the shape function of the FGM beams can be expressed using the superposition method by linearly combining the buckling modes as 𝑊(𝑋) = ∑∞ 𝑗=1,3,5… 𝛼𝑗 𝜓𝑗 (𝑋) ; (3-11) where 𝛼𝑗 are the coefficients that determine the contributions of the buckling modes to the shape function. The energy method is used to determine the weight coefficients 𝛼𝑗 in Eq. (3-11). In particular, the energy minimization has been applied to minimize the total energy of the postbuckled FGM beams between the bilateral constraints, since the total energy of the beams is, at any equilibrium states, the minimum. Due to the negligibility of the kinetic energy (i.e., quasi-static loading), the potential energy is considered. The total potential energy of the FGM beams is given as the summation of the bending strain energy 𝑢𝑏 and the energy resulting from compressive strain and external work 𝑢𝑠 . Because the beam is under the displacement control, the axial force 𝑝̂ induced by the gradually applied axial displacement and the midplane rotation due to buckling is: ̂ (𝑋) 2 d𝑤 𝐿 ( ) 𝐸(𝑋) d𝑥 𝑝̂ = 𝛥𝑏𝑡 [∫0 (𝐸(𝑋) − ) 𝑑𝑋]. (3-12) 2𝛥 The potential energy components can be obtained in terms of the overall variation in length due to the buckled configuration and the axial compressive displacement as 2 𝑏𝑡 3 𝐿 d2 𝑤 ̂ (𝑋) 𝑢𝑏 = ∫0 𝐸(𝑋) ( ) d𝑋 ; (3-13) 24 d𝑋 2 and 36 𝑆 𝑢𝑠 = 2 (𝛥𝑐 − 𝛥); (3-14) where 𝛥, 𝑆, and 𝛥𝑐 refer to the entire variation of the beam’s length, axial compressive force, and axial compressive deformation, respectively, which are defined in terms of the overall deformation 𝛥 as ̂ (𝑋) 2 d𝑤 𝐿 ( ) 𝐸(𝑋) d𝑥 𝑆= 𝛥𝑏𝑡 [∫0 (𝐸(𝑋) − ) 𝑑𝑋] ; (3-15) 2𝛥 𝑝̂ 𝐿 1 𝛥𝑐 = 𝑏𝑡 ∫0 d𝑥 ; (3-16) 𝐸(𝑥) and 1 ̂ (𝑋) 2 𝐿 d𝑤 𝛥𝑐 = 𝛥 − 2 ∫0 ( ) d𝑥. (3-17) d𝑥 Substituting Eqs. (3-15) to (3-17) into Eqs. (3-13) and (3-14), the total potential energy can be written as 37 2 𝑏𝑡 3 𝐿 d2 𝑤 ̂ (𝑋) 𝑝̂ ̂ (𝑋) 2 𝐿 d𝑤 𝛱 = 𝑢𝑏 + 𝑢𝑠 = ∫0 𝐸(𝑋) ( ) d𝑋 − 4 [∫0 ( ) d𝑋]. (3-18) 24 d𝑥 2 d𝑋 Figure 9. (a) Numerical simulations of the FGM beams defined as bi-layered lamination and the postbuckling shape configurations in mode 1, 3 and 5 (i.e., Φ1, Φ3 and Φ5). (b) Comparison of the normalized force-displacement relationship between the theoretical, numerical and experimental results for the PLA/ABS and PLA/nylon beams. Substituting Eqs. (3-5) and (3-11) into Eq. (3-18), the total potential energy with respect to the normalized deflection 𝑊(𝑋) can be expressed as 2 2 𝑏𝑡 3 𝐿 𝐿 𝑑2 𝜓𝑗 (𝑋) 𝑝̂ 𝑑𝜓𝑗 (𝑋) 2 𝛱 = ℎ2 [ 24 ∫0 𝐸(𝑥) [∑∞ 𝑗=1 𝛼𝑗 𝑑𝑋 2 ] 𝑑𝑋 − 4𝐿 [([∑∞𝑗=1 𝛼𝑗 𝑑𝑋 ] 𝑑𝑋) ]] ; (3-19) where, 𝑗 = 1, … … . , ∞ A total of 30 buckling modes are used to numerically solve the minimization of the total energy mainly due to the considerations of accuracy and computational cost. The 30-mode model offers the most accurate theoretical results comparing with the experimental results, especially for the postbuckling mode transitions (i.e., Φ3, Φ5 and Φ7). On the other hand, the computational cost of solving the 30-mode theoretical model is approximately two times of that for the 20-mode model. Thus, coefficients 𝛼𝑗 can be obtained as, 38 Min[𝛱 (𝛼𝑗 )], 𝑗 = 1, … , ∞ { . (3-20) 0 ≤ 𝑊(𝑋) ≤ 1 It is worthwhile to point out that the reported governing equation Eq. (3-3) is valid only for the buckling response of the FGM beams. This study uses the energy method to obtain the postbuckling behavior by minimizing the total potential energy of the postbuckled FGM beams between the bilateral constraints. Therefore, the energy minimization is calculated under 0 ≤ 𝑊(𝑋) ≤ 1, where 1 is the normalized constraints gap. Due to the complexity of the total potential energy (i.e., higher order expansions), numerical tools are used to obtain the postbuckling response of the beam assembly. In particular, the Nelder-Mead algorithm is deployed because of its well efficiency and accuracy in solving the energy minimization problem. The algorithm needs a set of n + 1 points for a function of n variables, to form polytope’s vertices in the n -dimensional space. The arrangement of the points are conducted in the form of 𝑓 (𝑥1 ) ≤ 𝑓 (𝑥2 ) ≤ . . . ≤ 𝑓 (𝑥𝑛+1 ). In order to replace the worst point 𝑥𝑛+1 , a new point is generated to the centroid of the polytope 1 consisting of the best 𝑛 points as 𝑐 = 𝑛 ∑𝑛𝑖=1 𝑥𝑖 . A trial point 𝑥𝑙 by reflecting the worst point through the centroid is expressed as 𝑥𝑙 = 𝑐 + 𝛾(𝑐 − 𝑥𝑛+1 ), where 𝛾 > 0 is a parameter. If the new point 𝑥𝑙 is neither a new worst n or best point, 𝑥𝑙 is replaced by 𝑥𝑛+1 . If the new generated point 𝑥𝑙 is better than the previous best point, the reflection can be moved forward to 𝑥𝑒 = 𝑐 + Ϛ(𝑥𝑙 − 𝑐), where Ϛ > 𝑙 is a parameter to expand the polytope. The Nelder–Mead algorithm accurately solves the energy minimization problem in this study. Detailed discussion can be found in [41]. 3.3.2 Numerical Modeling and Results Comparison The FE models of the bi-walled FGM beams are developed using Abaqus v16.1. The FGM beam and bilateral constraints are simulated using the shell elements (S4R), as shown in Figure 9(a). According to the 3D printed samples in the experiments, the FGM beams are partitioned into 6 39 parts and the bi-layer laminations [0/90] of PLA/ABS and PLA/nylon are defined in the material properties. Contact interaction was defined to take into account the conflict between the beam and constraints. Buckling and postbuckling calculations are conducted to obtain the postbuckling response. In particular, the buckling analysis is defined as the linear perturbation/buckle, and the postbuckling analysis is the dynamic implicit with Nlgeom. The dynamic loading is applied to the top edge of the FGM beams, and the theoretical, experimental and numerical results are obtained for the force-displacement relationships of the PLA/ABS and PLA/nylon samples. The geometric and material properties in Table 2 and Table 3 are used. The interaction between the beam and constraints was defined as “hard” contact without friction (i.e., frictionless) in Abaqus, mainly due to two reasons: 1) the axial displacement was relatively small comparing to the beam length, which means that the sliding of the contact region is negligible; and 2) lubrication was applied to the alumina constraints to reduce the friction in the experiments. Figure 9(b) compares the theoretical, experimental and numerical results of the force-displacement relationships for the PLA/ABS and PLA/nylon beams. Satisfactory agreements are obtained, which demonstrates the accuracy of the theoretical prediction. 3.4 Parametric study Figure 10 studies the influence of the material functions (i.e., the Young’s modulus in Eq. (3-2)) on the postbuckling response 40 Figure 10. (a) Effects of the material functions on the Young’s modulus of the FGM beams, and the (b) postbuckling resistance with respect to the material functions. (c) Effect of the volume fraction r on the Young’s modulus, and the (d) postbuckling resistance with respect to the volume fraction r. (e) Effect of the length L on the Young’s modulus, and the (f) postbuckling resistance with respect to different length L (gap, width and thickness of the FGM beams are fixed as 4 mm, 25 mm and 1.32 mm, respectively, for all cases). of the FGM beams. The volume fraction r and beam length L are particularly varied. Note that the FGM beams are investigated with the complete material A as PLA (i.e., 𝐸𝐴 = 3.47 GPa) and 41 complete material B as ABS (i.e., 𝐸𝐵 = 2.07 GPa). The gap, width and thickness are fixed as 4 mm, 25 mm and 1.32 mm, respectively. The results show that changing the volume fraction r or the distribution function increases the axial force. However, the stiffness (i.e., slope of the force-displacement relations) and the spacing between the mode transitions Φ3, Φ5 and Φ7 are not affected. Figure 10(a) and Figure 10(b) investigate the influences of changing the distribution functions of the Young’s modulus. The material functions are defined as sinusoidal, natural exponential, cosine squared and decimal logarithm functions, and the beam with and volume fraction are fixed as L = 185 mm and 𝑟 = 0.5, respectively. Figure 10(c) shows the effect of changing the volume fraction r on the young’s modulus if the power function was adopted, and Figure 10(d) shows the influence of r on the normalized force-displacement relations. The non-dimensional axial force at which the buckling transitions happen are drastically affected by varying the volume fraction r and the distribution functions. Figure 10(e) and Figure 10(f) study the influence of the beam’s length L on the postbuckling behavior (i.e., the volume fraction is 𝑟 = 1). It can be seen that changing the length affects the stiffness which is likely to postpone further buckling modes. This provides an approach to control over the location and magnitude of the postbuckling mode transitions. 3.5 Material Topology to Design the FGM Beams with Maneuverable Postbuckling Response The reported bi-walled FGM beams release the stored strain energy when postbuckling mode transitions occur. The FGM beams are expected to release more strain energy if the postbuckling response is more significant (i.e., drop-forces are more critical in buckling snap- throughs). As a consequence, the total drop-forces can be used to control the postbuckling response, which leads to the maneuverability over the structural instability of the bi-walled FGM 42 beams. We particularly use the total drop-forces to design the volume fraction r and the Young’s modulus of the softer material (i. e. , 𝐸𝐵 ) in the material functions of the FGM beams. The total drop-forces can be obtained as 𝛥𝑝̂ = ∑ 𝑝̂𝑖 , (𝑖 = 3,5,7); (3-21) where 𝑝̂𝑖 represents the drop-forces happened in buckling mode Φ3, Φ5 and Φ7, which are given in Eq. (3-12). Figure 11 presents the flowchart of the postbuckling analysis and material topology for the FGM beams. Following the former, material topology addresses the total drop-forces from the force-displacement relations of the FGM beams using the energy minimization. Force limit state ε is defined for the control criterion as Figure 11. Flowchart of the optimization of the FGM beams subjected to bilateral constraints, including the postbuckling analysis for structural instability and the material topology for material function optimization. 43 Figure 12. (a) Drop-forces prior and posterior to buckling mode transitions for the FGM beams with the power material function. Distributions of the drop-forces in buckling modes Φ3, Φ5, Φ7 and the summation with respect to E2 and r for the FGM beams with the material functions defined in (b) power, (c) cosine square, (d) decimal logarithm, (e) natural exponential, and (f) sinusoidal functions (𝐿 = 200 𝑚𝑚, 𝑏 = 25 𝑚𝑚, 𝑡 = 1.82 𝑚𝑚 and 𝐸𝐴 = 3.85 𝐺𝑃𝑎 for all cases). 𝛥𝑝̂ < 𝜀; (3-23) 44 where, 𝜀 represent the desired total drop-forces. In particular, if the drop-forces are larger than the limit state, the releasing of the stored strain energy is too significant, and therefore, the topology algorithm is calculated again by updating the Young’s modulus 𝐸 and volume fraction 𝑟. The loop keeps running until Eq. (3-23) is met. Here, we calculate the material functions as the power, cosine squared, decimal logarithm, natural exponential, and sinusoidal functions, as defined in Eq. (3-2). Note that the length, width, thickness, and the Young’s modulus of the stiffer material are fixed as: 𝐿 = 200 mm, 𝑏 = 25 mm, 𝑡 = 1.82 mm and 𝐸𝐴 = 3.85 GPa , respectively. The optimization variables are the volume fraction 𝑟 and the Young’s modulus of the softer material 𝐸𝐵 , which are varied as 0 ≤ 𝑟 ≤ 1 and 1.65 ≤ 𝐸𝐵 ≤ 3.85 GPa, respectively. Figure 12 investigates the influences of the volume fraction 𝑟 and the Young’s modulus of the softer material 𝐸𝐵 on the releasing of the stored strain energy in terms of the drop-forces. Following Eq. (2), the material functions of the FGM beams are defined as Figure 12(b) power, Figure 12(c) cosine squared, Figure 12(d) decimal logarithm, Figure 12(e) natural exponential, and Figure 12(f) sinusoidal functions. The results show that varying 𝐸𝐵 and 𝑟 results in significant changes in the drop-forces at each buckling mode transition. Note that the cosine squared and decimal logarithm functions cause decreasing of the drop-forces when 𝐸𝐵 and 𝑟 are increased. Theses variation patterns are the opposite of the power, natural exponential and sinusoidal cases. Figure 13 compares the variation patterns for the material functions in the power, cosine squared, decimal logarithm, natural exponential, and sinusoidal functions. It can be seen that the cosine squared function is likely to be the most critical function to affect the releasing of the stored energy in the FGM beams, while the power, natural exponential and sinusoidal functions tend to affect the releasing at the same level. More interestingly, increasing 𝐸𝐵 and 𝑟 lead to the same 45 drop-forces for all the cases, which demonstrates that material functions can be used to affect the releasing of the stored energy, and eventually, tune the structural instability of the bi-walled FGM beams. According to the reported material topology, the material functions of the FGM beams can be controlled, and thus, optimized by programming 𝐸𝐵 and 𝑟 , which provides an effective approach to maneuver the structural instability for different applications. Figure 13. Total drop-forces posterior to buckling mode transitions for the FGM beams with power, cosine square, decimal logarithm, natural exponential, and sinusoidal functions. 3.6 Additional Studies Table 4 presents the values in the material distribution functions used to define the programmable Young’s modulus, and the obtained axial forces for the postbuckling snap-through in the 3rd, 5th and 7th modes in Section 4. 46 Table 4. Material distribution functions in the theoretical modelling and the axial force for the postbuckling snap-through in the 3rd, 5th and 7th modes (all in mm and N). 𝐹5 𝐹7 𝐸(𝑥) r L 𝐹3 (Φ3) (Φ5) (Φ7) Natural exp 0.5 185 453.22 1131.81 2341.46 Decimal log 0.5 185 329.37 822.44 1701.43 0.25 185 233 583.57 1207.27 0.5 185 252.14 629.63 1302.57 160 365.98 917.7 1700.48 185 275.21 687.21 1412.68 1 190 258.56 650.99 1352.13 Power 220 194.41 488.4 987.37 250 150.81 375.57 782.09 2 185 298.27 744.79 1540.80 4 185 316.72 790.85 1636.09 Sinusoidal 1 185 285.29 749.92 1470.31 Cosine2 1 185 306.24 804.95 1591.01 In this study, a total of 30 buckling modes are taken into account in the present study mainly due to the considerations of: 1) accuracy; and 2) computational cost. Figure 14 compares the theoretical force-displacement relations of the bi-walled isotropic beam with 10, 20 and 30 buckling modes and the experiments. The 30-mode model offers the most accurate theoretical results comparing with the experimental results, especially for the postbuckling mode transitions 47 (i.e., Φ 3, Φ 5 and Φ 7). On the other hand, the computational cost of solving the 30-mode theoretical model is approximately two times of that for the 20-mode model, as listed in Table 5. The personal computer used to solve the theoretical model consists of the Intel(R) Core (TM) i7- 4770T CPU @ 2.5 GHz processor, 8 GB installed memory (RAM) and 64-bit operating system. Figure 14. Comparison of the postbuckling responses of the bilaterally constrained beam obtained using the 10-, 20- and 30-mode theoretical model and experiments. Table 5. Buckling mode vs. computing time. Buckling mode Time (min) 10 158 20 638 30 1249 48 3.7 Conclusions This study reported on a design and optimization approach for bilaterally constrained FGM beams with programmed material functions aimed at controlling structural instabilities (i.e., postbuckling). Theoretical and numerical models were developed to characterize the postbuckling response, and the results were compared with experimental observations. Satisfactory agreements were obtained between the theoretical, numerical and experimental results. Parametric studies were carried out to demonstrate the tunability of the reported FGM beams with respect to the geometric and material parameters (i.e., beam length L, Young's modulus E and volume ratio r). The material topology was then used to investigate the variation of total potential energy (i.e., releasing of the stored energy) in buckling mode transitions. The drop-forces were used to demonstrate the influences of the Young's modulus and volume ratio on the variation of the potential energy. It was found that the cosine squared function is the most critical function to affect the energy release in the FGM beams, while the power, natural exponential and sinusoidal functions affected the release of energy at a similar level. As a consequence, it was concluded that material distributions functions can be controlled based on the release of the stored energy, which allows for well-defined controllability over the structural instability of the FGM beams. The reported bi-walled FGM beams offer a novel approach to program structural instability in functional composites for multipurpose applications. 49 Chapter 4: Maneuverable Postbuckling of Extensible Mechanical Metamaterials Using Functionally Graded Materials and Carbon Nanotubes 4.1 Introduction Mechanical metamaterials (MM) provide the possibility of improving the mechanical properties of bulk materials from the structural perspective. Utilizing various engineered microstructures [23], MM have been reported on their outstanding mechanical characteristics resulted in the critical geometric nonlinearity. Studies have been conducted on controlling and optimizing the mechanical response of MM [72], such as the enhanced tension [73], programmable bending and buckling, compression-induced twisting [31], ultra-light and ultra-stiff [25, 74], negative Poisson’s ratio [5], programmable unusual bulk modulus and mass density, and complete recovery from large strains in compression [75-77], etc. Exploring and programming the mechanical characteristics in a desirable manner, MM have been used in various innovative ap- plications, e.g., mechanical sensor [78], heat transferring [79, 80], energy absorption [10, 81], artificial muscles [82], drug delivery [83, 84], and soft robots [85]. However, MM have been facing the dilemma in obtaining the better controllability in the mechanical response (i.e., tunability of performance) and less complexity in the microstructures (i.e., feasibility of fabrication). To address such research dilemma, functional composite materials can be used in MM, such as functionally graded materials (FGMs) and carbon nanotubes (CNTs). FGMs have been reported as the functional composite materials that the material properties are varied in the length, width, or thickness directions [86]. Studies have been conducted to explore the principles and applications of FGMs [87-89]. FGMs have been applied in many fields, including the reduction of the in-plane and lateral thickness stresses, improvement of the residual stress distribution, increment of the fracture toughness, reduction of the stress strength factor, and 50 enhancement of the thermoelectric energy conversion properties [90, 91]. Studies have been conducted to investigate the instability behavior (i.e., buckling) of FGMs at the multiscale. The buckling behavior of the Bernoulli-Euler nano-beams was studied under different boundary conditions [92], and the size-dependent response of the functionally graded nano-beams was investigated under the uniformly distributed transverse load [93]. The nonlocal integral formulation was proposed for the nano-beams in the non-isothermal environment by adopting the proper Gibbs potential [94]. CNTs, on the contrary, have been developed by eventually mixing carbon nanotubes in the matrix [95], which have been used in different applications [96-98]. Studies have shown that the addition of a small percentage of nanotubes are likely to considerably increase the mechanical, electrical, or thermal properties of CNTs [99-102]. The integral elasticity model was proposed to characterize the mechanical behavior of the single-walled CNTs under axial force and kinematic boundary using the molecular dynamics simulations [103]. Although the functionalities of FGMs and CNTs have been investigated by tuning the material properties of the composites, recent research interests have been dedicated to exploring their utilities using the method of structural design (e.g., MM). As a consequence, MM made by FGMs (i.e., FGM-MM) and CNTs (i.e., CNT-MM) are expected to behave more effective maneuverability with less complicated microstructures. Here, we develop the microscale, extensible FGM-MM, and CNT- MM by applying the structural design method (i.e., MM) to the functional composite materials (i.e., FGMs and CNTs). The FGM-MM and CNT-MM are designed in the plate-like configurations, which are placed between the bilateral constraints to obtain the postbuckling response under the axial compression. The material models are developed to characterize the effective material properties of the reported engineered metastructures, and the theoretical models are reported to investigate the postbuckling response (i.e., buckling mode transitions and force- 51 displacement relations). Numerical simulations are carried out to vali-date the theoretical results and satisfactory agreements are obtained. The programmable mechanical behaviors of the FGM- MM and CNT-MM are investigated in terms of the microstructures and material functions. This work combines the structural method (i.e., corrugation) and material method (i.e., FGMs and CNTs) to obtain the more maneuverable postbuckling response, which is novel since studies have yet been carried out on tuning the mechanical performance of MM using FGMs or CNTs. The rest of the chapter can be drawn as: Section 2 develops the material models to characterize the effective materials of the microscale FGM-MM and CNT-MM. Section 3 reports the theoretical models to investigate the postbuckling response of the FGM-MM and CNT-MM. Section 4 conducts the numerical simulations to validate the theoretical results of the postbuckling behaviors. Section 5 carries out the parametric studies to unveil the influences of the structure designs in the MM and material functions in the FGMs and CNTs on the postbuckling performance of the FGM-MM and CNT-MM. Section 6 summarizes the main findings of this study. 4.2 Problem Formulation of the Microscale FGM-MM and CNT-MM In this study, we theoretically formulate the FGM-MM and CNT-MM by considering the structural characteristics of MM and the material characteristics of FGMs and CNTs. Figure 15 illustrates the designs of the maneuverable, extensible FGM-MM and CNT-MM designed using the structural method of MM and material method of CNTs and FGMs. Due to the slenderness of the plate-like structures (i.e., thickness ≪ width), the material properties of FGMs are assumed to be distributed following the material functions in the length direction 𝑥. In contrast, the material properties of CNT-MM are evenly distributed (i.e., the material functions are not varied with 𝑥). Table 6 compares the material functions between the FGM-MM and CNT-MM. 52 Table 6. Comparison of the material functions between the FGM-MM and CNT-MM. FGM-MM CNT-MM Young’s modulus 𝐸 𝐸FM = 𝑓1 (𝑥)* 𝐸CM = 𝑓2 (𝑉) Shear modulus 𝐺 𝐺FM = 𝑔1 (𝑥) 𝐺CM = 𝑔2 (𝑉) * 𝑓𝑖 and 𝑔𝑖 (𝑖=1, 2) refer to the material functions, 𝑉 indicates the volume fraction in CNTs, and the subscripts FM and CM indicate the FGM-MM and CNT-MM, respectively. 53 Figure 15. Principles of the maneuverable, extensible FGM-MM and CNT-MM designed using the structural method of MM and material method of FGMs and CNTs. 4.3 Effective Materials Properties of the FGM-MM and CNT-MM 4.3.1 Structural Characteristics of MM Figure 16(a) demonstrates the postbuckling deformation for the extensible CNT-MM and FGM-MM subjected to bilateral constraints. According to [104], the effective stress-strain relations of plate MM can be written as 54 𝜌 𝟐 𝐸MM = 𝜁MM 𝐸 ( 𝜌MM ) ; (4-1) hb where 𝐸 is the Young’s modulus of the material in MM, and 𝜁MM , 𝜌MM and 𝜌hb are, respectively, the weight factor accounts for the influence of the cylindrical webbing patterns on the MM, density of the MM, and density of a hollow beam that has the equivalent length 𝐿, width W𝑏 , height H, thickness 𝑡MM and mass 𝑀 as the MM. The volume of the MM in FGM-MM and CNT-MM, and the volume of the hollow beam can be obtained as [105], therefore the effective Young’s modulus of the plate MM can be rewritten 2 2(𝐷+𝐺)2 (𝑊𝑏 +𝐻−2𝑡MM ) 𝐸MM = 𝜁MM [ ] 𝐸; (4-2) 𝑊𝑏 ((𝐷+𝐺)2 +𝜋𝐷𝐻) where 𝐷 and 𝐺 are the diameter and the distance between two cells in the cylindrical configuration. 4.3.2 Material Characteristics of FGMs and CNTs Assuming postbuckling deformation are resulted in the normal stress among the longitudinal direction (x-axis) 𝜎𝑥 and couple stress 𝑚𝑥𝑦 , the stress-strain relations of the FGM-MM and CNT- MM can be written as [106], 𝜎𝑥 𝑄 𝜀𝑥 𝜎 = {𝑚 } = 𝑄𝜀 = [ 11 ] {𝜒 } . (4-3) 𝑥𝑦 𝑄̂44 𝑥𝑦 Since the Poisson’s ratio can be neglected for the slender structures in this study [107], the components of the 𝑄 matrix can be reduced to [106]: 𝑄11 = 𝐸11 {̂ , (4-4) 𝑄44 = 𝐺13 𝑙 2 where 𝐸11 , 𝐺13 and 𝑙 are the Young’s modulus in CNT direction, shear modulus in 1-3 plan and 55 material length scale factor, respectively. Figure 16. (a) Demonstration of the postbuckling deformation for the extensible CNT-MM and FGM-MM subjected to bilateral constraints. (b) Deformation analysis of an arbitrary segment in the postbuckled CNT-MM and FGM-MM. First, let us define the effective material properties of the FGM-MM. The coefficients in Eq. (4-4) can be written as 1 𝑄11 = 𝐸(𝑥) and 𝐺13 = 2 𝐸(𝑥); (4-5) where 𝐸(𝑥) is the material functions of the FGMs. In this study, we particularly consider two material functions for FGMs as 𝐸(𝑥),1 = 𝛼(1 + 𝛽𝑥)𝛾 { 𝛽𝑥 ; (4-6) 𝐸(𝑥),2 = 𝛼𝑒 𝛾 where α, 𝛽 and 𝛾 are the factors used to define the Young’s modulus. Taking Eqs. (4-2) and (4-6) into Eq. (4-5), the applied coefficients in Eq. (4-4) can be written as 56 2 2(𝐷+𝐺)2 (𝑊𝑏 +𝐻−2𝑡MM ) 𝑄11,FM = 𝜁MM [ ] 𝐸(𝑥) 𝑊𝑏 ((𝐷+𝐺)2 +𝜋𝐷𝐻) { 2 ; (4-7) 1 2(𝐷+𝐺)2 (𝑊𝑏 +𝐻−2𝑡MM ) 𝑄̂44,FM = 𝜁 [ ] 𝐸(𝑥)𝑙 2 2 MM 𝑊𝑏 ((𝐷+𝐺)2 +𝜋𝐷𝐻) where the subscript FM refers to the FGM-MM. Next, let us determine the effective materials properties of the CNT-MM. The Young’s modulus of CNTs can be estimated by the Mori-Tanaka scheme or the rule of mixtures. To take into account the small scale effects, the CNTs efficiency parameters (𝜂1 , 𝜂3 ) are introduced in the extended version of this rule and evaluated by matching the effective properties of CNT-reinforced composites determined from the molecular dynamics simulations with those from the rule of mixture [108]. The Young's modulus and shear modulus of CNTs can be expressed as [109, 110]. 𝐴 𝐸11 = 𝜂1 𝑟𝑉𝐵 𝐸11 + (1 − 𝑟𝑉𝐵 )𝐸 𝐵 { 𝐴 𝐺𝐵 𝜂 𝐺13 ; (4-8) 𝐺13 = (𝑟𝐺𝐵3 +𝐺 𝐴 )𝑉 13 𝐵 𝐴 𝐴 where 𝐸11 , 𝐺13 , 𝐸 𝐵 , 𝐺 𝐵 are the CNTs’ Young’s modulus, shear modulus and matrix’s Young’s 𝑉 modulus and shear modulus, respectively, and 𝑟 = 𝑉𝐴 describes the variation in the volume 𝐵 fraction. Note that 𝐴 and 𝐵 refer to the CNTs and matrix, respectively. Substituting Eqs. (4-8) and (4-2) into Eq. (4-4), the effective material properties of the CNT-MM can be obtained by combining the Young’s modulus of CNTs with MM as 2 2(𝐷+𝐺)2 (𝑊𝑏 +𝐻−2𝑡MM ) 𝐴 𝑄11,CM = 𝜁MM [ ] (𝜂1 𝑟𝑉𝐵 𝐸11 + (1 − 𝑟𝑉𝐵 )𝐸 𝐵 ) 𝑊𝑏 ((𝐷+𝐺)2 +𝜋𝐷𝐻) { 2 𝐴 𝐺𝐵 ; (4-9) 2(𝐷+𝐺)2 (𝑊𝑏 +𝐻−2𝑡MM ) 𝜂3 𝐺12 𝑄̂44,CM = 𝜁MM [ 𝑊 ((𝐷+𝐺) 2 +𝜋𝐷𝐻) ] 𝐵 (𝑟𝐺 +𝐺12 𝐴 )𝑉 𝑙 2 𝑏 𝐵 where the subscript CM refers to the CNT-MM. 57 4.4 Postbuckling Analysis of the FGM/CNT-MM Subjected to Bilateral Constraints 4.4.1 General Governing Equations According to the modified couple stress theory, the strain energy density can be written as a function of the strain (i.e., conjugated with stress) and curvature (i.e., conjugated with couple stress) as 1 𝑈 = 2 ∭(𝜎𝜀 + 𝑚𝜒)d𝑉 ; (4-10) where 𝜎, 𝜀, 𝑚 and 𝜒 are the Cauchy stress tensor, conventional strain tensor, deviatoric couple stress tensor and symmetric curvature tensor, respectively. Based on the Euler-Bernoulli beam theory as shown in Figure 16(b), the displacement components of the FGM-MM and CNT-MM can be obtained as 𝑢𝑥 = −𝑢(𝑥, 𝜏) − 𝑧𝜙(𝑥, 𝜏) 𝑢𝑦 = 0 𝑢𝑧 = 𝑤(𝑥, 𝑡) ; (4-11) 𝜕𝜙 { 𝜒𝑥𝑦 = − 𝜕𝑥 𝜕𝑤 where 𝜙 = , 𝜏 and 𝑡 are the rotation angle of the neutral axis, time term and the thickness, 𝜕𝑥 respectively. 𝑢 and 𝑤 are the axial displacement and transverse deflection, respectively. Thus, the nonzero components of the strain curvature are obtained as 𝜕𝑢 𝜕2 𝑤 𝜀𝑥 − −𝑧 2 𝜀 = {𝜒 } = { 𝜕𝑥 𝜕2𝑤 𝜕𝑥 }. (4-12) 𝑥𝑦 − 𝜕𝑥 2 According to the Hamilton’s principle, the total energy can be express as 𝑇1 𝑇1 ∫𝑇2 ẟ𝛺d𝜏 = ∫𝑇2 (ẟ𝐾 − ẟ𝑈 + ẟ𝑊)d𝜏 = 0; (4-13) 58 where ẟ𝑈, ẟ𝑊 and ẟ𝐾 are the virtual potential energy, virtual external work, and virtual kinetic energy, respectively, which can be evaluated as 𝐿′ 𝑏 𝑡 ẟ𝑈 = ∫0 ∫0 [∫0 𝜎 𝑇 ẟ𝜀 + 𝑚𝑇 ẟ𝜒d𝑧] d𝑦d𝑥 ; (4-14) 𝐿′ 𝜕𝜙 𝑑𝑢 𝑥 = 𝐿′ ẟ𝑊 = 𝑃 ∫0 [sin 𝜙 𝜕𝑥 ẟ𝑢 + sin 𝜙 (1 − 𝑑𝑥 ) ẟ𝜙]d𝑥 + 𝑃 cos 𝜙 ẟ𝑢 | ; (4-15) 𝑥=0 and 𝐿′ 𝑏 𝑡 𝜕𝑢 𝜕𝑢 𝜕2𝑤 𝜕2𝑤 𝜕𝑤 𝜕𝑤 ẟ𝐾 = ∫0 ∫0 [∫0 𝜌 ( ) ẟ ( ) + 𝑧 2 ( )ẟ( ) + ( ) ẟ ( ) d𝑧] d𝑦d𝑥 . (4-16) 𝜕𝜏 𝜕𝜏 𝜕𝑥𝜕𝜏 𝜕𝑥𝜕𝜏 𝜕𝜏 𝜕𝜏 Note that 𝑃 and 𝐿′ are the axial force and beam length after deformation, respectively. The force, moments, unit mass, and moment of inertia of the FGM-MM and CNT-MM can be written as 𝑏 𝑡 𝑧2 {𝑁𝑥 , 𝑀𝑥 , 𝑀𝑥𝑦 , 𝑚, 𝐼} = ∫0 ∫0 {𝜎𝑥 , 𝑧𝜎𝑥 , 𝑚𝑥𝑦 , 𝜌, 𝑦 } d𝑧d𝑦. (4-17) Substituting Eq. (4-17) into Eq. (4-12), we have 𝜕𝑢 2 𝑁𝑥 = 𝑄̅11 𝜕𝑥 − 𝐽11 ̅ 𝜕 𝑤2 𝜕𝑥 2 𝑀𝑥 = 𝐽11 ̅ 𝜕𝑢 − 𝐼11 ̅ 𝜕 𝑤2 ; (4-18) 𝜕𝑥 𝜕𝑥 𝜕2 𝑤 { 𝑀𝑥𝑦 = −𝑄̿44 𝜕𝑥 2 where 𝑄̅11, 𝐽11 ̅ , 𝐼11̅ and 𝑄 can be expressed as 44 59 𝑡 𝑏 𝑄̅11 = ∫0 ∫2𝑡 𝑄11 d𝑧d𝑦 − 2 𝑡 𝑏 ̅ = ∫ ∫ 𝑡 𝑧𝑄11 d𝑧d𝑦 𝐽11 2 0 − 2 𝑡 . (4-19) ̅ = ∫𝑏 ∫ 𝑡 (𝑧)2 𝑄11 d𝑧d𝑦 𝐼11 2 0 − 2 𝑡 𝑏 𝑄44 = ∫0 ∫−2𝑡 𝑄̂44 d𝑧d𝑦 { 2 Taking Eqs. (4-14) to (4-17) into Eq. (4-13) and integrating by parts, we obtain 𝑇1 𝑇 2 𝐿′ 𝜕𝜙 𝜕𝑁𝑥 𝜕2𝑢 d𝑢 𝜕2 𝑤 𝜕2 𝑀𝑥 𝜕2 𝑀𝑥𝑦 ∫𝑇2 ẟ𝛺d𝜏 = ∫𝑇1 ∫0 [(𝑃 sin 𝜙 𝜕𝑥 − −𝑚 2 ) ẟ𝑢 + (−𝑃 (1 − ) + + + 𝜕𝑥 𝜕𝑡 d𝑥 𝜕𝑥 2 𝜕𝑥 2 𝜕𝑥 2 2 𝜕2𝑤 𝜕2𝑤 𝑇 𝜕𝑢 d𝑢 𝑚𝐼 (𝜕𝑥𝜕𝜏) − 𝑚 𝜕𝜏2 ) ẟ𝑤] d𝑥d𝜏 + ∫𝑇 2 [(𝑃 cos 𝜙 + 𝑁𝑥 + 𝑚 𝜕𝜏 ) ẟ𝑢 + (𝑃 sin 𝜙 (1 − d𝑥 ) − 1 𝜕𝑀𝑥 𝜕𝑀𝑥𝑦 𝜕𝑤 𝜕3𝑤 𝜕2 𝑤 𝜕𝑤 𝑥 = 𝐿′ − +𝑚 − 𝑚𝐼 𝜕𝑥𝜕𝜏2) ẟ𝑤 + (𝑀𝑥 + 𝑀𝑥𝑦 + 𝐼 𝜕𝑥𝜕𝜏) ẟ ( 𝜕𝑥 )] | d𝜏 = 0 ; (4-20) 𝜕𝑥 𝜕𝑥 𝜕𝜏 𝑥=0 which leads to the governing equation: 𝜕𝜙 𝜕𝑁𝑥 𝜕2𝑢 𝑃 sin 𝜙 𝜕𝑥 − − 𝑚 𝜕𝑡 2 = 0 𝜕𝑥 { 2 ; (4-21) d𝑢 𝜕2 𝑤 𝜕2 𝑀𝑥 𝜕2 𝑀𝑥𝑦 𝜕2𝑤 𝜕2𝑤 −𝑃 (1 − d𝑥 ) 𝜕𝑥 2 + 𝜕𝑥 2 + 𝜕𝑥 2 + 𝑚𝐼 (𝜕𝑥𝜕𝜏) − 𝑚 𝜕𝜏2 = 0 and the boundary conditions 𝜕𝑢 𝑃 cos 𝜙 + 𝑁𝑥 + 𝑚 𝜕𝜏 = 0 d𝑢 𝜕𝑀𝑥 𝜕𝑀𝑥𝑦 𝜕𝑤 𝜕3𝑤 𝑥 = 𝐿′ 𝑃 sin 𝜙 (1 − d𝑥 ) − − + 𝑚 − 𝑚𝐼 = 0 | . (4-22) 𝜕𝑥 𝜕𝑥 𝜕𝜏 𝜕𝑥𝜕𝜏2 𝑥=0 2 𝜕 𝑤 { 𝑀𝑥 + 𝑀𝑥𝑦 + 𝐼 𝜕𝑥𝜕𝜏 = 0 } Substituting Eqs. (4-18) and (4-19) into Eqs. (4-21) and (4-22), the governing equation and boundary conditions are obtained as 60 𝜕𝜙 𝜕2𝑢 𝜕3𝑤 𝜕2𝑢 𝑃 sin 𝜙 𝜕𝑥 − 𝑄11 𝜕𝑥 2 + 𝐽11 𝜕𝑥 3 − 𝑚 𝜕𝜏2 = 0 { 2 ; (4-23) d𝑢 𝜕2 𝑤 𝜕3𝑢 𝜕4𝑤 𝜕2𝑤 𝜕2𝑤 𝑃 (1 − d𝑥 ) 𝜕𝑥 2 − 𝐽11 𝜕𝑥 3 + (𝐼11 + 𝑄44 ) 𝜕𝑥 4 − 𝑚𝐼 (𝜕𝑥𝜕𝜏) + 𝑚 𝜕𝜏2 = 0 and 𝜕𝑢 𝜕2𝑤 𝜕𝑢 𝑃 cos 𝜙 + 𝑄11 𝜕𝑥 − 𝐽11 𝜕𝑥 2 + 𝑚 𝜕𝜏 = 0 d𝑢 𝜕2𝑢 𝜕3𝑤 𝜕𝑤 𝜕3𝑤 𝑥 = 𝐿′ 𝑃 sin 𝜙 (1 − d𝑥 ) − 𝐽11 𝜕𝑥 2 + (𝐼11 + 𝑄44 ) 𝜕𝑥 3 + 𝑚 𝜕𝜏 − 𝑚𝐼 𝜕𝑥𝜕𝜏2 = 0 | ; (4-24) 𝑥=0 𝜕𝑢 𝜕2𝑤 𝜕2𝑤 { 𝐽 ̅ 11 − (𝐼 ̅ 11 + 𝑄̿ 44 ) 2 + 𝐼 =0 } 𝜕𝑥 𝜕𝑥 𝜕𝑥𝜕𝜏 𝑑𝑢 where 𝛬 = 𝑑𝑥 is the extensible factor resulted in the cylindrical corrugation of the FGM-MM and CNT-MM. It is worthwhile to point out that Eqs. (4-23) and (4-24) are the general governing equations for the extensible FGM-MM and CNT-MM. Neglecting the extensible terms (i.e., 𝛬 = 0), Eq. (4- 23) and (4-24) are reduced to: 𝜕𝜙 𝜕2𝑢 𝜕3𝑤 𝜕2𝑢 𝑃 sin 𝜙 𝜕𝑥 − 𝑄11 𝜕𝑥 2 + 𝐽11 𝜕𝑥 3 − 𝑚 𝜕𝜏2 = 0 { 2 ; (4-25) 𝜕3𝑢 𝜕4𝑤 𝜕2𝑤 𝜕2𝑤 𝜕2𝑤 −𝐽11 𝜕𝑥 3 + (𝐼11 + 𝑄44 ) 𝜕𝑥 4 + 𝑃 𝜕𝑥 2 − 𝑚𝐼 (𝜕𝑥𝜕𝜏) + 𝑚 𝜕𝜏2 = 0 and 𝜕𝑢 𝜕2𝑤 𝜕𝑢 𝑃 cos 𝜙 + 𝑄11 𝜕𝑥 − 𝐽11 𝜕𝑥 2 + 𝑚 𝜕𝜏 = 0 𝜕2𝑢 𝜕3𝑤 𝜕𝑤 𝜕𝑤 𝜕3𝑤 𝑥 = 𝐿′ −𝐽11 + (𝐼11 + 𝑄44 ) +𝑃 +𝑚 − 𝑚𝐼 =0 | . (4-26) 𝜕𝑥 2 𝜕𝑥 3 𝜕𝑥 𝜕𝜏 𝜕𝑥𝜕𝜏2 𝑥=0 𝜕2𝑤 𝜕2𝑤 { ̅ 𝜕𝑢 − (𝐼11 𝐽11 ̅ + 𝑄̿44 ) 2 + 𝐼 =0 } 𝜕𝑥 𝜕𝑥 𝜕𝑥𝜕𝜏 Further neglecting the time terms, the governing equations are the same as the static micro- composite beams model proposed in [111]. 61 In order to compare with the micro-isotropic beams in the existing study, the coefficients for the FGMs and CNTs materials are reduced to isotropic as 𝑄̅11 = 𝐸𝑏𝑡 ̅ =0 𝐽11 𝐸𝑏𝑡 3 . (4-27) ̅ = 𝐼11 12 {𝑄44 = 𝑏𝑡𝐺𝑙 2 Taking Eq. (4-27) into Eq. (4-25), the governing equation for the micro-isotropic beams is 2 𝜕4𝑤 𝜕2𝑤 𝜕2𝑤 𝜕2𝑤 (𝐸𝐼 + 𝐺𝐴𝑙 2 ) + 𝑃 𝜕𝑥 2 − 𝑚𝐼 (𝜕𝑥𝜕𝜏) + 𝑚 𝜕𝜏2 = 0. (4-28) 𝜕𝑥 4 2 𝜕2 𝑤 Omitting the higher-order rotation term (𝜕𝑥𝜕𝜏) , Eq. (4-28) can be simplified to 𝜕4𝑤 𝜕2𝑤 𝜕2𝑤 𝐵 𝜕𝑥 4 + 𝑃 𝜕𝑥 2 + 𝑚 = 0; (4-29) 𝜕𝜏2 where 𝐵 = 𝐸𝐼 + 𝐺𝐴𝑙 2 is the bending stiffness. Eq. (4-29) is the same as the governing equation of the micro-isotropic beams in [112]. Further dropping the length scale factor 𝑙 and time related 𝜕2𝑤 term , the governing equation for classical Euler-Bernoulli beams is obtained. 𝜕𝜏2 4.4.2 Postbuckling Analysis of the FGM-MM and CNT-MM To solve the governing equation given in Eq. (4-23), the following assumptions are taken into account in this section: 1) The FGM-MM and CNT-MM are assumed to be symmetric about the neutral axis, and therefore, 𝐽11 = 0. 2) The axial force is gradually applied to the end, which lead to quasi-static postbuckling response 𝜕2𝑤 𝜕2 𝑤 of the FGM-MM and CNT-MM, and therefore, 𝜕𝑥𝜕𝜏 = 𝑚 𝜕𝜏2 = 0. 62 According to the assumptions, the governing equation in Eq. (4-23) can be reduced to 𝜕4 𝑤 𝜕2 𝑤 (𝐼11 + 𝑄44 ) 𝜕𝑥 4 + (𝑃 − 𝑃𝛬) 𝜕𝑥 2 = 0; (4-30) and the boundary conditions are 𝑥 = 𝐿′ d𝑤(𝑥) 𝑥 = 𝐿′ 𝑤(𝑥) | = d𝑥 | = 0; (4-31) 𝑥=0 𝑥=0 𝜕2 𝑤 where 𝑃𝛬 𝜕𝑥 2 is the extensibility term caused by the cylindrical corrugation. The extensible factor can be expressed as [107]. 𝑃 𝛬 = 𝑄̅ cos 𝜙 ; (4-32) 11 where 𝜙 is the rotation angle of the FGM-MM and CNT-MM resulted in the postbuckling deformation. For small deformations, 𝜙 ≈ 0, and cos 𝜙 ≈ 1. 4.4.2.1 Postbuckling Analysis of the FGM-MM The material factors of the FGM-MM are given as 𝑡 2 𝐸FGM−MM (𝑥) 𝐸FGM−MM (𝑥)𝑙2 {𝑄11 , 𝐽11 , 𝐼11 , 𝑄44 } = 𝑏𝑡 {𝐸FM (𝑥), 0, , }; (4-33) FM 12 2 and therefore, the governing equation for the postbuckling behavior of the FGM-MM can be rewritten as 1 𝜕4𝑤 𝜕2𝑤 [𝐸FM (𝑥)𝐼 + 2 𝐸FM (𝑥)𝐴𝑙 2 ] 𝜕𝑥 4 + (𝑃 − 𝑃𝛬) 𝜕𝑥 2 = 0. (4-34) 𝑥 𝑤(𝑋𝐿′ ) Introducing the non-dimensional variables 𝑋 = 𝐿′ and 𝑊(𝑋) = , Eq. (4-34) can be ℎ normalized as 63 𝜕4𝑊 𝜕2𝑊 + 𝑁(𝑋) 𝜕𝑋 2 = 0; (4-35) 𝜕𝑋 4 where the correlated normalized axial load 𝑁(𝑋) is obtained as 2 2𝐿′ (𝑃−𝑃𝛬) 𝑁(𝑋) = (2𝐼+𝐴𝑙2 )𝐸 . (4-36) FM Taking Eq. (4-7) into Eq. (4-36), 𝑁(𝑋) can be rewritten as 2 2 (𝑃−𝑃𝛬)𝑏 2 ((𝐷+𝐺)2 +𝜋𝐷𝐻) 𝐿′ 𝑁(𝑋) = 2𝜁 2 )(𝐷+𝐺)4 (𝑊 2 𝛼(1+𝛽𝑋)𝛾 . (4-37) MM (2𝐼+𝐴𝑙 𝑏 +𝐻−2𝑡MM ) According to [21], the second-order general solution of Eq. (4-35) can be written as d2 𝑊(𝑋) = 𝑘1 𝑆1 (𝑋) + 𝑘2 𝑆2 (𝑋); (4-38) d𝑋 2 where the integral functions are expressed as 1 2 𝑆1 (𝑋) = √1 + 𝛽𝑋𝑱 1 (𝑛 (1 + 𝛽𝑋)4−2𝛾 ) 2−𝛾 2−𝛾 1 ; (4-39) 2 𝑆 (𝑋) = √1 + 𝛽𝑋𝑱 1 (𝑛 2−𝛾 (1 + 𝛽𝑋) 4−2𝛾 ) { 2 𝛾−2 in which 𝑱 indicates the first Bessel functions, and 𝑛 is 𝑏[(𝐷+𝐺)2 +𝜋𝐷𝐻]𝐿′ 2(𝑃−𝑃𝛬) 𝑛 = 2𝛽(𝐷+𝐺)2(𝑊 √(2𝐼+𝐴𝑙2 )𝜁 . (4-40) 𝑏 +𝐻−2𝑡MM ) MM 𝛼 Integrating Eq. (4-39), we obtain 1 d𝑊(𝑋) 𝑋 2 𝑋 2 = 𝑘1 ∫0 √1 + 𝛽𝜉𝑱 1 (𝑛 (1 + 𝛽𝜉)4−2𝛾 ) d𝜉 + 𝑘2 ∫0 √1 + 𝛽𝜉𝑱 1 (𝑛 (1 + d𝑋 2−𝛾 2−𝛾 𝛾−2 2−𝛾 1 𝛽𝜉)4−2𝛾 ) d𝜉 + 𝑘3 ; (4-41) and 64 1 𝑋 𝑡 2 𝑊(𝑋) = 𝑘1 ∫0 ∫0 √1 + 𝛼𝜉𝑱 1 (𝑛 (1 + 𝛽𝜉)4−2𝛾 ) d𝜉d𝑡 + 2−𝛾 2−𝛾 1 𝑋 𝑡 2 𝑘2 ∫0 ∫0 √1 + 𝛽𝜉𝑱 1 (𝑛 (1 + 𝛽𝜉)4−2𝛾 ) d𝜉d𝑡 + 𝑘3 𝑋 + 𝑘4 . (4-42) 𝛾−2 2−𝛾 where 𝑘𝑖 (𝑖 = 1, … ,4) are the integral constants. Taking Eqs. (4-41) and (4-42) into Eq. (4-31), we have 𝑘3 = 0 𝑘4 = 0 1 1 1 2 1 2 𝑘1 ∫0 √1 + 𝛽𝜉𝑱 1 (𝑛 (1 + 𝛽𝜉)4−2𝛾 ) d𝜉 + 𝑘2 ∫0 √1 + 𝛽𝜉𝑱 1 (𝑛 (1 + 𝛽𝜉)4−2𝛾 ) d𝜉 + 𝑘3 = 0 . 2−𝛾 2−𝛾 𝛾−2 2−𝛾 1 1 1 𝑡 2 1 𝑡 2 𝑘 ∫ ∫ √1 + 𝛼𝜉𝑱 1 (𝑛 (1 + 𝛽𝜉) 4−2𝛾 ) d𝜉d𝑡 + 𝑘2 ∫0 ∫0 √1 + 𝛽𝜉𝑱 1 (𝑛 (1 + 𝛽𝜉) 4−2𝛾 ) d𝜉d𝑡 + 𝑘3 + 𝑘4 = 0 { 1 0 0 2−𝛾 2−𝛾 𝛾−2 2−𝛾 (4-43) Eq. (4-43) represents the eigenvalue problem for 𝑛𝑖 , which has the characteristic equation as 1 1 1 2 1 𝑡 2 ∫0 √1 + 𝛽𝜉𝑱 1 (𝑛 (1 + 𝛽𝜉)4−2𝛾 ) d𝜉 ⋅ ∫0 ∫0 √1 + 𝛽𝜉𝑱 1 (𝑛 (1 + 𝛽𝜉)4−2𝛾 ) d𝜉d𝑡 − 2−𝛾 2−𝛾 𝛾−2 2−𝛾 1 1 1 2 1 𝑡 2 ∫0 √1 + 𝛽𝜉𝑱 1 (𝑛 (1 + 𝛽𝜉)4−2𝛾 ) d𝜉 ⋅ ∫0 ∫0 √1 + 𝛽𝜉𝑱 1 (𝑛 (1 + 𝛽𝜉)4−2𝛾 ) d𝜉d𝑡 = 0. 𝛾−2 2−𝛾 2−𝛾 2−𝛾 (4-44) A set of value of 𝑛𝑖 can be obtained by numerically solving Eq. (4-44). Assuming 𝛼 = 2300, 𝛽 = 1.2 and 𝛾 = 0.3 and using the first material function in Eq. (4-6), the general form of the buckling shape function is obtained as 𝑋 𝑡 𝑊(𝑋) = ∑∞ 𝑗=1 𝐴𝑗 ∫0 ∫0 √1 + 1.2𝜉[𝑱0.588 (𝑛𝑗 1.176(1 + 1.2𝜉) 0.294 )+ 𝑘𝑗 𝑱−0.588 (𝑛𝑗 1.176(1 + 1.2𝜉)0.294 )]d𝜉d𝑡 ; (4-45) 65 where 𝐴𝑗 are the coefficients to determine the contributions of the buckling modes to the postbuckling shape functions of the FGM-MM. 𝑛𝑗 and 𝑘𝑗 are determined as 𝑛𝑗 = 1.78𝜋, 2.55𝜋, 3.56𝜋, 4.38𝜋 …, and 𝑘𝑗 = 1.80, −14.10, 2.03, 0.79 … for 𝑗 = 1, 2, 3, 4 …. In the same manner, assuming 𝛼 = 2300, 𝛽 = 1 , and 𝛾 = 2 and using the second material function in Eq. (4-6), the postbuckling shape functions of the FGM-MM can be obtained as 𝑋 𝑡 𝑊(𝑋) = ∑∞ 𝑗=1 𝐴𝑗 ∫0 ∫0 𝑱0 (2𝑛𝑗 𝑒 −0.5𝜉 ) − 𝑘𝑗 𝒀0 (2𝑛𝑗 𝑒 −0.5𝜉 )d𝜉d𝑡 ; (4-46) where 𝒀 indicates the second Bessel functions, and 𝑛𝑗 and 𝑘𝑗 are obtained as 𝑛𝑗 = 2.58𝜋, 3.56𝜋, 5.10𝜋, 6.21𝜋 …, and 𝑘𝑗 = −0.81, −0.06, −0.36, 0.86 … for 𝑗 = 1, 2, 3, 4 … 4.4.2.2 Postbuckling Analysis of the CNT-MM In this section, the postbuckling behavior of the CNT-MM is analyzed. The governing equation can be expressed as 𝜕4𝑊 𝜕2𝑊 + 𝑁2 = 0. (4-47) 𝜕𝑋 4 𝜕𝑋 2 where the correlated normalized axial load and coefficients of the CNT-MM are 2 𝑃(1−𝛬)𝐿′ 𝑁2 = ; (4-48) 𝐼11 CM +𝑄44 CM and 𝑡 2 𝑄11,CM {𝑄11 , 𝐽11 , 𝐼11 , 𝑄44 } = 𝑏𝑡 {𝑄11,CM , 0, , 𝑄̂44,CM 𝑙2 }. (4-49) CM 12 The general solution to Eq. (4-47) is: d2 𝑊(𝑋) = 𝑘1 sin 𝑁𝑋 + 𝑘2 cos 𝑁𝑋 ; (4-50) d𝑋 2 66 where 𝑘1 and 𝑘2 are the constants. Integrating Eq. (4-50), we have d𝑊(𝑋) cos 𝑁𝑋 sin 𝑁𝑋 = −𝑘1 + 𝑘2 + 𝑘3 ; (4-51) d𝑋 𝑁 𝑁 and sin 𝑁𝑋 cos 𝑁𝑋 𝑊(𝑋) = −𝑘1 − 𝑘2 + 𝑘3 𝑋 + 𝑘4 ; (4-52) 𝑁2 𝑁2 where 𝑘𝑖 (𝑖 = 1,2,3,4) indicate the integral constants. Using the boundary conditions in Eq. (4- 31), the relations of the constants can be determined as 𝑘 − 𝑁1 + 𝑘3 = 0 cos 𝑁 sin 𝑁𝑋 −𝑘1 + 𝑘2 + 𝑘3 = 0 𝑁 𝑁 𝑘 . (4-53) − 𝑁22 + 𝑘4 = 0 sin 𝑁 cos 𝑁 {−𝑘1 𝑁2 − 𝑘2 𝑁2 + 𝑘3 + 𝑘4 = 0 Expressing 𝑘1 , 𝑘2 and 𝑘3 in terms of 𝑘4 , the eigenvalues and buckling shape functions 𝑦𝑗 of the CNT-MM can be numerically obtained. Thereafter, the general form of the postbuckling shape function are obtained as 𝑁 sin 𝑁 sin 𝑁𝑋 𝑦(𝑋) = 1 − cos 𝑁𝑋 + 1−cos 𝑁 ( − 𝑋). (4-54) 𝑁 Eq. (4-54) can be expressed in terms of the symmetric and asymmetric terms as [46]. 𝑦𝑗 (𝑋) = 1 − cos 𝑁𝑗 𝑋 } 𝑗 = 1,3,5, … (4-55) 𝑁𝑗 = (𝑗 + 1)𝜋 and 2 sin 𝑁𝑗 𝑋 𝑦𝑗 (𝑋) = 1 − 2𝑋 − cos 𝑁𝑗 𝑋 + 𝑁𝑗 } 𝑗 = 2,4,6, … (4-56) 𝑁𝑗 = 2.86𝜋, 4.92𝜋, 6.94𝜋, 8.95𝜋, … 67 Using the superposition theory, the postbuckling shape functions of the CNT-MM can be obtained as 2 sin 𝑁𝑗 𝑋 𝑊(𝑋) = ∑∞ ∞ 𝑗=1,3,5… 𝐴𝑗 (1 − cos 𝑁𝑗 𝑋) + ∑𝑗=2,4,6… 𝐴𝑗 (1 − 2𝑋 − cos 𝑁𝑗 𝑋 + ); 𝑁𝑗 (4-57) where 𝐴𝑗 are the coefficients to determine the contributions of the buckling modes to the postbuckling shape functions. 4.4.3 Energy Method to Solve the Postbuckling Response In this section, the postbuckling behaviors of the FGM-MM and CNT-MM are resolved by determining weight coefficients 𝐴𝑗 using the energy method. In particular, 𝐴𝑗 are obtained by minimizing the total potential energy of the postbuckled FGM-MM and CNT-MM. The total potential energy of the FGM-MM and CNT-MM consists of the bending strain energy Ub , compressive strain energy Uc and external work Up . The total strain energy yields 2 1 𝐿′ 𝜕2𝑤 1 U = Ub + Uc = 2 ∫0 (𝐼11 + 𝑄44 ) ( 𝜕𝑥 2 ) d𝑥 + 2 𝑃∆𝑐 ; (4-58) and the external work caused by the axial force 𝑃 is 1 Up = − 2 𝑃; (4-59) where the longitudinal displacement ∆, axial shortening ∆𝑐 and variation of length ∆𝑏 resulted in the buckling-induced midplane rotation are ∆= ∆𝑐 + ∆𝑏 ∆𝑐 = 𝐿 − 𝐿′ { . (4-60) 1 𝐿′ 𝜕𝑤 2 ∆𝑏 = ∫ 2 0 ( 𝜕𝑥 ) d𝑥 68 As a consequence, the total potential energy 𝚷 can be written as 𝚷 = Ub + Uc + Up . (4-61) The normalized length after deformation can be written as ℎ2 1 𝜕𝑊 2 𝐿′ = 𝐿 − ∆𝑐 = 𝐿 − ∆ + 2𝐿′ ∫0 ( 𝜕𝑋 ) d𝑋 . (4-62) Considering an arbitrary segment in the FGM-MM and CNT-MM, the shortening of the element might be written as 𝑃d𝑥 d∆𝑐 = . (4-63) 𝑄11 Integrating Eq. (4-63) and substituting into Eq. (4-60), the axial force can be written as 1 ℎ2 1 𝜕𝑊 2 𝑃=𝐸 [∆ − 2𝐿′ ∫0 ( 𝜕𝑋 ) d𝑋 ] ; (4-64) 𝑚 𝐿 where 𝐸𝑚 is written as 1 1 𝐸𝑚 = ∫0 d𝑋. (4-65) 𝑄11 Substituting Eq. (4-64) into Eqs. (4-58) and (4-59), the normalized potential energies and axial force are expressed, in terms of the normalized axial displacement 𝑑, as 2 1 2 1 𝜕2 𝑦𝑗 (𝑋) 𝑉𝑏 = 2𝐼 [∑∞ 𝑗=1 𝐴𝑗 ∫0 (𝐼11 + 𝑄44 ) ( ) d𝑋] ; (4-66) 𝑚 𝜕𝑋 2 2 2 𝑘ℎ2 1 1 𝜕𝑦 (𝑋) 𝑉𝑐 = 2𝐸 [𝑑 − ∑∞ 𝐴 2 ∫0 ( 𝑗 ) d𝑋] ; (4-67) 𝑚 𝐼𝑚 2 𝑗=1 𝑗 𝜕𝑋 2 𝑘ℎ2 𝑑 2 1 𝜕𝑦𝑗 (𝑋) 𝑉𝑃 = − 2𝐸 [𝑑2 − 2 ∑∞𝑗=1 𝐴𝑗 ∫0 ( ) d𝑋 ] ; (4-68) 𝑚 𝐼𝑚 𝜕𝑋 69 and 𝑘ℎ 2 1 2 1 𝜕𝑦𝑗 (𝑋) 2 𝑃̃ = 𝐸 𝐼 [𝑑 − 2 ∑∞ 𝑗=1 𝐴𝑗 ∫0 ( 𝜕𝑋 ) d𝑋 ] ; (4-69) 𝑚 𝑚 where the variables and normalization factors are introduced as 2 3 3 3 ∆𝐿′ 𝐿′ 𝑡 𝑃𝐿′ 𝑈𝑏 𝐿′ 𝑈𝑐 𝐿′ 𝑈𝑃 𝐿′ 1 𝑑= , 𝑘 = 𝐿 , 𝑄 = ℎ, 𝑃̃ = , 𝑉𝑏 = ,𝑉 = , 𝑉𝑃 = and 𝐼𝑚 = ∫0 𝐼11 d𝑋 ℎ2 𝐼𝑚 𝐼𝑚 ℎ2 𝑐 𝐼𝑚 ℎ2 𝐼𝑚 ℎ2 (4-70) Substituting the postbuckling shape functions in Eqs. (4-45), (4-46) and (4-57) into Eqs. (4-66), (4-67) and (4-68), respectively, the normalized total potential energy of the FGM-MM is obtained as 2 1 𝐴𝑙2 2 1 𝜕2 𝑦𝑗 (𝑋) 6𝑘 𝚷FM = (2𝐼 + 4𝐼 ) [∑∞ 𝑗=1 𝐴𝑗 ∫0 𝐸FM (𝑋) ( ) d𝑋] + 𝑄2 [𝑑 − 𝑚 𝑚 𝜕𝑋 2 2 2 2 1 1 𝜕𝑦𝑗 (𝑋) 6𝑘 𝑑 1 𝜕𝑦𝑗 (𝑋) ∑∞ 𝐴 2 ∫0 ( ) d𝑋] − 𝑄2 [𝑑2 − 2 ∑∞ 2 𝑗=1 𝐴𝑗 ∫0 ( ) d𝑋] ; (4-71) 2 𝑗=1 𝑗 𝜕𝑋 𝜕𝑋 and the normalized total potential energy of the CNT-MM is 1 1 𝑄44 CM 2 2 𝚷CM = 2 ∫0 (1 + ) (∑∞ ∞ 𝑗=1,3,5… 𝐴𝑗 𝑁𝑗 cos 𝑁𝑗 𝑋 + ∑𝑗=2,4,6… 𝐴𝑗 (𝑁𝑗 cos 𝑁𝑗 𝑋 − 𝐼11 CM 2 𝑄̅11 CM ℎ2 𝑘 1 1 2𝑁𝑗 sin 𝑁𝑗 𝑋)) d𝑋 + [𝑑 − 2 ∫0 (∑∞ ∞ 𝑗=1,3,5… 𝐴𝑗 𝑁𝑗 sin 𝑁𝑗 𝑋 + ∑𝑗=2,4,6… 𝐴𝑗 (2 + 2𝐼11 CM 2 2 𝑄̅11 CM ℎ2 𝑘 𝑑 1 𝑁𝑗 sin 𝑁𝑗 𝑋 + 2 cos 𝑁𝑗 𝑋)) d𝑋] − [𝑑 2 − 2 ∫0 (∑∞ 𝑗=1,3,5… 𝐴𝑗 𝑁𝑗 sin 𝑁𝑗 𝑋 + 2𝐼11 CM 2 ∑∞ 𝑗=2,4,6… 𝐴𝑗 (2 + 𝑁𝑗 sin 𝑁𝑗 𝑋 + 2 cos 𝑁𝑗 𝑋)) d𝑋]. (4-72) 70 Minimizing total potential energies in Eqs. (4-71) and (4-72) with respect to the deflection 𝑊(𝑋) between bilateral constraints, the weight coefficients 𝐴𝑗 of the postbuckling shape functions can be determined for the FGM-MM and CNT-MM. Min[𝚷(𝐴𝑗 )] { . (4-73) 0 ≤ 𝑊(𝑋) ≤ 1 Given the complexity of the total potential energies, numerical approach is used to solve the energy minimization in Eq. (4-73). In particular, for sake of improve calculation efficiency while ensure accuracy, the first 20 modes were used to express the general buckling shape function, then the finite total potential energy expression were obtained. By using the built-in functions in Mathematica, the minimum potential energies with different axial displacement can be acquired. Figure 17 illustrates the flowchart of numerically solving the postbuckling response of FGM-MM and CNT-MM using the energy minimization. 71 Figure 17. Flowchart of numerically solving the postbuckling response of FGM-MM and CNT-MM using the energy minimization. 4.5 Numerical Modelling and Validation of the FGM-MM and CNT-MM In this section, numerical simulations are conducted to validate the presented theoretical model. The FGM-MM and CNT-MM are developed to obtain the postbuckling response in Abaqus R2016x. Linear perturbation buckle algorithm is used in the buckling analysis and dynamic implicit with Nlgeom is applied in the postbuckling analysis. Buckling imperfection is accounted the FE model by modifying the input files [113]. User-defined subroutine USDFLD is developed to characterize the material properties of the FGM-MM and CNT-MM. 4.5.1 Numerical Setup and Modelling Figure 18(a) displays the meshed FE models of the FGM-MM and CNT-MM subjected to the bilateral confinements. One end of the FGM-MM and CNT-MM (i.e., 𝑥 = 0) is fixed and the other 72 end can slide in the loading direction. To investigate the effect of the bilateral constraints on the postbuckling response of the FGM-MM and CNT-MM, contact interaction is considered in the FE models. In particular, the interaction properties are defined as “hard” contact that allows separation after frictionless contact for the normal and tangential behaviors are. In the FE models, linearly increasing loading is gradually applied to the edge of the beam ends (𝑥 = 𝐿) as axial displacement 𝑑, which is defined as 𝑑 = 𝐾𝑇; (4-74) where the amplitude and time period are 𝐾 = 60 μm and 𝑇 = 10 s, respectively. Figure 18(b) presents the postbuckling shape configurations between the CNT-MM, FGM1- MM and FGM2-MM. In particular, the postbuckled beams are demonstrated from the first buckling mode (Φ1), flatten first buckling mode (Φ1 flatten), third buckling mode (Φ3), flatten third buckling mode (Φ3 flatten), and fifth buckling mode. It can be seen that the psotbuckling shapes are significantly affected by the material functions of the FGMs and CNTs, especially for the beams in the higher buckling modes. Note that FGM1-MM denote the plate-like MM designed by the first type of FGMs and FGM2-MM denote the plates designed by the second FGMs. Following Eq. (4-6), the first and second types of FGMs are given as 73 Figure 18. (a) Meshed CNT-MM and FGM-MM with cylindrical corrugation under the axial force and bilateral constraints. (b) Postbuckling shape configurations of the CNT-MM, FGM1-MM and FGM2-MM from the first buckling mode (Φ1), flatten first buckling mode (Φ1 flatten), third buckling mode (Φ3), flatten third buckling mode (Φ3 flatten), and fifth buckling mode (In this study, FGM1-MM denote the MM designed with the first type of FGMs, and FGM2-MM are designed by the second FGMs). FGM1: 𝐸M (1 + 1.2𝑋)𝑏 { ; (4-75) FGM2: 𝐸𝑀 𝑒 𝑏𝑋 where 𝑏 is varied from 0.3 to 0.9. The geometric and material properties of the FE models are summarized in Table 7 [95]. 4.5.2 Comparison of the Theoretical and Numerical Models Figure 19 compares the postbuckling response of the CNT-MM and FGM-MM between the theoretical and numerical results. The geometric and materials in Table 7 are used to obtain the results. Figure 19(a) presents the comparison for the CNT-MM in terms of the force-displacement 74 relations, and the postbuckling shape configurations are provided in the buckling mode transitions (i.e., Φ1- Φ3 and Φ3- Φ5). It can be seen that when buckling mode transitions happen, the force- displacement curves behave sharply dropping in the axial force due to the displacement-control loading conditions. The theoretical model accurately captures the snap-through of the postbuckled CNT-MM in the force-displacement relations and deformed shape configurations. Figure 19(b) compares the postbuckling behaviors of FGM-MM between the theoretical and numerical results and good agreements are obtained for the force-displacement response and postbuckling shape configurations, which demonstrates the accuracy of the developed theoretical model in predicting the postbuckling performance of the architected beams. The theoretical and numerical results show clear differences of the postbuckling response results in the material properties, i.e., same geometries in MM corrugations but different materials in CNTs and FGMs, which indicate the possibility of using the material properties (e.g., material functions for CNTs and FGMs) to harness the posbuckling response. 75 Table 7. Geometric and material properties, element size and type, and loading condition of the FE models for the CNT-MM and FGM-MM. Length 𝐿 (mm) 2 Overall Width 𝑊𝑏 (mm) 0.4 Thickness 𝑡MM (μm) 10 Geometry Height 𝐻 (μm) 20 property Constraints net gap h (μm) 30 Corrugation Diameter 𝐷 (μm) 60 Rib width 𝐺 (μm) 20 Young’s modulus 𝐸11 (TPa) 5.6466 Shear modulus 𝐺13 (TPa) 1.9445 CNTs Volume fraction 𝑉CNT 0.028 Efficiency parameters 𝜂1 0.0058 Material Efficiency parameters 𝜂3 0.642 property Function 𝐸(𝑋),1 𝐸M (1 + 1.2𝑋)0.9 FGMs Function 𝐸(𝑋),2 𝐸M 𝑒 0.3𝑋 Matrix Young’s modulus 𝐸M (GPa) 2000 Length scale factor 𝑙 𝑡MM /10 FGM-MM and CNT-MM 0.01 Size 𝑙 (mm) Element size and type Bilateral constraints 0.05 Type S4R 76 Figure 19. Comparison of the force-displacement relations and postbuckling shape configurations at buckling mode transitions between the theoretical and numerical results for the (a) CNT-MM and (b) FGM1-MM and FGM2-MM. 4.5.3 Comparison of the Theoretical and Numerical Results with the Existing Study The CNT-MM and FGM-MM models are reduced to the uncorrugated, isotropic beams under bi-walls (i.e., without the corrugation in MM and composition in FGMs and CNTs). The reduced theoretical and numerical results are then compared with the postbuckling response reported in the existing study [63]. Table 8 shows the geometric and material factors used to compare the theoretical and numerical results with the theoretical results in the literature. 77 Table 8. Geometric and material properties used in the comparison between the presented and existing studies. Length L 250 Geometric Width 𝑊b 30 Properties Thickness 𝑡𝑀𝑀 2.3 (mm) Constraints gap H 4 Material properties Young’s modulus E (MPa) 2300 Figure 20. Comparison of the force-displacement relations between the reduced CNT-MM and FGM-MM models in this study, and the theoretical results in the existing study [63]. Figure 20 presents the comparison of the postbuckling response for the uncorrugated, isotropic beams under bilateral constraints between the theoretical and numerical models in this study and the theoretical results in the literature [63]. The developed CNT-MM and FGM-MM models are simplified to compare the force-displacement relations with the existing study, and satisfactory agreements are obtained. 78 4.6 Maneuverability of the Postbuckling Response for the FGM-MM and CNT-MM In this section, the presented theoretical models are used to investigate the influence of the material functions in FGMs and CNTs and the corrugation in MM aon the postbuckling response (i.e., force-displacement relations and postbuckling shape configurations) of the CNT-MM and FGM-MM. The material influence is studied in terms of the volume fraction 𝑉 for CNTs (see Eq. (8)) and index 𝑏 for FGMs (see Eq. (6)), and the geometric influence is particularly investigated with respect to the diameter 𝐷 and height 𝐻 of the corrugation in MM (see Figure 20). 4.6.1 Material Influence of FGMs and CNTs The effect of the material functions on the postbuckling response is presented on the CNT- MM and FGM-MM. The beam length 𝐿, width 𝑊b , thickness 𝑡, constraints net gap ℎ0 , and matrix Young’s modulus 𝐸M are fixed as 2 mm, 0.4 mm, 10 μm, 30 μm, and 130 GPa, respectively. Figure 21(a) shows the effect of the volume fraction 𝑉 on the force-displacement relations of the CNT- MM. Since different volumes of CNTs have different efficient parameters [95], the effective bending stiffness of the CNT-MM is critically changed, especially due to the fact that CNTs have extremely high stiffness-to-weight ratio. Figure 21(b) and Figure 21(c) display the effect of the index 𝑏 on the force-displacement curves for the FGM1-MM and FGM2-MM, respectively, using the FGMs defined in Eq. (4-75). Although the index has shown influence on the effective bending stiffness, less significance 79 Figure 21. Influence of the material functions on the postbuckling response of the CNT-MM, FGM1-MM and FGM2-MM. Force-displacement relations of the (a) CNT-MM affected by the volume fraction 𝑉, (b) FGM1-MM affected by the index 𝑏, and (c) FGM2-MM affected by the index 𝑏. Postbuckling shape configurations of the (d) CNT-MM in Φ1 for 𝑉 = 2.8%, Φ3 for 𝑉 = 12%, and Φ5 for 𝑉 = 28%, and the (e) FGM1-MM and (f) FGM2-MM in Φ1 for 𝑏 = 0.3, Φ3 for 𝑏 = 0.6, and Φ5 for 𝑏 = 0.9 (𝐿 = 2 𝑚𝑚, 𝑊𝑏 = 0.4 𝑚𝑚, 𝑡𝑀𝑀 = 10 𝜇𝑚, ℎ0 = 30 𝜇𝑚, and 𝐸𝑀 = 130 𝐺𝑃𝑎). is observed. Figure 21(d) shows the postbuckling shape configurations in the first buckling model (Φ1) for CNT-MM with 𝑉 = 2.8%, the third buckling model (Φ3) for 𝑉 = 12%, and the fifth buckling mode (Φ5) for 𝑉 = 28%. The presented buckling snap-through events are highlighted in Figure 21(a). Figure 21(e) and Figure 21(f) present the postbuckling shape configurations for the FGM1-MM and FGM2-MM, respectively, in the first buckling model (Φ1) for 𝑏 = 0.3, the third buckling model (Φ3) for 𝑏 = 0.6, and the fifth buckling model (Φ5) for 𝑏 = 0.9. 4.6.2 Geometric Influence of the Cylindrical Corrugation in MM Figure 22 depicts the effect of the cylindrical corrugation on the postbuckling response of the CNT-MM and FGM-MM. Figure 22(a) indicates the influence of the diameter-to-height ratio (i.e., 80 𝐷 = 0.2 − 5) of the corrugation on the force-displacement relations of the CNT-MM. Figure 22(b) 𝐻 𝐷 and Figure 22(c) show the effects of 𝐻 on the postbuckling response of the FGM1-MM and FGM2- 𝐷 MM, respectively. Note that 𝐻 are obtained from the cases of 20 μm-to-100 μm, 20 μm-to-20 μm and 100 μm-to-20 μm. Next, the material and geometric influences are investigated and compared on the postbuckling performance of the CNT-MM and FGM-MM. Figure 23(a) compares the force-displacement 𝐷 relations for the CNT-MM with 1) diameter-to-height ratio is 𝐻 = 0.2 and volume fraction is 𝑉 = 𝐷 𝐷 2.8% , 2) = 1 and 𝑉 = 12% , and 3) = 5 and 𝑉 = 28% . Figure 23(b) and Figure 23(c) 𝐻 𝐻 compares the force-displacement relations for the FGM1-MM and FGM2-MM, respectively, with 𝐷 𝐷 𝐷 1) = 0.2 and 𝑏 = 0.3, 2) = 1 and 𝑏 = 0.6, and 3) = 5 and 𝑏 = 0.9. According to the 𝐻 𝐻 𝐻 comparison, the volume fraction of CNTs proved more significant 81 Figure 22. Influence of the cylindrical corrugation on the postbuckling response of the CNT-MM, FGM1- MM and FGM2-MM. Force-displacement relations of the (a) CNT-MM, (b) FGM1-MM and(c) FGM2-MM 𝐷 affected by the diameter-to-height ratio . Postbuckling shape configurations of the (d) CNT-MM, (e) 𝐻 𝐷 𝐷 𝐷 𝐷 FGM1-MM and (f) FGM2-MM in Φ1 for 𝐻 = 0.2, Φ3 for 𝐻 = 1, and Φ5 for 𝐻 = 5. (𝐻 are obtained from the cases of 20 𝜇𝑚-to-100 𝜇𝑚, 20 𝜇𝑚-to-20 𝜇𝑚 and 100 𝜇𝑚-to-20 𝜇𝑚). influence on the CNT-MM, since the effective bending stiffness is not decreased critically when 𝐷 is reduced from 1 to 0.2 and the 𝑉 is increased from 2.8% to 12%, However, the FGM-MM are 𝐻 𝐷 found to be more sensitive to 𝐻 rather than the index b. Different geometric and material functions lead to different postbuckling response for the CNT-MM and FGM-MM. Figure 23(d) presents the distributions of the first four buckling mode coefficients (𝐴1 , 𝐴3 , 𝐴5 and 𝐴7 ) for the CNT-MM 𝐷 with 𝐻 = 0.2 and 𝑏 = 0.3. Figure 23(e) and Figure 23(f) show the distributions of 𝐴1 , 𝐴3 , 𝐴5 and 𝐷 𝐴7 for the FGM1-MM and FGM2-MM, respectively, with 𝐻 = 0.2 and 𝑏 = 0.3. It can be seen that certain mode coefficient dominants the postbuckling response in the stable buckling modes (i.e., 82 Φ1, Φ3 and Φ5), and significant coefficient fluctuations are obtained in buckling snap-throughs (i.e., Φ1-Φ3, Φ3-Φ5 and Φ5-Φ7). Figure 24 shows the variation of the larger axial force during the buckling mode transitions (i.e., 𝐹Φ3 for the Φ1-Φ3 snap-through) for the (a) CNT-MM, (b) FGM1-MM and (c) FGM2-MM. Different corrugation factors, 𝑉 and 𝑏 are used to investigate the influences of the material and 2 2(𝐷+𝐺)2 (𝑊+𝐻−2𝑡MM ) geometric functions. In particular, the corrugation factor is defined as [ ] based 𝑊((𝐷+𝐺)2 +𝜋𝐷𝐻) on Eq. (2). It can be seen that higher corrugation factors (i.e., more significant corrugation) and larger volume fraction and index b (i.e., more obvious functional materials) tend to provide higher controllability on the buckling mode transitions. In addition, CNTs are likely to be more maneuverable over the FGMs on the postbuckling response of the bilaterally constrained CNT- MM and FGM-MM. 83 Figure 23. Comparison of the material and geometric influences on the force-displacement relations for 𝐷 𝐷 𝐷 the (a) CNT-MM with = 0.2 and 𝑉 = 2.8%, = 1 and 𝑉 = 12%, and = 5 and 𝑉 = 28%, and (b) 𝐻 𝐻 𝐻 𝐷 𝐷 𝐷 FGM1-MM and (c) FGM2-MM with = 0.2 and 𝑏 = 0.3, = 1 and 𝑏 = 0.6, and = 5 and 𝑏 = 0.9. 𝐻 𝐻 𝐻 𝐷 Distributions of 𝐴1 , 𝐴3 , 𝐴5 and 𝐴7 for the (d) CNT-MM with = 0.2 and 𝑏 = 0.3, and the (e) FGM1-MM 𝐻 𝐷 and (f) FGM2-MM with 𝐻 = 0.2 and 𝑏 = 0.3. 84 Figure 24. Variation of the maximum axial force in buckling mode transitions for the (a) CNT-MM, (b) FGM1-MM and (c) FGM2-MM. 4.7 Additional Studies Figure 25(a) demonstrates the FGM-MM and CNT-MM designed with cylindrical corrugation, and Figure 25(b) shows the hollow beam that has the same length, width, height and thickness as the FGM-MM and CNT-MM. The volume of the MM in FGM-MM and CNT-MM, and the volume of the hollow beam can be obtained as [105]. 𝜋𝐷𝐻 𝑉MM = 𝐿𝑊𝑡MM (1 + (𝐷+𝐺)2) { . 𝑉hb = 2𝐿𝑡MM (𝑊 + 𝐻 − 2𝑡MM ) (A1) 85 Taking Eq. (A1) into Eq. (4-1) leads to Eq. (4-2). In addition, the weight factor 𝜁MM of the cylindrical corrugation for the FGM-MM and CNT- MM is numerically calibrated using Abaqus R2016x. The mean value is 0.43 using the parameters in Table 2. 4.7.1 User-Defined Subroutine USDFLD for the FGM-MM In this study, the material functions of the FGMs in the numerical modelling are defined using the field variables. The user-defined subroutine for the FGM-MM is given as, Figure 25. Demonstration of the (a) FGM-MM and CNT-MM designed with cylindrical corrugation, and the (b) hollow beam that has the same length, width, height and thickness as the FGM-MM and CNT-MM. 4.8 Conclusions In this study, we investigated the postbuckling response (i.e., force-displacement relations and postbuckling shape configurations) of the extensible mechanical metamaterials at the microscale designed by functionally graded materials (FGM-MM) and carbon nanotubes (CNT-MM). Theoretical models were developed to investigate the postbuckling response of the FGM-MM and CNT-MM subjected to bilateral constraints and solved using the energy method. Numerical 86 simulations were conducted to validate the theoretical results and good agreements were obtained. The theoretical models were used to investigate the maneuverability of the FGM-MM and CNT- MM with respect to the material properties (i.e., FGMs and CNTs) and geometric properties (i.e., corrugated microstructures). The findings showed that more significant corrugation in MM and more critical composition in FGMs and CNTs provided higher controllability on the buckling mode transitions, and CNTs were likely to be more maneuverable over the FGMs on the postbuckling response of the bilaterally constrained CNT-MM and FGM-MM. The reported FGM- MM and CNT-MM provide a novel direction of programming mechanical response of artificial materials. 87 Chapter 5: Structural Instability for the Design and Control of Computational Materials and Structures 5.1 Introduction Numerical computations, such as solving differential equations, are pervasive and well-known in scientific research and engineering areas [114, 115], as are many other tasks that involve simulation and optimization. The basic characteristics of any system can be represented by a set of equations and the prediction of the system’s behavior is their solution [116]. These differential equations arise in mathematical modeling of many real life application [117, 118], such as biosciences [119, 120]; human pupil reflex [121]; physics [122, 123]; weather prediction [124]; engineering [125] biology [126]; economy [127]. However, complicated mathematical approaches are required to achieve accurate and exact solutions of PDEs [128]. Moreover, understanding how to obtain solutions from differential equations, and their application to describe the physical processes in some cases is still challenging [129]. The latest developments in technology provides many different analytical and numerical methods solving tools [130]. However, these technologies are computer-based and still suffer severe drawbacks [131], it becomes a matter of interest to introduce other disruptive ways and approaches that will save cost, efforts for quicker results while improving the overall performance. The revival of analog computing might be the recourse for filling the gap. Different studies have been conducted to propose different analog computing system, for instance, an acoustic analog computing system was proposed, based on the concept on simulated transmitted pressure, as an ordinary differential equation solver [132]. On other hand, a silicon microring resonator that relies on the speed and wide band of the photonic integrated circuit, was used as an alternative solution to solve ordinary differential equation [133]. Moreover, another approach in which 88 electronic circuits were defined by periodic networks of repeated identical cells, was introduced and successfully obtained the finite difference approximation (i.e., locally) of the partial differential equation [134]. Also, recently a complete memristor-based hardware and software system with a high precise performance was proposed to solve static and time-evolving differential equation problems [135]. Structural Instability of slender elements (i.e., postbuckling) has been designed in advanced structures to meet multiple functional purposes in recent studies, such as autonomous sensing systems, medical devices and bioengineering [136-138]. Taking the advantage of the programmable, postbuckling-induced instability of the bi-walled beam, we propose a solution for analog computing. To overcome the issues associated with the performance of these smart applications, different studies that proposed different strategies (material strategies, geometric strategies, lateral confinement strategies) was conducted in order to achieve to better control over the postbuckling response of the bi-walled beams [139-141]. Our approach is based on bi-walled beam system capable of, in which deformed shapes are defined as solutions for the differential equation under study. This bi-walled beam must have the properties of being elastically deformed then return back to its original shape, when the force is released. Basically, for an arbitrary force to a bi-walled anisotropic non-uniform beam associated with a prescribed differential equation operator, the solution of such an equation is generated as a complex-valued output mode transition. By exploiting beam-walls interactions under quasi-static load our bi-walled beam-based analog computer may provide a route to achieve chip-scale, fast, and integrable computing elements. 5.2 Conceptual Design A conceptual representation of our idea for solving 4th order nonlinear differential equations using the postbuckling behavior of a non-uniform anisotropic bi-walled beam can be described by 89 characteristically considering the differential equation as an inverse problem, and the postbuckling mechanism can be used to generate the solution. We present a partial differential equation (PDE) solver for a class of higher-order nonlinear PDEs, as in Eq. (1). An appropriate analog approach was introduced to reach solutions that are consistent with the target equations in this research. Two applications (e.g., trigonometric functions, exponential and polynomial functions) and a parametric study were presented to verify the effectiveness and compatibility of the proposed technique. Figure 26 schematically illustrates the post-buckling process and the design principle of the PDE solver. 𝜕4 𝑓(𝑥) 𝜕2 𝑓(𝑥) 𝑔(𝑥, 𝑦, 𝑧) + 𝐶𝑜𝑛𝑠𝑡𝑎𝑛𝑡 = 0; (5-1) 𝜕𝑥 4 𝜕𝑥 2 𝜕2 𝑓(𝑥) where, represents the unknown solution of 4th order nonlinear differential equation. 𝜕𝑥 2 Figure 26. Postbuckling shape configurations in mode 1st, 3rd, 5th and 7th. 90 5.3 Theoretical formulation for PDE solver The prementioned fourth order partial differential equation, is represented by a bi-walled beam consists of the length 𝐿, width 𝑏, thickness 𝑡, cross section area 𝐴 = 𝑏𝑡, moment of inertia 𝐼 = 𝑏𝑡 3 , as shown in Figure 27. The beam has a fixed-fixed boundary condition and positioned between 12 two bilateral fixed frictionless walls along the length of the beam, with the net gap between the lateral constraints ℎ. The net gap can be simply defined as the difference between the total gap and the thickness of the beam ℎ = ℎ𝑜 − 𝑡. It worth to point out that the beam is modeled based on ℎ small deformation assumptions (since the ratio of the beam’s gap to length is relatively small 𝐿 ≪ 1). Figure 27. PDE solver arbitrary cross-section. 5.4 Trigonometric Functions Let’s assume that we have such 4th order partial differential equation: 𝑎03 2𝜋𝑥 2𝜋𝑥 𝜕4 𝑓(𝑥) 𝜕2 𝑓(𝑥) [(𝑎1 + |𝑎1 − 𝑎2 | 𝑐𝑜𝑠 ( 𝑎 )) (𝑎4 + 2𝑎5 𝑠𝑖𝑛 ( 𝑎 ))] + 𝑎6 = 0; (5-2) 12 3 3 𝜕𝑥 4 𝜕𝑥 2 91 To solve the above equation analogy, we are going to the build a beam with specific geometric and material properties in which its transverse displacement, while deflecting quasi-statically between frictionless bilateral constraints, will represent the solution of the partial differential equation. Since quasi-static assumptions are adopted, therefore the potential energy will represent the total energy and the net between the bending strain energy u_b, compressive strain energy and the energy of external work u_(p-c), is equal to the total potential energy. The total potential energy Ω can be written as 𝛺 = 𝑢𝑏 − 𝑢𝑝−𝑐 . (5-3) The energy components can be obtained as 2 1 𝑎 2𝜋𝑥 2𝜋𝑥 𝜕2 𝑓(𝑥) 𝑢𝑏 = 12 ∫0 3 [𝑎03 (𝑥) (𝑎4 (𝑥) + 2𝑎5 𝑠𝑖𝑛 ( 𝑎 )) (𝑎1 + |𝑎1 − 𝑎2 | 𝑐𝑜𝑠 ( 𝑎 )) ( ) ] d𝑥 ; 3 3 𝜕𝑥 2 (5-4) and 𝑎6 ∆𝑝−𝑐 𝑢𝑝−𝑐 = ; (5-5) 2 where 𝛥𝑝−𝑐 refers to the net deformation (net difference between the overall variation in the length of the beam and the axial compressive deformation) which can be obtained as 1 𝑎 𝜕𝑓(𝑥) 2 ∆𝑝−𝑐 = 2 ∫0 3 ( ) 𝑑𝑥. (5-6) 𝜕𝑥 Let’s assume that Ϛ represents the minimum energy path that PDEs solver will follow (rule of physics), Eq. (5-1) can be re-written as, 92 1 𝑎 2𝜋𝑥 Ϛ = min(𝛺) = 12 ∫0 3 𝑎03 (𝑥) (𝑎1 + |𝑎1 − 𝑎2 | 𝑐𝑜𝑠 ( 𝑎 )) (𝑎4 (𝑥) + 3 2 2 2𝜋𝑥 𝜕2 𝑓(𝑥) 𝑎6 𝑎3 𝜕𝑓(𝑥) 2𝑎5 𝑠𝑖𝑛 ( 𝑎 )) ( ) d𝑥 − ∫0 ( ) 𝑑𝑥 = 0. (5-7) 3 𝜕𝑥 2 2 𝜕𝑥 𝜕𝑓𝑏 (𝑥) First, let’s consider the bending strain energy; the solution 𝑓𝑏 (𝑥) and its first and second 𝜕𝑥 𝜕2 𝑓𝑏 (𝑥) derivatives can be expressed as, 𝜕𝑥 2 𝜕𝑢𝑏 (𝑥) 𝜕2 𝑓𝑏 (𝑥) 12 𝜕𝑥 =√ (5-8) 𝜕𝑥 2 2𝜋𝑥 2𝜋𝑥 𝑎03 (𝑎1 +|𝑎1 −𝑎2 | cos( ))(𝑎4 +2𝑎5 𝑠𝑖𝑛( )) 𝑎3 𝑎3 𝜕2 𝑓𝑏 (𝑥) and by double integrating the , 𝑓𝑏 (𝑥) resulted from bending can be obtained as, 𝜕𝑥 2 𝑎 𝑎 𝜕2 𝑓𝑏 (𝑥) 𝑓𝑏 (𝑥) = ∫0 3 ∫0 3 d𝑥 (5-9) 𝜕𝑥 2 Now, by considering both the compressive strain energy and the energy of external work, the 𝜕𝑓𝑐−𝑝 (𝑥) first derivative can be expressed as, 𝜕𝑥 𝜕𝑓𝑝−𝑐 (𝑥) 2 𝜕𝑢𝑝−𝑐 = √ (5-10) 𝜕𝑥 √𝑎6 𝜕𝑥 and the transverse displacements resulted from compressive and external work can be expressed as, 3 2 −1 𝑓𝑝−𝑐 (𝑥) = 3 4 √(𝜕𝑢𝑝−𝑐 ) (𝜕 𝑢𝑝−𝑐 ) + 𝑘3 (5-11) √𝑎6 𝜕𝑥 𝜕𝑥 2 where 𝑘𝑎 (𝑎 = 1,2,3) are the unknown integration constants that can be determine using the boundary conditions that describe the lateral confinement (clamped-clamped) can be described as, 93 𝑓𝑏 (0) = 𝑓𝑏 (1) = 𝑓𝑝−𝑐 (0) = 𝑓𝑝−𝑐 (1) = 0; {𝑑𝑓𝑏(𝑥) 𝑑𝑓𝑏 (𝑥) 𝑑𝑓𝑝−𝑐 (𝑥) 𝑑𝑓𝑝−𝑐 (𝑥) (5-12) |𝑥=0 = |𝑥=1 = |𝑥=0 = |𝑥=1 = 0. 𝑑𝑥 𝑑𝑥 𝑑𝑥 𝑑𝑥 The total 𝑓(𝑥) is the net of the transverse displacements resulted from bending, compressive and external work, 𝑓(𝑥) = ∑ 𝑓𝑏 (𝑥) − 𝑓𝑝−𝑐 (𝑥) ; (5-13) According to Salem et al. [139], the 4th order partial differential equation that govern the postbuckling behavior of beam, under quasi-static increase loading, constrained between frictionless bilateral walls can be expressed as, 𝜕4𝑤 ̂ (𝑥) 𝑝̂ 𝜕2 𝑤 ̂ (𝑥) + 𝐸𝐼 = 0; (5-14) 𝜕𝑥 4 𝜕𝑥 2 Thus, in order to solve the 4th order partial differential equation Eq. (5-2), the material distribution should be assumed to follow the cosine distribution as 𝐸(𝑥) = 𝑎1 + |𝑎1 − 2𝜋𝑥 𝑎2 | cos ( 𝑎 ) and the geometric variation should be assumed to follow the sinusoidal case as 3 𝑎03 (𝑥) 2𝜋𝑥 𝐼(𝑥) = (𝑎4 (𝑥) + 2𝑎5 𝑠𝑖𝑛 ( 𝑎 )) . Therefore, we need to build a walled beam with the 12 3 properties presented in Table 7 to obtain the solution of the equation under consideration. 94 Table 9. PDE-solver (trigonometric functions), parameter definition. 4th order PDE PDE-solver Definition (of the beam) 𝑎0 𝑡 Thickness of the beam 𝑎1 𝐸𝐴 Young’s modulus of stiffer material 𝑎2 𝐸𝐵 Young’s modulus of softer material 𝑎3 𝑙 Length of the beam 𝑎4 𝑏 Width of the beam 𝑎5 𝐴 Sine amplitude 𝑎6 𝑝̂ Axial compressive force 5.5 Exponential and Polynomial Functions Following the previous approach, here we are trying to solve another 4th order partial differential equation: 𝑎03 |𝑎1 −𝑎2 | 𝑥 𝜕4 𝑓(𝑥) 𝜕2 𝑓(𝑥) [(𝑎1 + 𝑥𝑎 ) (𝑎4 + (𝑎5 − 𝑎4 ) (𝑎 ))] + 𝑎6 = 0; (5-15) 12 ( 7) 3 𝜕𝑥 4 𝜕𝑥 2 𝑒 𝑎3 The minimum energy path Ϛ can be re-written as, 2 1 𝑎 |𝑎1 −𝑎2 | 𝑥 𝜕2 𝑓(𝑥) Ϛ = min(𝛺) = 12 ∫0 3 𝑎03 (𝑥) (𝑎1 + 𝑥𝑎 ) (𝑎4 + (𝑎5 − 𝑎4 ) (𝑎 )) ( ) d𝑥 − ( 7) 3 𝜕𝑥 2 𝑒 𝑎3 𝑎6 𝑎3 𝜕𝑓(𝑥) 2 2 ∫0 ( 𝜕𝑥 ) 𝑑𝑥 = 0. (5-16) Thus, in order to solve the 4th order partial differential equation Eq. (5-15), the material |𝑎1 −𝑎2 | distribution should be assumed to follow the distribution as 𝐸(𝑥) = 𝑎1 + 𝑥 and the geometric ( ) 𝑒 𝑎3 95 𝑎03 (𝑥) 𝑥 variation should be assumed as 𝐼(𝑥) = (𝑎4 + (𝑎5 − 𝑎4 ) (𝑎 )). Therefore, we need to build 12 3 a walled beam with the properties presented in Table 10 to obtain the solution of the differential equation under consideration. Table 10. PDE-solver (exponential and polynomial), parameter definition. 4th order PDE PDE-solver Definition (of the beam) 𝑎0 𝑡 Thickness of the beam 𝑎1 𝐸𝐴 Young’s modulus of stiffer material 𝑎2 𝐸𝐵 Young’s modulus of softer material 𝑎3 𝑙 Length of the beam 𝑎4 𝑏𝑡𝑜𝑝 Top width of the beam 𝑎5 𝑏𝑏𝑜𝑡𝑡𝑜𝑚 Bottom width of the beam 𝑎6 𝐿 Length of the beam 5.6 PDE solver Fabrication and Experimental Validation In this study, the PDE solver were fabricated using 3D printing technique to satisfy both low cost and well flexibility [142]. Figure 28 illustrates the principle and printing process of the Ultimaker printer. The polymer filaments were fed into the 3D printer through a hot extruder and then heated to glass transition temperature, to assure the filament be smoothly extruded from the brass nozzles [143, 144]. The printer deposits melted thermoplastic in a layer-by-layer format (with thickness tolerance of 0.06 mm) to create the PDE-solver designed in AutoCAD 2019 and Siemens NX. In 3D printing, the extrusion head rate was set to be 25 mm/s, and the printing bed temperature 96 was varied from 0-85 oC (based on the melting temperatures of the feeding filaments) and the temperatures of the nozzles were set at 200-245 oC. Furthermore, the standby temperatures for was fixed as 175 oC. The fan speeds used to print PLA was set as 40-100%. Figure 28. 3D printer for the fabrication of the PDE solver. 5.7 Experimental Validation To ensure the high level of tunability and controllability, the 3D printed beam was placed between two flat and rigid bilateral constraints with clamped-clamped boundary conditions. The constraints were made from aluminum with frictionless surface to omit the influence of friction on the behavior of the PDE-solver. Gradually increasing axial force was applied to the top of the PDE- solver by MTS Flextest 40 and a loading frame unit model370. The loading stage was displacement-control and was limited to a maximum total shorting of 6 mm with a loading period of 40 seconds. To validate the proposed model, the trigonometric function case was first considered. A bi-walled beam with a thickness of 𝑎0 = 1.2 𝑚𝑚, young’s modulus (for both the 97 soft and stiff materials) of 𝑎1 = 𝑎2 = 3.5 𝐺𝑃𝑎, length of 𝑎3 = 185 𝑚𝑚, width of 𝑎4 = 25 𝑚𝑚 and has a sine amplitude of 𝑎5 = 0.5 𝑚𝑚 . In the second set of experiment, we tested the exponential and polynomial functions case by considering a beam with thickness of 𝑎0 = 2.76 𝑚𝑚, young’s modulus of stiff material of 𝑎1 = 3.5 𝐺𝑃𝑎, young’s modulus of soft material of 𝑎2 = 0.9 𝐺𝑃𝑎, length of 𝑎3 = 170 𝑚𝑚, width of 𝑎4 = 25 𝑚𝑚 and volume fraction of 𝑎7 = 0.76 𝑚𝑚. The solution for the first set (Case A) can be represented as, 98 0.181539(1 − Cos[2𝜋𝑋]) + 0.10201190066875167(1 − Cos[4𝜋𝑋]) + 0.11127544904624068(1 − Cos[6𝜋𝑋]) + 0.024966657231540405(1 − Cos[8𝜋𝑋]) + 0.010977076443737004(1 − Cos[10𝜋𝑋]) − 0.0036583452682027937(1 − Cos[12𝜋𝑋]) − 0.0019335898911041057(1 − Cos[14𝜋𝑋]) − 0.007487237382895869(1 − Cos[16𝜋𝑋]) − 0.0013705099605957258(1 − Cos[18𝜋𝑋]) + 0.00034608282770272846(1 − Cos[20𝜋𝑋]) − 0.00035149035271551254(1 − Cos[22𝜋𝑋]) + 0.00019157727488008534(1 − Cos[24𝜋𝑋]) − 0.0008648467277794851(1 − Cos[26𝜋𝑋]) − 0.00037732327103937795(1 − Cos[28𝜋𝑋]) + 0.00004081026325742306(1 − Cos[30𝜋𝑋]) + 0.17029917874405862(1 − 2𝑋 − Cos[8.986839944858962𝑋] + 0.22254763768705213Sin[8.986839944858962𝑋]) − 0.4212368072332286(1 − 2𝑋 − Cos[15.450352670354603𝑋] + 0.12944688336063062Sin[15.450352670354603𝑋]) + 0.00609891052812938(1 − 2𝑋 − Cos[21.808307882689626𝑋] + 0.09170816969195042Sin[21.808307882689626𝑋]) + 0.009546616684841008(1 − 2𝑋 − Cos[28.132333894365882𝑋] + 0.07109257296283349Sin[28.132333894365882𝑋]) + 0.0009529577231114289(1 − 2𝑋 − Cos[34.44159442057026𝑋] + 0.058069320937287935Sin[34.44159442057026𝑋]) − 0.011899721568971356(1 − 2𝑋 − Cos[40.74268680587531𝑋] + 0.0490885642748428Sin[40.74268680587531𝑋]) + 0.01159708156921189(1 − 2𝑋 − Cos[47.03875264293461𝑋] 99 + 0.04251813425372382Sin[47.03875264293461𝑋]) + 0.0029629018855069737(1 − 2𝑋 − Cos[53.33199104660569𝑋] + 0.037500943819109296Sin[53.33199104660569𝑋]) + 0.001096948627932418(1 − 2𝑋 − Cos[59.62334449468461𝑋] + 0.03354390829548817Sin[59.62334449468461𝑋]) + 0.0006509484167195499(1 − 2𝑋 − Cos[65.91281298717136𝑋] + 0.030343114022295798Sin[65.91281298717136𝑋]) − 0.0008729809646230986(1 − 2𝑋 − Cos[72.20228147965814𝑋] + 0.02769995572180733Sin[72.20228147965814𝑋]) + 0.0006570504807220779(1 − 2𝑋 − Cos[78.48860837949131𝑋] + 0.025481404770613697Sin[78.48860837949131𝑋]) + 0.00039430615912742367(1 − 2𝑋 − Cos[84.77587775712057𝑋] + 0.023591616541322268Sin[84.77587775712057𝑋]) + 0.00019831880723003418(1 − 2𝑋 − Cos[91.06220465695374𝑋] + 0.021963008775532374Sin[91.06220465695374𝑋]) + 0.00017450091658007073(1 − 2𝑋 − Cos[97.34821739752155𝑋] + 0.020544803525605383Sin[97.34821739752155𝑋]) While the solution for the 1st mode transition in the second set (Case B) can be represented as, 100 0.5136794559191887(1 − Cos[2𝜋𝑋]) − 0.0000506642263922924(1 − Cos[4𝜋𝑋]) − 0.011919205492912258(1 − Cos[6𝜋𝑋]) − 0.000005986220498192588(1 − Cos[8𝜋𝑋]) − 0.0012529910689168886(1 − Cos[10𝜋𝑋]) − 0.000002407715704313205(1 − Cos[12𝜋𝑋]) − 0.0003092981214378512(1 − Cos[14𝜋𝑋]) − 0.000001302192335233192(1 − Cos[16𝜋𝑋]) − 0.00011054681524310559(1 − Cos[18𝜋𝑋]) − 8.116778710742132 × 10−7 (1 − Cos[20𝜋𝑋]) − 0.000048831852312373294(1 − Cos[22𝜋𝑋]) − 5.507803749402448 × 10−7 (1 − Cos[24𝜋𝑋]) − 0.00002476190266275978(1 − Cos[26𝜋𝑋]) − 3.959943009497743 × 10−7 (1 − Cos[28𝜋𝑋]) − 0.000013839623206363896(1 − Cos[30𝜋𝑋]) + 0.0067440092045444(1 − 2𝑋 − Cos[8.986839944858962𝑋] + 0.22254763768705213Sin[8.986839944858962𝑋]) + 0.0002521984227809577(1 − 2𝑋 − Cos[15.450352670354603𝑋] + 0.12944688336063062Sin[15.450352670354603𝑋]) − 0.00008459169489831279(1 − 2𝑋 − Cos[21.808307882689626𝑋] + 0.09170816969195042Sin[21.808307882689626𝑋]) + 0.00003025128126188748(1 − 2𝑋 − Cos[28.132333894365882𝑋] + 0.07109257296283349Sin[28.132333894365882𝑋]) − 0.000017866910615388942(1 − 2𝑋 − Cos[34.44159442057026𝑋] 101 + 0.058069320937287935Sin[34.44159442057026𝑋]) + 0.000009411843208285503(1 − 2𝑋 − Cos[40.74268680587531𝑋] + 0.0490885642748428Sin[40.74268680587531𝑋]) − 0.000006580834272156238(1 − 2𝑋 − Cos[47.03875264293461𝑋] + 0.04251813425372382Sin[47.03875264293461𝑋]) + 0.000004024549436129806(1 − 2𝑋 − Cos[53.33199104660569𝑋] + 0.037500943819109296Sin[53.33199104660569𝑋]) − 0.000003045591031604344(1 − 2𝑋 − Cos[59.62334449468461𝑋] + 0.03354390829548817Sin[59.62334449468461𝑋]) + 0.000002079244403972599(1 − 2𝑋 − Cos[65.91281298717136𝑋] + 0.030343114022295798Sin[65.91281298717136𝑋]) − 0.000001587101564243122(1 − 2𝑋 − Cos[72.20228147965814𝑋] + 0.02769995572180733Sin[72.20228147965814𝑋]) + 0.000001175843876178375(1 − 2𝑋 − Cos[78.48860837949131𝑋] + 0.025481404770613697Sin[78.48860837949131𝑋]) − 9.756578315261839 × 10−7 (1 − 2𝑋 − Cos[84.77587775712057𝑋] + 0.023591616541322268Sin[84.77587775712057𝑋]) + 7.300999721994216 × 10−7 (1 − 2𝑋 − Cos[91.06220465695374𝑋] + 0.021963008775532374Sin[91.06220465695374𝑋]) − 6.179859254044011 × 10−7 (1 − 2𝑋 − Cos[97.34821739752155𝑋] + 0.020544803525605383Sin[97.34821739752155𝑋]) and the solution for the 3rd mode transition in the second set (Case B) can be represented as, 102 5.459239694979734 × 10−8 (1 − Cos[2𝜋𝑋]) + 0.5022508118783651(1 − Cos[4𝜋𝑋]) − 6.30433557144345 × 10−8 (1 − Cos[6𝜋𝑋]) − 4.410371209712353 × 10−10 (1 − Cos[8𝜋𝑋]) − 4.054170783345386 × 10−9 (1 − Cos[10𝜋𝑋]) − 0.001962345115767776(1 − Cos[12𝜋𝑋]) − 1.031181542896072 × 10−9 (1 − Cos[14𝜋𝑋]) + 4.603705323998832 × 10−11 (1 − Cos[16𝜋𝑋]) − 2.814254007003814 × 10−10 (1 − Cos[18𝜋𝑋]) − 0.0002300873594734294(1 − Cos[20𝜋𝑋]) − 1.25989059538196 × 10−11 (1 − Cos[22𝜋𝑋]) − 1.264615821391044 × 10−10 (1 − Cos[24𝜋𝑋]) − 1.171493730657464 × 10−10 (1 − Cos[26𝜋𝑋]) − 0.000058361557561085415(1 − Cos[28𝜋𝑋]) − 7.095178792867994 × 10−11 (1 − Cos[30𝜋𝑋]) + 0.00041978114509530186(1 − 2𝑋 − Cos[8.986839944858962𝑋] + 0.22254763768705213Sin[8.986839944858962𝑋]) − 0.00046107460179950255(1 − 2𝑋 − Cos[15.450352670354603𝑋] + 0.12944688336063062Sin[15.450352670354603𝑋]) + 0.000007752188889690367(1 − 2𝑋 − Cos[21.808307882689626𝑋] + 0.09170816969195042Sin[21.808307882689626𝑋]) + 7.150297094650262 × 10−7 (1 − 2𝑋 − Cos[28.132333894365882𝑋] + 0.07109257296283349Sin[28.132333894365882𝑋]) − 0.000002945815934153597(1 − 2𝑋 − Cos[34.44159442057026𝑋] 103 + 0.058069320937287935Sin[34.44159442057026𝑋]) − 0.00000135116569878663(1 − 2𝑋 − Cos[40.74268680587531𝑋] + 0.0490885642748428Sin[40.74268680587531𝑋]) − 8.79691913311177 × 10−8 (1 − 2𝑋 − Cos[47.03875264293461𝑋] + 0.04251813425372382Sin[47.03875264293461𝑋]) − 4.46109071003682 × 10−8 (1 − 2𝑋 − Cos[53.33199104660569𝑋] + 0.037500943819109296Sin[53.33199104660569𝑋]) − 2.137233094143693 × 10−7 (1 − 2𝑋 − Cos[59.62334449468461𝑋] + 0.03354390829548817Sin[59.62334449468461𝑋]) − 1.922464600075722 × 10−7 (1 − 2𝑋 − Cos[65.91281298717136𝑋] + 0.030343114022295798Sin[65.91281298717136𝑋]) + 4.547821896314855 × 10−7 (1 − 2𝑋 − Cos[72.20228147965814𝑋] + 0.02769995572180733Sin[72.20228147965814𝑋]) − 6.443382503319151 × 10−8 (1 − 2𝑋 − Cos[78.48860837949131𝑋] + 0.025481404770613697Sin[78.48860837949131𝑋]) − 6.411825893513957 × 10−8 (1 − 2𝑋 − Cos[84.77587775712057𝑋] + 0.023591616541322268Sin[84.77587775712057𝑋]) − 6.954231304079054 × 10−8 (1 − 2𝑋 − Cos[91.06220465695374𝑋] + 0.021963008775532374Sin[91.06220465695374𝑋]) − 3.808032485636905 × 10−9 (1 − 2𝑋 − Cos[97.34821739752155𝑋] + 0.020544803525605383Sin[97.34821739752155𝑋]) and the solution for the 5th mode transition in the second set (Case B) can be represented as, 104 0.0025403607736504486(1 − Cos[2𝜋𝑋]) − 0.003769348141774792(1 − Cos[4𝜋𝑋]) + 0.5044865571485785(1 − Cos[6𝜋𝑋]) + 0.004006866875739551(1 − Cos[8𝜋𝑋]) + 0.002102686319344609(1 − Cos[10𝜋𝑋]) + 0.000006571735023804506(1 − Cos[12𝜋𝑋]) − 0.00040710038866607293(1 − Cos[14𝜋𝑋]) + 0.00002461810044278898(1 − Cos[16𝜋𝑋]) − 0.003295597954262204(1 − Cos[18𝜋𝑋]) + 0.000026822546098006285(1 − Cos[20𝜋𝑋]) + 0.0000832969476807227(1 − Cos[22𝜋𝑋]) + 0.000001132782605060111(1 − Cos[24𝜋𝑋]) − 0.00006158028826421117(1 − Cos[26𝜋𝑋]) + 6.276303333789091 × 10−7 (1 − Cos[28𝜋𝑋]) − 0.00036814342403603575(1 − Cos[30𝜋𝑋]) + 0.0034113173187788(1 − 2𝑋 − Cos[8.986839944858962𝑋] + 0.22254763768705213Sin[8.986839944858962𝑋]) − 0.01044503115586656(1 − 2𝑋 − Cos[15.450352670354603𝑋] + 0.12944688336063062Sin[15.450352670354603𝑋]) + 0.010428692570783752(1 − 2𝑋 − Cos[21.808307882689626𝑋] + 0.09170816969195042Sin[21.808307882689626𝑋]) − 0.003294260495508135(1 − 2𝑋 − Cos[28.132333894365882𝑋] + 0.07109257296283349Sin[28.132333894365882𝑋]) + 0.0002891712304945036(1 − 2𝑋 − Cos[34.44159442057026𝑋] + 0.058069320937287935Sin[34.44159442057026𝑋]) − 0.00016449470978471652(1 − 2𝑋 − Cos[40.74268680587531𝑋] 105 + 0.0490885642748428Sin[40.74268680587531𝑋]) + 0.00014014267307558932(1 − 2𝑋 − Cos[47.03875264293461𝑋] + 0.04251813425372382Sin[47.03875264293461𝑋]) + 0.00023158546366680712(1 − 2𝑋 − Cos[53.33199104660569𝑋] + 0.037500943819109296Sin[53.33199104660569𝑋]) − 0.0000936789832390618(1 − 2𝑋 − Cos[59.62334449468461𝑋] + 0.03354390829548817Sin[59.62334449468461𝑋]) − 0.00007089675115829527(1 − 2𝑋 − Cos[65.91281298717136𝑋] + 0.030343114022295798Sin[65.91281298717136𝑋]) + 0.000014497342787543317(1 − 2𝑋 − Cos[72.20228147965814𝑋] + 0.02769995572180733Sin[72.20228147965814𝑋]) − 0.000020922135506748997(1 − 2𝑋 − Cos[78.48860837949131𝑋] + 0.025481404770613697Sin[78.48860837949131𝑋]) + 0.000025860859351685054(1 − 2𝑋 − Cos[84.77587775712057𝑋] + 0.023591616541322268Sin[84.77587775712057𝑋]) + 0.0000363747703382087(1 − 2𝑋 − Cos[91.06220465695374𝑋] + 0.021963008775532374Sin[91.06220465695374𝑋]) − 0.000020743049617062955(1 − 2𝑋 − Cos[97.34821739752155𝑋] + 0.020544803525605383Sin[97.34821739752155𝑋]) Figure 29 displays the experimental and theoretical F-D curves and the normalized deflected shape for the trigonometric and exponential and polynomial functions under consideration up to the fifth mode transition. As shown, changing the material (Young’s modulus) and geometric 106 (shape configuration) properties led to well-defined tenability for the postbuckling behavior of the proposed bilateral constrained beam, which resulted in better controllability over the overall mechanism of the PDE-solver. Case B shows the approximate solutions for the fourth order PDE at different loading stages. Each equilibrium configuration (snap-through event) represents the minimum energy path (i.e., minimum potential energy) the PDE solver will follow to reach to the stable configuration. Figure 29. Comparison of the force-displacement relationship between the theoretical and experimental results for the material and geometric variation cases. 5.8 Conclusions Due to the possible advantages of analog computing, we have tried to add other features to this type of computational methods. We were able to develop a novel approach to solve a class of higher-order partial differential equations. The comparisons between theoretical model and the solution obtained from the experimental testing, indicates that the behavior of the bilateral constrained beam (solution of the 4th order PDE) can be predicted. Thus, the use of structural instability in modeling PDEs can be suitable alternative for the FEM approaches and other types of solving PDEs. It’s worth mentioning that the mechanism followed in our proposed analytical 107 tool to solve PDE didn’t exceed few seconds compared to the time required (hours) to solve similar equation using numerical operations (Mathematica software). We observed that our postbuckling- based solver is fast, simple and capable to solve such form of PDEs, which could broaden the horizon for computer-free computing approaches. 108 Chapter 6: Conclusions and Future Work 6.1 Conclusions This study examined the potential deployment of structural instability (postbuckling) in computation (i.e., solving high order partial differential equations). To achieve this objective, the main research contributions of the dissertation can be summarized as follows, 6.1.1 Post-Buckling Analysis of Bi-walled Non-Uniform Beams Using Small Deformation Assumptions Theoretical models are developed to investigate the mechanical responses of the laterally confined systems subjected to a quasi-static axial force. Nonlinear basic equations as of the microbeam are derived using the basis of modified couple stress theory and Euler-Bernoulli beam theory. To the achieve the desired level of controllability over the beam behavior, it is of necessity to accurately control the mode transitions and the travelling distances between them. Since, mode transitions’ locations can highly affect the development of the higher order modes (5th and 7th), this study proposed theoretical models that account for different shape configurations (e.g., sinusoidal, linear) and examined the effect of geometric parameters (e.g., width, length, thickness) on the postbuckling response. 6.1.2 Post-Buckling Analysis of Bi-walled Functionally Graded Material Beams Using Small Deformation Assumptions This study developed the mechanical response (i.e., postbuckling behavior under the bilateral constraints) of deformable functionally graded materials beams. The mathematical model was first presented to theoretically characterize the effect of material properties (i.e., effective young’s modulus) of the deformable FGMs bi-walled beam and then the energy minimization approach was employed to determine the minimum energy paths of the anisotropic constrained beam, along 109 with the corresponding mode transitions. In addition, topology algorithm was developed to optimize the stored strain energy by optimizing the young’s modulus and volume fraction between different materials. 6.1.3 Design and Control of Partial Differential Equations-Solver This study investigated the potential application of structural instabilities in single microbeams for computational capabilities into structural systems, in order to mimic the performance of analog partial differential equations solver using the mechanical deformation in structural elements. It indicated that the response of bi-walled beam can effectively solve the 4th order nonlinear partial differential equations (that represent different systems) analytically but was insufficient in solving the full range of loading (e.g., dynamic loading) since kinetic energy is negligible. To validate the application of the proposed analog computing approach, experiments were conducted covering both material and geometric variation cases. 6.2 Future Work 6.2.1 Optimization of the PDE-solver Parameters to Account for Different PDEs The conduct studies in this work provide the essential tools for the useable implementation of the proposed mechanism in solving 4th order nonlinear partial differential equations. The boundary conditions and the postbuckling governing equation restrict the problem under consideration to the 4th order presented form. Therefore, other boundary conditions, structural elements (cylinders, plates) and confinement environment (flexible, irregular walls) should be studied such that the proposed analog computing mechanism covers wider applications. 6.2.2 Optimization of the Algorithm to Include Different Loading Conditions In this work, the effect of the kinetic energy of the bi-walled beam was not included. It would be of interest to expand the theoretical model and account for the dynamic response of the 110 microbeam which could open different avenues for new applications that can be formed in PDE forms. Also, changing the loading type (i.e., applying electrical or magnetic field) would be a good path to study. 111 REFERENCES 112 REFERENCES 1. Timoshenko, S.P. and J.M. Gere, Theory of elastic stability. 2009: Courier Corporation. 2. Gao, X. and H. Ma, Topology optimization of continuum structures under buckling constraints. Computers & Structures, 2015. 157: p. 142-152. 3. Dunning, P.D., et al., Level‐set topology optimization with many linear buckling constraints using an efficient and robust eigensolver. International Journal for Numerical Methods in Engineering, 2016. 107(12): p. 1029-1053. 4. Wang, Y., et al., Buckling and postbuckling of dielectric composite beam reinforced with Graphene Platelets (GPLs). Aerospace Science and Technology, 2019. 91: p. 208-218. 5. Bertoldi, K., et al., Negative Poisson's ratio behavior induced by an elastic instability. Advanced materials, 2010. 22(3): p. 361-366. 6. Lin, C., et al., Nanocardboard as a nanoscale analog of hollow sandwich plates. Nature communications, 2018. 9(1): p. 1-8. 7. Dou, C., Y.-L. Pi, and W. Gao, Shear resistance and post-buckling behavior of corrugated panels in steel plate shear walls. Thin-Walled Structures, 2018. 131: p. 816-826. 8. Hu, N. and R. Burgueño, Buckling-induced smart applications: recent advances and trends. Smart Materials and Structures, 2015. 24(6): p. 063001. 9. Zhou, H., et al., Circumferential buckling and postbuckling analysis of thin films integrated on a soft cylindrical substrate with surface relief structures. Extreme Mechanics Letters, 2020: p. 100624. 10. Chen, Q., X. Zhang, and B. Zhu, Design of buckling-induced mechanical metamaterials for energy absorption using topology optimization. Structural and Multidisciplinary Optimization, 2018. 58(4): p. 1395-1410. 11. Stamatelos, D., G. Labeas, and K. Tserpes, Analytical calculation of local buckling and post- buckling behavior of isotropic and orthotropic stiffened panels. Thin-Walled Structures, 2011. 49(3): p. 422-430. 12. Ferreira, P.S. and F.B. Virtuoso, Semi-analytical models for the post-buckling analysis and ultimate strength prediction of isotropic and orthotropic plates under uniaxial compression with the unloaded edges free from stresses. Thin-Walled Structures, 2014. 82: p. 82-94. 13. Hu, N. and R. Burgueño, Harnessing seeded geometric imperfection to design cylindrical shells with tunable elastic postbuckling behavior. Journal of Applied Mechanics, 2017. 84(1). 14. Jiao, P., et al., Enhancement of quasi-static strain energy harvesters using non-uniform cross- section post-buckled beams. Smart Materials and Structures, 2017. 26(8): p. 085045. 113 15. Asadi, H. and M. Aghdam, Large amplitude vibration and post-buckling analysis of variable cross- section composite beams on nonlinear elastic foundation. International Journal of Mechanical Sciences, 2014. 79: p. 47-55. 16. Jiao, P., W. Borchani, and N. Lajnef, Large deformation solutions to static and dynamic instabilities of post-buckled beam systems. International Journal of Solids and Structures, 2017. 128(1): p. 85-98. 17. Genna, F. and G. Bregoli, Small amplitude elastic buckling of a beam under monotonic axial loading, with frictionless contact against movable rigid surfaces. Journal of Mechanics of Materials and Structures, 2014. 9(4): p. 441-463. 18. Katz, S. and S. Givli, The post-buckling behavior of a beam constrained by springy walls. Journal of the Mechanics and Physics of Solids, 2015. 78: p. 443-466. 19. Chen, J.-S. and Z.-S. Wen, Deformation and vibration of a buckled beam constrained by springy walls. European Journal of Mechanics-A/Solids, 2019. 77: p. 103791. 20. Borchani, W., et al., Control of postbuckling mode transitions using assemblies of axially loaded bilaterally constrained beams. Journal of Engineering Mechanics, 2017. 143(10): p. 04017116. 21. Li, Q., Buckling of multi-step non-uniform beams with elastically restrained boundary conditions. Journal of Constructional Steel Research, 2001. 57(7): p. 753-777. 22. Li, Q., H. Cao, and G. Li, Stability analysis of a bar with multi-segments of varying cross-section. Computers & structures, 1994. 53(5): p. 1085-1089. 23. Yu, X., et al., Mechanical metamaterials associated with stiffness, rigidity and compressibility: A brief review. Progress in Materials Science, 2018. 94: p. 114-173. 24. Schaedler, T.A., et al., Ultralight metallic microlattices. Science, 2011. 334(6058): p. 962-965. 25. Zheng, X., et al., Ultralight, ultrastiff mechanical metamaterials. Science, 2014. 344(6190): p. 1373-1377. 26. Meza, L.R., S. Das, and J.R. Greer, Strong, lightweight, and recoverable three-dimensional ceramic nanolattices. Science, 2014. 345(6202): p. 1322-1326. 27. Bauer, J., et al., Approaching theoretical strength in glassy carbon nanolattices. Nature materials, 2016. 15(4): p. 438-443. 28. Grima, J.N., R. Gatt, and P.S. Farrugia, On the properties of auxetic meta‐tetrachiral structures. physica status solidi (b), 2008. 245(3): p. 511-520. 29. Silverberg, J.L., et al., Using origami design principles to fold reprogrammable mechanical metamaterials. science, 2014. 345(6197): p. 647-650. 30. Xu, S., et al., Design of lattice structures with controlled anisotropy. Materials & Design, 2016. 93: p. 443-447. 31. Frenzel, T., M. Kadic, and M. Wegener, Three-dimensional mechanical metamaterials with a twist. Science, 2017. 358(6366): p. 1072-1074. 114 32. Chakraborty, A., S. Gopalakrishnan, and J. Reddy, A new beam finite element for the analysis of functionally graded materials. International journal of mechanical sciences, 2003. 45(3): p. 519- 539. 33. Shahba, A. and S. Rajasekaran, Free vibration and stability of tapered Euler–Bernoulli beams made of axially functionally graded materials. Applied Mathematical Modelling, 2012. 36(7): p. 3094-3111. 34. Bıˆrsan, M., et al., Mechanical behavior of sandwich composite beams made of foams and functionally graded materials. International Journal of Solids and Structures, 2013. 50(3-4): p. 519- 530. 35. Ke, L.-L., et al., Nonlinear free vibration of size-dependent functionally graded microbeams. International Journal of Engineering Science, 2012. 50(1): p. 256-267. 36. Li, L. and Y. Hu, Nonlinear bending and free vibration analyses of nonlocal strain gradient beams made of functionally graded material. International Journal of Engineering Science, 2016. 107: p. 77-97. 37. Thai, H.-T. and T.P. Vo, Bending and free vibration of functionally graded beams using various higher-order shear deformation beam theories. International Journal of Mechanical Sciences, 2012. 62(1): p. 57-66. 38. Wattanasakulpong, N., B.G. Prusty, and D.W. Kelly, Thermal buckling and elastic vibration of third-order shear deformable functionally graded beams. International Journal of Mechanical Sciences, 2011. 53(9): p. 734-743. 39. Trabelsi, S., et al., Thermal post-buckling analysis of functionally graded material structures using a modified FSDT. International Journal of Mechanical Sciences, 2018. 144: p. 74-89. 40. Evci, C. and M. Gülgeç, Functionally graded hollow cylinder under pressure and thermal loading: Effect of material parameters on stress and temperature distributions. International Journal of Engineering Science, 2018. 123: p. 92-108. 41. Zhang, Y., et al., Thermal shock resistance of functionally graded materials with mixed-mode cracks. International Journal of Solids and Structures, 2019. 164: p. 202-211. 42. Frikha, A., et al., A new higher order C0 mixed beam element for FGM beams analysis. Composites Part B: Engineering, 2016. 106: p. 181-189. 43. Zghal, S. and F. Dammak, Vibrational behavior of beams made of functionally graded materials by using a mixed formulation. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 2020: p. 0954406220916533. 44. Zghal, S., D. Ataoui, and F. Dammak, Static bending analysis of beams made of functionally graded porous materials. Mechanics Based Design of Structures and Machines, 2020: p. 1-18. 45. Chen, D., J. Yang, and S. Kitipornchai, Free and forced vibrations of shear deformable functionally graded porous beams. International journal of mechanical sciences, 2016. 108: p. 14-22. 115 46. Chen, M., et al., Isogeometric three-dimensional vibration of variable thickness parallelogram plates with in-plane functionally graded porous materials. International Journal of Mechanical Sciences, 2020. 169: p. 105304. 47. Nguyen, L.B., et al., An isogeometric Bézier finite element method for vibration analysis of functionally graded piezoelectric material porous plates. International Journal of Mechanical Sciences, 2019. 157: p. 165-183. 48. Nikbakht, S., S. Kamarian, and M. Shakeri, A review on optimization of composite structures Part II: Functionally graded materials. Composite Structures, 2019. 214: p. 83-102. 49. Ghatage, P.S., V.R. Kar, and P.E. Sudhagar, On the numerical modelling and analysis of multi- directional functionally graded composite structures: A review. Composite Structures, 2020. 236: p. 111837. 50. Xu, X.-J. and J.-M. Meng, A model for functionally graded materials. Composites Part B: Engineering, 2018. 145: p. 70-80. 51. Zhang, Y.-W., et al., Supersonic aerodynamic piezoelectric energy harvesting performance of functionally graded beams. Composite Structures, 2020. 233: p. 111537. 52. Bouhamed, A., et al., Homogenization of elasto-plastic functionally graded material based on representative volume element: application to incremental forming process. International Journal of Mechanical Sciences, 2019. 160: p. 412-420. 53. Nejad, M.Z. and P. Fatehi, Exact elasto-plastic analysis of rotating thick-walled cylindrical pressure vessels made of functionally graded materials. International Journal of Engineering Science, 2015. 86: p. 26-43. 54. Şimşek, M., Bi-directional functionally graded materials (BDFGMs) for free and forced vibration of Timoshenko beams with various boundary conditions. Composite Structures, 2015. 133: p. 968- 978. 55. Wang, Z.-h., et al., Free vibration of two-directional functionally graded beams. Composite structures, 2016. 135: p. 191-198. 56. Liu, B., H. Chen, and W. Cao, A novel method for tailoring elasticity distributions of functionally graded porous materials. International Journal of Mechanical Sciences, 2019. 157: p. 457-470. 57. Leseduarte, M. and R. Quintanilla, Decay rates of Saint-Venant type for functionally graded heat- conducting materials. International Journal of Engineering Science, 2019. 139: p. 24-41. 58. Challis, V., A. Cramer, and A. Roberts, An optimised family of anisotropic microstructures with application to functionally graded materials. International Journal of Solids and Structures, 2019. 171: p. 17-29. 59. Jiao, P., et al., A new solution of measuring thermal response of prestressed concrete bridge girders for structural health monitoring. Measurement Science and Technology, 2017. 28(8): p. 085005. 60. Liu, Z., et al., Hierarchically buckled sheath-core fibers for superelastic electronics, sensors, and muscles. Science, 2015. 349(6246): p. 400-404. 116 61. Chai, H., Contact buckling and postbuckling of thin rectangular plates. Journal of the Mechanics and Physics of Solids, 2001. 49(2): p. 209-230. 62. Chen, X., et al., Static and dynamic analysis of the postbuckling of bi-directional functionally graded material microbeams. International Journal of Mechanical Sciences, 2019. 151: p. 424-443. 63. Borchani, W., N. Lajnef, and R. Burgueño, Energy method solution for the postbuckling response of an axially loaded bilaterally constrained beam. Mechanics Research Communications, 2015. 70: p. 114-119. 64. Nejad, M.Z., A. Hadi, and A. Rastgoo, Buckling analysis of arbitrary two-directional functionally graded Euler–Bernoulli nano-beams based on nonlocal elasticity theory. International Journal of Engineering Science, 2016. 103: p. 1-10. 65. Jiao, P., et al., Static and dynamic post-buckling analyses of irregularly constrained beams under the small and large deformation assumptions. International Journal of Mechanical Sciences, 2017. 124: p. 203-215. 66. Jiao, P., W. Borchani, and N. Lajnef, Large deformation solutions to post-buckled beams confined by movable and flexible constraints: A static and dynamic analysis. International Journal of Solids and Structures, 2017. 128: p. 85-98. 67. Gnanasekaran, K., et al., 3D printing of CNT-and graphene-based conductive polymer nanocomposites by fused deposition modeling. Applied materials today, 2017. 9: p. 21-28. 68. Isakov, D., et al., 3D printed anisotropic dielectric composite with meta-material features. Materials & Design, 2016. 93: p. 423-430. 69. Tian, X., et al., Interface and performance of 3D printed continuous carbon fiber reinforced PLA composites. Composites Part A: Applied Science and Manufacturing, 2016. 88: p. 198-205. 70. Galantucci, L.M., et al., Analysis of dimensional performance for a 3D open-source printer based on fused deposition modeling technique. Procedia Cirp, 2015. 28: p. 82-87. 71. Garland, A. and G. Fadel, Design and manufacturing functionally gradient material objects with an off the shelf three-dimensional printer: challenges and solutions. Journal of Mechanical Design, 2015. 137(11). 72. Del Vescovo, D. and I. Giorgio, Dynamic problems for metamaterials: review of existing models and ideas for further research. International Journal of Engineering Science, 2014. 80: p. 153-172. 73. Rafsanjani, A., A. Akbarzadeh, and D. Pasini, Snapping mechanical metamaterials under tension. Advanced Materials, 2015. 27(39): p. 5931-5935. 74. Han, S.C., J.W. Lee, and K. Kang, A new type of low density material: Shellular. Advanced Materials, 2015. 27(37): p. 5506-5511. 75. Bauer, J., et al., Nanolattices: an emerging class of mechanical metamaterials. Advanced Materials, 2017. 29(40): p. 1701850. 117 76. Mirzaali, M., et al., Rational design of soft mechanical metamaterials: Independent tailoring of elastic properties with randomness. Applied Physics Letters, 2017. 111(5): p. 051903. 77. Janbaz, S., M. McGuinness, and A.A. Zadpoor, Multimaterial control of instability in soft mechanical metamaterials. Physical Review Applied, 2018. 9(6): p. 064013. 78. Jiao, P., et al., Mechanical metamaterial piezoelectric nanogenerator (MM-PENG): Design principle, modeling and performance. Materials & Design, 2020. 187: p. 108214. 79. Basu, S. and M. Francoeur, Penetration depth in near-field radiative heat transfer between metamaterials. Applied Physics Letters, 2011. 99(14): p. 143107. 80. Joulain, K., J. Drevillon, and P. Ben-Abdallah, Noncontact heat transfer between two metamaterials. Physical Review B, 2010. 81(16): p. 165119. 81. Mohsenizadeh, M., et al., Additively-manufactured lightweight Metamaterials for energy absorption. Materials & Design, 2018. 139: p. 521-530. 82. Li, S., et al., Fluid-driven origami-inspired artificial muscles. Proceedings of the National academy of Sciences, 2017. 114(50): p. 13132-13137. 83. Lee, J.B., et al., A mechanical metamaterial made from a DNA hydrogel. Nature nanotechnology, 2012. 7(12): p. 816-820. 84. Jin, Y., Engineering plasmonic gold nanostructures and metamaterials for biosensing and nanomedicine. Advanced Materials, 2012. 24(38): p. 5153-5165. 85. Rafsanjani, A., K. Bertoldi, and A.R. Studart, Programming soft robots with flexible mechanical metamaterials. arXiv preprint arXiv:1906.00306, 2019. 86. Taati, E., Analytical solutions for the size dependent buckling and postbuckling behavior of functionally graded micro-plates. International Journal of Engineering Science, 2016. 100: p. 45- 60. 87. Kieback, B., A. Neubrand, and H. Riedel, Processing techniques for functionally graded materials. Materials Science and Engineering: A, 2003. 362(1-2): p. 81-106. 88. Srivastava, M., et al., Design and processing of functionally graded material: review and current status of research, in 3D Printing and additive manufacturing technologies. 2019, Springer. p. 243- 255. 89. Birman, V. and L.W. Byrd, Modeling and analysis of functionally graded materials and structures. 2007. 90. Mueller, E., et al., Functionally graded materials for sensor and energy applications. Materials Science and Engineering: A, 2003. 362(1-2): p. 17-39. 91. Miyamoto, Y., et al., Functionally graded materials: design, processing and applications. Vol. 5. 2013: Springer Science & Business Media. 92. Barretta, R., et al., Buckling loads of nano-beams in stress-driven nonlocal elasticity. Mechanics of Advanced Materials and Structures, 2020. 27(11): p. 869-875. 118 93. Apuzzo, A., et al., Nonlocal strain gradient exact solutions for functionally graded inflected nano- beams. Composites Part B: Engineering, 2019. 164: p. 667-674. 94. Barretta, R., M. Čanađija, and F.M. de Sciarra, Nonlocal integral thermoelasticity: A thermodynamic framework for functionally graded beams. Composite Structures, 2019. 225: p. 111104. 95. Shen, H.-S. and Y. Xiang, Postbuckling of nanotube-reinforced composite cylindrical shells under combined axial and radial mechanical loads in thermal environment. Composites Part B: Engineering, 2013. 52: p. 311-322. 96. Spinks, G.M., et al., Mechanical properties of chitosan/CNT microfibers obtained with improved dispersion. Sensors and Actuators B: Chemical, 2006. 115(2): p. 678-684. 97. Cheng, Q., et al., High mechanical performance composite conductor: multi‐walled carbon nanotube sheet/bismaleimide nanocomposites. Advanced Functional Materials, 2009. 19(20): p. 3219-3225. 98. Esawi, A.M., et al., Effect of carbon nanotube (CNT) content on the mechanical properties of CNT- reinforced aluminium composites. Composites Science and Technology, 2010. 70(16): p. 2237- 2241. 99. Song, Y.S. and J.R. Youn, Modeling of effective elastic properties for polymer based carbon nanotube composites. Polymer, 2006. 47(5): p. 1741-1748. 100. Zhu, H., et al., Nacre-like composite films with a conductive interconnected network consisting of graphene oxide, polyvinyl alcohol and single-walled carbon nanotubes. Materials & Design, 2019. 175: p. 107783. 101. Wuite, J. and S. Adali, Deflection and stress behaviour of nanocomposite reinforced beams using a multiscale analysis. Composite Structures, 2005. 71(3-4): p. 388-396. 102. Acierno, S., et al., Experimental evaluations and modeling of the tensile behavior of polypropylene/single-walled carbon nanotubes fibers. Composite Structures, 2017. 174: p. 12-18. 103. Barretta, R., M. Čanadija, and F. Marotti de Sciarra, Modified nonlocal strain gradient elasticity for nano-rods and application to carbon nanotubes. Applied Sciences, 2019. 9(3): p. 514. 104. Gibson, L.J. and M.F. Ashby, Cellular solids: structure and properties. 1999: Cambridge university press. 105. Jiao, P. and A.H. Alavi, Buckling analysis of graphene-reinforced mechanical metamaterial beams with periodic webbing patterns. International Journal of Engineering Science, 2018. 131: p. 1-18. 106. Mohammad-Abadi, M. and A. Daneshmehr, Modified couple stress theory applied to dynamic analysis of composite laminated beams by considering different beam theories. International Journal of Engineering Science, 2015. 87: p. 83-102. 107. Zhang, H., et al., Extensible Beam-like Metastructures at the Microscale: Theoretical and Modified Hencky Bar-chain Modelling. International Journal of Mechanical Sciences, 2020: p. 105636. 119 108. Ansari, R., et al., Analytical solution for nonlinear postbuckling of functionally graded carbon nanotube-reinforced composite shells with piezoelectric layers. Composites Part B: Engineering, 2016. 90: p. 267-277. 109. Shen, H.-S., Postbuckling of nanotube-reinforced composite cylindrical shells in thermal environments, Part I: Axially-loaded shells. Composite Structures, 2011. 93(8): p. 2096-2108. 110. Shen, H.-S., Nonlinear bending of functionally graded carbon nanotube-reinforced composite plates in thermal environments. Composite Structures, 2009. 91(1): p. 9-19. 111. Abadi, M.M. and A. Daneshmehr, An investigation of modified couple stress theory in buckling analysis of micro composite laminated Euler–Bernoulli and Timoshenko beams. International Journal of Engineering Science, 2014. 75: p. 40-53. 112. Xia, W., L. Wang, and L. Yin, Nonlinear non-classical microscale beams: static bending, postbuckling and free vibration. International Journal of Engineering Science, 2010. 48(12): p. 2044-2053. 113. Jiao, P., et al., Micro-composite films constrained by irregularly bilateral walls: a size-dependent post-buckling analysis. Composite Structures, 2018. 195: p. 219-231. 114. Palmer, T., Modelling: Build imprecise supercomputers. Nature, 2015. 526(7571): p. 32-33. 115. Altrock, P.M., L.L. Liu, and F. Michor, The mathematics of cancer: integrating quantitative models. Nature Reviews Cancer, 2015. 15(12): p. 730-745. 116. Estakhri, N.M., B. Edwards, and N. Engheta, Inverse-designed metastructures that solve equations. Science, 2019. 363(6433): p. 1333-1338. 117. Corduneanu, C., Integral equations and applications. 1991: Cambridge University Press. 118. Vaid, M.K. and G. Arora, Solution of second order singular perturbed delay differential equation using trigonometric B-spline. International Journal of Mathematical, Engineering and Management Sciences, 2019. 4(2): p. 349-360. 119. Bocharov, G.A. and F.A. Rihan, Numerical modelling in biosciences using delay differential equations. Journal of Computational and Applied Mathematics, 2000. 125(1-2): p. 183-199. 120. Rihan, F.A. Delay Differential Equations in Biosciences: Parameter estimation and sensitivity analysis. in Recent Advances in Applied Mathematics and Computational Methods: Proceedings of the 2013 International Conference on Applied Mathematics and Computational Methods (Venice, Italy September 2013). 2013. 121. Longtin, A. and J.G. Milton, Complex oscillations in the human pupil light reflex with “mixed” and delayed feedback. Mathematical Biosciences, 1988. 90(1-2): p. 183-199. 122. Sommerfeld, A., Partial differential equations in physics. 1949: Academic press. 123. McCloskey, D.N., History, differential equations, and the problem of narration. History and theory, 1991. 30(1): p. 21-36. 120 124. Bauer, P., A. Thorpe, and G. Brunet, The quiet revolution of numerical weather prediction. Nature, 2015. 525(7567): p. 47-55. 125. DiPerna, R.J. and P.-L. Lions, Ordinary differential equations, transport theory and Sobolev spaces. Inventiones mathematicae, 1989. 98(3): p. 511-547. 126. Nelson, P.W. and A.S. Perelson, Mathematical analysis of delay differential equation models of HIV-1 infection. Mathematical biosciences, 2002. 179(1): p. 73-94. 127. Achdou, Y., et al., Partial differential equation models in macroeconomics. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2014. 372(2028): p. 20130397. 128. El-Ajou, A., et al., Smooth expansion to solve high-order linear conformable fractional PDEs via residual power series method: Applications to physical and engineering equations. Ain Shams Engineering Journal, 2020. 11(4): p. 1243-1254. 129. Zhukovsky, K., Operational solution for some types of second order differential equations and for relevant physical problems. Journal of Mathematical Analysis and Applications, 2017. 446(1): p. 628-647. 130. Chapra, S.C. and R.P. Canale, Numerical methods for engineers. 2010: Boston: McGraw-Hill Higher Education. 131. Caulfield, H.J. and S. Dolev, Why future supercomputing requires optics. Nature Photonics, 2010. 4(5): p. 261. 132. Zuo, S., et al., Acoustic analog computing system based on labyrinthine metasurfaces. Scientific reports, 2018. 8(1): p. 1-8. 133. Yang, T., et al., All-optical differential equation solver with constant-coefficient tunable based on a single microring resonator. Scientific reports, 2014. 4(1): p. 1-6. 134. Ratier, N. Analog computing of partial differential equations. in 2012 6th International Conference on Sciences of Electronics, Technologies of Information and Telecommunications (SETIT). 2012. IEEE. 135. Zidan, M.A., et al., A general memristor-based partial differential equation solver. Nature Electronics, 2018. 1(7): p. 411-420. 136. Pronk, S., P.L. Geissler, and D.A. Fletcher, Limits of filopodium stability. Physical review letters, 2008. 100(25): p. 258102. 137. Hirose, H., et al., Right ventricular rupture and tamponade caused by malposition of the Avalon cannula for venovenous extracorporeal membrane oxygenation. Journal of cardiothoracic surgery, 2012. 7(1): p. 1-4. 138. Andò, B., et al., Toward a Self-Powered Vibration Sensor: The Signal Processing Strategy. Energies, 2021. 14(3): p. 754. 121 139. Salem, T., et al., Functionally graded materials beams subjected to bilateral constraints: Structural instability and material topology. International Journal of Mechanical Sciences, 2021. 194: p. 106218. 140. Salem, T., et al., Maneuverable postbuckling of extensible mechanical metamaterials using functionally graded materials and carbon nanotubes. Thin-Walled Structures, 2021. 159: p. 107264. 141. Salem, T., et al., Programmable assembly of bi-walled nonuniform beams: Concept, modeling and performance. Mechanics Research Communications, 2020. 110: p. 103624. 142. Kampker, A., et al., Review on machine designs of material extrusion based additive manufacturing (AM) systems-Status-Quo and potential analysis for future AM systems. Procedia CIRP, 2019. 81: p. 815-819. 143. Chia, H.N. and B.M. Wu, Recent advances in 3D printing of biomaterials. Journal of biological engineering, 2015. 9(1): p. 1-14. 144. Mohammed, J.S., Applications of 3D printing technologies in oceanography. Methods in Oceanography, 2016. 17: p. 97-117. 122