## 3 Finite Element Expressions  
The difference between the Cosserat and classical FEAs are the strain transformation matrix, \(\mathbf{B}\), stiffness matrix, \(\mathbf{D}\), and degrees of freedom at each node. In the Cos-FEA formulation, each node is associated with three displacements or three rotational degrees of freedom. The vector of the \(i\)-th node degrees of freedom is defined as:  
\[\left[\mathbf{n}^{c}\mathbf{\omega}^{c}\right]_{i}=\left[\imath_{1}^{c},\ \imath_{2}^{c},\ \imath_{3}^{c},\ \omega_{1}^{c},\ \omega_{2}^{c},\ \omega_{3}^{c} \right]_{i}^{T} \tag{5}\]  
Using a notation similar to Voigt notation, the second-order strain and curvature tensors can be expressed in the following vectorial form:  
\[\left\{ \begin{array}{l}\mathbf{\epsilon}^{c}=\left[\varepsilon_{11}^{c} \varepsilon_{22}^{c}\varepsilon_{33}^{c}\varepsilon_{12}^{c}\varepsilon_{21} ^{c}\varepsilon_{32}^{c}\varepsilon_{23}^{c}\varepsilon_{13}^{c}\right]^{T} \\ \mathbf{\kappa}^{c}=\left[\kappa_{13}^{c}\kappa_{31}^{c}\right]^{T}\end{array} \right. \tag{6}\]  
Finally, by employing Finite Element Method (FEM) discretization techniques along with the interpolation function, \(\mathbf{N}\), the strain and curvature fields can be interpolated in relation to the nodal degrees of freedom vectors \(\mathbf{u}_{i}^{c}\) and \(\mathbf{\omega}_{i}^{c}\) as follows:  
\[\left\{ \begin{array}{l}\mathbf{\varepsilon}\\ \mathbf{\kappa}\end{array} \right\}=\mathbf{B}_{i}\left\{\begin{array}{l}\mathbf{u}_{i}^{ c}\\ \mathbf{\omega}_  
with  
\[B_{i1}= \begin{bmatrix}N_{i,1}&0&0\\ 0&N_{i,2}&0\\ 0&0&N_{i,3}\\ 0&N_{i,1}&0\\ N_{i,2}&0&0\\ 0&N_{i,3}&0\\ 0&0&N_{i,2}\\ 0&0&N_{i,1}\end{bmatrix} ,\ \textbf{B}_{i2}= \begin{bmatrix}0&0&0\\ 0&0&0\\ 0&0&0\\ 0&0&-N_{i}\\ 0&0&N_{i}\\ N_{i}&0&0\\ -N_{i}&0&0\\ 0&N_{i}&0\end{bmatrix} ,\ \text{and}\ \textbf{B}_{i3}= \begin{bmatrix}N_{i,3}&0\\ 0&0\\ 0&N_{i,1}\end{bmatrix} ^{T} \tag{9}\]  
where \(\textbf{N}_{i}\) represents the shape function associated with the \(\boldsymbol{i}\)-th node. Additionally, \(\textbf{N}_{i}\) is utilized for interpolating both the displacement and rotation fields. In the context of small deformation analysis, the internal force \(\boldsymbol{F}^{\text{int}}\) is defined as:  
\[\boldsymbol{F}^{\text{int}}=\int_{V}\left(\boldsymbol{B}_{i1}^{T}\boldsymbol{ \sigma}+\boldsymbol{B}_{i2}^{T}\boldsymbol{\sigma}+\boldsymbol{B}_{i2}^{T} \boldsymbol{m}\right)dv \tag{10}\]  
Finally, in the interlayer material stiffness analysis, the stiffness matrix \(\boldsymbol{K}_{ij}\) is defined as:  
\[\boldsymbol{K}_{ij}\boldsymbol{=}\boldsymbol{B}_{i}^{T}\boldsymbol{D} \boldsymbol{B}_{j} \tag{11}\]  
where \(\boldsymbol{D}\) is a block symmetric matrix that establishes the relationship between the stress and couple stress measures and their work-conjugate measures, namely strains and curvatures, through appropriate constitutive laws \(\boldsymbol{D}_{1}\) and \(\boldsymbol{D}_{2}\). In the local coordinate system, indicated by the hat symbol, \(\check{\ }\), the constitutive laws are expressed as:  
\[ \begin{split}\left\{\begin{matrix}\hat{\sigma}\\ \hat{\boldsymbol{m}}\end{matrix}\right\}=\hat{\boldsymbol{D}}\begin{pmatrix} \hat{\boldsymbol{\epsilon}}\\ \hat{\boldsymbol{\kappa}}\end{pmatrix}=\begin{bmatrix}\hat{\boldsymbol{D}}_{ 1}&\boldsymbol{\theta}_{8\times 2}\\ \boldsymbol{\theta}_{2\times 8}&\hat{\boldsymbol{D}}_{2}\end{bmatrix} \begin{bmatrix}\hat{\boldsymbol{\epsilon}}\\ \hat{\boldsymbol{\kappa}}\end{split}  \tag{12}\]  
In three-dimensional analysis, in order to obtain the \(\boldsymbol{D}_{1}\) and \(\boldsymbol{D}_{2}\) in the global coordinate system, a transformation matrix \(\boldsymbol{L}\) should be applied to both \(\hat{\boldsymbol{D}}_{1}\) and \(\hat{\boldsymbol{D}}_{2}\). It is expressed in the following form:  
\[\boldsymbol{D}= \begin{bmatrix}\boldsymbol{D}_{1}&\boldsymbol{\theta}_{8 \times 2}\\ \boldsymbol{\theta}_{2\times 8}&\boldsymbol{D}_{2}\end{bmatrix} = \boldsymbol{L}^{T}\hat{\boldsymbol{D}}\boldsymbol{L}=