A modified analytical model to study the sensing performance of a flexible capacitive tactile sensor array

For analytical modeling, the sensing unit is simplified into a three-layer plate structure: lower PET, middle solid PDMS and upper PET films, as shown in figure 3(a). The thicknesses of the three layers are 25 μm, 20 μm and 25 μm, respectively. As shown in figure 3(a), an xyz coordinate system is defined with the origin point located at the center of the structure and z-axis perpendicular to the upper PET film. According to [25], the top PDMS bump can be replaced by applying a uniform pressure q₀ on the upper PET film with an area of 3000  ×  3000 μm at the center. The coordinates of the plate structure are normalized as $X=2x/4000,$ $Y=2y/4000,$ and $Z=2z/70$ with $X,Y,Z\in \left[-1,1\right].$

Zoom In Zoom Out Reset image size

Since the Young's modulus of the PET (4000 MPa) is much larger than that of the PDMS (0.55 MPa), the deformation of the PET film is much smaller than the PDMS film. This will produce strong nonlinear relation between the deformation and coordinates of the z-axis. To achieve accurate deformation prediction of the sensing unit, a nonlinear function which describes the deformation tendency along the lower PET, PDMS and upper PET films of the three-layer plate structure in the z-axis will be obtained by analyzing the deformation of the PET and PDMS films under the same pressure as detailed in appendix A. The function is obtained as $f(Z)=[322/(1+{{\text{e}}^{-7Z-0.03}})-0.3],$ which is shown in equation (A.6). Using the nonlinear function, the general model of the three-layer structure can be expressed as

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} u=\left(1-{{X}^{2}}\right)\left(1-{{Y}^{2}}\right)\left(\frac{322}{1+{{\text{e}}^{-7Z-0.03}}}-0.3\right)\underset{i=1}{\overset{I}{\sum}} \underset{j=1}{\overset{J}{\sum}} \underset{k=1}{\overset{K}{\sum}} C_{ijk}^{u}{{P}_{i}}{{P}_{j}}{{P}_{k}} \\ v=\left(1-{{X}^{2}}\right)\left(1-{{Y}^{2}}\right)\left(\frac{322}{1+{{\text{e}}^{-7Z-0.03}}}-0.3\right)\underset{i=1}{\overset{I}{\sum}} \underset{j=1}{\overset{J}{\sum}} \underset{k=1}{\overset{K}{\sum}} C_{ijk}^{v}{{P}_{i}}{{P}_{j}}{{P}_{k}} \\ w=\left(1-{{X}^{2}}\right)\left(1-{{Y}^{2}}\right)\left(\frac{322}{1+{{\text{e}}^{-7Z-0.03}}}-0.3\right)\underset{i=1}{\overset{I}{\sum}} \underset{j=1}{\overset{J}{\sum}} \underset{k=1}{\overset{K}{\sum}} C_{ijk}^{w}{{P}_{i}}{{P}_{j}}{{P}_{k}} \\ \end{array}\right.\end{eqnarray} \tag{ 1 }$$

where u, v, w are the displacement components in the x-, y- and z-axis, respectively; $C_{ijk}^{u},C_{ijk}^{v},C_{ijk}^{w}$ are the coefficients to be solved; i, j and k are the serial index and I, J and K are the maximal number of the serial index; ${{P}_{i}},{{P}_{j}},{{P}_{k}}$ are Chebyshev polynimials [26].

The relations between the strain and displacement are expressed in normalized coordinate as

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} {{\varepsilon}_{x}}=\frac{\partial u}{\partial x}=\frac{2}{0.004}\cdot \frac{\partial u}{\partial X} \\ {{\varepsilon}_{y}}=\frac{\partial v}{\partial y}=\frac{2}{0.004}\cdot \frac{\partial v}{\partial Y} \\ {{\varepsilon}_{z}}=\frac{\partial w}{\partial z}=\frac{2}{0.00007}\cdot \frac{\partial w}{\partial Z} \\ {{\gamma}_{xy}}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}=\frac{2}{0.004}\cdot \frac{\partial u}{\partial Y}+\frac{2}{0.004}\cdot \frac{\partial v}{\partial X} \\ {{\gamma}_{xz}}=\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}=\frac{2}{0.00007}\cdot \frac{\partial u}{\partial Z}+\frac{2}{0.004}\cdot \frac{\partial w}{\partial X} \\ {{\gamma}_{yz}}=\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}=\frac{2}{0.00007}\cdot \frac{\partial v}{\partial Z}+\frac{2}{0.004}\cdot \frac{\partial w}{\partial Y} \\ \end{array}\right.\end{eqnarray} \tag{ 2 }$$

The material properties of the sensing unit are obtained by stretching and compression tests. The Young's modulus of the PDMS with 10:1 and 20:1 are measured as 2.03 MPa and 0.55 MPa, respectively, while the PET is 4000 MPa. Accordingly, the Poisson's ratio of the PDMS with 10:1, 20:1 and PET are 0.49, 0.49 and 0.4, respectively. The Lame's first parameter λ and second parameter G [27] can be calculated as

$$\begin{eqnarray}&&\lambda =\frac{E\nu}{\left(1+\nu \right)\left(1-2\nu \right)}=\left\{\begin{array}{*{35}{l}} 5.71\times {{10}^{9}},\,\,\,\,\,\,\,\,\, Z\in \left[-1,-10/35\right] \\ 1.85\times {{10}^{7}},\,\,\,\,\,\,\,\,\, Z\in \left(-10/35,10/35\right] \\ 5.71\times {{10}^{9}},\,\,\,\,\,\,\,\,\, Z\in \left(10/35,1\right] \\ \end{array}\right.\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&G=\frac{E}{2\left(1+\nu \right)}=\left\{\begin{array}{*{35}{l}} 1.43\times {{10}^{9}},\,\,\,\,\,\,\,\,\, Z\in \left[-1,-10/35\right] \\ 1.85\times {{10}^{5}},\,\,\,\,\,\,\,\,\, Z\in \left(-10/35,10/35\right] \\ 1.43\times {{10}^{9}},\,\,\,\,\,\,\,\,\, Z\in \left(10/35,1\right] \\ \end{array}\right.\end{eqnarray} \tag{ 3b }$$

where E and v are the Young's modulus and Poisson's ratio, respectively; $Z\in \left[-1,-10/35\right],$ $Z\in \left[-10/35,10/35\right],$ and $Z\in \left[10/35,1\right]$ represent the coordinate ranges in the z-axis of the lower PET, PDMS and upper PET films, respectively.

Based on the Ritz method [28], the energy of the plate structure is minimized by setting the partial derivative of the energy with respect to the undetermined coefficients $C_{ijk}^{u},$ $C_{ijk}^{v}$ and $C_{ijk}^{w}$ equal to zero, so we have

$$\begin{eqnarray}&&\begin{array}{*{20}{l}} {\partial \left\{ {\frac{1}{2}\int_{ - 1}^1 {\int_{ - 1}^1 {\int_{ - 1}^1 {\left[ {\left(\lambda + 2G\right)\left(\varepsilon _x^2 + \varepsilon _y^2 + \varepsilon _z^2\right)} \right.} } } } \right.}\\ {\,\,\,\,\,\,\, + 2\lambda \left({\varepsilon _x}{\varepsilon _y} + {\varepsilon _x}{\varepsilon _z} + {\varepsilon _y}{\varepsilon _z}\right)}\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left. { + 8.75 \times {{10}^{ - 12}}G\left(\gamma _{xy}^2 + \gamma _{xz}^2 + \gamma _{yz}^2\right)} \right]{\rm{d}}X{\rm{d}}Y{\rm{d}}Z}\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left. { - 4 \times {{10}^{ - 6}}\int_{ - 0.75}^{0.75} {\int_{ - 0.75}^{0.75} {{q_0}w\left(X,Y,1\right){\rm{d}}X{\rm{d}}Y} } } \right\}/\partial C_{ijk}^\theta = 0} \end{array}\end{eqnarray} \tag{ 4 }$$

where $\theta =u,v,w.$

According to equations (1)–(4), the governing equation of this three-layer plate structure can be expressed as

$$\begin{eqnarray}\left[\begin{array}{c} \left[{{K}^{uu}}\right] & \left[{{K}^{uv}}\right] & \left[{{K}^{uw}}\right] \\ \left[{{K}^{vu}}\right] & \left[{{K}^{vv}}\right] & \left[{{K}^{vw}}\right] \\ \left[{{K}^{wu}}\right] & \left[{{K}^{wv}}\right] & \left[{{K}^{ww}}\right] \\ \end{array}\right]\left[\begin{array}{c} \left[{{C}^{u}}\right] \\ \left[{{C}^{v}}\right] \\ \left[{{C}^{w}}\right] \\ \end{array}\right]=\left[\begin{array}{c} 0 \\ 0 \\ \left[Q\right] \\ \end{array}\right]\end{eqnarray} \tag{ 5 }$$

where ${{K}^{\alpha \beta}}$ ($\alpha ,\ \beta =u,\ v,\ w$) and Q are the stiffness sub-matrices and the loading vector, respectively. ${{C}^{\theta}}$ ($\theta =u,\ v,\ w$) are the column vectors containing undetermined coefficients of the general model in equation (1). The sub-matrices, loading vectors and column vectors in equation (5) are calculated in appendix B1.

By substituting I = J = K = 20, q₀ = 100 kPa into equations (1)–(4) and solving equation (5), the undetermined coefficients in column vectors ${{C}^{\theta}}$ ($\theta =u,\ v,\ w$) are obtained and the general model in equation (1) is solved. Theoretically, setting a larger value of I, J and K will increase the accuracy because more terms of functions are utilized [28]. However, a practical limit to the number of terms that are utilized always exists because of larger truncation error when computing the stiff-matrices in equation (5) and it may generate sparse matrices that cannot be solved. These factors will reduce the accuracy and stability of the model. To accurately predict the deformation without increasing the values of I, J and K, a modified model will be developed by studying the central, edge and corner regions separately.

The three-layer plate structure can be divided into nine regions: one central, four edges, and four corners, as in figure 3(a). Due to the symmetrical location of the four edge and corner regions, one edge and one corner region together with the central region are selected for consideration. The selected edge region locates at the negative side of the x-axis and the selected corner region locates at both negative sides of x- and y-axes, as shown in figures 3(c) and (d). The dimensions of the central, edge and corner regions are 2000  ×  2000  ×  70 μm, 2000  ×  1000  ×  70 μm and 1000  ×  1000  ×  70 μm, respectively. The distance between the edge of the central region and the loading region is 500 μm, which is seven times larger than the thickness of the plate (70 μm). According to the Saint–Venant theory [27], it can be inferred that the pressure in the corners and edges of the loading region will not deform the central region. Similarly, pressure in the corners of the loading region will not deform the edge region. For abbreviation, three local coordinate systems (${{x}_{1}}{{y}_{1}}{{z}_{1}},$ ${{x}_{2}}{{y}_{2}}{{z}_{2}}$ and ${{x}_{3}}{{y}_{3}}{{z}_{3}}$) are created for these three regions. In each local coordinate system, the origin point locates at the center and the z-axis is perpendicular to the upper PET film, as shown in figures 3(b)–(d).

In the central region, a uniform pressure q₀ is applied on the upper PET film with a loading area of 2000  ×  2000 μm, as shown in figure 3(b). The coordinates are normalized as ${{X}_{1}}=2{{x}_{1}}/2000,$ ${{Y}_{1}}=2{{y}_{1}}/2000,$ and ${{Z}_{1}}=2{{z}_{1}}/70$ with ${{X}_{1}},{{Y}_{1}},{{Z}_{1}}\in \left[-1,1\right].$ Since the uniform pressure is applied on the whole region surface, the deformation of the central region in the arbitrary plane parallel to the z₁-axis is zero in the x₁- and y₁-axis. The bottom of the lower PET film is fixed. Thus, the displacement field function of the central region can be expressed as

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} {{u}_{1}}=0 \\ {{v}_{1}}=0 \\ {{w}_{1}}=\left(\frac{322}{1+{{\text{e}}^{-7{{Z}_{1}}-0.03}}}-0.3\right)\underset{i=1}{\overset{I}{\sum}} \underset{j=1}{\overset{J}{\sum}} \underset{k=1}{\overset{K}{\sum}} C_{ijk}^{w}{{P}_{i}}{{P}_{j}}{{P}_{k}} \\ \end{array}\right.\end{eqnarray} \tag{ 6 }$$

where $C_{ijk}^{w}$ are the coefficients to be calculated.

When uniform normal force is applied, the strain of the central region is

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} {{\varepsilon}_{{{z}_{1}}}}=\frac{\partial w}{\partial {{z}_{1}}}=\frac{2}{0.00007}\cdot \frac{\partial {{w}_{1}}}{\partial {{Z}_{1}}} \\ {{\gamma}_{{{x}_{1}}{{z}_{1}}}}=\frac{\partial w}{\partial {{x}_{1}}}=\frac{2}{0.002}\cdot \frac{\partial {{w}_{1}}}{\partial {{X}_{1}}} \\ {{\gamma}_{{{y}_{1}}{{z}_{1}}}}=\frac{\partial w}{\partial {{y}_{1}}}=\frac{2}{0.002}\cdot \frac{\partial {{w}_{1}}}{\partial {{Y}_{1}}} \\ \end{array}\right.\end{eqnarray} \tag{ 7 }$$

Using the procedure described in the general model, total energy of the central region is minimized as

$$\begin{eqnarray}&&\begin{array}{l} \partial \left\{ {\frac{1}{2}\int_{ - 1}^1 {\int_{ - 1}^1 {\int_{ - 1}^1 {\left[ {\left( {\lambda + 2G} \right)\varepsilon _{{z_1}}^2 + G\left( {\gamma _{xz}^2 + \gamma _{yz}^2} \right)} \right]} } } } \right.\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\, \times 3.5 \times {10^{ - 11}}{\rm{d}}{X_1}{\rm{d}}{Y_1}{\rm{d}}{Z_1}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left. { - \int_{ - 1}^1 {\int_{ - 1}^1 {{q_0}w({X_1},{Y_1},1) \times {{10}^{ - 6}}{\rm{d}}{X_1}{\rm{d}}{Y_1}} } } \right\}/\partial C_{ijk}^w = 0 \end{array}\end{eqnarray} \tag{ 8 }$$

According to equations (3) and (6)–(8), the governing equation of the central region can be obtained as

$$\begin{eqnarray}&&{[^1}{K^{ww}}]{[^1}{C^w}] = {[^1}Q]\end{eqnarray} \tag{ 9 }$$

where $^{1}{{K}^{ww}}$ and $^{1}Q$ are the stiffness matrix and the loading vector, respectively. $^{1}{{C}^{w}}$ is the column vector containing undetermined coefficients of the model in equation (6). The matrix and vectors in equation (9) are described in appendix B2.

For the edge region, a uniform pressure q₀ is applied on the upper PET film with an area of 2000  ×  500 μm, as shown in figure 3(c). The coordinates can be normalized as ${{X}_{2}}=2{{x}_{2}}/1000,$ ${{Y}_{2}}=2{{y}_{2}}/2000,$ and ${{Z}_{2}}=2{{z}_{2}}/70$ with ${{X}_{2}},{{Y}_{2}},{{Z}_{2}}\in \left[-1,1\right].$ The bottom (${{Z}_{2}}=-1$) and left side (${{X}_{2}}=-1$) of the edge region are fixed. The right side (${{X}_{2}}=1$), which is connected to the central region, has the same boundary condition as the central region. Due to the uniform pressure applied along the y₂-axis being constant, the deformation in arbitrary plane perpendicular to the y₂-axis is zero in the y₂-axis. Thus, the boundary conditions are such that ${{u}_{2}},{{v}_{2}},{{w}_{2}}=0$ when ${{X}_{2}},{{Z}_{2}}=-1,$ ${{u}_{2}},{{v}_{2}}=0$ when ${{X}_{2}}=1,$ and ${{v}_{2}}=0$ when ${{Y}_{2}}=\pm\ 1.$ Then, the displacement field function of the edge region can be expressed as

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} {{u}_{2}}=\left(1-X_{2}^{2}\right)\left(\frac{322}{1+{{\text{e}}^{-7{{Z}_{2}}-0.03}}}-0.3\right)\underset{i=1}{\overset{I}{\sum}} \underset{j=1}{\overset{J}{\sum}} \underset{k=1}{\overset{K}{\sum}} C_{ijk}^{u}{{P}_{i}}{{P}_{j}}{{P}_{k}} \\ {{v}_{2}}=0 \\ {{w}_{2}}=\left(1+{{X}_{2}}\right)\left(\frac{322}{1+{{\text{e}}^{-7{{Z}_{2}}-0.03}}}-0.3\right)\underset{i=1}{\overset{I}{\sum}} \underset{j=1}{\overset{J}{\sum}} \underset{k=1}{\overset{K}{\sum}} C_{ijk}^{w}{{P}_{i}}{{P}_{j}}{{P}_{k}} \\ \end{array}\right.\end{eqnarray} \tag{ 10 }$$

where $C_{ijk}^{u}$ and $C_{ijk}^{w}$ are the undetermined coefficients.

According to equation (10), the strain field of the edge region is given as

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} {{\varepsilon}_{{{x}_{2}}}}=\frac{\partial u}{\partial {{x}_{2}}}=\frac{2}{0.001}\cdot \frac{\partial u}{\partial {{X}_{2}}} \\ {{\varepsilon}_{{{z}_{2}}}}=\frac{\partial w}{\partial {{z}_{2}}}=\frac{2}{0.00007}\cdot \frac{\partial w}{\partial {{Z}_{2}}} \\ {{\gamma}_{{{x}_{2}}{{y}_{2}}}}=\frac{\partial u}{\partial {{y}_{2}}}=\frac{2}{0.002}\cdot \frac{\partial u}{\partial {{Y}_{2}}} \\ {{\gamma}_{{{x}_{2}}{{z}_{2}}}}=\frac{\partial u}{\partial {{z}_{2}}}+\frac{\partial w}{\partial {{x}_{2}}}=\frac{2}{0.00007}\cdot \frac{\partial u}{\partial {{Z}_{2}}}+\frac{2}{0.001}\cdot \frac{\partial w}{\partial {{X}_{2}}} \\ {{\gamma}_{{{y}_{2}}{{z}_{2}}}}=\frac{\partial w}{\partial {{y}_{2}}}=\frac{2}{0.002}\cdot \frac{\partial w}{\partial {{Y}_{2}}} \\ \end{array}\right.\end{eqnarray} \tag{ 11 }$$

Similarly, the following equations minimize the energy of the edge region as

$$\begin{eqnarray}&&\begin{array}{l} \partial \left\{ {\frac{1}{2}\int_{ - 1}^1 {\int_{ - 1}^1 {\int_{ - 1}^1 {\left[ {\left( {\lambda + 2G} \right)\left( {\varepsilon _{{x_2}}^2 + \varepsilon _{{z_2}}^2} \right) + 2\lambda {\varepsilon _{{x_2}}}{\varepsilon _{{z_2}}}} \right.} } } } \right.\\ \,\,\,\,\,\,\,\,\,\,\,\left. { + 1.75 \times {{10}^{ - 11}}G\left( {\gamma _{{x_2}{y_2}}^2 + \gamma _{{x_2}{z_2}}^2 + \gamma _{{y_2}{z_2}}^2} \right)} \right]{\rm{d}}{X_2}{\rm{d}}{Y_2}{\rm{d}}{Z_2}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left. { - \int_{ - 1}^1 {\left[ {5 \times {{10}^{ - 8}}\int_0^1 {{q_0}w\left( {{X_2},{Y_2},1} \right){\rm{d}}{X_2}} } \right]{\rm{d}}{Y_2}} } \right\}/\partial C_{ijk}^\theta = 0 \end{array}\end{eqnarray} \tag{ 12 }$$

where $\theta =u,\ w.$

Based on equations (3) and (10)–(12), the following governing equation of the edge region is

$$\begin{eqnarray}\left[ {\begin{array}{*{20}{c}} {{[^2}{K^{uu}}]}&{{[^2}{K^{uw}}]}\\ {{[^2}{K^{wu}}]}&{{[^2}{K^{ww}}]} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{[^2}{C^u}]}\\ {{[^2}{C^w}]} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 0\\ {{[^2}Q]} \end{array}} \right]\end{eqnarray} \tag{ 13 }$$

where $^{2}{{K}^{\alpha \beta}}$ ($\alpha ,\ \beta =u,\ w$) and $^{2}Q$ are the stiffness sub-matrices and the loading vector, respectively. $^{2}{{C}^{\theta}}\ \left(\theta =u,\ w\right)$ are the column vectors containing undetermined coefficients of the model in equation (10). The sub-matrices and vectors mentioned above are calculated in appendix B2.

For the corner region, a uniform pressure q₀ is applied with an area of 500  ×  500 μm, as shown in figure 3(d). The coordinates are normalized as ${{X}_{3}}=2{{x}_{3}}/1000,\ {{Y}_{3}}=2{{y}_{3}}/1000\ \text{and}\ {{Z}_{3}}=2{{z}_{3}}/70$ with ${{X}_{3}},{{Y}_{3}},{{Z}_{3}}\in \left[-1,1\right].$ The bottom (${{Z}_{3}}=-1$), left side (${{X}_{3}}=-1$) and back side (${{Y}_{3}}=-1$) of the corner region are fixed. The right side (${{X}_{3}}=1$) and front side (${{Y}_{3}}=1$) are connected to the edge region. Thus, the boundary conditions of the corner region are such that ${{u}_{3}}=0,{{v}_{3}}=0,{{w}_{3}}=0$ when ${{X}_{3}},{{Y}_{3}},{{Z}_{3}}=-1,\ {{u}_{3}}=0$ when ${{X}_{3}}=1$ and ${{v}_{3}}=0$ when ${{Y}_{3}}=1.$ The displacement field function of the corner region is

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} {{u}_{3}}=\left(1-{{X}_{3}}^{2}\right)\left(1+{{Y}_{3}}\right)\left(\frac{322}{1+{{\text{e}}^{-7{{Z}_{3}}-0.03}}}-0.3\right)\underset{i=1}{\overset{I}{\sum}} \underset{j=1}{\overset{J}{\sum}} \underset{k=1}{\overset{K}{\sum}} C_{ijk}^{u}{{P}_{i}}{{P}_{j}}{{P}_{k}} \\ {{v}_{3}}=\left(1+{{X}_{3}}\right)\left(1-{{Y}_{3}}^{2}\right)\left(\frac{322}{1+{{\text{e}}^{-7{{Z}_{3}}-0.03}}}-0.3\right)\underset{i=1}{\overset{I}{\sum}} \underset{j=1}{\overset{J}{\sum}} \underset{k=1}{\overset{K}{\sum}} C_{ijk}^{v}{{P}_{i}}{{P}_{j}}{{P}_{k}} \\ {{w}_{3}}=\left(1+{{X}_{3}}\right)\left(1+{{Y}_{3}}\right)\left(\frac{322}{1+{{\text{e}}^{-7{{Z}_{3}}-0.03}}}-0.3\right)\underset{i=1}{\overset{I}{\sum}} \underset{j=1}{\overset{J}{\sum}} \underset{k=1}{\overset{K}{\sum}} C_{ijk}^{w}{{P}_{i}}{{P}_{j}}{{P}_{k}} \\ \end{array}\right.\end{eqnarray} \tag{ 14 }$$

where $C_{ijk}^{u},$, $C_{ijk}^{v}\ \text{and}\ C_{ijk}^{w}$ are the undermined coefficients.

The strain fields of the corner region are 3D functions that can be calculated as

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} {{\varepsilon}_{{{x}_{3}}}}=\frac{\partial u}{\partial {{x}_{3}}}=\frac{2}{0.001}\cdot \frac{\partial u}{\partial {{X}_{3}}} \\ {{\varepsilon}_{{{y}_{3}}}}=\frac{\partial v}{\partial {{y}_{3}}}=\frac{2}{0.001}\cdot \frac{\partial v}{\partial {{Y}_{3}}} \\ {{\varepsilon}_{{{z}_{3}}}}=\frac{\partial w}{\partial {{z}_{3}}}=\frac{2}{0.00007}\cdot \frac{\partial w}{\partial {{Z}_{3}}} \\ {{\gamma}_{{{x}_{3}}{{y}_{3}}}}=\frac{\partial u}{\partial {{y}_{3}}}+\frac{\partial v}{\partial {{x}_{3}}}=\frac{2}{0.001}\cdot \frac{\partial u}{\partial {{Y}_{3}}}+\frac{2}{0.001}\cdot \frac{\partial v}{\partial {{X}_{3}}} \\ {{\gamma}_{{{x}_{3}}{{z}_{3}}}}=\frac{\partial u}{\partial {{z}_{3}}}+\frac{\partial w}{\partial {{x}_{3}}}=\frac{2}{0.00007}\cdot \frac{\partial u}{\partial {{Z}_{3}}}+\frac{2}{0.001}\cdot \frac{\partial w}{\partial {{X}_{3}}} \\ {{\gamma}_{{{y}_{3}}{{z}_{3}}}}=\frac{\partial v}{\partial {{z}_{3}}}+\frac{\partial w}{\partial {{y}_{3}}}=\frac{2}{0.00007}\cdot \frac{\partial v}{\partial {{Z}_{3}}}+\frac{2}{0.001}\cdot \frac{\partial w}{\partial {{Y}_{3}}} \\ \end{array}\right.\end{eqnarray} \tag{ 15 }$$

Using the same procedure as described in the central and edge regions, the total energy of the corner region is minimized as

$$\begin{eqnarray}&&\begin{array}{l} \partial \left\{ {\frac{1}{2}\int_{ - 1}^1 {\int_{ - 1}^1 {\int_{ - 1}^1 {\left[ {\left( {\lambda + 2G} \right)\left( {\varepsilon _{{x_3}}^2 + \varepsilon _{{y_3}}^2 + \varepsilon _{{z_3}}^2} \right)} \right.} } } } \right.\\ \,\,\,\,\,\,\, + 2\lambda \left( {{\varepsilon _{{x_3}}}{\varepsilon _{{y_3}}} + {\varepsilon _{{x_3}}}{\varepsilon _{{z_3}}} + {\varepsilon _{{y_3}}}{\varepsilon _{{z_3}}}} \right)\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left. { + 8.75 \times {{10}^{ - 12}}G\left( {\gamma _{{x_3}{y_3}}^2 + \gamma _{{x_3}{z_3}}^2 + \gamma _{{y_3}{z_3}}^2} \right)} \right]{\rm{d}}{X_3}{\rm{d}}{Y_3}{\rm{d}}{Z_3}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left. { - 2.5 \times {{10}^{ - 7}}\int_0^1 {\int_0^1 {{q_0}w\left( {{X_3},{Y_3},1} \right){\rm{d}}{X_3}{\rm{d}}{Y_3}} } } \right\}/\partial C_{ijk}^\theta = 0 \end{array}\end{eqnarray} \tag{ 16 }$$

where $\theta =u,\ v,\ w.$

According to equations (3) and (14)–(16), the governing equation of the corner region can be expressed as

$$\begin{eqnarray}\left[ {\begin{array}{*{20}{c}} {{[^3}{K^{uu}}]}&{{[^3}{K^{uv}}]}&{{[^3}{K^{uw}}]}\\ {{[^3}{K^{vu}}]}&{{[^3}{K^{vv}}]}&{{[^3}{K^{vw}}]}\\ {{[^3}{K^{wu}}]}&{{[^3}{K^{wv}}]}&{{[^3}{K^{ww}}]} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{[^3}{C^u}]}\\ {{[^3}{C^v}]}\\ {{[^3}{C^w}]} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 0\\ 0\\ {{[^3}Q]} \end{array}} \right]\end{eqnarray} \tag{ 17 }$$

where $^{3}{{K}^{\alpha \beta}}\ \left(\alpha ,\ \beta =u,\ v,\ w\right)$ and $^{3}Q$ are the stiffness sub-matrices and the loading vector, respectively. $^{3}{{C}^{\theta}}\ \left(\theta =u,\ v,\ w\right)$ are the column vectors containing undetermined coefficients of the model in equation (14). The sub-matrices and vectors in equation (17) are calculated in appendix B2.

By setting I = J = K = 20, q₀ = 100 kPa and solving the equations (9), (13) and (17), the displacement field functions of the central, edge and corner regions can be obtained. Thus, the deformation of the sensing units can be predicted by the combination of the displacement functions which are calculated in the central, edge and corner regions.

For the developed tactile sensor, the deformations of the lower and upper electrodes (on the lower and upper PET films) need to be studied. They determine the capacitance change of the sensing unit. Here, we utilize the modified analytical model developed in the above section to study the deformations of the upper and lower PET films under uniform normal force. Figures 4(a) and (b) show the calculated deformations of the bottom surface of the upper PET film and the top surface of the lower PET film, respectively. Two cross-sections P–P' (crossing central and edge regions) and Q–Q' (crossing edge and corner regions) are created with distances of 2 mm and 0.5 mm away from the side of sensing unit, respectively. Also, figures 5(a)–(b) show the deformations in the cross-sections P–P' and Q–Q'. On the bottom surface of the upper PET film (figure 4(a)), the maximum deformation is 0.22 μm, about 600 times larger than that of the top surface of lower PET film (0.35 nm in figure 4(b)). This result is expected because the upper PET film is laminated on a soft PDMS film, thus can be deformed easily when external force applied. Therefore, when calculating the capacitance between the upper and lower electrodes, the deformation of the lower electrode could be omitted, because it has little effect on the capacitance change. Comparing figures 5(a) with (b), the deformations in cross-section P–P' are larger than those in Q–Q', because the loading area on the central region is larger than that on the edge and corner regions, as shown in figures 3(b)–(d). In figures 5(a) and (b), slight waves are observed in the edge and corner regions. When force is applied, the PDMS film in the central region will be squeezed towards the edge and corner regions and generate small waves on the edges of the loading area. Results demonstrate that the modified analytical model can accurately describe the deformation of the sensing unit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Besides uniform normal force, uneven normal force and shear force usually occur on the sensing unit in actual applications. In this study, the uneven normal force and shear force, which can be expressed by 2D loading functions as in [25], are applied on the sensing unit. Then, the modified model is used to calculate the deformations of the sensing unit under uneven normal force and shear force. The calculation procedures are the same as those described in section 3.2 by replacing the uniform normal force q₀ with the uneven normal force and shear force loading functions, respectively. Figures 6(a) and 7(a) show the calculated deformations of the upper PET film under uneven normal force and shear force, respectively. For the applied uneven normal force, the upper PET film will generate uneven deformation. For comparison, the deformations at the cross-section P–P' of the upper PET film when applying uniform and uneven normal forces with the same magnitude are shown in figure 6(b). The applied uneven normal force generates much greater deformation at the central region. For the applied shear force, the upper PET film is stretched on the negative x-axis side and is compressed on the other side, as shown in figure 7(b). The deformation predictions by the modified model are the same as the results obtained in [25]. According to the calculated deformations in figures 6 and 7, the percentage capacitance changes under uneven normal force and shear force can be calculated as 1.06% and 0.04%, respectively, while the uniform normal force causes 0.98% capacitance change. Results demonstrate that the uneven normal force, compared with uniform normal force, will generate greater capacitance change and that the sensing unit is not sensitive to shear force.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 8(a)–(d) shows the fabrication process of the flexible capacitive tactile sensor array. To enhance sensitivity, a softer PDMS (Sylgard 184, A:B = 20:1 in weight) is used to make the dielectric layer, while the bump is made of harder PDMS (Sylgard 184, A:B = 10:1 in weight). For the fabrication of the electrode layer, a 25 μm-thick PET film is firstly adhered to a 3 inch silicon wafer, as shown in figure 8(a). Then, Cu electrodes are magnetron sputtered on the PET film. PDMS (20:1) is spin-coated on the electrodes at 7000 rpm for 40 s to form a 10 μm-thick PDMS film. The PDMS film is cured at room temperature for 24 h to generate the electrode layer. To form a capacitor, two electrode layers are fabricated with the PDMS surface activated by the oxygen plasma. Then, these two electrode layers are positioned and bonded, forming a capacitive layer, as shown in figure 8(b).

Zoom In Zoom Out Reset image size

Figure 8(c) shows the process of fabricating the PDMS bump. The liquid-state PDMS (10:1) is casted and cured at 80 °C for 3 h in a stainless steel mold, which is manufactured by electrical discharge machining (EDM), to replicate the shape of the bump and generate a PDMS mold. After peeling off, the PDMS mold is treated by the release agent (mixture of ethanol and detergent) for 30 min. Then, the liquid-state PDMS (10:1) is casted into the PDMS mold and laminated by the fabricated capacitive layer. The PDMS bump is cured at room temperature for 24 h and then peeled off from the PDMS mold together with the fabricated capacitive layer. The final capacitive tactile sensor array is shown in figure 8(d).

Figure 9(a) shows the fabricated capacitive tactile sensor array, which has 64 sensing units. The overall dimensions for the sensor array are 35  ×  35 mm. Figure 9(b) shows the close-up view of the sensing units; the spatial resolution of the sensor array is 4 mm. The sensor array features high flexibility, as shown in figure 9(c).

Zoom In Zoom Out Reset image size

To study the sensing performance of the sensor array, experiments have been conducted. Figure 10(a) shows the overall experimental setup. The sensor array is placed on an electronic balance, which has 0.01 g resolution and 2000 g measurement range. A round loading bar with diameter of 4 mm is fixed onto a three-axis positioner (Newport M460P) to apply forces on the PDMS bump of the sensing unit, as shown in figure 10(b). The applied force will be measured by the electronic balance and the corresponding capacitance change of the sensing unit is measured by a LCR meter (Applent AT2818), which has a resolution of 0.01 fF and measurement range of 1 F.

Zoom In Zoom Out Reset image size

Experiment procedures are as follows: (1) Initial capacitances of the sensor array are measured, without the loading bar contacting with the sensing unit; (2) For each sensing unit, 0 to 10 N force is applied by the loading bar and the capacitance changes are recorded by the LCR meter. For each sensing unit, three repeated tests have been conducted.

Figure 11(a) shows the measured initial capacitances of all the 64 sensing units. The average and standard derivation of the initial capacitances are 10.64  ±  1.03 pF. The minimum and maximum capacitances are 7.82 pF (row 2, column 3) and 12.63 pF (row 8, column 5), respectively. Figure 11(b) shows the measured percentage capacitance changes versus applied normal forces of the sensing unit (at row 3, column 6). C and C₀ are the capacitance induced by the applied force and initial capacitance, respectively. The upper and lower bounds of the error bars present the maximum and minimum of the measured percentage capacitance changes. A least squares linear curve fits the experiment data, as shown in figure 11(b). The slope of the fitted line (1.24%/N) represents the sensitivity of this sensing unit.

Zoom In Zoom Out Reset image size

By using the same method, the sensitivities of all the 64 sensing units are obtained, as shown in figure 12. The sensitivities vary from 1.19 to 1.26%/N and the average and standard derivation of the sensitivities are 1.24  ±  0.01%/N. Table 1 summarizes the measured initial capacitances and experimental fitted sensitivities of the 64 sensing units.

Zoom In Zoom Out Reset image size

To verify the developed models, the capacitance of the sensing unit is calculated using the general model (equation (1)) and the modified model (equations (6), (10) and (14)). The deformation of the bottom surface of the upper PET film can be calculated as

$$\begin{eqnarray}&&{{w}_{F}}\left(X,Y,Z\right)=\frac{F/{{S}_{F}}}{{{q}_{0}}}w\left(X,Y,Z\right)\end{eqnarray} \tag{ 18 }$$

where F is the applied normal force ranging from 0 to 10 N, ${{w}_{F}}\left(X,Y,Z\right)$ and $w\left(X,Y,Z\right)$ are the deformations induced by F and q₀, respectively; ${{S}_{F}}=3\times 3$ mm² and ${{q}_{0}}=100$ kPa.

Then, the capacitance of the sensing unit can be calculated as

$$\begin{eqnarray}&&C=\frac{\varepsilon}{4\pi k}\int_{-0.75}^{0.75} \int_{-0.75}^{0.75} \frac{1}{20\times {{10}^{-6}}+{{\omega}_{F}}\left(X,Y,10/35\right)}\text{d}X\text{d}Y\end{eqnarray} \tag{ 19 }$$

where $\varepsilon =2.68$ is relative permittivity of the PDMS film [29] and k = 9.0  ×  10⁹ is static electricity constant.

By setting ${{w}_{F}}\left(X,Y,Z\right)=0$ in equation (19), the initial capacitance of the sensing unit can be easily calculated as 10.66 pF, which is almost consistent with the measured average initial capacitance (10.64 pF). According to equations (18)–(19), the general model and modified model are utilized to calculate the capacitance change of the sensing unit and then compared to the experiment data, as shown in figure 11(b). It can be seen that the modified model can accurately predict the capacitance changes compared with the general model. For the sensing unit located at row 3, column 6, the predicted sensitivities by the general and modified models are 1.05%/N and 1.25%/N, respectively. The discrepancy between the predicted sensitivities by general model and experimental measurement are 15.3%, which is much larger than the discrepancy (0.8%) between the modified model and experiment. Results demonstrate that the modified model can accurately predict the sensitivity of the sensing unit at row 3, column 6. Also, the discrepancies of the sensitivities between the modified model predicted and experimental measurements are summarized in table 1. The maximum error is 4.88% which occurs at column 4, row 4. For all the 64 sensing units, the average error is only 0.56%. High accuracy demonstrates that the developed modified model can accurately predict the capacitance change and sensitivity of all the sensing units for normal force measurement.

For the tactile sensor array, the adjacent units will inevitably be affected when applying force on a sensing unit. We apply a 5.21 N normal force on the sensing unit U₄₄ (at row 4, column 4) and then measure the capacitance changes of this sensing unit and its adjacent units. Experiment measurement shows that the capacitance of U₄₄ increased from 9.68 to 10.27 pF, about 6.1% capacitance change. The capacitances of its eight adjacent units decrease in the range of 0.015%–0.019%. The result demonstrates that the interactions between adjacent units are small. The gap between two adjacent units is 1000 μm, about 18.2 times larger than the thickness of the upper PET and PDMS film (55 μm in total). Thus, according to the Saint–Venant theory [27], the loading force applied on one sensing unit will have little effect on its adjacent units.

The sensing performance of the sensor array mounted on a curved surface with 30 mm radius of curvature has been studied, as shown in figure 13. The average initial capacitances of the sensing units increase 0.35% compared to those on the flat surface. When the sensor array is bent onto a curved surface, the dielectric layer is compressed, which leads to the capacitance increase. Since the un-patterned PDMS film has been utilized as the dielectric layer, it will not be compressed too much during bending. Thus, the capacitance increase is small (0.35%). Then, we apply 0–10 N normal force on the sensing unit U₄₄ (at row 4, column 4) and the capacitance change is measured as 15.1% when the force is 10 N. The capacitance change calculated by the modified model is 12.5%. The discrepancy between the measured and calculated results is 2.6%. When the sensor has been mounted onto the curved surface, the bump surface will also be bent into a curved shape. The loading bar, which contacts with the curved bump surface, will induce uneven force on the sensing unit and thus generates larger capacitance increase, as described in section 3.3. As the radius of curvature of the curved surface decreases, the contact force on the sensing unit becomes more uneven, thus the accuracy of the modified model will be decreased. For accurate prediction of the capacitance change, the effects of the curved surface need to be considered and incorporated into the developed models, which has become our future work.

Zoom In Zoom Out Reset image size

By defining $K_{\overline{ijk}ijk}^{mn}\ \left(m,n=u,v,w\right)$ as the element in $\overline{ijk}$ th row and ijkth column, the elements of the stiffness sub-matrices in equation (5) are calculated as

$$\begin{eqnarray}&&\left\{ {\begin{array}{*{20}{l}} {K_{\overline{ijk}ijk}^{uu} = 3.5 \times {{10}^{ - 5}}D_{i\bar i}^{11}E_{j\bar j}^{00}F_{k\bar k}^{00} + 3.5 \times {{10}^{ - 5}}D_{i\bar i}^{00}E_{j\bar j}^{11}T_{k\bar k}^{00}}\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 0.114D_{i\bar i}^{00}E_{j\bar j}^{00}T_{k\bar k}^{11}}\\ {K_{\overline{ijk}ijk}^{uv} = 3.5 \times {{10}^{ - 5}}D_{i\bar i}^{10}E_{j\bar j}^{01}T_{k\bar k}^{00} + 3.5 \times {{10}^{ - 5}}D_{i\bar i}^{01}E_{j\bar j}^{10}S_{k\bar k}^{00}}\\ {K_{\overline{ijk}ijk}^{u\omega } = 0.002D_{i\bar i}^{10}E_{j\bar j}^{00}T_{k\bar k}^{01} + 0.002D_{i\bar i}^{01}E_{j\bar j}^{10}S_{k\bar k}^{10}}\\ {K_{\overline{ijk}ijk}^{vu} = 3.5 \times {{10}^{ - 5}}D_{i\bar i}^{01}E_{j\bar j}^{10}T_{k\bar k}^{00} + 3.5 \times {{10}^{ - 5}}D_{i\bar i}^{10}E_{j\bar j}^{01}S_{k\bar k}^{00}}\\ {K_{\overline{ijk}ijk}^{vv} = 3.5 \times {{10}^{ - 5}}D_{i\bar i}^{00}E_{j\bar j}^{11}F_{k\bar k}^{00} + 0.114D_{i\bar i}^{00}E_{j\bar j}^{00}T_{k\bar k}^{11}}\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 3.5 \times {{10}^{ - 5}}D_{i\bar i}^{11}E_{j\bar j}^{00}T_{k\bar k}^{00}}\\ {K_{\overline{ijk}ijk}^{v\omega } = 0.002D_{i\bar i}^{00}E_{j\bar j}^{10}T_{k\bar k}^{01} + 0.002D_{i\bar i}^{00}E_{j\bar j}^{01}S_{k\bar k}^{10}}\\ {K_{\overline{ijk}ijk}^{\omega u} = 0.002D_{i\bar i}^{01}E_{j\bar j}^{00}T_{k\bar k}^{10} + 0.002D_{i\bar i}^{10}E_{j\bar j}^{01}S_{k\bar k}^{01}}\\ {K_{\overline{ijk}ijk}^{\omega v} = 0.002D_{i\bar i}^{00}E_{j\bar j}^{01}T_{k\bar k}^{10} + 0.002D_{i\bar i}^{00}E_{j\bar j}^{10}S_{k\bar k}^{01}}\\ {K_{\overline{ijk}ijk}^{\omega \omega } = 0.114D_{i\bar i}^{00}E_{j\bar j}^{00}F_{k\bar k}^{11} + 3.5 \times {{10}^{ - 5}}D_{i\bar i}^{00}E_{j\bar j}^{11}T_{k\bar k}^{00}}\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 3.5 \times {{10}^{ - 5}}D_{i\bar i}^{11}E_{j\bar j}^{00}T_{k\bar k}^{00}} \end{array}} \right.\end{eqnarray} \tag{ B.1 }$$

in which

$$\begin{eqnarray}&&\left\{ {\begin{array}{*{20}{l}} {D_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{{{\rm{d}}^r}[(1 - {X^2}){P_\sigma }]}}{{{\rm{d}}{X^r}}} \cdot \frac{{{{\rm{d}}^s}[(1 - {X^2}){P_\theta }]}}{{{\rm{d}}{X^s}}}{\rm{d}}X} }\\ {E_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{{{\rm{d}}^r}[(1 - {Y^2}){P_\sigma }]}}{{{\rm{d}}{Y^r}}} \cdot \frac{{{{\rm{d}}^s}[(1 - {Y^2}){P_\theta }]}}{{{\rm{d}}{Y^s}}}{\rm{d}}Y} }\\ {F_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{E(1 - \nu )}}{{(1 + \nu )(1 - 2\nu )}} \cdot \frac{{{{\rm{d}}^r}[(1 + Z){P_\sigma }]}}{{{\rm{d}}{Z^r}}} \cdot \,\frac{{{{\rm{d}}^s}[(1 + Z){P_\theta }]}}{{{\rm{d}}{Z^s}}}{\rm{d}}Z} }\\ {S_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{E\nu }}{{(1 + \nu )(1 - 2\nu )}} \cdot \frac{{{{\rm{d}}^r}[(1 + Z){P_\sigma }]}}{{{\rm{d}}{Z^r}}} \cdot \,\frac{{{{\rm{d}}^s}[(1 + Z){P_\theta }]}}{{{\rm{d}}{Z^s}}}{\rm{d}}Z} }\\ {T_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{E}{{2(1 + \nu )}} \cdot \frac{{{{\rm{d}}^r}[(1 + Z){P_\sigma }]}}{{{\rm{d}}{Z^r}}} \cdot \frac{{{{\rm{d}}^s}[(1 + Z){P_\theta }]}}{{{\rm{d}}{Z^s}}}{\rm{d}}Z} } \end{array}} \right.\end{eqnarray} \tag{ B.2 }$$

where $r,s=0,1;\sigma =i,j,k;\theta =\bar{i},\bar{j},\bar{k}.$

The elements of the loading vector Q in equation (5) can be calculated as

$$\begin{eqnarray}&&{{Q}_{ijk}}=\frac{{{a}^{2}}{{q}_{0}}}{2}\int_{-0.75}^{0.75} \left(1-{{X}^{2}}\right){{P}_{i}}\text{d}X\cdot \int_{-0.75}^{0.75} \left(1-{{Y}^{2}}\right){{P}_{j}}\text{d}Y\cdot {{P}_{k}}(1)\end{eqnarray} \tag{ B.3 }$$

where ${{Q}_{ijk}}$ is the element in kth row of the matrix Q.

The explicit form of the undetermined vector ${{C}^{\theta}}\ \left(\theta =u,v,w\right)$ in equation (5) is:

$$\begin{eqnarray}&&\begin{array}{l} {C^\theta } = \left[ {C_{111}^\theta ,C_{112}^\theta ,\cdot\cdot\cdot,C_{11K}^\theta ,C_{121}^\theta ,C_{122}^\theta ,\cdot\cdot\cdot,} \right.\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{\left. {C_{12K}^\theta ,\cdot\cdot\cdot,C_{1JK}^\theta ,\cdot\cdot\cdot,C_{IJK}^\theta } \right]^T} \end{array}\end{eqnarray} \tag{ B.4 }$$

where $\theta =u,v,w.$

By defining $K_{\overline{ijk}ijk}^{{}}$ as the element in $\overline{ijk}$ th row and ijk th column, the elements of the stiffness matrix in equation (9) can be calculated as

$$\begin{eqnarray}&&\begin{array}{l} {K_{\overline{ijk}ijk}} = \frac{1}{{35}}\int_{ - 1}^1 {{P_i}{P_{\bar i}}{\rm{d}}{X_1}} \cdot \int_{ - 1}^1 {{P_j}{P_{\bar j}}{\rm{d}}{Y_1}} \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \cdot \int_{ - 1}^1 {\left(\lambda + 2G\right)\frac{{\partial \left(f{P_k}\right)}}{{\partial {Z_1}}} \cdot \frac{{\partial \left(f{P_{\bar k}}\right)}}{{\partial {Z_1}}}{\rm{d}}{Z_1}} \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 3.5 \times {10^{ - 5}}\int_{ - 1}^1 {\frac{{{\rm{d}}{P_i}}}{{{\rm{d}}{X_1}}} \cdot \frac{{{\rm{d}}{P_{\bar i}}}}{{{\rm{d}}{X_1}}}{\rm{d}}{X_1}} \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \cdot \int_{ - 1}^1 {{P_j}{P_{\bar j}}{\rm{d}}{Y_1}} \cdot \int_{ - 1}^1 {G{f^2}{P_k}{P_{\bar k}}{\rm{d}}{Z_1}} \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 3.5 \times {10^{ - 5}}\int_{ - 1}^1 {{P_i}{P_{\bar i}}{\rm{d}}{X_1}} \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \cdot \int_{ - 1}^1 {\frac{{{\rm{d}}{P_j}}}{{{\rm{d}}{Y_1}}} \cdot \frac{{{\rm{d}}{P_{\bar j}}}}{{{\rm{d}}{Y_1}}}{\rm{d}}{Y_1}} \cdot \int_{ - 1}^1 {G{f^2}{P_k}{P_{\bar k}}{\rm{d}}{Z_1}} \end{array}\end{eqnarray} \tag{ B.5 }$$

The elements of the loading vector $^{1}Q$ in equation (9) are calculated as

$$\begin{eqnarray}&&{{Q}_{k}}=3.22\times {{10}^{-4}}{{q}_{0}}{{P}_{k}}(1)\int_{-1}^{1} {{P}_{i}}\text{d}{{X}_{1}}\cdot \int_{-1}^{1} {{P}_{j}}\text{d}{{Y}_{1}}\end{eqnarray} \tag{ B.6 }$$

where ${{Q}_{k}}$ is the element in kth row of the matrix $^{1}Q.$

The explicit form of the undetermined vectors in equation (9) is

$$\begin{eqnarray}&&\begin{array}{*{20}{l}} {{}^1{C^w} = \left[ {C_{111}^w,C_{112}^w, \cdot \cdot \cdot ,C_{11K}^w,C_{121}^w,C_{122}^w, \cdot \cdot \cdot ,} \right.}\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\left. {C_{12K}^w, \cdot \cdot \cdot ,C_{1JK}^w, \cdot \cdot \cdot ,C_{IJK}^w} \right]}^T}} \end{array}\end{eqnarray} \tag{ B.7 }$$

By defining $K_{\overline{ijk}ijk}^{mn}\left(m,n=u,w\right)$ as the element in $\overline{ijk}$ th row and ijk th column, the elements of the stiffness matrices in equation (13) can be calculated as

$$\begin{eqnarray}&&\left\{ {\begin{array}{*{20}{l}} {K_{\overline{ijk}ijk}^{uu} = 7 \times {{10}^{ - 5}}\int_{ - 1}^1 {\frac{{{\rm{d}}(1 - {X_2}^2){P_i}}}{{{\rm{d}}{X_2}}} \cdot \frac{{d(1 - {X_2}^2){P_{\bar i}}}}{{{\rm{d}}{X_2}}}{\rm{d}}{X_2}} }\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \cdot \int_{ - 1}^1 {{P_j}{P_{\bar j}}{\rm{d}}{Y_2}} \cdot \int_{ - 1}^1 {(\lambda + 2G){f^2}{P_k}{P_{\bar k}}{\rm{d}}{Z_2}} }\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 1.75 \times {{10}^{ - 5}}\int_{ - 1}^1 {{{(1 - {X_2}^2)}^2}{P_i}{P_{\bar i}}{\rm{d}}{X_2}} \cdot \int_{ - 1}^1 {\frac{{{\rm{d}}{P_j}}}{{{\rm{d}}{Y_2}}} \cdot \frac{{{\rm{d}}{P_{\bar j}}}}{{{\rm{d}}{Y_2}}}{\rm{d}}{Y_2}} }\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \cdot \int_{ - 1}^1 {G{f^2}{P_k}{P_{\bar k}}{\rm{d}}{Z_2}} + \frac{1}{{70}}\int_{ - 1}^1 {{{(1 - {X_2}^2)}^2}{P_i}{P_{\bar i}}{\rm{d}}{X_2}} }\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \cdot \int_{ - 1}^1 {{P_j}{P_{\bar j}}{\rm{d}}{Y_2}} \cdot \int_{ - 1}^1 {G\frac{{{\rm{d}}{P_k}f}}{{{\rm{d}}{Z_2}}} \cdot \frac{{{\rm{d}}P{f_{\bar k}}}}{{{\rm{d}}{Z_2}}}{\rm{d}}{Z_2}} }\\ {K_{\overline{ijk}ijk}^{uw} = 0.001\int_{ - 1}^1 {(1 + {X_2}){P_i} \cdot \frac{{{\rm{d}}(1 - {X_2}^2){P_{\bar i}}}}{{{\rm{d}}{X_2}}}{\rm{d}}{X_2}} }\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \cdot \int_{ - 1}^1 {{P_j}{P_{\bar j}}{\rm{d}}{Y_2}} \cdot \int_{ - 1}^1 {\lambda \frac{{{\rm{d}}{P_k}f}}{{{\rm{d}}{Z_2}}} \cdot f{P_{\bar k}}{\rm{d}}{Z_2}} }\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 0.001\int_{ - 1}^1 {\frac{{{\rm{d}}(1 + {X_2}){P_i}}}{{{\rm{d}}{X_2}}}(1 - {X_2}^2){P_{\bar i}}{\rm{d}}{X_2}} }\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \cdot \int_{ - 1}^1 {{P_j}{P_{\bar j}}{\rm{d}}{Y_2}} \cdot \int_{ - 1}^1 {Gf{P_k} \cdot \frac{{{\rm{d}}f{P_{\bar k}}}}{{{\rm{d}}{Z_2}}}{\rm{d}}{Z_2}} }\\ {K_{\overline{ijk}ijk}^{wu} = K_{ijk\overline{ijk} }^{uw}}\\ {K_{\overline{ijk}ijk}^{ww} = \frac{1}{{70}}\int_{ - 1}^1 {{{(1 + {X_2})}^2}{P_i} \cdot {P_{\bar i}}{\rm{d}}{X_2}} \cdot \int_{ - 1}^1 {{P_j}{P_{\bar j}}{\rm{d}}{Y_2}} }\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \cdot \int_{ - 1}^1 {(\lambda + 2G)\frac{{{\rm{d}}f{P_k}}}{{{\rm{d}}{Z_2}}} \cdot \frac{{{\rm{d}}f{P_{\bar k}}}}{{{\rm{d}}{Z_2}}}{\rm{d}}{Z_2}} }\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 0.00007\int_{ - 1}^1 {\frac{{{\rm{d}}(1 + {X_2}){P_i}}}{{{\rm{d}}{X_2}}} \cdot \frac{{{\rm{d}}(1 + {X_2}){P_{\bar i}}}}{{{\rm{d}}{X_2}}}{\rm{d}}{X_2}} }\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \cdot \int_{ - 1}^1 {{P_j}{P_{\bar j}}{\rm{d}}{Y_2}} \cdot \int_{ - 1}^1 {G{f^2}{P_k}{P_{\bar k}}{\rm{d}}{Z_2}} }\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 1.75 \times {{10}^{ - 5}}\int_{ - 1}^1 {{{(1 + {X_2})}^2}{P_i}{P_{\bar i}}{\rm{d}}{X_2}} }\\ {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \cdot \int_{ - 1}^1 {\frac{{{\rm{d}}{P_j}}}{{{\rm{d}}{Y_2}}} \cdot \frac{{{\rm{d}}{P_{\bar j}}}}{{{\rm{d}}{Y_2}}}{\rm{d}}{Y_2}} \cdot \int_{ - 1}^1 {G{f^2}{P_k}{P_{\bar k}}{\rm{d}}{Z_2}} } \end{array}} \right.\end{eqnarray} \tag{ B.8 }$$

The elements of the loading vector $^{2}Q$ in equation (13) are calculated as

$$\begin{eqnarray}&&{{Q}_{ijk}}=1.61\times {{10}^{-4}}{{q}_{0}}{{P}_{k}}(1)\int_{0}^{1} \left(1+{{X}_{2}}\right){{P}_{i}}\text{d}{{X}_{2}}\cdot \int_{-1}^{1} {{P}_{j}}\text{d}{{Y}_{2}}\end{eqnarray} \tag{ B.9 }$$

where ${{Q}_{ijk}}$ is the ijkth element of the matrix $^{2}Q.$

The explicit form of the undetermined vectors in equation (13) is:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}} {{}^2{C^\theta } = \left[ {C_{111}^\theta ,C_{112}^\theta , \cdot \cdot \cdot ,C_{11K}^\theta ,C_{121}^\theta ,C_{122}^\theta , \cdot \cdot \cdot ,} \right.}\\ {{{\left. {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,C_{12K}^\theta , \cdot \cdot \cdot ,C_{1JK}^\theta , \cdot \cdot \cdot ,C_{IJK}^\theta } \right]}^T}} \end{array}\end{eqnarray} \tag{ B.10 }$$

where $\theta =u,w.$

By defining $K_{\overline{ijk}ijk}^{mn}\left(m,n=u,v,w\right)$ as the element in $\overline{ijk}$ th row and ijkth column, the elements of the stiffness matrices in equation (17) can be calculated as

$$\begin{eqnarray}&&\left\{ \begin{array}{l} K_{\overline{ijk}ijk}^{uu} = 3.5 \times {10^{ - 5}}{ \times ^2}D_{i\bar i}^{11}{ \cdot ^1}E_{j\bar j}^{00} \cdot F_{k\bar k}^{00}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 3.5 \times {10^{ - 5}}{ \times ^2}D_{i\bar i}^{00}{ \cdot ^1}E_{j\bar j}^{11} \cdot T_{k\bar k}^{00}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 7.14 \times {10^{ - 3}}{ \times ^2}D_{i\bar i}^{00}{ \cdot ^1}E_{j\bar j}^{00} \cdot T_{k\bar k}^{11}\\ K_{\overline{ijk}ijk}^{uv} = 3.5 \times {10^{ - 5}}{ \times ^3}D_{i\bar i}^{01}{ \cdot ^4}E_{j\bar j}^{10} \cdot S_{k\bar k}^{00}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 3.5 \times {10^{ - 5}}{ \times ^3}D_{i\bar i}^{10}{ \cdot ^4}E_{j\bar j}^{01} \cdot T_{k\bar k}^{00}\\ K_{\overline{ijk}ijk}^{uw} = 5 \times {10^{ - 4}}{ \times ^3}D_{i\bar i}^{01}{ \cdot ^1}E_{j\bar j}^{00} \cdot S_{k\bar k}^{10}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 5 \times {10^{ - 4}}{ \times ^3}D_{i\bar i}^{10}{ \cdot ^3}E_{j\bar j}^{00} \cdot T_{k\bar k}^{01}\\ K_{\overline{ijk}ijk}^{vu} = K_{ijk\overline{ijk} }^{uv}\\ K_{\overline{ijk}ijk}^{vv} = 3.5 \times {10^{ - 5}}{ \times ^1}D_{i\bar i}^{00}{ \cdot ^2}E_{j\bar j}^{11} \cdot F_{k\bar k}^{00}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 3.5 \times {10^{ - 5}}{ \times ^1}D_{i\bar i}^{11}{ \cdot ^2}E_{j\bar j}^{00} \cdot T_{k\bar k}^{00}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 7.14 \times {10^{ - 3}}{ \times ^1}D_{i\bar i}^{00}{ \cdot ^2}E_{j\bar j}^{00} \cdot T_{k\bar k}^{11}\\ K_{\overline{ijk}ijk}^{vw} = 5 \times {10^{ - 4}}{ \times ^1}D_{i\bar i}^{00}{ \cdot ^3}E_{j\bar j}^{01} \cdot S_{k\bar k}^{10}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 5 \times {10^{ - 4}}{ \times ^1}D_{i\bar i}^{00}{ \cdot ^3}E_{j\bar j}^{10} \cdot T_{k\bar k}^{01}\\ K_{\overline{ijk}ijk}^{wu} = K_{ijk\overline{ijk} }^{uw}\\ K_{\overline{ijk}ijk}^{wv} = K_{ijk\overline{ijk} }^{vw}\\ K_{\overline{ijk}ijk}^{ww} = 7.14 \times {10^{ - 3}}{ \times ^1}D_{i\bar i}^{00}{ \cdot ^1}E_{j\bar j}^{00} \cdot F_{k\bar k}^{11}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 3.5 \times {10^{ - 5}}{ \times ^1}D_{i\bar i}^{11}{ \cdot ^1}E_{j\bar j}^{00} \cdot T_{k\bar k}^{00}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 3.5 \times {10^{ - 5}}{ \times ^1}D_{i\bar i}^{00}{ \cdot ^1}E_{j\bar j}^{11} \cdot T_{k\bar k}^{00} \end{array} \right.\end{eqnarray} \tag{ B.11 }$$

in which

$$\begin{eqnarray}&&\left\{ \begin{array}{l} ^1D_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{{{\rm{d}}^r}[(1 + {X_3}){P_\sigma }]}}{{{\rm{d}}{X_3}^r}} \cdot \frac{{{{\rm{d}}^s}[(1 + {X_3}){P_\theta }]}}{{{\rm{d}}{X_3}^s}}{\rm{d}}{X_3}} \\ ^2D_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{{{\rm{d}}^r}[(1 - {X_3}^2){P_\sigma }]}}{{{\rm{d}}{X_3}^r}} \cdot \frac{{{{\rm{d}}^s}[(1 - {X_3}^2){P_\theta }]}}{{{\rm{d}}{X_3}^s}}{\rm{d}}{X_3}} \\ ^3D_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{{{\rm{d}}^r}[(1 + {X_3}){P_\sigma }]}}{{{\rm{d}}{X_3}^r}} \cdot \frac{{{{\rm{d}}^s}[(1 - {X_3}^2){P_\theta }]}}{{{\rm{d}}{X_3}^s}}{\rm{d}}{X_3}} \\ ^4D_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{{{\rm{d}}^r}[(1 - {X_3}^2){P_\sigma }]}}{{{\rm{d}}{X_3}^r}} \cdot \frac{{{{\rm{d}}^s}[(1 + {X_3}){P_\theta }]}}{{{\rm{d}}{X_3}^s}}{\rm{d}}{X_3}} \\ ^1E_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{{{\rm{d}}^r}[(1 + {Y_3}){P_\sigma }]}}{{d{Y_3}^r}} \cdot \frac{{{{\rm{d}}^s}[(1 + {Y_3}){P_\theta }]}}{{{\rm{d}}{Y_3}^s}}{\rm{d}}{Y_3}} \\ ^2E_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{{{\rm{d}}^r}[(1 - {Y_3}^2){P_\sigma }]}}{{{\rm{d}}{Y_3}^r}} \cdot \frac{{{{\rm{d}}^s}[(1 - {Y_3}^2){P_\theta }]}}{{{\rm{d}}{Y_3}^s}}{\rm{d}}{Y_3}} \\ ^3E_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{{{\rm{d}}^r}[(1 + {Y_3}){P_\sigma }]}}{{{\rm{d}}{Y_3}^r}} \cdot \frac{{{{\rm{d}}^s}[(1 - {Y_3}^2){P_\theta }]}}{{{\rm{d}}{Y_3}^s}}{\rm{d}}{Y_3}} \\ ^4E_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{{{\rm{d}}^r}[(1 - {Y_3}^2){P_\sigma }]}}{{{\rm{d}}{Y_3}^r}} \cdot \frac{{{{\rm{d}}^s}[(1 + {Y_3}){P_\theta }]}}{{{\rm{d}}{Y_3}^s}}{\rm{d}}{Y_3}} \\ F_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{E(1 - \nu )}}{{(1 + \nu )(1 - 2\nu )}} \cdot \frac{{{{\rm{d}}^r}(f{P_\sigma })}}{{{\rm{d}}{Z^r}}} \cdot \,\frac{{{{\rm{d}}^s}(f{P_\theta })}}{{{\rm{d}}{Z^s}}}{\rm{d}}Z} \\ S_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{{E\nu }}{{(1 + \nu )(1 - 2\nu )}} \cdot \frac{{{{\rm{d}}^r}(f{P_\sigma })}}{{{\rm{d}}{Z^r}}} \cdot \,\frac{{{{\rm{d}}^s}(f{P_\theta })}}{{{\rm{d}}{Z^s}}}{\rm{d}}Z} \\ T_{\sigma \theta }^{rs} = \int_{ - 1}^1 {\frac{E}{{2(1 + \nu )}} \cdot \frac{{{{\rm{d}}^r}(f{P_\sigma })}}{{{\rm{d}}{Z^r}}} \cdot \frac{{{{\rm{d}}^s}(f{P_\theta })}}{{{\rm{d}}{Z^s}}}{\rm{d}}Z} \end{array} \right.\end{eqnarray} \tag{ B.12 }$$

where $r,s=0,1;\sigma =i,j,k;\theta =\bar{i},\bar{j},\bar{k}.$

The elements of the loading vector $^{3}Q$ in equation (17) are calculated as

$$\begin{eqnarray}&&{{Q}_{ijk}}=8.05\times {{10}^{-5}}{{P}_{k}}(1)\int_{0}^{1} \int_{0}^{1} {{q}_{0}}\left(1+{{X}_{3}}\right){{P}_{i}}\left(1+{{Y}_{3}}\right){{P}_{j}}\text{d}{{X}_{3}}\text{d}{{Y}_{3}}\end{eqnarray} \tag{ B.13 }$$

where ${{Q}_{ijk}}$ is the ijkth element of the matrix $^{3}Q.$

The explicit form of the undetermined vector in equation (17) is

$$\begin{eqnarray}&&\begin{array}{*{20}{l}} {{}^3{C^\theta } = \left[ {C_{111}^\theta ,C_{112}^\theta , \cdot \cdot \cdot ,C_{11K}^\theta ,C_{121}^\theta ,C_{122}^\theta , \cdot \cdot \cdot ,} \right.}\\ {{{\left. {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,C_{12K}^\theta , \cdot \cdot \cdot ,C_{1JK}^\theta , \cdot \cdot \cdot ,C_{IJK}^\theta } \right]}^T}} \end{array}\end{eqnarray} \tag{ B.14 }$$

where $\theta =u,v,w.$