In this work, a macroscopic model of a commercial lithium-ion cell in pouch format is used. The cell weighs around 0.9 kg and has dimensions of 260 × 216 × 7.8 mm (Figure 1); Kovachev et al. (2019) carried out a detailed analysis of the cell. The structure essentially consists of a stack of 42 separators, 21 cathodes and 22 anodes, which is enclosed, in the pouch. The anodes are copper foils coated on both sides with graphite. The core of the cathodes is made of aluminium and the active material on both sides is LiNiMnCoO₂ (NMC). The porous separators, positioned between the electrodes, consist mainly of polypropylene. Due to the porous structure, the separator is permeated with the electrolyte, which enables the flow of ions. The pouch, the anode, the cathode and the separator have a thickness of 190 μm, 140 μm, 170 µm and 20 μm, respectively.

This cell was modelled in the software LS-DYNA version 9.3.1. The basic concept is shown in Figure 2A. It is based on the model described by Raffler et al. (2017) where beam elements also play a central role (Raffler et al., 2017). However, instead of a cylindrical cell, a cell in pouch format is modelled here. The basic concept can be divided into four essential components, which are highlighted in Figure 2A. The model was calibrated to fit the experimental data described by Raffler et al. (2022); Schmid et al. (2022). The macroscopic model is shown in Figure 2B.

The outer layers ① represent the pouch shell. This is modelled by shell elements. An isotropic elastic-plastic material model (*MAT_PLASTICITY_COMPRESSION_TENSION) is used for this. The required parameters can be determined by tensile tests of samples of the pouch. As with the other components, samples of 15 × 5 mm are tested at 20 and 600 mm/min. The samples are taken in different directions to assess the anisotropy in addition to the strain-rate influence. All samples were saturated with substitute electrolyte to test them as realistic as possible. For statistical validation, each test configuration was repeated 5 times. In order to increase flexibility during calibration, a non-symmetrical material behaviour is used (different tensile and compressive behaviour). Thus, for example, the buckling of the shell under compressive load can also be represented in a simple way. The parameters for the compression behaviour are calibrated by a quasi-static three-point bending test of the cell (Figure 3A). The velocity of the impactor is 1 mm/s, as in all quasi-static cell tests.

The quasi-static behaviour in thickness direction is reflected by beam elements ②. To avoid decreasing time step (explicit simulation) upon transverse compression, a discrete element formulation is used. The required load curves for the non-linear behaviour are calibrated by a quasi-static indentation test with cylindrical impactor (Figure 3B).

The middle layer, employing shell elements, ③ replicates the in-plane behaviour of the jelly stack. The same isotropic elastic-plastic material model was used for this as for the outer layers. For the parameterisation of the tensile behaviour, the results of the tensile tests of the individual components are used. According to Eq. 1, the total stress σ ε M L results from the sum of the three components ( C = 3 ), taking into account the cross-section of the middle layer S M L . S c is the cross-section of the components, n c the number of components in the jelly stack and σ ε c the tensile behaviour of the corresponding component. σ ε M L S M L = ∑ c = 1 C S c n c σ ε c (1)

The last characteristic component of this modelling approach are the solid elements ④, which are used to replicate strain-rate dependency and incompressibility of the electrolyte. Cells show a considerable strain-rate dependency, both w.r.t. Initial stiffness and displacement to failure, particularly upon out-of-plane compression (Tancogne-Dejean et al., 2022). The solid elements provide the dynamic over-stress, i.e., the difference between the dynamic and quasi-static behaviour. For this purpose, a strain-rate dependent honeycomb material model (*MAT_MODIFIED_HONEYCOMB) is used. The same test configuration as shown in Figure 3B is used for the calibration. However, the impactor velocity is 3 m/s, which is significantly higher than the quasi-static level (1 mm/s). (Remark: Clearly, this could have been modelled through rate-dependent dampening of the beam elements, too. However, it turned out, that the stable time-step chosen by the solver drops massively once dampening is enabled in the non-linear 6 DoF beam material model).

To summarize, the macroscopic model is calibrated through tensile tests of the four components anode, cathode, separator and pouch. In addition, the data from two quasi-static and one dynamic cell tests are used to calibrate the model.

The quasi-static experiments were performed with a hydraulic press. The system had a maximum force of 420 kN and allowed displacements of up to 400 mm. The possible load velocities were between 0.5 and 6 mm/s. To measure the indentation force, a load cell type K 500 kN from GTM in Bickenbach, Germany, was utilized. For displacement measurement, a linear glass scale with a 1 µm resolution and an accuracy of ± 0.01 mm/m was used. The dynamic experiments were carried out with the sled test rig. Ellersdorfer et al. described its construction and function in detail (Ellersdorfer et al., 2023).

For the prediction of the time history curves, as well as for the contours, the Proper Generalized Decomposition–PGD–method is used (Chinesta et al., 2013), in particular two of its non-intrusive variants, the sparse PGD–sPGD–and the ANOVA-PGD, both presented in (Sancarlos et al., 2021). This choice is mostly motivated by two reasons. The first reason is the capability of these regression tools to address the high-dimensionality (i.e., a high number of parameters involved in the model) at the low-data limit (i.e., a small number of data points to build the model). Indeed, when dealing with complex models and problems, sample data points (numerical simulations) may be expensive to gather, so that we would like to reduce their number as much as possible. Second, non-intrusive PGD techniques have demonstrated good efficiency in industrial settings involving various engineering applications (Champaney et al., 2022a; el FallakiIdrissi et al., 2022; Champaney et al., 2022b; Chinesta et al., 2022; Champaney et al., 2021; Pasquale et al., 2022). Such techniques are developed by ESI Group Chair teams at ENSAM ParisTech and industrialized by the ESI Group Research and Innovation teams (Scientific Department) through the Hybrid-Twin AdMoRe platform, whose workflow and roadmap can be found in (Champaney et al., 2021; Chinesta et al., 2022).

The non-intrusive PGD techniques (Sancarlos et al., 2021; Champaney et al., 2022a; Pasquale et al., 2022) are nonlinear regression methods inspired by the separation of variables, which is at the heart of the PGD (Chinesta et al., 2013). This means that a generic function of space x , time t and N p parameters μ is approximated as a sum of M products of functions depending only on space X , time T and parameters M , respectively (Eq. 2). u x , t , μ 1 , μ 2 , μ 3 , … , μ N p ≈ ∑ i = 1 M T i t X i x ∏ j = 1 Q = N p M i j μ j (2)

Even though in this work we are interested in the prediction of curves and contours (vector-valued functions), for the sake of simplicity, we will present the algorithm in the simplest form of a scalar quantity prediction. The reader will find all the generalizations concerning vector-valued quantities in (Champaney et al., 2022a).

Moreover, it shall be noticed that in this work we aim at building a parametric model, which also involves geometric parameters, such as the diameter of the impactor and its orientation. This feature is handled directly by AdMoRe platform. To be able to interpolate the snapshots, all the contours are mapped onto a common reference mesh where the regression is performed (usually the one corresponding to the last run of the DoE). This mapping is executed as a data pre-processing prior to the execution of all the Model Order Reduction algorithms.

An indentation test with a cylindrical impactor, as shown in Figure 3B, is used as an application. Here, however, a dynamic case is examined. Four different parameters are varied. All four are related to the impactor but can be divided into geometric and non-geometric. The geometric parameters include the diameter and orientation of the impactor. For orientation, the impactor is positioned over the centre of the cell. The rotation occurs around the w-axis. The initial position refers to the configuration shown in Figure 3B. Because of the change in diameter, the position in the w-direction is adjusted such that there is always an initial distance of 0.2 mm between the cell and the impactor. The non-geometric parameters are the initial velocity and the mass of the impactor. The value ranges of the individual parameters are listed in Table 1.

To generate a data basis, 100 simulations with different parameter combinations are carried out. The corresponding combinations were determined using Latin Hypercube Sampling (LHS) for the first regression model based on the sPGD. This sampling method is an alternative to the classical random distribution (Monte Carlo Strategy). For the second regression method (ANOVA-PGD) another enhanced LHS sampling strategy is used. For each parameter, a nominal, which corresponds to mean value, is used. Those nominal values of all parameters are listed in Table 1. For the first run, all the parameters are set to the nominal value. The following runs are then defined by changing only one parameter while the others remain at the nominal value, until all the line segments around the nominal point are discretized, creating a sort of multi-dimensional cross in the parametric domain. These runs will be used by the model to study the effect of each parameter on the output independently of the others. Additional runs are generated from an LHS to allow the method to capture interactions between the parameters.

All runs are carried out with the same simulation duration. From all simulations, both the visual result in the form of the contour plot and the time history curves are available. In addition to impactor force over penetration depth, the kinetic energy of the impactor over time is also considered. For the prediction of the force-displacement curve the sPGD method is used. The kinetic energy over time as well as the contour plot is predicted using ANOVA-PGD.

Finally, the data basis is divided into training data and test data in a ratio of 80/20. The training data is used to create the meta-model. The quality can then be checked with the help of the test data, as these are not used for training.

Before the data can be used for the training of the meta-model, some physics-informed data pre-processing can help to improve regression results and is sometimes necessary. This mainly concerns the regression on time history curves since those data do not have the same end-points or may exhibit clearly some localizations (Champaney et al., 2022a). An example of the former are the force-displacement curves, while the latter is observed in the case of the kinetic energy curves due to the classical evolution of this quantity (localization at shifted minimum points).

A curve alignment and dilatation strategy based on Gappy POD is proposed for the pre-processing of the force-displacement curves, while the kinetic energy evolution curves are aligned with respect to the minimum point, as suggested in (Champaney et al., 2022a). Since the maximum points of the force-displacement curves are not uniform, these curves have different lengths. To achieve a better prediction, all curves are extended to the range d max . This parameter is the maximum value of all m simulations with respect to their penetrations d (Eq. 12). d max = max m d m (12)

Thus, a prediction of two aspects with respect to the force-displacement curve is required: the path of the curve and its end point. To extend all curves, missing data points must be added. These points should be chosen in such a way that the course fit the overall characteristic. The Gappy POD method is used for this.

This method is based on the Proper Orthogonal Decomposition approach. The approach assumes that a high-dimensional system can be approximated by a linear combination of low-dimensional basis functions (eigenmodes). This does not require a-priori knowledge of the system itself, as it is a method of data analysis. This is shown in Eq. 13, where f x , t is the collected data and f ∼ x , t its approximation. The data was measured at N points. The approximation is calculated as a finite sum of the temporal coefficients α i and the spatial modes φ i   x . Therefore, only those P eigenmodes with the biggest energetic proportion are used (Liang et al., 2002). f x , t k ≈ f ∼ x , t k = ∑ i = 1 P α i k φ i   x   k ∈ 1 , N (13)

There are different methods for calculating the coefficients and modes. The snapshot method is very well known (Sirovich, 1987). As shown in Eq. 14, the spatial modes φ i   x consist of a linear combination of the temporal coefficients α i and the snapshots f x , t k . φ i   x = ∑ i = 1 N α i k f x , t k (14)

Thus, the temporal coefficients α i are the last remaining unknowns. These can be calculated by solving the eigenvalue problem shown in Eq. 15 Therein K is the correlation matrix according to Eq. 16 and Λ is a diagonal matrix consisting of the eigenvalues λ which are ordered in descending order of energy content. K α = Λ α (15) K = f   f T N (16)

The Gappy Proper Orthogonal Decomposition (GPOD) is an extension of the classical POD method. It is used for data preparation. Data gaps can be filled, incorrect entries can be replaced and data sets can be completed. Several GPOD methods are well known in literature, such as the Everson Sirovich Method (Everson and Sirovich, 1995), Gunes Method (Gunes et al., 2006) or Field Smoothness Method (Raben et al., 2012). In this work, the Gunes Method is used to bring all force-displacement curves to the length d max .

This is an addition of missing data points, which is done iteratively in several steps. In the first step of this method, all missing entries are replaced by an initial guess. This is shown in Eq. 17. If there is a valid entry ( x ∈ x a ) for any value t k it is taken over. Otherwise ( x ∈ x m ), the mean value f x of all data records at this value t k is taken as the initial value. f ¯ x , t k n = 0 = f x , t k   x ∈ x a f x   x ∈ x m (17)

In the next step, a POD analysis is performed on this prepared data set, see Eq. 18 By solving the problem as described in Eqs 14–16, a new approximated data set f ¯ ∼ x , t k n is computed. f ¯ ∼ x , t k n = ∑ i = 1 P α i k φ i   x (18)

The newly calculated values are now taken over for those places where the data set is to be extended ( x ∈ x m ), see Eq. 19 This is now the starting point for the next iteration n + 1 . f ¯ x , t k n = f x , t k   x ∈ x a f ¯ ∼ x , t k n   x ∈ x m (19)

This routine is carried out until convergence occurs. To assess this, the parameter ε is introduced, see Eq. 20 This describes the change between the individual iteration steps. The Frobenius norm is used for this. The iteration is carried out until the value falls below a limit value. ε = ( ∑ f ¯ x , t n − f ¯ x , t n − 1 2 ) 1 / 2 (20)

This procedure is used to bring all n force-displacement curves to the same final length. This is necessary because they have been capped at their respective maximum penetration point, once kinetic energy is zero and has been fully converted to internal energy. For the kinetic energy over time this is not needed.

For anode, cathode and pouch, both anisotropy and strain-rate dependency were extremely low. As a result, both effects were neglected in the modelling for the components mentioned. Accordingly, only the mean values of all curves are shown for the results in Figures 4A–C.

With separators, both anisotropy and strain-rate effect cannot be neglected. The corresponding mean value curves are shown in Figure 4D. However, if Eq. 1 is considered and the ratio of the component thicknesses is added, it is noticeable that the separator has a minor influence on the behaviour of the entire jelly stack. This justifies the assumption of an isotropic material for the middle layer.

In addition to the component data, data from cell experiments are required for calibration. The comparison between the experiments and the simulation of the three-point bending test is shown in Figure 4E. The maximum impactor travel thereby is 30 mm. The resulting force is about 1.2 kN.

The comparison of experiment and simulation of the indentation test with cylindrical impactor is shown for both the dynamic (3 m/s) and the quasi-static (1 mm/s) case in Figure 4F. The maximum force in the quasi-static test is around 220 kN. The impactor penetrates the cell by about 2.4 mm. In the dynamic case, a penetration of about 1.6 mm is already sufficient to cause the cell to fail.

To create the database, 100 simulations were carried out with different parameter combinations. The results are shown in Figure 5. The left subplot shows the contact force between impactor and cell over the penetration depth. The maximum force reached around 340 kN. The respective parameter combination also gives the largest penetration depth of around 2.8 mm. It can be seen that all force-displacement curves have a similar characteristic shape. The force increases up to the maximum and then the contact force decreases.

The right subplot of the figure shows the kinetic energy of the impactor over time. The time range that was simulated was the same for all 100 runs. The termination time was 5 m. However, the initial value varies between the individual parameter combinations. It lies in the range of 4.6 and 465 kNmm. Furthermore, it can be seen that the assumption is confirmed that all curves with a minimum value have the same characteristic shape. Please note that different sampling strategies were used for the results shown in Figure 5.

In order to improve the meta-model quality, as suggested in (Champaney et al., 2022a), a suitable data pre-processing based on curves alignment was performed. The force-displacement curves in Figure 5A were first clipped at the maximum point as shown in Figure 6A. Then, within the data set, the curve showing maximum displacement (abscissa) is taken as reference and all the others are filled up to this displacement using the Gappy POD strategy. In such a way, we obtain the aligned curves in Figure 6C. A value of 10 − 2 was chosen as the limit value for the convergence control in Gappy POD algorithm, see Eq. 20 In Figure 6D it is clearly visible that the algorithm leads to convergence. To reach the limit, 609 iterations were required.

For what concerns the kinetic energy curves in Figure 5B, the alignment and expansion with respect to the minimum point was performed, following the procedure explained in (Champaney et al., 2022a). This leads to the curves of Figure 6B.

Of course, after the regression is performed on the pre-processed data, the original curves are reconstructed through inverse mapping (Champaney et al., 2022a).

The database consists of 100 runs for both methods, sPGD and ANOVA-PGD. Of these, 80 are used for the training phase. The remaining 20 runs are available for validation. Here, five comparisons between HD model and meta-model are discussed as examples. This is the case for the prediction of the force-displacement curve, the kinetic energy over time and the contour plot. These combinations were picked randomly out of the validation data of the sPGD DoE. The corresponding parameter combinations are listed in Table 2.