Introduction

The expansion of the hydrogen economy offers an excellent opportunity to reduce pollutant emissions and establish hydrogen as a secondary energy carrier in the context of the energy turnaround [1]. Electrolyzers and fuel cell stacks, which are required as electrochemical energy converters, include numerous bipolar plates (BPP) and bipolar half plates (BPHP) for media distribution as essential components. For passenger car applications, a single currently available PEMFC stack (model NM12, 123 kW, ElringKlinger AG) requires the production of approximately 720 BPHP. With increasing demand for stacks the need for individual BPPs rises in correspondingly higher numbers. Therefore, their high-volume production has to be ensured to make fuel cell technology available to the mass market [2]. Metallic BPHPs can be formed by several technological approaches, which differ in the basic working principle and the achievable production rates [3].

Hydrofoming is an established process for forming BPHPs, which uses fluid pressure of up to 250 N/mm2 to form the sheet into the die [4]. The technological advantages are good surface quality, comparatively low flatness deviation of the formed BPHP, and high shape accuracy of the channel cross sections [5]. In addition, only one forming tool is required, eliminating the need to align multiple dies to each other with high accuracy requirements [6]. A limiting aspect is the high closing force required for sealing (e. g. 30,000 kN for an exemplary BPHP with a projected pressure area of 800 cm2), which has to be applied by special presses. In addition, auxiliary devices are required to generate the forming pressure. In perspective, passive high-pressure sheet metal forming can be used with suitable systems to generate the forming pressure directly in the forming tool by the closing motion of the press without external pressure devices [7]. This offers high potential to increase the production rate, which is currently estimated at 7 BPHP/min (without parallelization) for conventional hydroforming process variants [3, 6].

Hollow embossing is another established forming process for the production of BPHPs [8]. Similar to classic deep drawing, the semi-finished part is clamped and formed between two rigid tools (punch and die), thus causing thinning and an additional blankholder can limit the lateral draw-in of the blank in order to reduce wrinkling and flatness deviation [3]. The die and punch are aligned by guides in the forming tool, which requires high manufacturing accuracy of the tool parts and guiding elements [6]. Current research includes multi-stage forming [9] and temperature-controlled process variants [10]. To reduce tooling costs, a rubber pad can be used instead of a rigid counter die in rubber pad forming [11, 12]. However, even using hollow embossing forming processes, only limited production rates (< 60 BPHP/minute) can be realized [3, 6].

In perspective, production of single BHPH by single point incremental forming can be useful for prototyping as low-cost tools are used [13]. Current research and development also focuses on electromagnetic forming of metallic BPHPs [14] to take advantage of the positive effect of higher strain rates on the forming process [15].

The potential for mass production of BPHPs by higher productivity and lower production costs has been a major driver in recent years for the development of rolling processes that can also be used to incrementally form channel structures in thin strip [3]. Several basic research papers have been published on this topic: The technological approach of conventional roll forming according to Zhang et al. [16] and Abeyrathna et al. [17] is limited to continuous longitudinal channels due to the multistage forming process by several rollers and therefore it is not relevant for functional flow field designs. Roll forming of transversal channel structures for BPHPs was published by Nikam et al. [18] and further investigated by Huang et al. [19] demonstrating the possibility to produce transversal channels by incremental roll forming. However, these approaches are also limited to continuous channel geometries.

The forming of discontinuous channel structures by hollow embossing rolling was first described by Bauer et al. [20] showing the potential for incremental forming of typical BPHP flow field designs for industrial applications. Further research activities provided initial approaches for forming process optimization based on FE simulations [21]. However, continuous further development of hollow embossing rolling is required to enable future series production of BPHP in common sizes and featuring typical flow field designs [6].

In this context, a key means is the modeling and simulation of the forming process to obtain deep process understanding as an essential requirement for optimization. The main challenges here arise from the incremental forming character and the resulting long calculation times. Furthermore, delicate channel structures within sub-millimeter scale and thin metal strips (0.05 – 0.10 mm) require extremely fine meshing of the strip section and tools resulting in high element numbers. In general, the finite element (FE) model must accurately predict different scales from the macroscopic plate flatness in the millimeter range to the local channel forming in the scale of some tens of a millimeter up to the strip thinning in the range of a few micrometers. Impacts on the frictional influence and material properties can also be expected due to the existing scaling effects [22].

The previous work of Bauer et al. [20, 21] provided the basis for simulation-based process design and optimization. However, the FE models using simufact forming have only been validated on simple quasi-continuous longitudinal channel designs of small size. The applicability of the modeling approaches to typical BPHP designs with characteristic discontinuities is currently unknown. Another disadvantage of the modeling approaches published so far is the use of volume elements for meshing the strip, which leads to long computation times for larger plate designs (100 × 50 mm) [23]. It is not clear whether it can be used for simulation-based process optimization considering standard industrial plate sizes in the range of 400 × 150 mm.

This investigation aims to analyze the potential of a shell-based FE model using LS-DYNA® to enable technological process optimization with acceptable computing times. For this purpose, the FE model is validated on an experimental basis of a representative forming demonstrator with characteristic geometrical plate design discontinuities of stainless steel 1.4404 (thickness \({t}_{0}\) = 0.1 mm). The evaluation comprises both the accuracy of roll geometry and the influence of roll system elasticity on forming-specific results.

Modeling of hollow embossing rolling

Material characterization and modeling

The austenitic stainless steel X2CrNiMo17 12 2 (1.4404) with a sheet thickness of \({t}_{0}\) = 0.1 mm (uncoated) and an average grain size of 18 µm (ASTM E112, line section method) was investigated as a typical representative material used in proton-exchange membrane fuel cell applications. Quasi-static tensile tests according to DIN EN ISO 6892–1 were performed using a ZwickRoell Z020 uniaxial testing machine (ZwickRoell AG, Ulm, Germany, Fig. 1 a) to determine basic mechanical properties and flow curves for different orientations to the rolling direction (RD 0°, 45° and 90°) at a strain rate of \({\dot{\varepsilon }}_{pl}\) = 0.001 s−1. The strain in the specimen was evaluated using an optical system (Aramis, Carl Zeiss GOM Metrology GmbH), because the small sheet thickness prevented the use of tactile sensors.

Fig. 1
figure 1

Tensile test set up: (a) quasi-static tensile testing, (b) tensile testing at higher strain rates

It is known from preliminary studies on the simulation of hollow embossing rolling that the strain rates can increase locally during channel formation. For this reason, additional tensile tests for RD 0° were performed at increased strain rates of \({\dot{\varepsilon }}_{pl}\) = 0.1, 5.0 and 20 s−1 using a ZwickRoell HTM 16020 high speed testing machine with optical strain measurement (Aramis, Carl Zeiss GOM Metrology GmbH with Photron FASTCAM SA1.1 high speed cameras, Fig. 1 b). Buckling of the thin, wire-eroded specimens due to the long free specimen length was avoided by applying support sheets on both sides. The material investigated is characterized by significant work hardening (\(UTS\)\(YS\)) and high uniform elongations \({A}_{G}\) (Fig. 2), which allows the experimental determination of true stress-true strain curves up to high strain levels (\({\varepsilon }_{pl}\) = 0.42 – 0.47, Fig. 3 a). At quasi-static strain rates, the material shows only minor differences in mechanical properties per orientation to the rolling direction, with only the \(r\)-values showing a stronger dependence on orientation to the rolling direction. A significant strain rate dependence of the mechanical properties can be observed particularly considering the yield strength. The ultimate strength is lower at \({\dot{\varepsilon }}_{pl}\) = 0.1 s−1 compared to the quasi-static reference but it increases again at even higher strain rates. The reason for this was identified by thermography. At higher strain rates the heating of the tensile specimen, which causes a local maximum temperature of 74 °C at \({\dot{\varepsilon }}_{pl}\) = 0.1 s−1, leads to less strain hardening. This can also be seen comparing the true stress-true strain curves at higher strain rates from \({\varepsilon }_{pl}\) > 0.15 (Fig. 3 b). In contrast, the \(r\)-values for RD 0° (\({r}_{0}\)) remain relatively stable at higher strain rates (Fig. 2 c).

Fig. 2
figure 2

Averaged tensile properties of 1.4404 for different strain rates and orientations to the rolling direction: (a) yield strength (\(YS\)), (b) ultimate tensile strength (\(UTS\)), (c) \(r\)-values

Fig. 3
figure 3

Averaged true stress-true strain curves and approximations for: (a) quasi-static strain rates and varying orientations to the rolling direction (RD), (b) varying strain rates

True stress-true strain curves were obtained by averaging three tensile test results for each parameter set (Fig. 3 a) and approximating the resulting mean curve according to Hockett-Sherby [24] (Eq. 1, Table 1) for RD 0°. The averaged true stress-true strain curves for RD 45° and RD 90° are only slightly below the yield stress level of RD 0° and show similar strain hardening behavior. Therefore, the experimental data points of the true stress-true strain curves for RD 45° and RD 90° are approximated by scaling the RD 0° approximation with the identified initial yield stresses \({\sigma }_{45}\) and \({\sigma }_{90}\) (Table 2) relative to \({\sigma }_{0}\). The derived true stress-true strain curves at higher strain rates for RD 0° (Fig. 3 b) were further approximated by a modified Cowper-Symonds approach [25] (Eq. 2, Table 1) using only the range up to \(\varepsilon\) ≤ 0.15 due to the sample heating problem described. Additive strain rate modeling according to Eq. 2 resulted in smaller approximation errors than the standard multiplicative variant.

Table 1 Parameters describing the RD 0° true stress-true strain approximation and strain rate sensitive modeling
Table 2 Parameters for yield curve approach
$${\sigma }_{HS}={\sigma }_{sat}-\left({\sigma }_{sat}-{\sigma }_{ini}\right)\bullet {e}^{-a\bullet {{\varepsilon }_{pl}}^{p}}$$
(1)
$$\sigma ={\sigma }_{HS}+{C\left(\frac{{\dot{\varepsilon }}_{pl}}{{\dot{\varepsilon }}_{pl,ref}}\right)}^{1/m}$$
(2)

Extended material testing for biaxial tension was performed by bulge tests in a test setup according to DIN EN ISO 16808 (Fig. 4 a) with optical strain measurement. The biaxial true stress-true strain approximation (Fig. 4 b) was determined by referencing the identified uniaxial approximation \({\sigma }_{HS}({\varepsilon }_{pl})\) and considering the ratio \({\sigma }_{b}/{\sigma }_{0}\) (Table 2), which was determined by least squares minimization of the averaged data from the bulge test.

Fig. 4
figure 4

Bulge testing: (a) test set up, (b) derived biaxial true stress-true strain curve

Modeling the elastic–plastic material behavior for varying, multiaxial loading conditions requires the selection of a suitable yield approach. Different yield approaches were compared with respect to their approximation quality using the values determined experimentally (Table 2). Established approaches of lower complexity, such as v. Mises or Hill90 [26], show deficits in the exact reproduction of the yield stress or \(r\)-value anisotropy and partly fail to meet the specification of the biaxial initial yield stress \({\sigma }_{b}\) (Fig. 5). For this reason, a more complex yield approach YLD2000-2d according to Barlat [27] was chosen with a yield exponent \(M\) = 6, which considers all determined material properties. Due to the comparable hardening behavior of the true stress-true strain curves for RD 0°, 45°, and 90°, as well as the nearly constant \(r\)-values at higher strain rates, isotropic hardening is assumed based on the true stress-true strain curve in RD 0°.

Fig. 5
figure 5

Selection of the yield locus approach based on the approximation quality of the measured characteristics

The selected yield approach was implemented in the FE simulation software LS-DYNA® R12.1 using the material law "MAT_BARLAT_YLD2000" (MAT133) with YLD2000-2d. The yield curve approximation, including the strain rate sensitivity, was stored in tabular form as yield curves for discrete strain rates. Isotropic strain hardening was assumed.

Bipolar plate demonstrator and model setup

Based on common flow field designs and cross-section geometries [6], a miniaturized BPP geometry was developed as a forming demonstrator for simulation studies (Fig. 6 a). This design advances previously published test geometries for validating experimental data with FE simulation models by considering geometric discontinuities typically occurring in real-size plate designs for automotive applications, such as changing channel orientations (longitudinal and transversal), channel arcs and step-shaped transitions.

Fig. 6
figure 6

Boundary conditions for model development: (a) BPP demonstrator, (b) process realization concept

The design of the FE model for hollow embossing rolling results from the process concept shown in Fig. 6 b. The strip first passes through feeding rollers, which use a braking torque to build up restraining tension in the strip. This restraining tension is referred as \({\sigma }_{R}\) in the subsequent discussion. After incremental forming of the channel structure in the second rolling stage, guiding rollers ensure that the strip is safely carried away. The model for the explicit forming simulation was set up in the FE software LS-DYNA® R12.1 mpp (Fig. 7 a). The geometries of the forming rollers (Ø105) were implemented as rigid shell meshes rotating around a stationary turning point at a rotational speed of 20 rpm. A target rolling gap corresponding to the initial strip thickness of \({t}_{0}\) = 0.10 mm was set between the roller surfaces \({s}_{w0}\). Taking advantage of symmetry, only half of the 30 mm wide strip was meshed in the forming area with a uniform element edge length of 0.075 mm (in total 110,000 elements).

Fig. 7
figure 7

Modeling of hollow embossing rolling: (a) forming simulation model, (b) boundary conditions for tool and strip meshing

The strip discretization was realized by using established shell element types for forming simulation with LS-DYNA (Belytschko-Lin-Tsay shells) with three integration points over the shell thickness. Additionally, two other rigid bodies were integrated into the model, thus implementing the translational constraint of the strip section synchronously to the roller rotation (holder) and applying the restraining tension \({\sigma }_{R}\) to the strip on the inlet side using force control (counterhold). The numerical setting parameters for mass scaling, time acceleration, and contact definition had been evaluated within separate numerical studies and were selected with respect to a short computation time and at the same time ensuring that they cause only insignificant changes of the resulting variables relevant for forming [28]. In particular, the required mesh fineness of the strip segment was examined in comparative simulations in order to find a compromise between computation time and model error. The meshing of the strip should therefore allow the smallest tool radius (here \(R\) 0.25) to be represented with five elements, which corresponds to an initial element edge length of 0.075 mm. In this context, the use of the reduced integrated element formulation was also tested, with the result that no significant changes occur in the present case compared to higher order approaches. The use of only 3 integration points over the shell thickness was also justified compared to a larger number of integration points (5 – 9). Following studies on scaling effects [29], the first approach assumed a conservative Coulomb friction model with a static friction coefficient of \({\mu }_{s}\) = 0.15.

Applying these boundary conditions, the computation of this model takes 8.7 h with 32 CPUs for a modelled strip size of ~ 10 cm2 (Table 3). Compared to the published state of the art usually using volume elements, the presented modeling strategy allows an 11-times finer discretization in the strip plane with a lower relative computational time, which illustrates the potential of the chosen shell-based approach. However, the presented performance comparison must be interpreted cautiously, as it is based on different FE systems, boundary conditions and plate designs.

Table 3 Comparison of the calculation efficiency to the state of the art

Analysis of the demonstrator model

Depending on the active forming zone, the incremental forming process alternates between continuous (longitudinal channels) and discontinuous forming conditions (transversal channels). The latter lead to alternating strip movement in the \(z\)-direction on the inlet side due to alternating die and punch side contact if the strip is fed straightly. In the investigated forming process, initial contact with the die side forming roller (top) occurs at 7° before the rolling axis plane is reached (Fig. 8, step 1). Subsequently, the progressive closing movement of the rollers results in a channel forming with locally increasing thinning (step 2), which is particularly noticeable within the channel flanks.

Fig. 8
figure 8

Analysis of the incremental forming behavior using the example of stepwise forming of a channel cross-section segment along the channel flank

When the forming roller axis plane is reached (step 3), the intended rolling gap (\({s}_{w0}\) = \({t}_{0}\) = 0.10 mm) is achieved, and the forming process is completed locally. The proceeding rotational movement finally causes the opening of the rollers and releases the formed strip section. In the present study, the forming simulation is always followed by implicitly springback calculation of the strip, which may have been cut free previously, the subsequent deposition of the demonstrator plate (gravity step) on a plane surface, and a clamping simulation for the final evaluation of the demonstrator plate flatness deviation.

From the analysis of the strain states by the major and minor plastic strains \({\varepsilon }_{pl,1}\) and \({\varepsilon }_{pl,2}\), it is evident that plane-strain conditions dominate in the middle shell layer in the region of the flow field (Fig. 9 a and b). This characteristic strain state is superimposed by a bending deformation, which occurs primarily in the radius and flank regions and results from the differences in true plastic strains between the upper and lower shell layers. Longitudinal channels show more significant true plastic strains \({\varepsilon }_{pl,1,mid}\) than transversal channels of identical cross-section (meander zone). The reason for this effect is that due to the simultaneous rolling of several longitudinal channels no local draw-in in transversal direction of the strip (\(y\)-direction) can occur, while in the case of transversal channels, the freely fed strip allows a certain amount of draw-in and the true plastic strains are consequently lower. Geometric discontinuities such as channel deflections, arcs and step-shaped transitions initiate biaxial strain states. Depending on the strip restraining tension \({\sigma }_{R}\), the part flange is characterized by areas with negative minor strains due to edge draw-in and a varying tendency to wrinkling. The strain paths are largely linear, but due to the bending effects occurring at the radius transitions and in the channel flanks, a conventional evaluation of the forming limits using a forming limit curve is restricted.

Fig. 9
figure 9

Analysis of strain states and strain rates: (a) qualitative strain states, coloring corresponds to the position at the forming limit diagram, (b) forming limit diagram distribution for different shell layers, (c) evolution of strain rates at characteristic points

The incremental nature of the process results in a short forming time with rapidly increasing and decreasing strain rates and a discontinuous dynamic history (Fig. 9 c). Thus, biaxial strain state regions show maximum strain rates up to \({\dot{\varepsilon }}_{pl}\) = 25 s−1 (region C). When forming continuous channel geometries, strain rates are highest in the radius and flank regions (\({\dot{\varepsilon }}_{pl}\) = 18 s−1, regions E – G). Bending effects also cause different strain rates in the individual shell layers. Due to these dynamic properties, the process is sensitive to strain rate dependent material behavior, whereby the strain rates shown are relative to the boundary conditions of the strip forming speed (here 157 mm/s with 20 rpm).

Experimental setup and procedure

For model validation, experimental rolling tests were carried out on a Bühler EW 105 × 100 strip rolling mill (Fig. 10 a) applying the boundary conditions described in chapter 2. The die and punch side roller surfaces (∅105 mm) derived from the upper and lower BPP demonstrator surface geometries were implemented in the forming rollers as milled inserts considering a strip thickness gap of 0.1 mm in between. Additional elements for axial (tongue and groove combination) and tangential alignment (gear profile) were implemented in the inserts to synchronize rolling (Fig. 10 b). Guide rollers were integrated on both sides for lateral strip guidance and to provide strip restraining tension \({\sigma }_{R}\) by weights on the inlet side. The height of the top roller was adjusted manually with a hand wheel thus ensuring secure clamping after contact with the rollers despite of any weight used on the inlet side. In the rolling tests, the restraining tension \({\sigma }_{R}\) on the inlet side was gradually increased from \({\sigma }_{R}\) = 0 N/mm2 to 35 N/mm2, 68 N/mm2, and 101 N/mm2. This was done to verify experimentally whether the restraining tension proposed by Bauer et al. [20] has a positive influence on the forming process and wrinkling tendency.

Fig. 10
figure 10

Experimental setup: (a) rolling machine with strip restraining tension \({\sigma }_{R}\), (b) milled roller inserts for BPP demonstrator forming

For the validation of the FE model, knowledge about manufacturing inaccuracies is important. Therefore, the roller inserts were measured using a coordinate measuring machine ZEISS PRISMO7S-ACC and the actual rolling gap \({\overline{s} }_{w0}\) between the punch and die side insert was determined. Section-wise analysis of the geometrical deviation from the nominal state of CAD showed that the actual gap is significantly smaller than the initially planned target gap due to manufacturing inaccuracies (Fig. 11). Instead of a target value of 100 µm, the average actual rolling gap \({\overline{s} }_{w0,actual}\) only amounts to 78 µm in the measured section areas. The deviations of the cross-section analysis were mapped onto the FE meshes of the rollers to integrate this gap change into the simulation model. The existing rolling gap in the simulation model was thus significantly smaller than the initial strip thickness, which leads to higher forming forces.

Fig. 11
figure 11

Average initial rolling gap \({\overline{s} }_{w0}\) based on section-measurements of the roller inserts

Results and discussion

Experimental results for model validation

Based on preliminary simulation studies, different restraining tensions \({\sigma }_{R}\) on the inlet side were selected for the experimental studies to reduce the wrinkling tendencies due to a decrease in the otherwise unimpeded draw-in in the longitudinal direction of the strip primarily occurring during the rolling of transversal channels. The formed demonstrator plates confirm the positive effect of the restraining tension \({\sigma }_{R}\) (inlet side) on the wrinkling tendency (Fig. 12). The rolling tests were carried out without additional lubrication.

Fig. 12
figure 12

Demonstrator plates formed by hollow embossing rolling with increasing restraining tension \({\sigma }_{R}\) and decreasing wrinkling tendency

The formed demonstrator plates were scanned using a Keyence profilometer (VR-5000), with the demonstrator plates uniformly clamped in a magnetic fixture to generate a validation database. Then both the lateral strip draw-in in the \(y\)-direction and the channel depths in the sections illustrated in Fig. 11 were analyzed. Furthermore, the thickness distributions in these sections were detected by microscopy of embedded specimens at a measuring distance of 40 µm. Three separate parts were evaluated for each restraining tension \({\sigma }_{R}\) to statistically validate the results, including averaging and standard deviation analysis.

Analysis of roller elasticity to estimate rolling gap changes

As the rolling gap is smaller than the initial thickness of the strip, the roller bodies will be elastically displaced and deformed due to the increasing forming forces. A separate FE model was created in LS-DYNA to estimate the extent of forming roller elasticity. This model comprises loading the elastically deformable roller body by a rigid plate of the same size as the strip (Fig. 13 a). Here, the bearings were assumed to be ideally rigid, which provides a conservative evaluation of bending behavior. The rollers are represented as tetrahedral meshes based on the assumption of purely elastic material behavior (using *MAT_ELASTIC). As a result, the \(z\)-displacement of the roller along the strip width of 15 mm was almost uniform here (\(\Delta z\) < 1 µm, Fig. 13 b). Under these conditions, the existing bending effects were negligible, and the roller elasticity could be approximated in the model as a constant increase of the rolling gap (rolling gap change \({\Delta s}_{w}\)).

Fig. 13
figure 13

Analysis of roller elasticity: (a) FE model, (b) distribution of \(z\) displacement, (c) rolling gap change due to elasticity of the forming rollers

The maximum \(z\)-displacement was evaluated as a function of \({F}_{z}\), which is equivalent to a rolling force dependent rolling gap change \({\Delta s}_{w}\) (Fig. 13 c). This means that for an assumed rolling force of 25 kN (maximum rolling force in the initial simulation model), a rolling gap change \({\Delta s}_{w}\) of 24 µm (corresponds to 5% of the channel depth) must be expected and from a technological point of view this is a magnitude relevant for the forming simulation. Thus, the model validation must take into consideration the change in the rolling gap due to elastic roller displacement.

Model validation by comparison to experimental results with different assumed constant rolling gap changes

In order to obtain the absolute deviation from the experimental reference, the simulation results were compared to the experimental results (average values from three parts). All restraining tension variants \({\sigma }_{R}\) of the experimental test series were simulated and compared in the following analysis. The model deviation was investigated as a function of an additional rolling gap \({\Delta s}_{w}\) (assumed to be constant) as a first approximate representation of the elastic \(z\) displacement of the rollers. Thus, the existing rolling gap \({s}_{w}\) is the sum of the measured initial rolling gap \({\overline{s} }_{w0}\) (Fig. 11) and the assumed additional rolling gap change \({\Delta s}_{w}\).

Figure 14 shows the averaged result deviation of the strip thickness \(\Delta t\) in the evaluation sections (Fig. 11). Total deviation decreases with greater \({\Delta s}_{w}\) up to an average deviation of \(\Delta t\) < 4 µm (~ 4% thinning error) at \({\Delta s}_{w}\) = 10 µm. The most considerable deviations occur in the channel flanks, while the deviations in the channel head and foot areas and the flange of the part become very small, i. e. approx. 2 µm. For \({\Delta s}_{w}\) = 0 µm, the shell-based model predicts significantly too much thinning due to the tight rolling gap, especially in the channel head and foot areas (\(\Delta t\)= 8 µm). An average absolute deviation of the predicted strip thickness of 4 µm over all sections and test series with varying restraining tension \({\sigma }_{R}\) can be evaluated as a satisfactory result considering the average standard deviation of \(\overline{\sigma }\) = 2.7 µm of the experimental evaluation.

Fig. 14
figure 14

Mean thickness deviation \(\Delta t\) of simulation model in the sections for all restraining tensions: (a) model deviation as a function of an assumed constant rolling gap change \({\Delta s}_{w}\), (b) exemplary result comparison for section A-A

Evaluating the lateral edge draw-in \({y}_{E}\) shows only a slight change of the results at different \({\Delta s}_{w}\) with minimal average deviations \({\Delta y}_{E}\) for \({\Delta s}_{w}\) = 5 µm (Fig. 15 a). Thus, the average model deviation of the edge draw-in is approx. 10% with a reference value in the experiments of \({y}_{E}\) = 0.32 mm in the central area of the plate. The differences in results increase from \({\Delta s}_{w}\) > 5 µm as the simulation predicts more severe draw-in than the measured values in the experiments.

Fig. 15
figure 15

Mean edge draw-in deviation \({\Delta y}_{E}\) of simulation model for all restraining tensions: (a) model deviation as a function of an assumed constant rolling gap change \({\Delta s}_{w}\), (b) exemplary result comparison for \({\sigma }_{R}\) = 101 N/mm.2

The evaluation of the channel depth deviation \(\Delta h\) presented in Fig. 16 shows the smallest deviations for \({\Delta s}_{w}\) = 5 to 15 µm. Thus, approximate consideration of the rolling gap variation by \({\Delta s}_{w}\) again leads to significant improvement of the model accuracy. Additionally, the mean model deviation at \({\Delta s}_{w}\) = 5 µm of 4 µm is minimal considering the average experimental standard deviation of 3 µm.

Fig. 16
figure 16

Mean channel depth deviation \(\Delta h\) of the simulation model in the sections for all restraining tensions: (a) model deviation as a function of an assumed constant rolling gap change \({\Delta s}_{w}\), (b) exemplary result comparison for \({\sigma }_{R}\) = 35 N/mm.2

The assumed constant increase of the rolling gap \({\Delta s}_{w}\) analog to the elastic displacement of the rollers is associated with a significant decrease of the simulation’s rolling forces. Figure 17 a shows an example for \({\sigma }_{R}\) = 68 N/mm2, demonstrating that the rolling force is not constant due to geometrical discontinuities and that the decrease of the rolling forces varies with increasing \({\Delta s}_{w}\). On average, \({\Delta s}_{w}\) = 10 µm (~ 2% of the maximum channel depth of 0.5 mm) for the actual rolling gap \({\overline{s} }_{w0}\) of 78 µm causes a decrease in the rolling force to 45% to the reference value (\({\Delta s}_{w}\)= 0 µm).

Fig. 17
figure 17

Correlation between rolling force and rolling gap change: (a) development of rolling force with increasing rolling gap and reference to simulation with target rolling gap, (b) approximate estimation of a state of equilibrium between rolling force and elastic displacement of the rollers

The rolling forces using the original target rolling gap of \({s}_{w0}\) = 100 µm are also illustrated in Fig. 17 a for comparison. The significantly lower rolling forces clearly show the importance of considering the actual rolling gap due to manufacturing inaccuracies. Comparative simulations with \({s}_{w0}\) = 100 µm lead to considerably more significant model deviations regarding the edge draw-in (\(\Delta {y}_{E}\) = 0.12 mm) and the channel depths (\(\Delta h\) = 10 µm). Therefore, a rolling gap change from \({s}_{w0}\) = 100 µm to 78 µm is relevant concerning model deviation. When the rolling gap is increased, the gap-dependent decrease in rolling force interacts with the resulting reaction force of the roller displacement, thus leading to an approximate state of equilibrium between these variables (Fig. 17 b). Considering the average rolling forces, this state of equilibrium is established as a function of strip restraining tension \({\sigma }_{R}\) at an average rolling gap increase \({\Delta s}_{w}\) of 9 – 11 µm. Rolling gap changes can be expected in a range of 8 – 12 µm due to the variation of the rolling forces by process time. These values correspond to the minimum deviations of the result comparisons, which finally confirms the correlation between the reduced rolling gap (\({\overline{s} }_{w0}\) = 78 µm) due to the manufacturing process and the resulting elastic \(z\)-displacement of the rollers due to the rolling forces.

Introduction of an optimized model approach with spring-controlled rollers

The existing model variant with fixed rotation axes (Fig. 7, reference model "RM") cannot represent a load-dependent change of the rolling gap \({\Delta s}_{w}\) due to roller body elasticity. Based on the relevance for model accuracy proven in the previous section, the reference model was optimized by releasing the displacement degrees of freedom of the rotation axis in \(z\)-direction (abbreviated as “OM”, Fig. 18 a). Additionally implemented spring elements with a spring stiffness of 2.07 kN/µm in the restrained axis represent the elasticity of the roller bodies (reference: overall model, Fig. 13), so that a load-dependent rolling gap change \({\Delta s}_{w}\) can occur during the incremental forming (Fig. 18 b). The gap variation is limited to constant values along the strip width disregarding any potential bending effect due to rigid surface meshes. Due to the explicit forming simulation, the use of damping elements is still necessary to prevent the system from oscillating. While the lower restrained axis is assumed to be stationary, the upper restrained axis allows implementing a force- or displacement-controlled movement, which is then transmitted via the spring and damper elements to the rotation axes and thus to the rigid roller meshes. Thus, the optimized model approach enables modeling of the strip clamping to be represented by a force boundary condition (clamping force \({F}_{z,c}\)) in addition to the load-dependent rolling gap change \({\Delta s}_{w}\). In the experiment, application of a clamping force via manual movement of the top forming roller was necessary in order to apply the reaction forces of the introduced restraining tension \({\sigma }_{R}\) and thus to avoid slippage between the translational strip movement and the forming roller rotation. In this context, the transversal channels of the rollers represent a form fit supporting strip clamping and friction forces between the strip and the roller surface counteract the restraining tension as well. The flat strip areas before and after the channel zone are more critical, because here the applied reaction forces must be high enough to absorb the strip restraining tension on the inlet side and at the same time low enough to prevent any plastic deformation of the strip. Separate simulation-based investigations on equivalent models with different roller diameters indicated a linear correlation between the applied clamping force in \(z\)-direction \({F}_{z,c}\) and the transmissible restraining tension on the inlet side \({\sigma }_{R}\). Using Coulomb's law of friction, the following relationship can be derived considering contact on both sides with the top and bottom forming roller:

Fig. 18
figure 18

Optimized model approach: (a) forming simulation model with spring-controlled rollers and clamping force application, (b) analysis of rolling force and rolling gap change during hollow embossing rolling

$${\sigma }_{R}=2\bullet \frac{{F}_{z,c}}{b\bullet {t}_{0}}\bullet {\mu }_{s}$$
(3)

For an exemplary restraining tension \({\sigma }_{R}\) of 68 N/mm2, this results in a required clamping force \({F}_{z,c}\) of 0.8 kN considering a safety factor of 1.2 (\({\mu }_{s}\) = 0.15; \(b\) = 30 mm; \({t}_{0}\) = 0.1 mm). For the shell models considered here, clamping force initiation results in contact penetrations causing negative rolling gap changes \({\Delta s}_{w}\) in the range of < -2 µm (Fig. 18 b).

Starting with channel forming, the rolling force increases from the clamping force level. Due to the load-dependent rolling gap change \({\Delta s}_{w}\) of up to 13 µm, the rolling forces are on average 55% lower than in the reference model (RM) without rolling gap variation (\({\Delta s}_{w}\) = 0 µm). Thus, significantly more thinning in the channel head and foot areas and in the channel flanks is predicted at RM due to the reduced rolling gap (\({\overline{s} }_{w}\) = 78 µm). These thinning effects are significantly reduced by using the optimized model approach OM (on average \(\Delta \overline{t }\) =  + 5 µm). Consideration of the load-dependent rolling gap change \({\Delta s}_{w}\) also affects the predicted edge draw-in, which is on average 0.08 mm more in the middle demonstrator area compared to the simulation using RM. This corresponds to a relative change in draw-in results of 14% (Fig. 19 a).

Fig. 19
figure 19

Comparison of the optimized model (OM) with the reference model (RM): (a) \(y\)-displacement and edge draw-in, (b) \(x\)-displacement and longitudinal draw-in, (c) true effective plastic strain, (d) deviation from the nominal geometry after springback considering clamping

Another advantage of the optimized model approach is that the holder on the outlet side, which was moved synchronously with the roller rotation in RM, is no longer required and any possible influence of this idealized constraint on the simulation results is eliminated. Instead, as on the inlet side, a restraining tension can now be introduced on the outlet side by implementing a force boundary condition \({F}_{h}\) corresponding to the experimental test boundary conditions (Fig. 18 a). In the example shown in Fig. 19 b, the use of the constrained holder at RM already results in 0.5 mm more strip movement \(\Delta x\) in the rolling direction compared to OM. In addition, the end of the strip is less drawn in, resulting in a total difference of \({\Delta x}_{e}\) = 0.36 mm in the longitudinal draw-in between the beginning and the end of the strip. As a result, RM leads to more stretching of the strip area on the outlet side due to the constraint of the holder, including plastic deformation of the strip area before the actual forming zone (Fig. 19 c). Based on the experimental results, this effect is unrealistic and illustrates that the simplified modeling of a constrained strip holder results in a relevant influence on the simulation results.

Using the optimized modeling with load-dependent rolling gap change, significant changes in the predicted rolling forces, rolling gaps, channel depths, and strip thicknesses of the channel cross sections, as well as different strip draw-in in the transversal and longitudinal directions, have been identified in comparison with RM. In summary, these changes also affect the resulting flatness deviation in the clamped state after springback calculation (Fig. 19 d). The example shows a significant change in the absolute flatness deviation from 0.28 mm (RM) to 1.15 mm (OM), which illustrates the influence of the optimized modeling on the simulation results. In this context, the small size of the demonstrator geometry should be explicitly mentioned. Larger, standard industrial plate sizes in the range of 150 × 400 mm would result in correspondingly higher forming forces and rolling gap changes, which in turn may influence the sensitivity of the results of other evaluation criteria.

Validation of the optimized model approach

Based on the evidence that the optimized model approaches lead to significant changes in results compared to the original reference model, the question of whether the optimized model approach with load-dependent rolling gap change and additional clamping force application does in fact have an accuracy advantage compared to the reference model is investigated subsequently. Similar to the procedure described in Figs. 14, 15 and 16, simulations corresponding to the experimental conditions considering the 4 discrete restraining tensions \({\sigma }_{R}\) were performed and the deviations of the numerical results from the experimental references in terms of strip thicknesses \(\left|\Delta t\right|\) channel depths \(\left|\Delta h\right|\) and edge draw-in \(\left|\Delta {y}_{e}\right|\) were analyzed. For the quantitative evaluation of the absolute result deviations, the permissible limit values \({\Delta t}_{limit}\) = 5 µm, \({\Delta h}_{limit}\) = 10 µm, and \({\Delta y}_{e,limit}\) = 0.04 mm are introduced.

Figure 20 shows that the optimized model approach OM leads to a significant reduction of the deviation between numerical and experimental results, assuming the same friction conditions (\({\mu }_{s}\) = 0.15) as in RM. The average rolling gap change with OM is 8.3 µm, which corresponds to the approximate rolling gap change determined in Fig. 17. With OM, the mean deviations of the strip thicknesses in the sections \(\left|\Delta t\right|\) and the predicted channel depth \(\left|\Delta h\right|\) are significantly lower and within the allowed limits. As already analyzed in Fig. 15, there is no significant change in the results of the edge draw-in.

Fig. 20
figure 20

Analysis of the absolute result deviation between simulation and experiment depending on the model variants and friction coefficients

In general, the simulation results with OM still show too much edge draw-in combined with too little thinning in the channel flanks compared to the experimental data, suggesting that the influence of friction has been underestimated so far. This is further supported by the scaling effects mentioned in the state of the art, which indicate a stronger frictional influence compared to conventional sheet metal forming simulation scales. By increasing the coefficient of friction, an improved agreement of the considered evaluation criteria can be observed with OM. In particular, the changes in the edge draw-in deviations are significant. The reference model (RM) also tends to show lower mean deviations from the experimental results, but individual deviations still remain above the level of the limits established, which confirms the deficiencies of the initial model approach.

With \({\mu }_{s}\) = 0.20, the optimized model approach meets all established limits for the first time. This demonstrates the achievable accuracy of the optimized model approach with load-dependent rolling gap change and clamping force application. The smallest deviations occur with \({\mu }_{s}\) = 0.25. Even larger friction coefficients did not prove advantageous due to the overestimation of the maximum thinning effects in the channel flanks.

Summary and conclusions

Providing a sufficient data base for model validation, hollow embossing was successfully implemented on a simple strip rolling mill. In contrast to many other demonstrator variants in various publications, the bipolar plate demonstrator used here is characterized by a combination of longitudinal and transversal channels as well as several geometric discontinuities. This ensures comparability of the characteristic flow fields with common plate designs of fuel cell drives already available. The hollow embossing rolling tests prove the positive influence of a restraining tension on the inlet side in favor of lower wrinkling tendencies.

Based on quasi-static and dynamic material testing for different orientations to the rolling direction, an anisotropic material model was parameterized, and the hollow embossing rolling process was modeled using a shell-based approach overcoming the disadvantages of volume-based strip meshing regarding computation times. Comparing the simulative and experimental results demonstrates that the shell-based simulation model achieves sufficient accuracy to represent the hollow embossing rolling process. Accordingly, it is justified to apply the established shell-based modeling approaches of conventional sheet metal forming simulation with the advantage of lower computation times for technological forming simulation and process optimization. An essential requirement for model deviation emerges considering the manufacturing-related existing rolling gaps in combination with the representation of the elasticity-induced rolling gap change. An average manufacturing-related rolling gap decrease in the range of ~ 20 µm is relevant. It emphasizes the importance of high accuracy in the manufacturing process of the roller bodies.

Elastic displacement of the rollers during hollow embossing rolling is predicted in a first simulation study with assumed constant rolling gaps at a range of approx. 8 – 12 µm here (~ 10% of the initial strip thickness) due to higher rolling forces by a decreased rolling gap. In further investigations an optimized model variant with spring-controlled forming rollers was introduced representing both a load-dependent rolling gap change and the force-controlled strip clamping. Comparisons with the experimental data shows that the optimized model achieves very good accuracy. This successfully demonstrated that the hollow embossing rolling process can be modeled using shell-based surface modeling with the advantage of lower computational times. The optimized model predicts an average rolling gap change of 8 µm, thus confirming the values estimated in advance. Despite the small width of the demonstrator geometry, this magnitude is significant for result deviation based on the comparison to experimental data.

When evaluating the findings on the importance of gap accuracy and roller elasticity, the boundary conditions of this study (small demonstrator plate with narrow width, narrower rolling gap due to manufacturing) must be considered. Future investigations should implement the spring-controlled model approach at larger bipolar plate geometries validating the suitability and accuracy of the forming simulation. The optimization of the hollow embossing rolling process based on these fundamentals is the subject of further research.