Numerical investigation of smart auxetic three-dimensional meta-structures based on shape memory polymers via topology optimization

Abstract

Today, the human being endeavors to manufacture devices and materials capable of doing something in an intelligent way. Shape memory polymers are a series of smart materials, capable of retrieving their original shape from a temporary form by applying external stimuli, for example, heat, electricity, magnetism, light, pH, and humidity. In this research, the behavior of temperature-sensitive shape memory polymer–based structures with positive and negative Poisson’s ratio has been analyzed. The purpose is the material design of smart structures with tunable Poisson’s ratio using topology optimization. In this study, a meta-structure is designed, which is made by a smart material. Not only does this structure have shape memory effects, but also it has negative Poisson’s ratio, which can be used in new sensors, actuators, and biomedical applications. After creation of the unit cell and the representative volume element and formation of final three-dimensional structure, finite element modeling is conducted based on a thermo-visco-hyperelastic constitutive model at large deformations. Examining the behavior of structures in tensile pre-strains of 20%, 10%, and 5%, it is observed that pre-strain has no considerable effect on Poisson’s ratio, but under compressive strain of 20%, it is concluded that the type of loading is effective on Poisson’s ratio and the results are different in tension and compression modes. Finally, the influence of temperature rate on the behavior of structures is inspected, and it is concluded that the more slowly the temperature changes, the more strain or shape recovery is accomplished at a specific temperature.

Keywords

Shape memory polymers metamaterials auxetic structures (negative Poisson’s ratio)topology optimization thermo-visco-hyperelastic behavior

1. Introduction

Shape memory polymers (SMPs) are a type of smart materials capable of recovering their original shape from a temporary shape under application of external stimuli, for example, heat, magnetic, electricity, pH, light, and moisture (Baniasadi et al., 2020; Hager et al., 2015; Meng and Li, 2013; Momeni et al., 2017; Xie et al., 2016). Low weight, low cost, high deformability, controllable recovery, biocompatibility, and biodegradability are advantages of SMPs over metallic shape memory materials (Ansari et al., 2019; Li and Wang, 2011; Liu et al., 2007; Roudbarian et al., 2019). On the contrary, pure SMPs have low deformation stiffness and low recovery stress. To overcome these limitations required in some applications, the study of shape memory polymer composite (SMPC) has attracted a great deal of interest (Abbasi-Shirsavar et al., 2019; Gall et al., 2000; Tandon et al., 2010; Yarali et al., 2019). After the introduction of polynorbornene as a smart polymer by CDF Chimie Company in France in 1984, these materials were the subject of many studies (Ratna and Karger-Kocsis, 2008). Tobushi et al. (1996) conducted researches on the shape fixity and recovery behavior of SMPs. They performed uniaxial experiments on polyurethane material. In 1997, based on experimental observations, Tobushi et al. presented a one-dimensional linear viscoelastic constitutive model in small strains. Also, for the first time, the concept of stored strain was introduced and used in modeling (Tobushi et al., 1997). In 2001, they proposed a nonlinear viscoelastic model based on their previous model, suitable for large deformations (Tobushi et al., 2001).

Liu et al. (2006) proposed a constitutive model through the introduction of some phenomenological concepts, that is, phase transition and stored strain for SMPs in 2006. It is one of the basic models in the constitutive modeling of SMPs so that several models presented after have their basics borrowed from Liu et al. (2006) model. Ghosh and Srinivasa (2011) developed a model for small strains in 2011, in which SMP is assumed to consist of a permanent and a temporary network. The permanent network is responsible for shape memory effect (SME) and retrieval properties, while the temporary network is responsible for shape fixation. Employing the thermodynamics of irreversible processes, Baghani et al. (2012) developed a three-dimensional (3D) constitutive model based on the Liu et al. model in small and large deformations. Taking into account Maxwell’s generalized model, Diani et al. (2012) presented a thermo-visco-hyperelastic model at large deformations, which by performing dynamic mechanical analysis (DMA), the coefficients and parameters of epoxy material were identified. In a similar work, Arrieta et al. (2014), proposed a 3D constitutive model to determine the response of SMPs at large deformations, which has some advantages over the other models presented so far that can accurately predict the behavior of SMPs in strain recovery (free recovery or shape recovery) as well as stress recovery (force recovery or constraint recovery).

Metamaterials are artificial and engineered structures with unique properties beyond natural materials which are used in myriad applications; these properties are due to structure and not a material composition (Ahmadivand et al., 2016; Chen et al., 2016; Nouh et al., 2016; Selimis and Farsari, 2017). They are usually arranged in repeating patterns and made from assemblies of multiple elements and materials like polymers and metals. They are extensively employed in the design of broadband light bending (Ni et al., 2012), sensors (Liu et al., 2012), and switches (Zheludev and Kivshar, 2012). Auxetics are structures with a negative Poisson’s ratio (NPR), that is, when stretched, they become thicker, and when compressed, they become thinner in the direction perpendicular to the applied force. They may have unique industrial applications, such as new sealants and cushions (Hassan et al., 2009; Tanaka et al., 2001; Wang et al., 2015; Xu and Li, 2011). Planar auxetic structures have the potential to influence a wide range of applications. Rossiter et al. (2014) showed that an SMP auxetic hexachiral two-dimensional (2D) structure could be tailored to prepare a tunable stiffness behavior in its fully deployed state by changing the angle of inter-hub connections. In another work, Jacobs et al. (2012) described the design and manufacturing of a deployable antenna for deep-space missions based on hybrid honeycomb truss made of a shape memory alloy (SMA). The deployable characteristics were enhanced by the equivalent auxetic behavior of the cellular configuration. Scarpa et al. (2003) illustrated various dynamic features of open cell compliant polyurethane foam with an auxetic behavior. Foam samples had been tested in a viscoelastic analyzer tensile test machine to determine Young’s modulus and loss factor for small dynamic strains. Shim et al. (2013) designed a class of 2D soft auxetic metamaterials whose architecture could be changed in response to an external stimulus. Clausen et al. (2015) worked on topology optimized architectures with programmable Poisson’s ratio over large deformations to design their structure with 3D printing and in order to fabricate the designs and validate the numerical behavior. They created architected materials with programmable Poisson’s ratios between negative and positive ranges over large deformations of up to 20% or more.

Since Bendsøe and Kikuchi introduced a general computational material distribution method, topology optimization has been growing both in academic research and industrial applications (Bendsøe, 1989; Bendsøe et al., 1988; Bendsøe and Sigmund, 1999). Over the past decade, several similar plans with slightly different attributes have been developed in the field of topology optimization. These attributes include density-based methods (Bendsøe, 1989; Zhou and Rozvany, 1991), evolutionary procedures (Huang and Xie, 2009; Xie and Steven, 1993; Zhu et al., 2007), and level set methods (Allaire et al., 2004; Sethian and Wiegmann, 2000; Wang et al., 2003), which their goal is to find optimal structural topologies or to optimize the layout of the material under the constraints and boundary conditions (BCs). Bendsøe (2003) presented a density-based method that was associated with simplification assumptions and known as the solid isotropic material with penalization (SIMP) method.

In this study, we present a numerical investigation on smart auxetic 3D meta-structures made from SMPs, constructed through topology optimization and the SIMP method. In section 1, an introduction to SMPs and the relevant literature review was presented. The modeling, which involves SMP constitutive equations, is described in section 2. Geometrical modeling along with topology optimization is discussed in section 3. In section 4, details of the finite element analysis (FEA) are reported. Results and discussion in section 5 include comparing positive Poisson’s ratio (PPR) and NPR structures, the effect of pre-strains, loading type, and temperature rate on NPR structure. Finally, in section 6, a summary is given, and some concluding remarks are presented.

2. A brief review of SMP constitutive equations

A typical programming path of an SMP material is shown in Figure 1, which consists of four steps: (1) loading at above the glass transition temperature $T_{g}$ , where SMP is in the rubbery phase; (2) cooling the structure below $T_{g}$ while the external strain has been preserved (shape fixation), which SMP enters the glassy phase; (3) unloading while SMP remains in the glassy phase and, (4) heating (strain recovery), so that SMP returns to the rubbery phase, and the original shape is recovered, and (4′) heating (stress recovery) at constrained strain, where SMP comes back to the rubbery phase, but due to the external constraint a recovery force is generated.

Figure 1.

Stress–strain–temperature diagram depicting the thermomechanical behavior of SMP in strain and stress recovery conditions.

Commonly, there are two approaches for modeling SMPs: (1) viscoelasticity-based approach, for example, Tobushi et al. (2001) and Diani et al. (2012) models; (2) phase transition approach, for example, Liu et al. (2006) and Baghani et al. (2012) models. In this study, we employ the thermo-visco-hyperelastic model proposed by Arrieta et al. (2014) for SMPs at large deformations, which is successful in the reproduction of both the strain and stress recovery stages. This model contains three main parts: viscoelastic, hyperelastic, and thermal branches, as illustrated in Figure 2.

Figure 2.

Schematic representation of the rheological constitutive model proposed by Arrieta et al. (2014) for SMPs.

As has already been stated, the total deformation gradient tensor $F$ can be decomposed into the following parts (Holzapfel, 2002)

F = F^{m} F^{t h} (= F^{v e} F^{t h} = F^{h e} F^{t h})

(1)

where $F^{m}$ , $F^{ve}$ , $F^{he}$ , and $F^{th}$ represent mechanical, viscoelastic, hyperelastic, and thermal parts of the deformation gradient tensor, respectively. At finite strains, the elastic and viscous parts can be described by a hyperelastic model and Prony series, respectively. The total Cauchy stress tensor $σ$ can be written as follows (Brinson and Brinson, 2008)

σ (ε, t) = \int_{0}^{t} g (t - ξ) \frac{d σ^{h e}}{d ξ} d ξ

(2)

where $ε$ stands for the total strain tensor, $σ^{he}$ denotes the hyperelastic part of the stress tensor which shows the strain-dependent part, and $g (t)$ represents time-dependent part coefficients of the stress. In the following, the methods to calculate $g (t)$ and $σ^{he}$ are presented.

The viscoelastic part is assumed as the Maxwell–Wiechert model (Brinson and Brinson, 2008), also known as the generalized Maxwell model. The general form of the constitutive equations and the stress–strain relation is given as

\sum_{k = 0}^{n} p_{k} \frac{\partial^{k} σ^{ve}}{\partial t^{k}} = \sum_{k = 0}^{m} q_{k} \frac{\partial^{k} ε^{ve}}{\partial t^{k}}

(3)

where $σ^{ve}$ and $ε^{ve}$ are the viscoelastic part of the stress and the strain tensor, respectively; t denotes the time; $p_{k}$ and $q_{k}$ are constant coefficients in terms of the elastic modulus E and shear modulus $μ$ . Given that the dynamic loading is applied to sample in DMA test, by defining the complex relaxation modulus, the viscoelastic true strain and the Cauchy stress are expressed as

ε^{ve} (t) = {ε_{0}}^{ve} \exp (i ω t)

(4)

σ^{ve} (t) = {ε_{0}}^{ve} E^{*} (i ω) \exp (i ω t)

(5)

in which $ω$ is angular frequency, and $E^{*}$ is complex modulus defined as

E^{*} (i ω) = E' (ω) + iE ″ (ω)

(6)

where $E'$ and $E ″$ are storage modulus and loss modulus, respectively. Now substituting equations (4) to (6) into equation (3), and also employing the relationship between the time-dependent shear modulus $g_{i}$ and the elastic modulus, the following relations are found

E' (ω) = E_{\infty} + \sum_{i = 1}^{n} \frac{E_{i} τ_{i}^{2} ω^{2}}{1 + τ_{i}^{2} ω^{2}}

(7)

E ″ (ω) = \sum_{i = 1}^{n} \frac{E_{i} τ_{i} ω}{1 + τ_{i}^{2} ω^{2}}

(8)

g_{i} = \frac{E_{i}}{E_{\infty} + \sum_{i = 1}^{n} E_{i}}

(9)

so that $E_{\infty}$ stands for long-term elastic modulus, $τ$ is relaxation time, and n represents the number of rheological branches assumed in Maxwell–Wiechert model, which chosen to be equal to 20 in this work according to Arrieta et al. (2014).

For the hyperelastic part, which is strain-dependent and plays a significant role in large deformations, we consider Neo-Hookean strain energy function, which can be written as

W = C_{10} ({\bar{I}}_{1}^{he} - 3) + \frac{1}{D_{1}} {(J^{he} - 1)}^{2}

(10)

where $C_{10}$ and $D_{1}$ are model coefficients, ${\bar{I}}_{1}^{he}$ represents the first invariant of deviatoric left Cauchy–Green deformation tensor, and $J^{he}$ is the hyperelastic volume ratio ( $J^{he} = J / J^{th}$ , in which J and $J^{th}$ are total volume ratio and thermal volume ratio, respectively). Thus, the following equation for the strain-dependent part of the stress is computed

σ^{he} = \frac{\partial W}{\partial J^{he}} I + \frac{2}{J^{he}} dev (\frac{\partial W}{\partial {\bar{I}}_{1}^{he}} {\bar{B}}^{he})

(11)

so that dev expresses the deviatoric part of the argument, and ${\bar{B}}^{he}$ represents left Cauchy–Green deformation tensor. Considering a nearly incompressible response, a quasi-static test is performed at high temperature to find $C_{10}$ and hyperelastic material parameters.

For the thermal part of the model, considering an isothermal process, thermal part of the deformation gradient tensor can be written as

F^{th} = {(J^{th})}^{1 / 3} I \approx (1 + ε^{th}) I = [1 + \int_{T_{0}}^{T} α d T] I

(12)

where $T_{0}$ represents thermal expansion reference temperature, and $α$ is the thermal expansion coefficient, which depends on the temperature in this work. Using equations (1) and (12), ${\bar{B}}^{he}$ can be calculated and $σ^{he}$ is determined. Material parameters for the proposed SMP network are listed in Table 1. It is noted that the thermal expansion parameters are reported for $1$ and $5 ° C / \min$ temperature rates in this Table. $k_{i}$ denotes the time-dependent bulk modulus associated with relaxation time. All $k_{1} - k_{20}$ are assumed to be zero due to elastic dilatation assumption related to stress analysis in the viscoelasticity theory (Brinson and Brinson, 2008).

Table 1.

Thermo-visco-hyperelastic material parameters for the proposed SMP network adopted from Arrieta et al. (2014).

Viscoelastic material parameters of SMP
$g_{1} - g_{10} (-)$	1.8706E−01	1.5840E−01	1.3761E−01	1.2378E−01	1.1283E−01	9.8684E−02	7.5819E−02	4.7521E−02	2.4477E−02	1.1651E−02
$g_{11} - g_{20} (-)$	5.8277E−03	3.1949E−03	1.8477E−03	1.0527E−03	5.6387E−04	2.8586E−04	1.4378E−04	7.4217E−05	3.9113E−05	2.0746E−05
$τ_{1} - τ_{10} (s)$	1.26E−06	2.94E−06	6.87E−06	1.60E−05	3.75E−05	8.75E−05	2.04E−04	4.78E−04	1.12E−03	2.61E−03
$τ_{11} - τ_{20} (s)$	6.08E−03	1.42E−02	3.32E−02	7.75E−02	1.81E−01	4.23E−01	9.88E−01	2.31E + 00	5.39E + 00	1.26E + 01
Hyperelastic material parameters of SMP						WLF equation parameters of SMP
$C_{10} (Pa)$ —Instantaneous			109,730,000			WLF equation parameters of SMP
$D_{1} (P a^{- 1})$ —Instantaneous			1.745e−009			$T_{ref} (° C)$			80
General material parameter of SMP						$C_{1} (-)$			6.9
$ρ (kg / m^{3})$			1200			$C_{2} (° C)$			87.9
Thermal material parameters of SMP at $\overset{\cdot}{T} = 1 ° C / \min$
$α (1 / ° C)$		9.1718E−05	9.5553E−05	9.9354E−05	1.0520E−04				1.1222E−04	1.2083E−04
$T (° C)$		20	25	30	35				40	45
$α (1 / ° C)$		1.3032E−04	1.3899E−04	1.4697E−04	1.5258E−04				1.5652E−04	1.5943E−04
$T (° C)$		50	55	60	65				70	75
Thermal material parameters of SMP at $\overset{\cdot}{T} = 5 ° C / \min$
$α (1 / ° C)$		1.1741E−04	1.2162E−04	1.2589E−04	1.3229E−04				1.4066E−04	1.4997E−04
$T (° C)$		20	25	30	35				40	45
$α (1 / ° C)$		1.5868E−04	1.6658E−04	1.7326E−04	1.7811E−04				1.8170E−04	1.8425E−04
$T (° C)$		50	55	60	65				70	75

WLF: Williams–Landel–Ferry; SMP: shape memory polymer.

Now, to consider time and temperature dependency for non-isothermal processes in an SME path, we may use the time-temperature superposition principle (TTSP). In other words, changing in time scale, the effects of temperature variation can be taken into account. In order to compute the total stress, we can write

σ (ε, t, T) = \int_{0}^{t} g (t_{r} - ξ_{r}) \frac{d σ^{he}}{d ξ} d ξ

(13)

so that $t_{r}$ is the reduced time. For changing time scales, we can use the Williams–Landel–Ferry (WLF) equation (Williams et al., 1955)

Log (a_{T}) = \frac{- C_{1} (T (t) - T_{ref})}{C_{2} + (T (t) - T_{ref})}

(14)

ξ (t) = \int_{0}^{t} \frac{d ξ}{a_{T} (ξ)}

(15)

in which $a_{T}$ represents the shift factor, T is temperature, $T_{ref}$ expresses the reference temperature, and $C_{1}$ and $C_{2}$ are empirical constants. Thanks to the calculations of shift factors, the time-temperature superposition is identified. Considering 20 rheological branches of Maxwell–Wiechert model, we may calibrate the viscoelastic material parameters by fitting the 20 terms of equations (7) and (8) simultaneously on the master curve of experimental DMA data and identify the parameters (Arrieta et al., 2014).

To model the problem, we need to optimize the topology, construct the representative volume element (RVE), derive constitutive equations, and simulate it via FEA. The SMP examined in this study is an acrylate polymer network composition, which includes 90% molar mass of benzyl methacrylate (BMA) with 10% molar mass of poly(ethylene glycol) dimethacrylate (PEGDMA), and the glass transition temperature of this polymer network is $T_{g} = 48 ° C$ based on the Gordon–Taylor copolymer equation (Gordon and Taylor, 1952).

3. Geometrical modeling of an RVE

Currently, one of the most widely used techniques for topology optimization based on the finite element method is the SIMP method with optimality criteria algorithm (OCA), developed in the late 1980s. For the use of the SIMP-OCA method in the field of material design, the density of each element $ρ_{e}$ is defined as the amount of material to be placed within the volumetric space of that element. In a two-phase (solid-void) microstructural topology, $ρ_{e} = 0$ means the element is formed from the void (empty) and $ρ_{e} = 1$ indicates that the element is composed of the solid (filled) base material. For interpolation function, a relation exists between stiffness and density of elements. Young’s modulus of solid phase defined by this method is $E_{e} (ρ_{e}) = E_{min} + ρ_{e}^{p} (E_{0} - E_{min})$ , where $E_{0}$ is Young’s modulus of the base material, $E_{min}$ indicates Young’s modulus of void phase, taken to be a small positive number to avoid singularity problems in the numerical calculation of the stiffness tensor, and the exponent $p \geq 1$ is the parameter of penalization. The densities between zero and one are subjected to a penalization process, in order to eventually determine their status as solid and voids. An essential part of the SIMP method is to select the penalty factor to express the physical values of the problem (such as stiffness) as a function of continuous design variables. According to Bendsøe and Sigmund (1999), Hashin et al. (1963), Moshki et al. (2018), the range of the penalty factor is obtained as the following equation and is based on Poisson’s ratio of the base material $ν_{0}$ , which is equal to 0.435 for BMA in pure SMP

p \geq max {15 \frac{1 - ν_{0}}{7 - 5 ν_{0}}, \frac{3}{2} \frac{1 - ν_{0}}{1 - 2 ν_{0}}}

(16)

So far, for structures made of history-dependent materials, such as SMPs, topology optimization has not been carried out, since these categories of materials are structurally related to the loading path and thus practically cannot be optimized uniquely. Therefore, in this work, an appropriate approximation is employed to simplify the relations to achieve an approximate optimized geometry and topology. It is considered that the material behaves almost linearly in the loading step. The basis of the relationships, coding, and the creation of the topology in this research is based on Andreassen et al. (2014), Moshki et al. (2018), Sigmund (2007) to develop RVEs. Unit cells are generated using a MATLAB code. The problem-solving algorithm in this section consists of four main steps: (1) obtaining the stiffness matrix of each unit element using FEA, (2) homogenizing based on asymptotic method to get the effective elasticity matrix, (3) optimizing the topology using the SIMP method based on minimizing the function of elasticity matrix arrays, and (4) filtering the results and determining the structure. The asymptotic homogenization theory is implemented to determine the effective elasticity tensor of the optimized RVEs. The theory of asymptotic homogenization makes it possible to predict the effective properties of the periodic microstructures using FEA of the periodic unit cell. In fact, the obtained elasticity tensor of the optimized RVEs can be generalized to the overall properties due to applying the double-scale asymptotic expansion and adopting the assumption of periodicity. Therefore, the elasticity tensor of the periodic unit cell is similar to the elasticity tensor of the stacked microstructure for the periodic unit cell (Hassani et al., 1998a, 1998b). It is noteworthy to say that the applied filtration scheme produces efficient mesh-independent designs. In other words, the filtering approach eliminates the problem of mesh dependency for the optimized RVE (Sigmund, 1994; Sigmund and Petersson, 1998). Eventually, to decrease the computational cost in the current research, the minimum allowable size of the samples has been used for the topology optimization process and construction of the unit cell. In this work, Poisson’s ratio is chosen as the target function. The constraint is the volume fraction of the structure, and the BCs are considered periodically. The flowchart in Figure 3 depicts a rational explanation of the steps involved. For more details on the SIMP-OCA method, one may refer to Bendsøe (2003), Sigmund (2001), and Xia et al. (2015).

Figure 3.

The flowchart of the SIMP methodology for the optimal design of a unit cell.

Initially, as shown in Figure 4, softer materials (elements with lower volume fraction) are located at the corners of the unit cell (consist of one-eighth of a sphere with the radius equals to one-third of the unit-cell edge) to avoid uniformly distributed sensitivity field and to update the element’s densities (Moshki et al., 2018; Xia et al., 2015).

Figure 4.

The illustration of the initial guess of topology to generate unit cells.

According to the type of application, different Poisson’s ratios can be obtained by changing the volume fraction as a constraint and filtration parameters. Setting 0.75 as the volume fraction and performing the topology optimization procedure with the objective function which is defined based on the arrays minimization of the effective compliance tensor, then, considering the effective Poisson’s ratio (Xia et al., 2015), smart SMP structures with negative and PPR of −0.184 and +0.213 were produced, respectively. NPR and PPR values of −0.184 and +0.213 are the results of the numerical topology optimization algorithm and cannot be interpreted as predetermined or special values. The generated unit cells and their middle sections are illustrated in Figure 5. In the next section, we compare the pure SMP structure with PPR and NPR and examine the effect of different mechanical and thermal parameters on the results, as well.

Figure 5.

The unit cell of (a) the conventional structure (PPR) and (b) the auxetic structure (NPR). The middle sections of the unit cell of (c) the PPR structure and (d) the NPR structure.

4. FEA of SMP microstructure

Finite element simulation is employed to implement the described SMP model and to study the results of the topology optimization in ABAQUS. Then, implementing the discussed SMP constitutive model, we properly account for the thermomechanical behavior of SMP. The whole steps of the thermomechanical cycle are explained in Table 2. The standard dynamic implicit solver with backward Euler time integration method is used while the nonlinear geometry effects are considered. Regarding the periodic pattern of the main structure, we repeat the unit cell, in three main directions to get RVE size.

Table 2.

Thermomechanical programming of SMP.

Step	Condition	Temperature $(° C)$	Time intervals $(s)$ at $\overset{\cdot}{T} = 1 or 5 ° C / \min$
(1) Loading	Applying a uniaxial tensile or compressive displacement	65	50\| 20
(2) Cooling	Cooling the structure while maintaining the pre-strain	65 to 25	2400\| 480
(3) Unloading	Unloading occurs suddenly	25	50\| 20
(4) Heating	Heating the structure (strain or stress recovery)	25 to 65	2400\| 480

To construct the final 3D periodic cell (based on the asymptotic homogenization theory), cell dependency is examined for $2 \times 2 \times 2$ to $8 \times 8 \times 8$ cells, and results revealed that $6 \times 6 \times 6$ RVE cells could have enough accuracy. One-eighth of the whole structures is shown in Figure 6. Thus, the final structure built from the repetition of the unit cell contains $6 \times 6 \times 6$ cells. Here, for simplicity, BCs applied to only one sample unit cell (single element) are shown in Figure 7, while for the whole structure, the same BCs are assumed. According to Figure 7, in all simulation steps in the plane with normal vector (−1,0,0), BCs are $U_{1} = U R_{2} = U R_{3} = 0$ , in the plane (0, −1,0), BCs are $U_{2} = U R_{1} = U R_{3} = 0$ , and finally, in the plane (0,0, −1) BCs are $U_{3} = U R_{1} = U R_{2} = 0$ . In the plane (0,0,1), a tensile strain is applied along the normal vector in the loading step, while it is removed in the unloading step. Thermal BCs are explained in Table 2. The structures are meshed using structured eight-node linear hexahedral elements with reduced integration mode (C3D8R elements), as shown in Figure 8.

Figure 6.

The RVE structures (a) with PPR = +0.213 and (b) with NPR = −0.184.

Figure 7.

The boundary conditions applied to the sample unit cell (single element).

Figure 8.

The structured mesh for the (a) PPR RVE and (b) NPR RVE.

4.1. Determination of elements size

To find the appropriate size of elements and optimal mesh, the final structure should be analyzed with different mesh grids. As shown in Figure 9, the mesh dependency analysis for the final auxetic structure at the heating step of the strain recovery cycle at 20% pre-strain and $5 ° C / \min$ temperature rate is performed for four different meshes. Results for the stress recovery cycle are almost similar and follow the same pattern (the structure geometry, loading, BCs, and all of the steps except the heating step are the same). The total number of elements is doubled in each mesh, and it is revealed that the mesh grid in the first and the second stage does not have enough accuracy, but the third one is accurate enough. Therefore, the total number of elements in the final auxetic structure is 328,016, and so for the conventional structure is 289,456.

Figure 9.

The mesh dependency analysis for the auxetic structure: the strain versus time in the heating step of the strain recovery cycle.

5. Results and discussion

To validate the numerical implementation of the SMP constitutive model, Arrieta et al.’s (2014) experiments are used and a pure SMP under uniaxial loading is considered in the simulation. According to Figure 10, the model results are compared with those of experiments at the heating step of the strain recovery and stress recovery stages under 20% pre-strain at $1 ° C / \min$ heating rate, and it is concluded that the results are accurate enough to verify the model. According to the experimental data and numerical simulation of Arrieta et al. (2014), there are some differences between the experimental and numerical results which may be due to imperfections of the specimen, some empirical and human errors, and also using some assumptions to simplify the modeling of the thermomechanical behavior of SMP. As shown in Figure 10, in the strain recovery, the maximum error occurs between $T = 42 ° C$ and $T = 52 ° C$ . Also, in the stress recovery, the maximum error happens between $T = 27 ° C$ and $T = 35 ° C$ . There are two competing phenomena, which are phase transition and thermal expansion, that affect the recovery path. Thermal expansion dominates over phase transition at the beginning of the recovery and may cause such deviation from experimental results.

Figure 10.

Comparison between the model results and experiments in Arrieta et al. (2014) at the heating step under 20% pre-strain at $1 ° C / \min$ heating rate in (a) the strain recovery path and (b) the stress recovery path.

5.1. Comparing PPR and NPR structures

Here, the effect of different parameters on structures is investigated. It is noted that for all results, the xyz- or 123-coordinate system is used. Furthermore, thermomechanical programming of structures, including loading sequences, related temperatures, and time intervals, is according to Figure 1 and Table 2, which are determined by numbers 1 to 4. First, a comparison between the PPR and the NPR structures behavior under 20% tensile pre-strain and at $5 ° C / \min$ temperature rate is carried out. As mentioned before, the applied load is along the z-axis, and the structures have almost the same behavior along the x-axis and y-axis due to the symmetry. The strain along the z-axis, the strain along the x-axis, Poisson’s ratio, and the stress along the z-axis versus time are shown in Figure 11 in both the strain and stress recovery paths. Figure 11 reveals that the strain along the loading, that is, $ε_{33}$ , is almost the same for all the structures, but the strain perpendicular to the loading, that is, $ε_{11}$ , is different. The reason is that both structures are restricted along the loading due to external displacement; however, they are free perpendicular to the loading, and because of their different Poisson’s ratio, their behavior must be different in this direction. It is noted that Poisson’s ratio defined as $ν_{31} = - ε_{11} / ε_{33}$ is selected as the optimization objective to be different for the structures. The average Poisson’s ratio is calculated and compared to that of the topology optimization as shown in Table 3, and it is observed that all of the relative errors are less than 10%. The trend of the stress along the loading, that is, $σ_{33}$ , is the same for all cases, but values are different, arising from the nonlinearity in the geometries and structures, and the stiffness matrices.

Figure 11.

Comparing PPR and NPR structures with pure solid SMP in the strain recovery and stress recovery paths: (a) the strain along the loading $ε_{33}$ , (b) the strain perpendicular to the loading $ε_{11}$ , (c) Poisson’s ratio $ν_{31}$ , and (d) the stress along the loading $σ_{33}$ .

Table 3.

Results of the average Poisson’s ratio of FE simulation and topology optimization for PPR and NPR structures in the strain recovery and stress recovery.

Structures	Topology optimization	Strain recovery		Stress recovery
		FE simulation	Relative error (%)	FE simulation	Relative error (%)
NPR	−0.184	−0.166	9.78	−0.167	9.24
PPR	+0.213	+0.234	9.86	+0.233	9.39
Pure SMP	–	+0.475	–	+0.467	–

5.2. The effect of pre-strain on the response of NPR structure

From this section on, for the sake of briefness, the results are reported only for NPR structure. Here, the behavior of NPR structure under 20%, 10%, and 5% pre-strains at $5 ° C / \min$ temperature rate is scrutinized. The applied load and BCs are the same as before. $ε_{33}$ , $ε_{11}$ , $ν_{31}$ , and $σ_{33}$ versus time are depicted in Figure 12 in the strain recovery and stress recovery paths as well. This figure reveals that the amplitude of $ε_{33}$ , $ε_{11}$ , and $σ_{33}$ is declined with almost the same ratio by reducing the pre-strain. The more the amplitude of pre-strain is decreased, the more the amplitude of $ε_{33}$ , $ε_{11}$ , and $σ_{33}$ is reduced. It is clear that the pre-strain almost does not affect Poisson’s ratio, because as the amount of pre-strain along the loading is declined, the amount of strain perpendicular to the loading is also decreased proportionally, and this causes Poisson’s ratio to be unchanged during the path, while the pre-strain values are changing.

Figure 12.

The effect of pre-strain on NPR structure in the strain and stress recovery paths: (a) the strain along the loading $ε_{33}$ ,(b) the strain perpendicular to the loading $ε_{11}$ , (c) Poisson’s ratio $ν_{31}$ , and (d) the stress along the loading $σ_{33}$ .

5.3. The effect of loading type on NPR structure

In this section, a comparison between the behavior of NPR structure under 20% compressive loading and 20% tensile loading at $5 ° C / \min$ temperature rate is conducted. The applied load and BCs are the same as before. $ε_{33}$ , $ε_{11}$ , $ν_{31}$ , and $σ_{33}$ versus time are depicted in Figure 13 in the strain recovery as well as stress recovery. As one may observe from Figure 13(a), $ε_{33}$ for the compressive loading is symmetric compared to the tensile one, but this symmetry does not occur in $ε_{11}$ as shown in Figure 13(b). The reason is that, in the cooling step, since the structure is not bounded in the x-direction and the temperature decreases simultaneously, the structure is contracted and $ε_{11}$ is reduced in both cases. In the heating step, the temperature reaches above $T_{g}$ , and $ε_{11}$ returns to zero in the strain recovery path, so the strain and the shape of the structure recover in both cases. Moreover, the stress and force are recovered in the stress recovery path. As a result, Poisson’s ratio behavior becomes different, as illustrated in Figure 13(c). According to Figure 13(d), in the cooling step, as the structure is bounded in the loading direction, by reducing the temperature, a positive stress is developed in the structure to counteract the decrease in the temperature and stands for non-contraction of the structure in this direction, while the other steps are like before. Thus, we can have smart SMP structures that reach a maximum and minimum volume in their temporary shape in tension and compression, respectively, which can be used for a variety of potential applications, especially in biomedical devices.

Figure 13.

Effect of loading type on NPR structure in the strain and stress recovery paths: (a) the strain along the loading $ε_{33}$ , (b) the strain perpendicular to the loading $ε_{11}$ , (c) Poisson’s ratio $ν_{31}$ , and (d) the stress along the loading $σ_{33}$ .

5.4. The effect of temperature rate on NPR structure

Finally, to study the effect of temperature rate, the behavior of NPR structure under 20% pre-strain at $5$ and $1 ° C / \min$ temperature rates is studied. As it is noted, thermal expansion coefficients depend on temperature and temperature rate, according to Arrieta et al. (2014), and data for thermal expansion for other temperature rates were not available. Again, the applied load and BCs are the same as before. $ε_{33}$ and $ε_{11}$ versus the temperature, $ν_{31}$ versus the normalized time, and $σ_{33}$ versus the temperature are shown in Figure 14 in both the strain and stress recovery paths. According to Figure 14(a), in the heating step of the strain recovery, at lower temperature rates, the structure has a lower strain at a certain temperature, and the diagram is shifted to the left, because a higher strain recovery is achieved at a certain temperature. This event occurs in Figure 14(b), as well. As illustrated in Figure 14(c), the temperature rate does not affect Poisson’s ratio, since the ratio of the strain perpendicular to the loading to the strain along the loading is unchanged. As shown in Figure 14(d), in the cooling step, as the temperature changes more slowly, the structure has lower stress at a certain temperature, because it senses the temperature changes faster, so smaller stress is developed in the structure. All these discussions stand correct for the stress recovery. For example, in Figure 14(d), in the heating step, as the temperature rises more slowly, less amount of stress is applied to the structure.

Figure 14.

Effect of temperature rate on NPR structure in the strain and stress recovery paths: (a) the strain along the loading $ε_{33}$ , (b) the strain perpendicular to the loading $ε_{11}$ , (c) Poisson’s ratio $ν_{31}$ , and (d) the stress along the loading $σ_{33}$ .

6. Summary and conclusion

In this research, the material design and the FEA of a smart auxetic 3D meta-structures based on SMPs were numerically examined. These shape memory metamaterials can be used in biomedical devices and structures, in which two rational designing and SME have been coupled. It means the auxetic or conventional properties of structures can be controlled via an external applied thermal source. For instance, these structures could be used as bone repairing tools or deployable scaffold in the body. In this work, first, tunable PPR and NPR unit cells were created by optimizing the topology using the SIMP method. Then, RVEs and the final 3D structures were generated by repeating the unit cells. Thereafter, FE modeling was done based on the thermo-visco-hyperelastic constitutive model at large deformations for SMPs. Results of topology optimization and FE simulation were compared, and the relative differences were less than 10%. The behavior of PPR and NPR smart structures was studied and evaluated. Also, the effect of mechanical and thermal parameters was explored which could be summarized as follows. (1) The amplitude of the strain and the stress is lowered with almost the same ratio by reducing the applied pre-strain; however, pre-strain nearly does not affect Poisson’s ratio. (2) Different types of loading may result in a different strain, stress, and Poisson’s ratio due to the history-dependent thermomechanical behavior of the SMP-based structures. (3) The temperature rate influences the behavior of smart structures, and the more slowly the temperature changes, the more strain and shape recover at a specific temperature.

Footnotes

Acknowledgements

The authors are grateful for the research support of the Iran National Science Foundation (INSF).

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Majid Baniassadi

Mostafa Baghani

References

Abbasi-Shirsavar

Baghani

Taghavimehr

, et al. (2019) An experimental–numerical study on shape memory behavior of PU/PCL/ZnO ternary blend. Journal of Intelligent Material Systems and Structures 30: 116–126.

Ahmadivand

Karabiyik

Pala

(2016) Plasmonic photodetectors. Photodetectors 2016; 157–193.

Allaire

Jouve

Toader

A-MJ

(2004) Structural optimization using sensitivity analysis and a level-set method. Journal of Computational Physics 194: 363–393.

Andreassen

Lazarov

Sigmund

(2014) Design of manufacturable 3D extremal elastic microstructure. Mechanics of Materials 69: 1–10.

Ansari

Golzar

Baghani

, et al. (2019) Force recovery evaluation of thermo-induced shape-memory polymer stent: material, process and thermo-viscoelastic characterization. Smart Materials and Structures 28: 095022.

Arrieta

Diani

Gilormini

(2014) Experimental characterization and thermoviscoelastic modeling of strain and stress recoveries of an amorphous polymer network. Mechanics of Materials 68: 95–103.

Baghani

Naghdabadi

Arghavani

, et al. (2012) A thermodynamically-consistent 3 D constitutive model for shape memory polymers. International Journal of Plasticity 35: 13–30.

Baniasadi

Maleki-Bigdeli

M-A

Baghani

(2020) Force and multiple-shape-recovery in shape-memory-polymers under finite deformation torsion-extension. Smart Materials and Structures 29: 055011.

Bendsøe

(1989) Optimal shape design as a material distribution problem. Structural Optimization 1: 193–202.

10.

Bendsøe

(2003) Sigmund, Topology Optimization-Theory, Methods and Applications. Berlin: Springer.

11.

Bendsøe

Kikuchi

(1988) Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering 71: 197–224.

12.

Bendsøe

Sigmund

(1999) Material interpolation schemes in topology optimization. Archive of Applied Mechanics 69: 635–654.

13.

Brinson

(2008) Polymer engineering science and viscoelasticity, 2nd edition. New York: Springer.

14.

Chen

Huang

(2016) A design of active elastic metamaterials for control of flexural waves using the transformation method. Journal of Intelligent Material Systems and Structures 27: 1337–1347.

15.

Clausen

Wang

Jensen

, et al. (2015) Topology optimized architectures with programmable Poisson’s ratio over large deformations. Advanced Materials 27: 5523–5527.

16.

Diani

Gilormini

Frédy

, et al. (2012) Predicting thermal shape memory of crosslinked polymer networks from linear viscoelasticity. International Journal of Solids and Structures 49: 793–799.

17.

Gall

Mikulas

Munshi

, et al. (2000) Carbon fiber reinforced shape memory polymer composites. Journal of Intelligent Material Systems and Structures 11: 877–886.

18.

Ghosh

Srinivasa

(2011) A two-network thermomechanical model of a shape memory polymer. International Journal of Engineering Science 49: 823–838.

19.

Gordon

Taylor

(1952) Ideal copolymers and the second-order transitions of synthetic rubbers. I. Non-crystalline copolymers. Journal of Applied Chemistry 2: 493–500.

20.

Hager

Bode

Weber

, et al. (2015) Shape memory polymers: past, present and future developments. Progress in Polymer Science 49–50: 3–33.

21.

Hashin

Shtrikman

(1963) A variational approach to the theory of the elastic behaviour of multiphase materials. Journal of the Mechanics and Physics of Solids 11: 127–140.

22.

Hassan

Scarpa

Mohamed

(2009) In-plane tensile behavior of shape memory alloy honeycombs with positive and negative Poisson’s ratio. Journal of Intelligent Material Systems and Structures 20: 897–905.

23.

Hassani

Hinton

(1998a) A review of homogenization and topology optimization I—homogenization theory for media with periodic structure. Computers & Structures 69: 707–717.

24.

Hassani

Hinton

(1998b) A review of homogenization and topology opimization II—analytical and numerical solution of homogenization equations. Computers & Structures 69: 719–738.

25.

Holzapfel

(2002) Nonlinear solid mechanics: a continuum approach for engineering science. Meccanica 37: 489–490.

26.

Huang

Xie

(2009) Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials. Computational Mechanics 43: 393.

27.

Jacobs

Coconnier

DiMaio

, et al. (2012) Deployable auxetic shape memory alloy cellular antenna demonstrator: design, manufacturing and modal testing. Smart Materials and Structures 21: 075013.

28.

Wang

(2011) Studies on the shape frozen/recovery behaviors of SMP-based sandwich structure. Journal of Intelligent Material Systems and Structures 22: 1605–1612.

29.

Liu

Qin

Mather

(2007) Review of progress in shape-memory polymers. Journal of Materials Chemistry 17: 1543–1558.

30.

Liu

Gall

Dunn

, et al. (2006) Thermomechanics of shape memory polymers: uniaxial experiments and constitutive modeling. International Journal of Plasticity 22: 279–313.

31.

Liu

Xia

Sun

, et al. (2012) Visible transmission response of nanoscale complementary metamaterials for sensing applications. Nanotechnology 23: 275503.

32.

Meng

(2013) A review of stimuli-responsive shape memory polymer composites. Polymer 54: 2199–2221.

33.

Momeni

Liu

(2017) A review of 4D printing. Materials & Design 122: 42–79.

34.

Moshki

Ghazavizadeh

Atai

, et al. (2018) 3D-printable unit cell design for cubic and orthotropic porous microstructures using topology optimization based on optimality criteria algorithm. International Journal of Applied Mechanics 10: 1850060.

35.

Emani

Kildishev

, et al. (2012) Broadband light bending with plasmonic nanoantennas. Science 335: 427–427.

36.

Nouh

Aldraihem

Baz

(2016) Periodic metamaterial plates with smart tunable local resonators. Journal of Intelligent Material Systems and Structures 27: 1829–1845.

37.

Ratna

Karger-Kocsis

(2008) Recent advances in shape memory polymers and composites: a review. Journal of Materials Science 43: 254–269.

38.

Rossiter

Takashima

Scarpa

, et al. (2014) Shape memory polymer hexachiral auxetic structures with tunable stiffness. Smart Materials and Structures 23: 045007.

39.

Roudbarian

Baniasadi

Ansari

, et al. (2019) An experimental investigation on structural design of shape memory polymers. Smart Materials and Structures 28: 095017.

40.

Scarpa

Ciffo

Yates

(2003) Dynamic properties of high structural integrity auxetic open cell foam. Smart Materials and Structures 13: 49.

41.

Selimis

Farsari

(2017) Laser-Based 3D Printing and Surface Texturing. Available at: http://scitechconnect.elsevier.com/wp-content/uploads/2016/09/3-s2.0-B9780128035818091712-main.pdf

42.

Sethian

Wiegmann

(2000) Structural boundary design via level set and immersed interface methods. Journal of Computational Physics 163: 489–528.

43.

Shim

Shan

Košmrlj

, et al. (2013) Harnessing instabilities for design of soft reconfigurable auxetic/chiral materials. Soft Matter 9: 8198–8202.

44.

Sigmund

(1994) Design of material structures using topology optimization. PhD Thesis, Technical University of Denmark, Lyngby.

45.

Sigmund

(2001) A 99 line topology optimization code written in Matlab. Structural and Multidisciplinary Optimization 21: 120–127.

46.

Sigmund

(2007) Morphology-based black and white filters for topology optimization. Structural and Multidisciplinary Optimization 33: 401–424.

47.

Sigmund

Petersson

(1998) Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Structural Optimization 16: 68–75.

48.

Tanaka

Todoh

Oki

, et al. (2001) Material structure design of stress-dependent adaptive stiffness. Journal of Intelligent Material Systems and Structures 12: 229–234.

49.

Tandon

Goecke

Cable

, et al. (2010) Environmental durability of fabric-reinforced shape-memory polymer composites. Journal of Intelligent Material Systems and Structures 21: 1365–1381.

50.

Tobushi

Hara

Yamada

, et al. (1996) Thermomechanical properties in a thin film of shape memory polymer of polyurethane series. Smart Materials and Structures 5: 483–491.

51.

Tobushi

Hashimoto

Hayashi

, et al. (1997) Thermomechanical constitutive modeling in shape memory polymer of polyurethane series. Journal of Intelligent Material Systems and Structures 8: 711–718.

52.

Tobushi

Okumura

Hayashi

, et al. (2001) Thermomechanical constitutive model of shape memory polymer. Mechanics of Materials 33: 545–554.

53.

Wang

Chang

Y-H

Chen

, et al. (2015) Designable dual-material auxetic metamaterials using three-dimensional printing. Materials & Design 67: 159–164.

54.

Wang

Guo

(2003) A level set method for structural topology optimization. Computer Methods in Applied Mechanics and Engineering 192: 227–246.

55.

Williams

Landel

Ferry

(1955) The temperature dependence of relaxation mechanisms in amorphous polymers and other glass-forming liquids. Journal of the American Chemical Society 77: 3701–3707.

56.

Xia

Breitkopf

(2015) Design of materials using topology optimization and energy-based homogenization approach in Matlab. Structural and Multidisciplinary Optimization 52: 1229–1241.

57.

Xie

Huang

Leng

, et al. (2016) Thermoset shape memory polymers and their composites. Journal of Intelligent Material Systems and Structures 27: 2433–2455.

58.

Xie

Steven

(1993) A simple evolutionary procedure for structural optimization. Computers & Structures 49: 885–896.

59.

(2011) A shape memory polymer based syntactic foam with negative Poisson’s ratio. Materials Science and Engineering: A 528: 6804–6811.

60.

Yarali

Baniassadi

Baghani

(2019) Numerical homogenization of coiled carbon nanotube reinforced shape memory polymer nanocomposites. Smart Materials and Structures 28: 035026.

61.

Zheludev

Kivshar

YSJ

(2012) From metamaterials to metadevices. Nature Materials 11: 917.

62.

Zhou

Rozvany

(1991) The COC algorithm, Part II: topological, geometrical and generalized shape optimization. Computer Methods in Applied Mechanics and Engineering 89: 309–336.

63.

Zhu

Zhang

Qiu

(2007) Bi-directional evolutionary topology optimization using element replaceable method. Computational Mechanics 40: 97