Introduction

First, please watch the following video.

The structure shown in the movie uses a metamaterial (patent No. JP7229616) that has the function of deforming out-of-plane due to in-plane loading.
This article introduces the mechanism of this metamaterial and a simple analysis method.

References

• Imperial College London. Shell metamaterials. (Ref. Sep. 2023)
• Rob Hewson. Multiscale Actuated Shells Structures. 15th World Congress on Computational Mechanics, 2022

Basic structure

Honeycomb, Re-entrant honeycomb

All of the structures in the video in this article are honeycomb-based shapes. The honeycomb shape is widely used as a pattern of densely packed hexagonal cells that allows efficient use of space.
Normal honeycomb has a positive Poisson's ratio. However when the angle is changed, the Poisson's ratio becomes negative. This is called a re-entrant honeycomb which is one of the typical auxetic structures (structure with negative Poisson's ratio).

Let us try to calculate the Poisson's ratio when the angle is varied by actually using finite element analysis. Fix the central position of the diagonal beam so as not to change the size of the unit structure enclosed by the red line in the figure below, and vary the angle \(θ\). When \(h=1.6w, t_h=t_s=0.05w\), the macroscopic Poisson's ratios in each direction are as shown below. \(ν_{LT}, ν_{TL}\) are Poisson's ratios representing strain in the T-direction under L-directional loading and strain in the L-direction under T-directional loading, respectively. We can see that Poisson's ratio is positive when \(θ\) is positive, and Poisson's ratio is negative when \(θ\) is negative.

Varying the dimensions of the unit structure allows for a continuous change in macro-physical properties.

Deformation mechanism

The basic structure of metamaterials that deform out-of-plane due to in-plane loading have shapes that are smoothly connected so that one face is a honeycomb structure with a positive Poisson's ratio and the other face is a re-entrant honeycomb structure with a negative Poisson's ratio.
Under tension load, the honeycomb side shrinks in the orthogonal direction and the re-entrant honeycomb side expands, resulting in overall bending deformation.

Homogenization method

Homogenization methods are methods that convert the properties of materials and structures with complex microstructures into simpler equivalent homogeneous material properties. This makes it possible to understand the behavior of the material as a whole. The homogenization method also simplifies the analysis of materials with complex microstructures and reduces computational costs. The calculation of Poisson's ratio shown in the graph shown above can also be said to be homogenization in a two-dimensional plane.
In this case, since the model has a thin plate-like shape, the shell element is used to create the model. Assuming elastic deformation, the properties of the shell element can be expressed as follows

\(\begin{bmatrix} [N] \\ [M] \end{bmatrix}
=
\begin{bmatrix} 
[A] & [B] \\ [B] & [D]
\end{bmatrix}
\begin{bmatrix} 
[\varepsilon] \\ [\kappa]
\end{bmatrix}
=
\begin{bmatrix} 
a_{11} & a_{12} & a_{16} & b_{11} & b_{12} & b_{16} \\
a_{12} & a_{22} & a_{26} & b_{12} & b_{22} & b_{26} \\
a_{16} & a_{26} & a_{66} & b_{16} & b_{26} & b_{66} \\
b_{11} & b_{12} & b_{16} & d_{11} & d_{12} & d_{16} \\
b_{12} & b_{22} & b_{26} & d_{12} & d_{22} & d_{26} \\
b_{16} & b_{26} & b_{66} & d_{16} & d_{26} & d_{66} \\
\end{bmatrix}
\begin{bmatrix} 
\varepsilon_{x} \\ \varepsilon_{x}\\ \gamma_{xy}\\
\kappa_{x} \\ \kappa_{y} \\ \kappa_{xy}
\end{bmatrix}\)

where \([N]\) and \([M]\) are the in-plane loads and moments per unit length, and \([\varepsilon]\) and \([\kappa]\) are the in-plane strain and curvature at the neutral plane. \([A], [B], [D]\) are called the in-plane stiffness matrix, coupling matrix, and bending stiffness matrix, respectively.
Transverse shear deformation is also considered in this case.

\(\begin{bmatrix} [V]\end{bmatrix}
=
\begin{bmatrix} 
[K]
\end{bmatrix}
\begin{bmatrix} 
[\gamma]
\end{bmatrix}
=
\begin{bmatrix} 
a_{55} & a_{45}\\
a_{45} & a_{44}
\end{bmatrix}
\begin{bmatrix} 
\gamma_{zx}\\ \gamma_{yz}
\end{bmatrix}\)

\([V]\) represents the transverse shear force per unit length and \([\gamma]\) represents the transverse shear strain. For isotropic materials, the transverse shear stiffness matrix \([K]\) is a diagonal matrix whose components are the modulus of transverse elasticity \(G\) multiplied by the plate thickness \(t\).

\(\begin{bmatrix} 
[K]
\end{bmatrix}
=
\begin{bmatrix} 
Gt & 0\\
0 & Gt
\end{bmatrix}\)

Calculating these stiffness matrices enables a simplified simulation of their behavior on a macroscale.

Classical Laminate Theory

In this article, we will calculate the stiffness matrix using Classical Laminate Theory(CLT). CLT is frequently used for composite laminates and is a method for calculating the elastic properties (in-plane stiffness and bending stiffness) of a plate under plane stress conditions easily.
Given the in-plane strain \([\varepsilon]\) at the neutral plane and the curvature \([\kappa]\), the strain at the plane at \(z\) from the neutral plane can be given by

\(\begin{bmatrix}
\varepsilon_{x} \\ \varepsilon_{y} \\ \gamma_{xy}
\end{bmatrix}+
z\begin{bmatrix}
\kappa_{x} \\ \kappa_{y} \\ \kappa_{xy}
\end{bmatrix}\)

Considering that the integral value of the stress is balanced by the external force, the in-plane load \([N]\) and moment \([M]\) per unit length can be represented using the in-plane stiffness matrix \([Q]\) at position \(z\) as follows. In the following, it is assumed that the neutral plane and the center plane coincide.

\(N_i = 
\int_{-t/2}^{t/2}\sigma_{i}dz=
\left(\int_{-t/2}^{t/2}Q_{ij}dz\right)\varepsilon_{i}+
\left(\int_{-t/2}^{t/2}Q_{ij}zdz\right)\kappa_{i}\\
M_i = 
\int_{-t/2}^{t/2}\sigma_{i}zdz=
\left(\int_{-t/2}^{t/2}Q_{ij}zdz\right)\varepsilon_{i}+
\left(\int_{-t/2}^{t/2}Q_{ij}z^2dz\right)\kappa_{i}\)

Based on this, the in-plane stiffness matrix \([A]\), the coupling matrix \([B]\), and the bending stiffness matrix \([D]\) described above can be calculated respectively as follows

\(A_{ij} = \int_{-t/2}^{t/2}Q_{ij}dz =
\sum_{k=1}^{N}Q_{ij}^{(k)}(z_k-z_{k-1})\\
B_{ij} = \int_{-t/2}^{t/2}Q_{ij}zdz =
\frac{1}{2}\sum_{k=1}^{N}Q_{ij}^{(k)}(z_k^2-z_{k-1}^2)\\
D_{ij} = \int_{-t/2}^{t/2}Q_{ij}z^2dz =
\frac{1}{3}\sum_{k=1}^{N}Q_{ij}^{(k)}(z_k^3-z_{k-1}^3)\)

In this way, by dividing the geometry into multiple layers in the thickness direction and calculating the in-plane stiffness matrix for each geometry, it is possible to calculate the stiffness matrix of the shell element.

The transverse shear stiffness is calculated as the sum of those values for each layer.

\(K_{ij} = \sum_{k=1}^{N}K_{ij}^{(k)}\)

Calculation of in-plane stiffness

The in-plane stiffness of extruded shapes can be calculated by 2-dimensional finite element analysis. Since the geometry in this case exhibits orthotropic anisotropy, with Young's modulus and Poisson's ratio in each direction as \(E_L\), \(E_T\), \(ν_{LT}\), \(ν_{TL}\) and transverse modulus as \(G_{LT}\), the load-displacement relationship per unit length in plane stress state is represented as follows.

\(\begin{bmatrix} \sigma \end{bmatrix}
=
\begin{bmatrix} Q \end{bmatrix}
\begin{bmatrix} 
\varepsilon \end{bmatrix}\\
\begin{bmatrix} 
\sigma_L \\ \sigma_T \\ \tau_{LT}
\end{bmatrix}
=
\begin{bmatrix} 
\frac{E_L}{1-\nu_{LT}\nu_{TL}} &
\frac{\nu_{LT}E_T}{1-\nu_{LT}\nu_{TL}} & 0 \\
\frac{\nu_{TL}E_L}{1-\nu_{LT}\nu_{TL}} &
\frac{E_T}{1-\nu_{LT}\nu_{TL}} &
0\\
0&0&G_{LT}
\end{bmatrix}
\begin{bmatrix} 
\varepsilon_{L} \\ \varepsilon_{T}\\ \gamma_{LT}
\end{bmatrix}\)

Each component of the stiffness matrix is calculated from the reaction forces obtained by analyzing the unit structure under three forced displacement conditions: \(L\) direction, \(T\) direction, and shear. Symmetrical or anti-symmetrical boundaries were utilized and a quarter model was used in the calculations.

With \(h = 1.6w, t_h = t_s = 0.05w\), and varying the angle \(θ\), each component of the in-plane stiffness matrix is as follows

These results are used to calculate the properties of the metamaterial.

Analysis with shell elements

Now that the in-plane stiffness has been calculated, the stiffness matrices of the shell element are calculated and a finite element analysis is performed.

Element stiffness matrix

In this geometry, when the angles of the beam at \(z=-t/2, t/2\) are \(θ_1, θ_2\), the angle of the beam at thickness \(z\) is expressed as follows.

\(\theta = \arctan{\left(
\frac{t/2-z}{t}\tan\theta_1 + \frac{z-t/2}{t}\tan\theta_2\right)}\)

Set \(t=2 {\text{mm}}\) and \(θ_1=30°, θ_2=30°\). The unit structure is divided into 30 layers, and the stiffness matrix calculated by CLT is as follows

\(\begin{bmatrix} 
[A] & [B] \\ [B] & [D]
\end{bmatrix}
= \begin{bmatrix} 
0.435 \text{ N/mm} & 0.060 \text{ N/mm} & 0 & 0.032 \text{ N} & 0.419 \text{ N}& 0 \\
 & 1.850 \text{ N/mm} & 0 & 0.419 \text{ N}& 0.049\text{ N}& 0\\
 &  & 0.009\text{ N/mm} & 0 & 0 & 0.003\text{ N} \\
 & & & 0.221\text{ Nmm} & 0.032 \text{ Nmm}& 0\\
 & \text{sym.}& & & 0.477 \text{ Nmm}& 0\\
 & & & & & 0.004 \text{ Nmm}\\
\end{bmatrix}\\
\begin{bmatrix} 
[K]
\end{bmatrix}
=
\begin{bmatrix} 
1.552 \text{ N/mm} & 0\\
0 & 1.552\text{ N/mm}
\end{bmatrix}\)

The transverse shear stiffness of each layer was approximated by the area ratio.

\(a_{44} = a_{55} = Gt\times(Cross-sectional area of the unit structure)/wh\)

Also, using the compliance matrix \([S]\), which is obtained by the inverse of the stiffness matrix, we see that the \(X\) directional load \(N_x\) generates the curvature \(\kappa_y\) in the \(Y\) direction.

\(\begin{bmatrix} 
[\varepsilon] \\ [\kappa]
\end{bmatrix}=
\begin{bmatrix} [S] \end{bmatrix}
\begin{bmatrix} [N] \\ [M] \end{bmatrix}\\
\begin{bmatrix} 
\varepsilon_{x} \\ \varepsilon_{x}\\ \gamma_{xy}\\
\kappa_{x} \\ \kappa_{y} \\ \kappa_{xy}
\end{bmatrix}=
\begin{bmatrix} 
15.07\text{ mm/N}\\
-0.14\text{ mm/N}\\
 0 \\
0.03\text{ N}^{-1}\\
13.22\text{ N}^{-1}\\
0 \\
\end{bmatrix}
N_x\)

Finite element analysis

The analysis is performed on solid and shell elements and the results are compared. Material non-linearity is not considered in the calculations, but geometric non-linearity is.
When the edge is displaced in the L direction for the geometry shown in the following figure, the results for the reaction force and out-of-plane displacement are as follows

• Solid: reaction force \(9.38×10^{-2} \text N\), out-of-plane displacement \(13.14 \text {mm}\)
• Shell: reaction force \(8.24×10^{-2} \text N\), out-of-plane displacement \(13.54 \text {mm}\)
  
  
  The values obtained are close enough to give a rough picture of the behavior.
  The following factors can be considered as possible sources of error.
• The stiffness matrix of the shell element was calculated under the assumption that the neutral plane is at mid-thickness at all positions, but in reality the neutral plane does not coincide with the mid-thickness plane.
• The mesh accuracy in the solid element analysis is not sufficient.

The above is a simple case, but calculations can be performed with different boundary conditions and geometries. The metamaterial is used to create a butterfly-like shape in the figure below and check its behavior when the center is displaced. A 1/4 model was used in the calculation. Note that the analysis of solid elements could not be performed for this shape because of the huge calculation time required for convergence.

The approximate shape of the prototype under tensile load is similar to the results of the shell element analysis, indicating that the deformation behavior can be calculated.

Summary

Finite element analysis was conducted to understand the deformation behavior of metamaterial structures that have the function of deforming out-of-plane due to in-plane loading.

It was shown that homogenization with shell elements can be used to simulate the deformation behavior easily. It is also shown that the stiffness matrix of the shell element can be approximated from 2D analytical results by using classical lamination theory.