Numerical Simulation of the Forming Process of Veneer Laminates

Abstract

In automotive manufacturing, laminated veneer sheets are formed to have 3D geometries for the production of trim parts with wood surfaces. Nowadays, investigation of the formability requires extensive tests with prototype tools, due to the brittle, anisotropic and inhomogeneous material behaviors. The present paper provides numerical methods for the simulation of the forming process of veneers with non-woven backings. Therefore, a conventional forming process of an interior trim part surface is carried out experimentally and numerically, using veneer samples with different individual textures originating from the characteristic growth ring structure. Gray scale images of these samples are mapped to finite element models to account for the wood-specific structure. The forming simulation process comprises two steps, where a gravity simulation depicts the initial position of the blank sheets and the closing of the tool induces the material deformation. The virtual forming of the digital twins accurately reproduces the wrinkling behavior observed in experimental studies. Based on the proposed methods, the design process of manufacturing wood trim parts based on tedious prototype tooling can be replaced by a fully virtual forming process taking into account the individual growth-related properties of the veneer structure.

1. Introduction

Decorative automotive interior trim parts with wood surfaces are conventionally realized with thin veneer sheets, non-woven backings and plastic supporting structures. A forming process gives the 3D shape of the veneer laminate surface material. The forming behavior depends on the grain orientation, the material moisture content and the environmental temperature due to the heated forming tools. Additionally, local deformation, fracture and wrinkling varies with the type of the used veneer, e.g., burled or sliced veneer. Nowadays, the development of a stable forming process for series manufacturing is carried out with time- and cost-intensive tests with prototype tools. A feasible trim part geometry and suitable process parameters are derived in trial-and-error forming tests. Thus, the focus of the present and previous contributions [1,2] is the development and validation of a numerical methodology that replaces experimental trials.

The numerical simulation of the forming process requires adequate material modelling of the veneer laminate. In Ormarsson and Sandberg’s work [3], a numerical analysis of the forming process of birch wood veneer layers into a curved chair seat was carried out. The study contributes to the understanding of the influence of heat, pressure and fiber orientation on distortion and shape stability of the furniture structure. Opposed to the more homogeneous wood species of birch wood, experimental analysis of ash wood veneer laminates of Zerbst et al. [2] showed that the local deformation and fracture behavior of the material strongly depend on the periodic inhomogeneity caused by early wood (EW) and late wood (LW) zones in the veneer. In early wood, wide vessels and thin cell walls cause lower stiffness and strength compared to late wood. These structural differences contribute significantly to the variability of the material properties [4,5,6,7,8].

Stochastic modelling approaches depicting the uncertainty of wood material mechanics are presented in references [9,10,11,12]. However, the design of the veneer laminate forming process requires discrete consideration of the structural differences of early and late wood in areas of large deformation as the range of veneers with different textures and properties are used for trim part surfaces. A modelling approach based on the extended finite element method with repetitive early and late wood unit cells was presented in Lukacevic et al.’s work [13] in order to depict local material fracture. Different sources deal with the prediction of the effective strength of timber, considering the influence of knots and fiber deviations in finite element analysis, with the mapping of surface laser scans [14,15,16,17]. A mapping scheme was developed in Zerbst et al.’s work [1] to capture the geometric arrangement of early and late wood over a veneer sheet due to the veneer manufacturing process, e.g., slicing or rotary cut technology. Thereby, the color differences of early and late wood, as a result of significant differences in density, were used for the mapping of the two zones from gray scale images to a finite element mesh. Stochastic numerical analyses of tensile tests showed very good agreement of the distributions of strength and ultimate strain compared to experimentally obtained data and validated the approach. These results promised to capture the variation of material fracture for the suggested procedure.

In this work, the before-mentioned modelling approach is applied to establish a virtual forming process in order to speed up and optimize the development of wood surfaces for automotive trim parts. Therefore, blank sheets of ash wood veneer laminate with different textures are analyzed experimentally and numerically with finite element models, created with the gray scale mapping procedure. The objective is to provide a preprocessing chain in combination with a two-step forming simulation that completes the virtual design of the forming process of veneer laminate surfaces.

2. Materials and Methods

2.1. Material

Ash wood veneer (Fraxinus excelsior L.) was used for the characterization of the forming properties and the numerical simulation of the forming process. The veneers were laminated with a non-woven fabric. The lamination was used to reduce the bending stiffness by reducing the stiffer wood layer in the effective cross-section [18]. According to the manufacturer PWG VeneerBackings GmbH (Walkertshofen, Germany), the non-woven fabric consists mainly of long-fiber cellulose with smaller proportions of synthetic fibers [19]. The bonding of the non-woven was carried out via hydro entanglement with additional use of a latex binder. In the lamination process, the layers of veneer and non-woven fabric were bonded via a phenolic resin adhesive and then calibrated to a total thickness of 0.5 mm using a sanding machine. In the calibrated laminate, the thickness of the veneer layer was 0.2 mm. During lamination, the adhesive penetrated into the wood and the non-woven fabric. Hence, the mechanics of the individual layers could no longer be considered separately [2]. In the works of Gindl et al. [20] and Konnerth and Gindl [21], it was shown that wood cell walls have a significantly higher stiffness in the transition area of a bond. Therefore, the characterization and material modeling were carried out on the composite, which will be referred to as veneer laminate in the following. Under normal climatic conditions (20 °C/65% relative humidity), the moisture content of the veneer laminate was 7% with an average density of 780 kg/m³.

2.2. Experimental Forming Analysis

Forming of wood surfaces for a trim part component was performed and the results are compared with simulations. The real forming was carried out in the same way as in the series production process for trim parts. A wooden insert for the palm rest of the Mercedes-Benz X167 series was used. This component is relatively small. Thus, the forming is more dependent on the individual annual ring structure and therefore well suited for model validation. The forming tool is made of aluminum. It consists of a male punch and a female die for forming the topology and two cylindrical pins for fixing the blank, onto which the blank is fitted with holes (Figure 1).

Four different blanks of laminated veneer were used for the experimental and numerical trials. They have the same sample geometry as in the series production process (Figure 2). A CNC-controlled milling machine cut the samples. Ash wood veneers with different textures were used to analyze the influence of the growth ring structure on the forming result. This range corresponds to veneer samples such as those used typically in Mercedes-Benz vehicles. The grain direction was aligned parallel to the length direction of the vehicle for three specimens (Figure 2a–c), as this represents the more critical case. Additionally, the grain direction of one specimen was aligned perpendicularly (Figure 2d).

Samples were put into water at 20 °C for 1 h in order to generate a defined moisture content, as proposed before by Zerbst et al. [2]. The samples achieved a moisture content of around 52%, as determined by the specific weight of a tensile sample after water immersion and oven drying. The forming tool was heated to a constant temperature of 140 °C, as was also carried out in the series production process. The blanks were then placed on the locating pins. The punch was moved at high speed to a distance of 10 mm from the female die. By closing the forming tool at a constant speed of 1 mm/s, the blank was formed and a 3D wooden shell was created under deformation and plasticizing influence of moisture and heat [22].

2.3. Virtual Forming Process

2.3.1. The Digital Twin of the Blank

A method for a virtual forming process was developed. Finite element (FE) models of the blank samples used in the experimental forming were created. A gray scale mapping scheme as described by Zerbst et al. [2] was implemented in the mapping software Envyo^®. It accounts for early wood (EW) and late wood (LW) areas in the veneer and transfers them to the finite element model. The blank samples, depicted in Figure 2, were scanned individually with a conventional flatbed scanner, resulting in a colored image of each structure. Using the image manipulation program GIMP, the picture was transformed to a gray scale image and exported as ASCII-based portable gray map (pgm) with gray scale values between 0 (completely black) and 255 (completely white) for each pixel. To avoid the transfer of single peaks or artefacts, a Gaussian filter was used to flatten the gray scale values of the picture. As a target for the mapping procedure, an FE model of mainly fully integrated quadrilateral shell elements was created, assuming plane stress for the small ratio of thickness to length of the selected component. The chosen element length of 0.5 mm led to a total number of 34,591 elements for the blank model. After aligning source pixel mesh and target finite element mesh, a nearest neighbor search assigned the corresponding gray scale values to the target elements. In a second step, the elements are assigned to parts in the FE software LS-Dyna that correspond to EW and LW material models (see Figure 3). To obtain consistent results, the gray scale range used to assign EW and LW zones was adjusted for each specimen in a way that a measured width of the EW zones corresponded to the width assigned as EW zone in the finite element mesh. The parameters found for the four specimens used for evaluation purposes of the method are given in

Table 1. Gray scale ranges for the mapping of early wood (EW) and late wood (LW) zones of all samples.

Sample No.	Gray Scale Ranges
	LW	EW
${\mathbf{1}}^{\parallel }$	0–212	213–255
${\mathbf{2}}^{\parallel }$	0–205	206–255
${\mathbf{3}}^{\parallel }$	0–205	206–255
${\mathbf{4}}^{\perp }$	0–204	205–255

.

Considering the global mechanical behavior of the laminate, a material model with the name *MAT_LAMINATED_COMPOSITE_FABRIC (*MAT_058) implemented in the FE software LS-DYNA was chosen (Figure 4). *MAT_058 is commonly used to model the transversally anisotropic behavior of laminated composite fabrics [23,24,25] but has been used by for other procedures due to its capability to describe the non-symmetric tension/compression behavior of wood [26]. The advantage of *MAT_058 opposed to other material models is the possibility to consider a non-linear stress–strain relationship under tension and shear. This depicts the behavior of veneer laminates very well, as was shown in our previous paper [1].

The implementation described by Schweizerhof et al. [23] is based on the work proposed by Matzenmiller et al. [27]. In order to reduce efforts for material testing and FE modeling, the material was characterized as a single laminate, avoiding the separate considerations of the single layers. In general, the laminate behavior with respect to stiffness and strength was dominated by the veneer sheet, as stated by Zerbst et al. [2]. For the veneer laminate, material failure was identified for the directions parallel (X) and transverse to the grain orientation (Y), with a maximum stress criterion, under tensile (Index T) or compressive loading (Index C) or under in-plane shear loading ($S_C$). As described by Zerbst et al. [1], those failure criteria were calibrated for EW and LW material through reverse-engineering, using the optimization software LS-OPT. For the adjustment of the non-linear curve fitting, Young’s moduli ($E_{ij}$) and maximum strains (${\epsilon }_{ij,T}$) were identified under tensile loading for the longitudinal (1) and the transverse (2) material directions. Material parameters are provided in

Table 2. Calibrated material parameter in tensile loading direction for EW and LW obtained in [2].

Part	${\mathit{E}}_{11}$	${\mathit{E}}_{22}$	${\mathit{X}}_{\mathit{T}}$	${\mathit{Y}}_{\mathit{T}}$	${\mathit{\epsilon }}_{11,\mathit{T}}$	${\mathit{\epsilon }}_{22,\mathit{T}}$	${\mathit{\nu }}_{12}$	${\mathit{\nu }}_{21}$
	(MPa)	(MPa)	(MPa)	(MPa)	(-)	(-)	(-)	(-)
LW	4452	600	61	27	0.039	0.116	0.42	0.05660
EW	2000	136	30	12	0.056	0.207	0.42	0.02751

.

Picture frame tests of Zerbst et al. [2] were used to calibrate material behavior under shear loading conditions. The derived parameters for failure ($S_C$, ${\gamma }_{12,S}$) under pure in-plane shearing conditions and the corresponding input values for the stress–strain curve ($G_{12}$, ${\tau }_{12,A}$, ${\gamma }_{12,A}$) are provided in

Table 3. Calibrated material parameter for in-plane shear behavior for EW and LW.

Part	${\mathit{G}}_{12}$	${\mathit{\tau }}_{12\mathit{A}}$	${\mathit{\gamma }}_{12\mathit{A}}$	${\mathit{S}}_{\mathit{C}}$	${\mathit{\gamma }}_{12,\mathit{S}}$
	(MPa)	(MPa)	(-)	(MPa)	(-)
LW	800	9	0.064	15	0.212
EW	300	7	0.15	8	0.360

.

Since testing and measuring of material behavior under compression is challenging due to the tendency of the laminar structure to buckle out of the loading plane, these parameters were derived from literature data. Based on measurements of tension and compression behavior of ash wood and beech wood of Niemz and Sonderegger [28] and Ozyhar et al. [29], compressive failure parameters were derived as provided in

Table 4. Estimated material parameters for compressive loading for EW and LW.

Part	${\mathit{X}}_{\mathit{C}}$	${\mathit{Y}}_{\mathit{C}}$	${\mathit{\epsilon }}_{11,\mathit{C}}$	${\mathit{\epsilon }}_{22,\mathit{C}}$
	(MPa)	(MPa)	(-)	(-)
LW	24	18	0.026	0.078
EW	12	8	0.038	0.140

, assuming the same ratio of tension to compression strength and maximum strain. Thereby, the tension–compression strength ratio was chosen as 3:1 parallel to the fiber direction and 1.5:1 perpendicular to the fiber direction for EW:LW. Strain values for maximum strain were reduced by 30% with respect to the tensile loading direction.

In order to account for the remaining stresses after reaching the maximum stress limit prior to element erosion due to maximum strain, *MAT_58 provided the SLIMi values that allow this minimum stress to be defined depending on the loading direction. In this work, stresses were reduced to 50% of the maximum stress in the tensile direction and were kept at the maximum stress level for in-plane shear and compressive direction, as illustrated in Figure 4. The initial main material orientation was assigned using the AOPT—flag equal to 2, allowing for a direction definition in the global coordinate system. The global x-direction (samples 1-3) or y-direction (sample 4) was chosen, which are parallel or transverse to the vehicle travel direction (see Figure 2), respectively. To consider the non-linear shearing behavior, the failure surface flag FS was set equal to −1. Other parameters not listed here were left undefined or were implemented as default values, respectively.

2.3.2. Simulation Step 1: Initial Conditions

This first simulation step accounts for the initial position of the wood component due to gravity loading. The static implicit analysis solution scheme was chosen in order to save calculation time, neglecting inertia and oscillation effects. Due to the expected small displacement of the specimen and a low bending degree, this is a valid approach. As shown in Figure 5, the specimen was placed in the XY plane of the global coordinate system above the die, which is modeled as a rigid body. A total of 43,355 shell elements were used. The gravity loading was applied to the blank part nodes. Segment-based contact algorithms were used between the guiding pens and the specimen. Standard two-way contacts were used between the tool and the specimen. For all contact definitions, a friction coefficient of 0.15 was assumed, based on literature values from measurements of friction between steel surfaces and different wood species at low velocities [30]. The overall time considered in the simulation is 1.5 s. All simulations have been performed on a Linux-based x86 cluster, using 96 parallel CPUs with 5 Gb RAM for both simulation steps. The results presented have been obtained using LS-DYNA-version R11.1 double precision. Within the used environment, a simulation time of about 1 h and 20 min was needed to run this first step of the process simulation.

2.3.3. Simulation Step 2: Forming

In the second step, the forming simulation was performed using the estimated initial position from simulation step 1. To do so, the rigid punch model was included (Figure 6). The contact between the guiding pens and specimen was changed to a standard penetration check for nodes to a surface, whereas the contact between the specimen and the tools stayed the same. The simulation was driven by the displacement of the stamping tool, which is supposed to travel 31 mm. The specimen was expected to generate wrinkles, which would lead to instabilities when considering self-contact of the specimen. Therefore, penetrations of the wood specimen were tolerated to allow the stamping tool to travel the full length without creating unphysically high contact forces. Since the method described should allow the detection of wrinkles and not the exact alignment of wood fibers during the forming processes, neglecting the self-contact does not contradict this task.

A total time of 0.1 s was considered in the simulation, which is a 31-fold increase, compared to the true forming velocity, using the explicit solution algorithm with single precision. In the environment described above, the total simulation time was between 8 and 10 min.

3. Results

For a validation of the used methods with respect to the quality of the transformation of the true process into a virtual process, experimental and numerical results were compared. Under gravity loading in the first step of the virtual forming process, the blank shifts with the cantilever sight towards the die and slides with a transient oscillation moves into the balance position. The balance position of the simulations (Figure 7a) was the same as observed for the forming experiments (Figure 7b), where the blank was placed by hand. Due to the different swelling behaviors of the veneer and the non-woven layer in consequence of the water pretreatment, the real samples had a curvature from the beginning. This was not considered in the simulation model. However, the gravity simulation provided a good approximation of the initial blank position before the forming started. The determination of the starting conditions of the forming process becomes even more important for larger pieces of trim part surfaces. Those veneer laminate blanks would bend much more in its lying position on the die topology. The effect of a curvature due to the water pretreatment is expected to be lower.

The simulation of the forming process enables the observation of the forming history during the moving punch. This allows analysis of the development of critical deformations (Figure 8). Due to the convex geometry, there is tendency of material compression within the chosen part. Those compressive loads can either be distributed over the sample or localize into buckling and wrinkling modes, depending on the arrangement of the growth rings and the grain orientation. Consequently, the main critical areas are at the circumferential edges.

The experimental and numerical forming results are compared below. Characteristic damage phenomena are identified within six zones on the top view images of the physically formed samples (Figure 9a, Figure 10a, Figure 11a and Figure 12a). For further visualization of critical deformations, the simulation results are given in terms of effective strain $\widehat{\epsilon }$ (Figure 9b, Figure 10b, Figure 11b and Figure 12b). The effective direction of strains as well as their compression or tensile origin is smeared within this scalar measure, with the following equation:

$$\widehat{\epsilon }=\sqrt{\frac{2}{3}{\epsilon }_{ij}{\epsilon }_{ij}}\hspace{1em} \hspace{1em}\mathrm{where} i,j=1,2$$

(1)

In general, the numerical results provide a very good agreement with the experimental results. The characteristic deformation was captured for all samples with the orientation of the grain direction being parallel as well as perpendicular to the vehicle direction. Corresponding with the specific wood structure, all physical samples and their discrete finite-element representations show buckling and wrinkling behavior, which is mainly induced through the early wood zones. Here, the surrounding stiffer late wood zones thrust the weaker early wood zone. If the geometrical arrangement of the early wood zones is unfavorable, this leads to out-of-plane movement. For early wood zones with higher width, this results into the development of wrinkles, as can be seen in samples ${\mathbf{1}}^{\parallel }$ and ${\mathbf{2}}^{\parallel }$ (Figure 9 and Figure 10) in the zones 2 and 3. The constitution of sample ${\mathbf{3}}^{\parallel }$, with early wood zones of smaller and uniform width, results in slender compression lines due to forming (Figure 11). The buckling around the pin holes (zones 4 and 6) is also depicted very well by the FE models of samples ${\mathbf{1}}^{\parallel }$, ${\mathbf{2}}^{\parallel }$ and ${\mathbf{3}}^{\parallel }$ (Figure 9, Figure 10 and Figure 11).

For the orientation of the grain direction perpendicular to the vehicle direction of sample ${\mathbf{4}}^{\perp }$, the simulation shows a larger wrinkle in zone number 3 as well as compression lines running from the pin holes (zones 5 and 6). This is in accordance with the experiment (Figure 12). Some damage is hard to see on the images of the physically formed samples, especially if compression lines developed parallel to the grain direction (Figure 9, Figure 10 and Figure 11 zones 1 and 5, Figure 12, zones 2 and 3). Consequently, such failure is evaluated to be less critical for the assessment of a stable forming process under esthetic criteria.

A larger wrinkle occurred in the middle of zone number 5 for the simulations of all samples. In the experimental forming process, this wrinkle was smaller and not always concentrated to just one early wood zone, as opposed to model predictions. It has to be noted that the parameters of the material models were calibrated in the experiments after water immersion. In the forming process, increased temperature of the tools was used, which generated additional straining capacity due to the softening of wood cell wall components [31,32,33,34,35]. Thus, it was assumed that the simulations represent a more critical case with respect to the material formability. However, the simulations with the proposed modeling process captured the forming behavior of individual wood surfaces accurately.

4. Conclusions

In the present study, veneer sheets were formed physically and numerically into an automotive trim part geometry. Digital twins of the natural samples were created by the gray scale mapping method. The initial position of the blank sheets was considered with a gravity simulation in a good approximation. This first step of the virtual forming process is more important for larger parts, where higher bending over the die is expected. The two-step forming simulation depicted the deformation behavior of each sample with very high accuracy. These results clearly showed that the structural differences between early and late wood and their distribution over the blank sheet are the critical factors for the prediction of the forming behavior of ash wood veneer surfaces. The presented virtual forming process together with the introduced material modelling approach can be used for the development of a stable design of manufacturing processes of veneer forming. Optimizations of the fixation concept or adaptions of the geometry of the tools can be derived based on these simulations. Furthermore, the natural variation of veneer surface composites can be analyzed to make suggestions for the improvement of the material design, e.g., by modifications of the supporting layer of the non-woven backing.

The forming simulation was validated with ring porous ash wood veneer laminate. Apart from that, the virtual process chain can be used for structural-mechanical analysis of many other wood species and applications. High structural density variations due to the season-dependent growth of the tree can also be found for softwood species such as spruce and pine. These species are frequently used in civil engineering, where numerical methods are in high demand for the design of timber structures [36]. These methods can also help to analyze individual conditions in existing constructions.