A building with timber saves 1 ton of CO² per cubic meter and replaces other materials made of finite resources such as metals and minerals. Additionally, wood is a recyclable material, which means that if managed in the right way all parts of the building are easily demountable. As a result, a building made with timber is energy positive and able to store thousand kg of CO², supporting the European Union to meet the net-zero economy by 2050. The use of wood as a material of construction also reduces the overall quantity of material used during the construction process due to its lightweight as well as decreases building process time and noise pollution on the construction site [Buildup-EU].
For this reasons, accurate modelling of cross laminated timber plates is of high importance. Therefore this project's main objective was also to develop a FE element models of CLT plates by finding formulations with damping models which represent the behaviour of CLT plates appropriately. Moreover, we validate our results by comparing them to an existing ANSYS-model [Schneider et al.].
Introduction
Material behaviour of wood
As it is stated for example in [Winter], strictly taken wood only shows orthotropic behaviour in a cylindrical coordinate system aligned with the natural radial, tangential and longitudinal directions of a trunk (like it is shown on the left part of the figure below). However, as one can see on the right part of the figure, planks are not sawn uniformly from the trunk, the cylindrical (natural) reference frame having always a somewhat arbitrary orientation compared to the plane surfaces of the plank. Therefore, making the transformation of the material matrix from the cylindrical frame to the global Cartesian frame for each individual plank is infeasible. So, to simplify the problem it is assumed that the material is orthotropic in the plank’s Cartesian reference frame with the longitudinal direction of the cylindrical frame matching one of the Cartesian axis.
Orthotropy in cylindrical coordinate system [Winter]
CLT plates
A CLT plate stands for cross laminated timber plate and consists of several layers put together perpendicularly to each other with glue (as shown on the figure on the right). Besides its structural functions, the CLT also creates and separate indoor environments. Additionally, it is durable, light, strong and sustainable. Therefore, the whole building can be constructed using CLT walls and ceilings [Winter].
As mentioned in [Brandner] the CLT shows high mechanical properties. Indeed the side-boards used to make the solid plate when locked in-plane, reduces swelling and shrinkage. Moreover, the benefit of the alternating layer disposition is not only its resistance against exceptional loads (e.g. earthquake) but also the capability for pre-fabrication, which leads to clean and fast construction. Therefore, the low mass, stiffness and bearing capacity against inplane and out-of-plane stresses make it perfect for multi-storey buildings.
Furthermore, as stated by the company Sterling [Sterling Solutions], CLT plates were successfully deployed for creating temporary access roads and working platforms. The company managed to replace traditional hardwood bolted mats made from trees that need 80-100 years to regrow with CLT plates made of softwood (southern yellow pine), which only need 25-30 years to regrow, therefore making their business more sustainable without a decrease in the quality of their product.
Assembly of a CLT plate [Sterling Solutions]
Previous work
The project heavily relied on the work done on CLT plates at the Chair of Structural Mechanics and the Engineering Risk Analysis Group at TUM. Especially the paper "Bayesian system identification in linear structural dynamics with frequency transformed data using rational surrogate models" [Schneider et al.] served as a starting point for our work. In this paper both a measurement done at the University of Applied Sciences in Rosenheim was investigated and a finite element analysis of the CLT plate was condacted using Ansys shell and solid elements (SHELL281 and SOLID185). During the measurement the CLT plate - depicted on the figure on the right below - was hanging freely from two cables attached at the edge x = 0. The excitation of the system was done via an impedance hammer at the marked point ( x = 0.33 m and y = 0.475 m ). The obtained eigenfrequencies and modal damping ratios are given in the table in the lower left corner. During the project, this modal damping ratio data was used to determine the parameters (\alpha-values) for the Rayleigh and Caughey damping models. During the measurement the frequency response of the system was also recorded at multiple points of the CLT plate, noted on the geometry. These frequency responses then were compared with the solutions obtained using Ansys. The comparison of the results is shown in the upper left graphs of the figure below.
| Mode 1 | Mode 2 | Mode 3 | Mode 4 | Mode 5 | Mode 6 | |
|---|---|---|---|---|---|---|
| Eigenfrequency [Hz] | 27.9 | 64.4 | 82.4 | 87.6 | 104.2 | 153.8 |
| Modal damping ratio [-] | 0.067 | 0.02 | 0.015 | 0.009 | 0.031 | 0.014 |
The results and setup of the experiment [Schneider et al.]
Our task was to develope a plate (or shell) and a solid finite element in MATLAB and compare the results with the one obtained with the above mentioned Ansys elements validated with the mentioned measurement. Since the comparison of the measurement data and the Ansys elements was done below the 200 Hz frequency line we later also restricted our comparisons with the Ansys elements to this particular domain.
The material properties used in our finite element analysis of CLT plates was also taken from the study mentioned above and is collected in the table below.
| Material data for the Ansys models | ||
|---|---|---|
| Parameter | Shell model | Solid model |
$E_x \hspace{2pt} [Pa]$ | $1.061 \times 10^{10}$ | $1.1 \times 10^{10}$ |
$E_y \hspace{2pt} [Pa]$ | $6.5 \times 10^8$ | $3.1 \times 10^8$ |
$G_{xy} \hspace{2pt} [Pa]$ | $4.83 \times 10^8$ | $4.83 \times 10^8$ |
$G_{xz} \hspace{2pt} [Pa]$ | $1.725 \times 10^8$ | $6.9 \times 10^8$ |
$G_{yz} \hspace{2pt} [Pa]$ | $9.857 \times 10^7$ | $6.9 \times 10^7$ |
$\zeta \hspace{2pt} [-]$ | $2 \times 10^{-2}$ | $2 \times 10^{-2}$ |
The material data used in the Ansys model [Schneider et al.]
On the modelling approach in general
Outline
In this section a general overview of the process followed for modelling CLT plates will be given, which is the application of the general finite element methodology for this specific case.
When trying to capture a physical phenomenon with a mathematical model one must realize that accounting for every single aspect of the process at hand is impossible, therefore the engineer must decide about the aspects of the phenomenon most relevant for the intended use of the object and try to formalize only these in terms of equations, neglecting everything else in the process. The application of this information filtering is a crutial step in the modelling process and a lot of errors can be made therefore one most reason about every simplification that was made properly.
At the end of this first modelling step one usually ends up with partial differential equations ( in our case the Kirchhoff, Reissner-Mindlin and general Continuum equations ). These are in general impossible to solve analytically therefore numerical methods are needed. For these a proper discretization of the system has to be done. One way of doing this is to use the finite element method. However this discretization of course introduces additional errors into the solution therefore again care must be taken to properly balance between the accuracy of approximation and the computational cost for solving the discrete model.
Finally the system of equation obtained via the discretization has to be solved by numerical methods. These also introduce some additional numerical error to the finally obtained solution. In the end from the solution of the discrete model predictions for the physical phenomenon can be made, but one must always keep in mind the simplifications made during the modelling process.
The modelling process
Mass and stiffness
For deriving the mass and stiffness matrices of the discretized system we reffered to the lecture notes of Proffesor Müller [Müller]. First it is assumed that no damping is present in the system. In this case if the system is in dynamic equilibrium the principle of virtual work for a finite element (with domain \Omega and boundary \Gamma) can be stated as
| \delta W = - \int_{(\Omega)} \delta \mathbf{\epsilon}^T \mathbf{\sigma} d\Omega -\int_{(\Omega)} \delta \mathbf{u}^T \rho \ddot{\mathbf{u}} d\Omega + \int_{(\Omega)} \delta \mathbf{u}^T \mathbf{p} d\Omega + \int_{(\Gamma)} \delta \mathbf{u}^T \mathbf{t} d\Gamma = 0 |
By choosing a set of appropriate shapefunctions \mathbf{\Phi} (containing the shapefunctions columnwise) and time dependent vector of coefficients \mathbf{\alpha} , the quantities in the above expression for the total virtual work can be expressed in the following manner
| \mathbf{u} = \mathbf{\Phi} \mathbf{\alpha}, \\ \mathbf{\varepsilon} = \mathbf{B} \mathbf{\alpha}, \\ \mathbf{\sigma} = \mathbf{E} \mathbf{B} \mathbf{\alpha}, \\ \ddot{\mathbf{u}} = \mathbf{\Phi} \ddot{\mathbf{\alpha}}, |
where \mathbf{E} is the elasticity matrix containing the material parameters and \mathbf{B} is the matrix containing the derivatives of the shapefunctions. With this shapefunction-approximation the principal of virtual work now can be stated in the form below
| \delta W = - \delta \mathbf{\alpha}^T \int_{(\Omega)} \mathbf{B}^T \mathbf{E} \mathbf{B} d\Omega - \mathbf{\alpha}^T \int_{(\Omega)} \rho \mathbf{\Phi}^T \mathbf{\Phi} d\Omega \ddot{\mathbf{\alpha}} + \mathbf{\alpha}^T \int_{(\Omega)} \delta \mathbf{\Phi}^T \mathbf{p} d\Omega + \mathbf{\alpha}^T \int_{(\Gamma)} \mathbf{\Phi}^T \mathbf{t} d\Gamma = 0 |
The mass and stiffness matrices from the above equations are the following expressions:
| \mathbf{M} = \int_{(\Omega)} \rho \mathbf{\Phi}^T \mathbf{\Phi} d\Omega \\ \mathbf{K} = \int_{(\Omega)} \mathbf{B}^T \mathbf{E} \mathbf{B} d\Omega |
For evaluating these expressions, shapefunctions appropriate for the mathematical model in question need to be chosen and the integrals to be evaluated.
Damping
The expressions and figures in this section are taken from the work of Humar [Humar].
Damping is frequently related to the dissipation of energy during mechanical vibrations. In other words, it is connected with the irreversible transformation of mechanical energy into other types of energy, for instance, thermal energy . For our project, we chose two proportional damping models for our CLT plate system: the Rayleigh damping and Caughey damping.
Therefore, we began by the analyse of free vibration. If the damping distribution of our systems is denoted as a matrix [C], the equation of motion will be described by:
| [M] (\ddot{x}) + [C] (\dot{x}) + [K] (x) = (0) |
In the above equation, the matrix [C] cannot be diagonalized. Hence, the damping distribution that allows the diagonalization of matrix [C], as well as mass matrix [M] and stiffness matrix [K], is the proportional damping. This was one of the reasons we chose these Rayleigh damping and Caughey damping models for our project. The other reason is that this type of damping does not require concentrated damper elements and thus, it is a valid approach.
The Rayleigh damping indicates that the damping forces are proportional to mass and stiffness matrices, it can be written as:
| C = α_0M + α_1K |
Here, the α and β are two free parameters. Therefore, we can specify the damping ratio for any two modes (ξi and ξj). Then we obtain:
| Φ^T_i C ϕ_i = 2 ζ_iω_i = α_0 + α_1ω^2_i \\ Φ^T_j C ϕ_j = 2 ζ_jω_j = α_0 + α_1ω^2_j |
In matrix notation:
| \frac{1}{2}\begin{bmatrix} \frac{1}{ω_i}\ & ω_i \\ \frac{1}{ω_j}\ & ω_j \\ \end{bmatrix} \begin{bmatrix} \ α_0 \\ \ α_1 \\ \end{bmatrix} = \begin{bmatrix} \ ξ_i \\ \ ξ_j \\ \end{bmatrix} |
Damping ratio ξ now varies with the frequency ω according to the curve in the Figure on the right.
The variation of the damping ratio over the frequency domain [Humar]
The alpha values for the Rayleigh damping have been calculated from the first and the last modal damping ratios showed in the section “Previous work”. The objective was to try to keep more or less the damping constant over the range of the modal damping ratios.
However, when more than two damping ratios need to be determined, the damping matrix can be given by the Caughey damping:
| C = α_0M + α_1K + α_2KM^{-1}K |
Now, using the orthogonality condition we have:
| Φ^T_i C ϕ_i = 2 ζ_iω_i = α_0 + α_1ω^2_i + α_2ω^4_i \\ Φ^T_j C ϕ_j = 2 ζ_jω_j = α_0 + α_1ω^2_j + α_2ω^4_j \\ Φ^T_k C ϕ_k = 2 ζ_kω_k = α_0 + α_1ω^2_k + α_2ω^4_k |
In matrix form:
| \frac{1}{2}\begin{bmatrix} \frac{1}{ω_i}\ & ω_i & ω^3_i \\ \frac{1}{ω_j}\ & ω_j & ω^3_j \\ \frac{1}{ω_k}\ & ω_k & ω^3_k \\ \end{bmatrix} \begin{bmatrix} \ α_0 \\ \ α_1 \\ \ α_2 \\ \end{bmatrix} = \begin{bmatrix} \ ξ_i \\ \ ξ_j \\ \ ξ_k \\ \end{bmatrix} |
It is important to highlight that the coefficient matrix for the unknowns α is ill-conditionate since the terms in this matrix diverge by several orders of magnitude. Hence, only a few damping ratios can be considered (for our case 3)
In the graph on the right it is possible to see that the Caughey damping frequency is heavily dependent on the choice of the damping ratio values. With 3 damping ratios, the results were reliable but with 4 damping ratios, the numerical calculation failed as the green frequency line on the graph. This is because of the ill-conditioned matrix or some calculation error.
FRFs with Caughey damping coefficients determined from different modal damping ratios
The Plate and Shell elements
The element formulation
Plate in bending element
In our first approach for modelling CLT plates we have implemented and examined two 2 dimensional elements. First a discrete Kirchhoff quadrilateral plate in bending action element was created within the “bm_mfem” MATLAB framework. This element formulation is based on the work of Barrales [Barrales] who builds on the formulation proposed by Batoz and Tahar [Batoz et al.].
This formulation starts from an 8-node serendipity element based on the Reissner-Mindlin theory and enforces the Kirchhoff assumptions at the nodes. The base-element has 8 \times 3 = 24 degrees of freedom. For these then the Kirchhoff assumptions are prescribed, in each coordinate direction at the corner nodes (4 \times 2 = 8 equations) and among the boundary direction at the edge nodes (4 \times 1 = 4 equations). Therefore, adding these 12 constraining equations to the initial 24 degrees of freedom Reissner-Mindlin serendipity element formulation, the final finite element derived becomes a 12 degrees of freedom quadrilateral. These are the transverse translation and two rotational degrees of freedoms at each node as it is depicted on the figure on the right.
This element formulation is very advantageous, because due to the Reissner-Mindlin formulation, the shapefunctions used can be of lower order, but at the same time it does not show the transverse shear locking phenomena observable for thin plate geometries in a classical Reissner-Mindlin finite element formulation. This was tested also during the project and the element was found indeed locking free.
Plate in bending element [Barrales]
Furthermore, in our implementation the element was created as geometrically linear. The material model aiming to represent the behaviour of the CLT plate was chosen as being orthotropic with the material properties taken from the experiment mentioned at the beginning of this project report.
Shell element
The previously discussed plate in bending element was extended later to a shell by adding a membrane formulation to it. This additional model is also based on the work done by Barrales [Barrales] and the element is created as a quadrilateral with two translational and one drilling degree of freedom. The element geometry and degrees of freedom can be seen on the figure on the left next to the previously described plate in bending action model.
For the implementation of the element in our case the same considerations were taken into account as in the case of the plate element. The formulation is geometrically linear, and the material behaviour was chosen to be orthotropic in the Cartesian frame of the element.
With this formulation now the degrees of freedom were increased from 3 to 6 per node therefore substantially increasing the computational cost of determining the solution for the CLT plate vibration problem.
Comparison of the two elements
For comparison of the behaviour of the elements two finite element models were set up within the “bm_mfem” MATLAB framework. One using the plate in bending action element and the other the shell formulation. The models resembled the measurement setup with free boundary conditions at all the edges of the geometry.
With solving the equation systems obtained through the two finite element discretizations the modeshapes and the eigenfrequencies (undamped) were determined for each model. It was observed that the order of the modes is the same in the two cases and the shapes coincide as well.
To be able to assess the differences observed when comparing the eigenfrequencies of the models these were plotted up to the 25th for both cases. The figure can be seen on the right. The blue and orange curves show the eigenfrequencies for the plate and the shell case respectively (with the values given on the left vertical axis) while the black graph shows the calculated difference of these values referring to the right vertical axis.
One can observe that the results up to the 8th mode are almost identical (there is less than 1% difference between the two solutions) and the obtained values only start to diverge after the 13th eigenfrequency. Above this, the plate in bending action solution behaves much stiffer. One explanation for this can be that – unlike the shell element – it is not capable to represent in-plane displacements which for more complicated, highly curved modeshapes can become significant. This could cause a stiffer behaviour compared to the shell formulation.
Since we intend to compare our results with the Ansys solutions only in the same range where its performance was compared with the measurement results, eigenfrequencies above 200 Hz are out of interest. The 8th eigenfrequency is already 204 Hz, therefore no higher modes are needed for the comparison. Concerning that within this range the difference between the plate and shell models is less than 1 %, examining only the plate solutions is sufficient, which significantly decreases the computational cost of the numeric study.
Comparison of the eigenfrequencies of the solutions with the two different elements
Validation
Hinged beam validation case
Before comparing the performance of the newly implemented MATLAB element with the more sophisticated Ansys formulations used in the study [Schneider et al.] mentioned at the beginning of this report, it had to be checked whether the MATLAB code can deliver accurate solutions for a few basic analytically solvable cases.
The first of these cases was a Euler-Bernoulli beam hinged at both of its ends (see upper right geometry of the figure on the right). The analytical expression for the circular eigenfrequencies is derived e.g., in [Müller] and it looks like the following
| \omega_j = \left(\frac{j \pi}{l} \right) \sqrt{\frac{E I}{\mu}} |
The MATLAB model was created with the newly implemented elements as it is depicted in the lower right corner of the figure on the right. The graph on the left part of the right figure shows the calculated difference between the numerical results for the beam eigenfrequencies and the analytical ones obtained using the above stated formula. One can see that even at the 9th eigenfrequency the error stays below 3.5 %. Here it must be also stated that already the 4th beam eigenfrequency is out of the earlier stated frequency range of interest.
The hinged Euler-Bernuolli beam validation case [Müller]
Eigenfrequencies of rectangular plates [Müller]
Validation via solutions for rectangular plates with various boundary conditions
As a second analytical validation case, a rectangular plate solution was taken with two different boundary condition cases. For these two cases analytical solutions for the first two eigenfrequencies are available given that the material model is isotropic. These expressions are given in the table on the left and are taken from the lecture notes of Professor Müller [Müller].
The results of the analytical calculations were compared with the corresponding MATLAB finite element results. The obtained data is given in the table below. As it can be seen, the numerical results obtained for both eigenfrequencies in both cases are quite accurate.
| CSCS Plate | CCCC Plate | |||
|---|---|---|---|---|
f_1 \hspace{2pt} [Hz] | f_2 \hspace{2pt} [Hz] | f_1 \hspace{2pt} [Hz] | f_2 \hspace{2pt} [Hz] | |
| Analytical | 354.97 | 413.21 | 363.56 | 444.75 |
| MATLAB | 358.08 | 421.18 | 365.82 | 450.69 |
| Error | 1% | 2% | 1% | 1% |
Comparison with the Ansys shell model
Finally, the performance of the implemented plate in bending action MATLAB element was compared with the Ansys shell formulation used in the study [Schneider et al.]. Compared to the previous section with the analytical plate solutions the material model here now was already taken as orthotropic. First the effect of different boundary conditions was examined. After this again the measurement setup was recreated in the model space with free boundaries at all edges and the differences between the MATLAB and Ansys models were analysed in terms of eigenfrequencies and frequency response functions.
To study the behaviour of the two models with respect to changing boundary conditions three different cases were set up. A plate with all edges free, another with all hinged and finally a fully clamped example was investigated. The Ansys and MATLAB modeshapes determined are in the same order in each of the cases and match up qualitatively. The results of this analysis in terms of eigenfrequencies of the plate are given on the leftmost graph below. The blue curves represent the MATLAB solutions with the different boundary conditions, while the ones belonging to the Ansys results are depicted in orange. It was observed that the more the degrees of freedom at the boundaries were restricted (first hinged than clamped) the stiffer the MATLAB plate behaved compared to the Ansys solution. While in the “all free” case the difference at the upper bound of the frequency range of interest is “only” 16% in the “all clamped” case it becomes multiples of that. The MATLAB code was examined several times, but we could not find any bugs related to the boundary conditions, but the observed phenomena suggests that the element implementation should be handled with a great amount of scepticism.
Restricting the comparison to the case with free boundary conditions at all edges, some further observations to the eigenfrequencies of the system can be made. For this a more detailed figure was made which can be seen in the middle below. The blue and orange curves depict the eigenfrequencies of the MATLAB and Ansys models respectively and are referring to the left vertical axis, while the error rate plotted in black refers to the right vertical axis. As stated earlier the error rate at the upper bound of the range of interest is around 16%. But looking at the error rate curve (in black) it is apparent that the difference between the solutions does not increase in a gradual manner. According to our observations more complicated modeshapes tend to come with a higher rate of error. This phenomenon might be smoothened out by applying selective mesh refinement in case of the more complex modes.
Comparison of the MATLAB and Ansys plate eigenfrequencies with different boundary conditions
Comparison of the MATLAB and Ansys eigenfrequencies with all boundaries free
MATLAB and Ansys frequency response functions with different damping models
Moving on to the comparison of the frequency response functions (only the real part addressed here for brevity) at the point of excitation, first the behaviour of the Ansys solution with hysteretic (used in the study [Schneider et al.]) and Rayleigh damping was examined. Together with the frequency response of the MATLAB model these functions are plotted on the graph in the rightmost figure above. Looking at the curves of the two Ansys results one can see that except for the last response observable the Rayleigh response is damped more than the hysteretic. In case of the response at the first eigenfrequency the Rayleigh response got almost damped out entirely. As for the MATLAB response it obvious from the figure that it significantly differs from the other two solutions. This will be addressed, however, based on the leftmost figure below.
On the mentioned graph the Ansys response solution with Rayleigh damping is depicted alongside the MATLAB responses obtained with Rayleigh and Caughey damping, furthermore with modal superposition of the first 6 modes using all the measured 6 damping ratios. This last approach is determined based on the expression below
| h\left(\Omega\right) = \sum_{i=1}^{n}\frac{\mathbf{\Phi}_i\mathbf{\Phi}_i^T}{\left(\omega_i^2-\Omega^2\right) + 2 i \zeta_i \omega_i \Omega}, |
where \omega_i are the circular eigenfrequencies belonging to the modes \mathbf{\Phi}_i, \zeta_i are the modal damping ratios and \Omega is the excitation circular frequency.
Returning to the analysis of the above-mentioned figure, one can state that the MATLAB results show the stiffer behaviour (higher eigenfrequencies) on this graph too. Furthermore, the damping of the different responses differs from the Ansys solution as well. The already small response of the Ansys solution gets damped out completely in case of the MATLAB results. This insignificant (or inexistent) response at the first eigenfrequency is due to the fact that the excitation point where the response was determined is close to the “neutral” (not moving) line of the first eigenshape. (The excitation point is marked below with the red dot on the graph in the middle and the right.) On the contrary the response of the MATLAB solutions at the second eigenfrequency is way larger than in case of the Ansys result. Furthermore, looking at the responses at the highest depicted mode one can see that for the MATLAB solution there are two responses very close together while the Ansys solution only shows one. As a final observation it can be noted that the MATLAB response using Caughey damping is the closest to the modal superposition approach which uses all the modal damping ratios obtained from the measurement.
As a conclusion of this section, it is to be stated that the implemented plate in bending action MATLAB element show a behaviour inconsistent with the Ansys element that showed acceptably close results to the measurement data in the study [Schneider et al.]. Therefore, the implemented element can be termed unreliable and further thorough examination of the code and the validity of the underlaying assumption for this case would be needed.
Comparing different MATLAB FRFs with the Ansys response using Rayleigh damping
The first modeshape of the MATLAB plate with free boundary conditions prescribed at all edges
The solid element
The element formulation
In this project, we decided to go for two approaches to model the wooden cross-ply laminate using solid elements. The reason for this was to investigate which of the two methods more accurately represents the dynamic behavior of the structure and whether any of the two methods has any computational advantage over the other.
The Homogenized element approach
The first method was using a homogenized solid element that represents the three layers of the cross-ply laminate. The complete model of the CLT made up from these elements can be seen on the right.
Each one of these layers is modeled as an orthotropic material and thus has an Orthotropic Matrix as shown below:
[C]=\frac{1}{D}\begin{bmatrix} Ex(1-\frac{Ez}{Ey}\nu yz^{2}) & Ey(\nu xy+\frac{Ez}{Ey}\nu xz\nu yz) & Ez(\nu xz+\nu xy\nu yz) & 0 & 0 & 0\\ Ey(\nu xy+\frac{Ez}{Ey}\nu xz\nu yz) & Ey(1-\frac{Ez}{Ex}\nu xz^{2}) & Ez(\nu yz+\frac{Ey}{Ex}\nu xy\nu xz) & 0 & 0 & 0\\ Ez(\nu xz+\nu xy\nu yz) & Ez(\nu yz+\frac{Ey}{Ex}\nu xy\nu xz) & Ez(1-\frac{Ey}{Ex}\nu xy^{2}) & 0 & 0 & 0\\ 0 & 0 & 0 & Gyz & 0 & 0\\ 0 & 0 & 0 & 0 & Gxz & 0\\ 0 & 0 & 0 & 0 & 0 & Gxy \end{bmatrix}
with D = \frac{ExEyEz-\nu yz^{2}ExEz^{2}-\nu xy^{2}Ey^{2}Ez-2\nu xy\nu xz\nu yzEyEz^{2}-\nu xz^{2}EyEz^{2}}{ExEyEz}
CLT FE model with the homogenized elements
The top and the bottom layer are aligned with the global axis; however, the second layer required a transformation of coordinates as the ply is rotated by 90 degrees with respect to the global coordinates. Another way to do this would be to assign the opposite value of Youngs Moduli, Shear Moduli, as well as Poisson’s ratio for the second layer with respect to the longitudinal and transverse direction. In this project, the elasticity tensor of the second layer is transformed using the stress and strain transformation matrices. In order to carry out this transformation, two coordinate systems are considered: The local and the global coordinates. We will refer to the local coordinate as the unprimed coordinate (p,q,r), and the global as the primed coordinate (p’,q’,r’). The relationship between these two coordinate systems can be expressed using the 9 direction cosines (r11, r21, r31), (r12, r22, r32), (r13, r23, r33). The relations between these two coordinate systems can be also expressed in terms of three angles Qp, Qq, Qr. These angles express the rotations about the x, y, and z-axis.
r11 = cos (\varphi pp)
r21 = cos (\varphi qp)
r31 = cos (\varphi rp)
The angles are related to the cosines through the following equation:
Qq=Atan2(-r31,\sqrt{r11^{2} +r21^{2}} ))
When Qq=90
then Qr = 0 and Qp = Atan2(r12, r22)
Then the stress transformation matrix is calculated as the following:
cp=cos(Qp ) , sp=sin(Qp )
cq=cos(Qq),sq=sin(Qq )
cr=cos(Qr),sr=sin(Qr )
[T\sigma p] = \begin{bmatrix} 1& 0& 0& 0& 0&0 \\ 0& cp^{2}& sp^{2} &2cpsp & 0 & 0\\ 0& sp^{2}& cp^{2} & -2cpsp & 0 & 0\\ 0& -cpsp &cpsp & cp^{2}-sp^{2} & 0 & 0\\ 0& 0& 0 & 0 & cp &-sp \\ 0&0 & 0 & 0 & sp & cp \end{bmatrix}[T\sigma q] = \begin{bmatrix} cq^{2}& 0& sq^{2}& 0& 2cqsq&0 \\ 0& 1& 0 &0 & 0 & 0\\ sq^{2}&0& cq^{2} & 0 & -2cqsq& 0\\ 0& -0 &0 & cq & 0 & -sq\\ -cqsq& 0& cqsq & 0 & cq^{2}-sq^{2} &0 \\ 0&0 & 0 & sq & 0 & cq \end{bmatrix}[T\sigma r]=\begin{bmatrix} cr^{2} & sr^{2} & 0 & 0 & 0 & 2crsr\\ sr^{2} & cr^{2} & 0 & 0 & 0 & -2crsr\\ 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & cr & -sr & 0\\ 0 & 0 & 0 & sr & cr & 0\\ -crsr & crsr & 0 & 0 & 0 & cr^{2}-sr^{2} \end{bmatrix}
[Tσ]=[Tσp][Tσq][Tσr]
The strain transformation Matrix will have the exact same result as the stress transformation matrix due to the rotation by 90 degrees around the r axis (Our case the z-axis).
Taking the stress-strain relation and multiplying both sides by the stress transformation matrix yields the following:
[Tσ][σ]=[Tσ][C][ε]
This can be extended to:
[T\sigma ][\sigma ]=[T\sigma ][C] [T\varepsilon ]^{-1} [T\varepsilon ][\varepsilon ]
[\sigma ]'=[T\sigma ][C] [T\varepsilon ]^{-1} [\varepsilon ]'
And thus
[C]'=[T\sigma ][C] [T\varepsilon ]^{-1}
After All the Elasticity tensors of the three layers are aligned with the global coordinates, the corresponding elements of the three elasticity tensors are added together. Each contributes by the ratio of the height of the layer to the total height of the three layers. As all three layers have the same height, this ratio is 1/3 for all layers.
[Ce] = (1/3)[C1] + (1/3)[C2]'+(1/3)[C3]
The three-layered solid elements approach
The second method used was modeling the three layers with three separate solid elements. Each with an elasticity tensor with the proper orientation. As in the previous method, only the second layer requires a transformation of the elasticity tensor. This method definitely required more elements to model the same structure.
X Y Z
The element type and allowed degrees of freedom
The element we chose was a hexahedron 8 noded element with only translational degrees of freedom and no rotational degrees of freedom.
The shape functions used to formulate the stiffness matrix were tri-linear shape functions:
N1 = (1/8)(1 − ξ)(1 − η)(1 + μ)\\ N2 = (1/8)(1 − ξ)(1 − η)(1 − μ)\\ N3 = (1/8)(1 − ξ)(1 + η)(1 − μ)\\ N4 = (1/8)(1 − ξ)(1 + η)(1 + μ)\\ N5 = (1/8)(1 + ξ)(1 − η)(1 + μ)\\ N6 = (1/8)(1 + ξ)(1 − η)(1 − μ)\\ N7 = (1/8)(1 + ξ)(1 + η)(1 − μ)\\ N8 = (1/8)(1 + ξ)(1 + η)(1 + μ)
These shape functions were used to calculate the B matrix as well as the Jacobian and correspondingly the determinant of the jacobian. The global elasticity matrix of the element is also required to compute the stiffness matrix. Two Guass points were used to assure full integration during the calculation of the stiffness matrix. The outcome is a 24x24 stiffness matrix for each solid element.
Ke = \int [B]^{T}[C][B]dV
The mass matrix was calculated using the same shape functions used for the stiffness matrix. The determinant of the Jacobian and the density of the material are required in order to complete the integration and calculate the mass matrix.
Me = \int \rho dV
Provided that dV = detJdξdηdμ
X Y Z
Comparing with well known analytical solutions
Simply supported Euler Bernoulli Beam
To test the element's application limits, it was compared with the simply supported Euler Bernoulli beam in terms of mode shapes and eigenfrequencies. Although the limits of application of an Euler Bernoulli beam do not apply to our structure as it is too wide to be modeled as a Bernoulli beam, this comparison was done to test the limitations of using this kind of element. Because the eigenfrequencies were very different for both the model and the beam, eigenfrequencies with the same mode shapes were compared.
In order to ensure that the boundary conditions were as close as possible to a simply supported beam and to allow a rotation at the fixed supports, Only the bottom ends of the plate were fixed. As the top nodes are allowed to move freely, this simulates a rotation about the fixed nodes.
The Eigenfrequencies of a simply supported beam depend on its length, Youngs Modulus, weight per unit length, and moment of Inertia.
\omega j = (\frac{j\pi }{l})^{2}\sqrt{\frac{EI}{\mu }}
The eigenfrequencies of our model were found by solving the eigenvalue problem using the global mass and stiffness matrices. The rows and columns of the Mass and Stiffness matrices corresponding to the fixed degrees of freedom are removed before solving the eigenvalue problem. The eigenvalues are the eigenfrequencies squared.
det(-\omega^{2}[M]-[K])=0
The mode shapes were found by assuming a value for one of the entries of the eigenvectors and finding the ratio between the displacements for each eigenvalue
| Modes with the same shapes (M,B) | 1,1 | 4,2 | 10,3 | 15,4 |
|---|---|---|---|---|
| Model (Hz) | 236.88 | 652.7 | 1547.54 | 2478.34 |
| EB beam (Hz) | 128.1 | 512.38 | 1152.86 | 2049.52 |
| Error (%) | 84.24 | 27.40 | 36.58 | 20.92 |
As expected, the values were not close to each other due to the limitations of Euler Bernoulli beam theory, the lack of rotational degrees of freedom, as well as coarse FE mesh.
Simply Supported Kirchhoff Plate
The model was also compared with a simply supported Kirchhoff plate in terms of the mode shapes and eigenfrequencies. In order to allow rotation around the supports, only the bottom nodes were fixed.
Mode | 1 | 2 |
|---|---|---|
| Model (Hz) | 1355.64 | 1722.25 |
| Hinged (Hz) | 788.95 | 1172.85 |
| error (%) | 71.83 | 46.84 |
The Eigenfrequencies of a Kirchhoff plate depend on its length, ratio of length to width, depth, Youngs Modulus, Poisson's ratio, and weight per unit length.
\omega j = \beta i\frac{2\pi }{a^{2}}\sqrt{\frac{Ed^{3}}{12(1-\nu ^{2})\mu }}\\ with\, \gamma = a/b\\ \beta \, depends \, on \, \gamma \,and\,on\,the\,boundary\,conditions
The model failed to represent rotation accurately due to a lack of rotational degrees of freedom and a coarse FE mesh, therefore, the error is quite significant for the first couple of mode shapes.
Clamped Kirchhoff Plate
In order to test the model with clamped boundary conditions, it was compared with a clamped Kirchhoff plate. To properly apply the boundary conditions, the plate must be able to neither move nor rotate at the supports. Therefore, both top and bottom edges were fixed.
Mode | 1 | 2 |
|---|---|---|
| Model (Hz) | 2004.42 | 2364.68 |
| Clamped (Hz) | 1610.768 | 1970.521 |
| error (%) | 24.4388 | 20.00278 |
This case showed the closest eigenfrequencies as the boundary conditions more accurately match and the Kirchhoff plate theory applications.
Comparison with the Ansys solid model
In order to test our model, the eigenfrequencies, as well as the frequency response at the load point, were compared with the Ansys model which uses the experimentally verified element.
Eigenfrequencies
As shown in the table above, the unsupported case results were always closer to the Ansys results for both the single layer as well as the three-layer model. Furthermore, the three-layered model results more accurately represented the cross-ply wooden laminate as the eigenfrequencies were closer to the Ansys eigenfrequencies than those of the single-layered homogenized model. It can also be observed that the more fine the FE mesh was, the closer the results were to those obtained from Ansys. Both models appear stiffer than the Ansys model because they use less sophisticated elements with fewer nodes and lower order shape functions. Therefore, the finer the mesh, the less stiff the models behaved. In order to further test this observation, another run was performed on the clamped three-layered model with a very fine mesh.
The fine mesh gave the most accurate results that behaved the least stiff and were closest to the Ansys results. However, it took the longest time to run and was extremely costly computationally. This shows that the choice of the tri-linear shape functions was not adequate enough and that due to the large difference in the dimensions of the solid element, the results were directionally biased and could not accurately capture rotation. Because of the shape of the wooden plate, a coarse mesh in the x-y plane direction give elements that have a much larger length and width than depth.
Frequency response
The frequency response was compared between all three models at the point of load application. The point lies roughly at the beginning of the length and in the middle of the width of the plate. At some of the eigenfrequencies, the mode shape shows that this point should not show a large response and causes some peaks to be higher than others in the frequency response graph. The damping model used in all cases is Rayleigh damping. The values of alpha and beta (also referred to above as alpha0 and alpha1) are the same as the ones used for the Kirchhoff plate and shell models mentioned in the previous section.
The clamped case
3 Layer Model
1 Layer
Ansys
Because of the computational time involved, a coarse mesh was used to evaluate both models, and thus both act quite stiffer than the Ansys model. This causes the eigenfrequencies to be shifted to the right side of the spectrum.
Unsupported case
3 Layer Model
1 Layer
Ansys
The same behavior can be also observed in the unsupported case. It is also noticeable that the response at the first eigenfrequency is nearly completely damped out in both models due to the shape of the first eigenmode and the location of the response measurement location.
Conclusion
In terms of numerical calculation, it was undoubtedly evident that the Ansys model and choice of elements, degrees of freedom, and order of shape functions were the most accurate and no proposed model could outperform its performance. In terms of the proposed models, the Kirchhoff plate element was the most computationally efficient and was quite accurate and close to the Ansys model results up to the 8th eigenfrequency for the unsupported case. The shell model was slightly more accurate, however, the gain in accuracy does not outweigh the computational cost. The solid elements with their proposed formulation were shown to be a poor choice in this case as the results were less accurate than those of the Kirchhoff plate and shell elements. In order for the results to be better using the solid elements, the length, width, and depth of the chosen element should be very close to each other and the mesh used should be fine. Another way to improve the results obtained by the solid elements is by using higher-order shape functions. It was also concluded that poor elements tend to overexaggerate the stiffness of the structure which should be taken into consideration when taking safety factor engineering decisions.
It was more accurate to use the three-layered approach than homogenizing the elements in representing the layered aspect of the cross-ply wooden laminate. However, when modeling large systems the saving in computational time would outweigh the accuracy limitations.
The plate dimensions made it quite favorable for the plate and shell theories to apply and be of high accuracy in representing it. The reason being is that the length and width of the wooden structure are far greater in magnitude than the thickness. However, when dealing with a higher thickness using such elements can give inaccurate results.
The Rayleigh damping proved to better represent the true damping behavior of the structure than hysteric damping used in the previous study. Moreover, Cauphey damping proved to be the best for this application as it used more variables to determine the damping matrix and the results were closest to the experimental values measured in the previous study.
References
[Barrales] Barrales, F. R. (2012). Development of a nonlinear quadrilateral layered membrane element with drilling degrees of freedom and a nonlinear quadrilateral thin flat layered shell element for the modeling of reinforced concrete walls.
[Humar] Humar, J. L. (2012). Dynamics of Structures.
[Schneider et al.] Schneider et al (2020). Bayesian system identification in linear structural dynamics with frequency transformed data using rational surrogate models.
[Müller] Müller, G. (2020). Structural Dynamics (Lecture Notes).
[Bletzinger] Bletzinger, K. U. (2013). Theory of Plates (Lecture Notes).
[Duddeck] Duddeck, F. (2019). Computational Modelling of Composite Structures (Lecture Notes).
[Winter] Winter, C. K. (n.d.). Frequency-Dependent Modeling for the Prediction of the Sound Transmission in Timber Constructions.
[Batoz et al.] Batoz, J.-L. and Tahar, M. B. (1982). Evaluation of a new quadrilateral thin plate bending element.
[Brandner] Brandner, R. (2013) Production and Technology of Cross Laminated Timber (CLT): A state-of-the-art Report
[Buildup-EU] Buildup, EU (2021) Zero Carbon Buildings 2050