A Material Model for the Orthotropic and Viscous Behavior of Separators in Lithium-Ion Batteries under High Mechanical Loads

: The present study is focused on the development of a material model where the orthotropic-visco-elastic and orthotropic-visco-plastic mechanical behavior of a polymeric material is considered. The increasing need to reduce the climate-damaging exhaust gases in the automotive industry leads to an increasing usage of electric powered drive systems using Lithium-ion (Li-ion) batteries. For the safety and crashworthiness investigations, a deeper understanding of the mechanical behavior under high and dynamic loads is needed. In order to prevent internal short circuits and thermal runaways within a Li-ion battery, the separator plays a crucial role. Based on results of material tests, a novel material model for ﬁnite element analysis (FEA) is developed using the explicit solver Altair Radioss. Based on this model, the visco-elastic-orthotropic, as well as the visco-plastic-orthotropic, behavior until failure can be modeled. Finally, a FE simulation model of the separator material is performed, using the results of different tensile tests conducted at three different velocities, 0.1 mm · s − 1 , 1.0 mm · s − 1 and 10.0 mm · s − 1 and different orientations of the specimen. The purpose is to predict the anisotropic, rate-dependent stiffness behavior of separator materials in order to improve FE simulations of the mechanical behavior of batteries and therefore reduce the development time of electrically powered vehicles and consumer goods. The present novel material model in combination with a well-suited failure criterion, which considers the different states of stress and anisotropic-visco-dependent failure limits, can be applied for crashworthiness FE analysis. The model succeeded in predicting anisotropic, visco-elastic orthotropic and visco-plastic orthotropic stiffness behavior up to failure.


Introduction
Increasing efficiency, defined as energy consumption relative to the distance travelled, is considered a core development goal of the automotive industry. Customer demands for vehicles with larger mobility range in combination with legal requirements for climate protection by limiting CO 2 emissions are the motivation for increasing the efficiency of electric vehicles. In addition to the ecological motivation, economic factors are a major incentive for vehicle manufacturers to focus on the reduction in the development time of their vehicles by using high fidelity computational simulation tools such as finite element analysis. Product development to raise battery performance of portable electronic devices such as laptops or mobile phones also benefits from an increased accuracy of FEA results [1,2]. To achieve this goal, more precise FE models are required to describe the mechanical behavior of all involved materials under high mechanical loads. Recent developments, especially in the

Separator Materials and Methods
The investigated commercial polymeric separator was provided to us from a battery manufacturing company. The exact chemical composition of the materials used in these samples was not specified. The separators had the porosity of 36-46% [17].
In this study, the mechanical properties of the specimen were studied by conducting tensile tests on small strips cut out in 2 perpendicular spatial directions corresponding to machined direction (MD) and transversal direction (TD) of the separator samples. Figure  2 shows the experimental setup which was used for performing the tensile tests. Each specimen with dimensions of 12 × 60 × 0.025 mm was glued at each end to metal specimen holders which led to a free test length of 36 mm. The specimens were marked with a gray pattern to enable digital image correlation (DIC), as shown in Figure 3a,b. The 3D ARAMIS system from GOM (GOM GmbH, Schmitzstraße 2, 38122 Braunschweig, Germany) was used for analysis of the strain field, shown in Figure 3c,d. Each test was repeated 5 times to ensure repeatability. Three different tension velocities were used for testing: 0.1 mm·s −1 , 1 mm·s −1 and 10 mm·s −1 corresponding to strain rates 0.002778 s −1 , 0.02778 s −1 , and 0.2778 s −1 .

Separator Materials and Methods
The investigated commercial polymeric separator was provided to us from a battery manufacturing company. The exact chemical composition of the materials used in these samples was not specified. The separators had the porosity of 36-46% [17].
In this study, the mechanical properties of the specimen were studied by conducting tensile tests on small strips cut out in 2 perpendicular spatial directions corresponding to machined direction (MD) and transversal direction (TD) of the separator samples. Figure 2 shows the experimental setup which was used for performing the tensile tests. Each specimen with dimensions of 12 × 60 × 0.025 mm was glued at each end to metal specimen holders which led to a free test length of 36 mm. The specimens were marked with a gray pattern to enable digital image correlation (DIC), as shown in Figure 3a,b. The 3D ARAMIS system from GOM (GOM GmbH, Schmitzstraße 2, 38122 Braunschweig, Germany) was used for analysis of the strain field, shown in Figure 3c,d. Each test was repeated 5 times to ensure repeatability. Three different tension velocities were used for testing: 0.1 mm·s −1 , 1 mm·s −1 and 10 mm·s −1 corresponding to strain rates 0.002778 s −1 , 0.02778 s −1 , and 0.2778 s −1 . Figure 3 shows one specimen before and after testing. The recording frequency for the force and displacement measurements ranged from 2 Hz for the slowest tension velocity up to 40 Hz for the fastest ones.
The measured elongation in MD direction is in the range of 10-19 mm, which results in 27-53% total strain. Figure 4 shows the engineering stress vs. engineering strain curves for the three tested velocities in MD. The Young's moduli are computed using the initial slope. Within this work, the statistical R-squared value (coefficient of determination) is used for the fitted regression line, starting with the first 2 measuring points. When the adjusted R-squared value drops below a defined threshold, we use the previous fit line as the Young's modulus. After failure, the specimen retained almost all its extension, which indicates a low relaxation. In conclusion, the material exhibits high plasticity and low elasticity behavior. Localization of strain occurs close to the center and proceeds to one of the fixed ends, where final rupture occurs.  Figure 3 shows one specimen before and after testing. The recording frequency for the force and displacement measurements ranged from 2 Hz for the slowest tension velocity up to 40 Hz for the fastest ones.  The measured elongation in MD direction is in the range of 10-19 mm, which results in 27-53% total strain. Figure 4 shows the engineering stress vs. engineering strain curves for the three tested velocities in MD. The Young's moduli are computed using the initial slope. Within this work, the statistical R-squared value (coefficient of determination) is used for the fitted regression line, starting with the first 2 measuring points. When the adjusted R-squared value drops below a defined threshold, we use the previous fit line as the Young's modulus. After failure, the specimen retained almost all its extension, which indicates a low relaxation. In conclusion, the material exhibits high plasticity and low Figure 3. Set-up of prepared specimens in MD before the test (a) and at the end of the test, after failure (b). Specimen using DIC measurement at the beginning of the test (c) and at the end of the test, briefly before failure (d), showing the inhomogeneous strain field distribution withing the specimen.  The measured elongation in TD direction was in the range of 40-53 mm, which results in 111-147% total strain. Figure 5 shows the engineering stress vs. strain curves for all three tested velocities in TD. The Young's moduli are computed using the initial slope as shown in Figure 6, which shows a cutout (indicated with the red-dotted rectangle in Figure 5).  The measured elongation in TD direction was in the range of 40-53 mm, which results in 111-147% total strain. Figure 5 shows the engineering stress vs. strain curves for all three tested velocities in TD. The Young's moduli are computed using the initial slope as shown in Figure 6, which shows a cutout (indicated with the red-dotted rectangle in Figure 5).  The measured elongation in TD direction was in the range of 40-53 mm, which results in 111-147% total strain. Figure 5 shows the engineering stress vs. strain curves for all three tested velocities in TD. The Young's moduli are computed using the initial slope as shown in Figure 6, which shows a cutout (indicated with the red-dotted rectangle in Figure 5).    Table 1 shows the estimated Young's moduli over the nominal strain rates (based on the used tension velocities). The values for the diagonal direction DD (= 45°) were obtained from Sahraei [17].

Viscosity in Elastic Region
Hereinafter the denotation of a matrix is defined as ∎ , and of a vector as ∎ .
The viscosity effect is easily treated by explicit solvers because of the usage the equation of motion and solving for unknown values without the need of inverting the entire model:

Viscosity in Elastic Region
Hereinafter the denotation of a matrix is defined as [ ], and of a vector as { }. The viscosity effect is easily treated by explicit solvers because of the usage the equation of motion and solving for unknown values without the need of inverting the entire model: [M] .. where: [M] = Mass matrix (kg), [C] = Damping matrix (kg·m·s −1 ), ..
x n = Velocity vector (m·s −1 ), Due to the usage of the small strain formulation, where the stress increment is integrated at each time step and considering a stiffness that depends on the local strain rate in individual time steps, the integration schema-which solves the equation of motion in the time domain-leads to high oscillations of the total strain rates. Figure 7a shows considerable strain rate oscillations within the specimen during loading in the elastic domain. To reduce these oscillations, several approaches are used by commercial solvers [18][19][20].

Plasticity
In the novel developed visco-elastic-plastic model, the true-stress vs. true-plas strain curves must be derived from the tested engineering-stress vs. engineering-str curves. For the calculation of the plasticity curves in the MD direction, the approach ba on Swift-Voce interpolation, described by Bulla [16], is used. However, the TD behav is showing a much higher plasticity range beyond 100% strain. Figure 8 shows o representative engineering stress vs. engineering strain curve close to 120% total str before the material fractures (green curve). For usage in an elasto-plastic material mod the parametrization for the engineering stress vs. engineering strain curve must modified further. The true stress vs. true strain curve ( Figure 8: red curve) is derived a then the elastic portion of the strain is subtracted from the total strain, using the lin strain-stress relation corresponding to the strain rate, dependent Young's modu ( Figure 6: dashed lines).
Due to the need for values beyond the curve, derived from testing, it is necessary extrapolate the true stress vs. true plastic strain curve. Based on the work of G'Sell et [22] who tested such high ductile polymer materials, the parameter identification carried out using the modified equation:  The approach chosen here is to apply an exponential moving average filter [21] to the strain rate as a time series filter that removes high frequency oscillations: where: ε n = Actual (not filtered) strain rate (s −1 ), α EMA = Exponential moving filter value (−).
With the α EMA parameter it is possible to control the response of the filter. If α EMA = 0, no filtering is applied. Values close to 1 apply a strong low-pass filter. Figure 7b shows the result of applying the strain rate filtering. It leads to a smoother strain rate distribution and more realistic behavior within the FE model. The strain rate filtering plays an important role in explicit time integration.

Plasticity
In the novel developed visco-elastic-plastic model, the true-stress vs. true-plasticstrain curves must be derived from the tested engineering-stress vs. engineering-strain curves. For the calculation of the plasticity curves in the MD direction, the approach based on Swift-Voce interpolation, described by Bulla [16], is used. However, the TD behavior is showing a much higher plasticity range beyond 100% strain. Figure 8 shows one representative engineering stress vs. engineering strain curve close to 120% total strain before the material fractures (green curve). For usage in an elasto-plastic material model, the parametrization for the engineering stress vs. engineering strain curve must be modified further. The true stress vs. true strain curve (Figure 8: red curve) is derived and then the elastic portion of the strain is subtracted from the total strain, using the linear strain-stress relation corresponding to the strain rate, dependent Young's modulus ( Figure 6: dashed lines).
= hardening coefficient (MPa); = hardening plastic strain coefficient (−) = 2nd hardening plastic strain coefficient (−); = 3rd hardening plastic strain coefficient (−); . = equivalent plastic strain (−). Table 2 shows the identified parameter used for the hardening curve, based on G'Sell equation.  Figure 9 shows very good agreement between the true stress vs. true plastic strain curves obtained from real test data (black curve) and the calculated hardening curve based on G'Sell's equation.  Due to the need for values beyond the curve, derived from testing, it is necessary to extrapolate the true stress vs. true plastic strain curve. Based on the work of G'Sell et al. [22] who tested such high ductile polymer materials, the parameter identification is carried out using the modified equation: where: A = initial yield stress parameter (MPa); B = hardening coefficient (MPa); C = hardening plastic strain coefficient (−) D = 2nd hardening plastic strain coefficient (−); F = 3rd hardening plastic strain coefficient (−); ε pl. = equivalent plastic strain (−). Table 2 shows the identified parameter used for the hardening curve, based on G'Sell equation.  Figure 9 shows very good agreement between the true stress vs. true plastic strain curves obtained from real test data (black curve) and the calculated hardening curve based on G'Sell's equation. Within the further FEA analysis, the G'Sell-based hardening curve will be used for the TD behavior. Figure 10 shows parametrized hardening curves resulting from test data for different strain rates.

Stability Investigation
In the explicit FEA, the numerical stability plays an important role, since the numerical method is unstable if the timestep is above a critical timestep [23,24]. Another instability results from the material behavior. This material instability is investigated further since the tested separator material shows a large strain in the TD direction.
In real tests with metallic materials, the engineering stress vs. engineering strain curve shows monotonic increase in stress with the increase in strain, starting with the slope of the Young's modulus. When approaching plastic deformation, the stress still increases, but with decreasing slope, until a maximum stress value is reached. Then, the measured stress starts to decrease. At this point of maximal stress, localization occurs, where a certain area continues plastic deformation without an increase in tension force. Within the further FEA analysis, the G'Sell-based hardening curve will be used for the TD behavior. Figure 10 shows parametrized hardening curves resulting from test data for different strain rates. Within the further FEA analysis, the G'Sell-based hardening curve will be used for the TD behavior. Figure 10 shows parametrized hardening curves resulting from test data for different strain rates.

Stability Investigation
In the explicit FEA, the numerical stability plays an important role, since the numerical method is unstable if the timestep is above a critical timestep [23,24]. Another instability results from the material behavior. This material instability is investigated further since the tested separator material shows a large strain in the TD direction.
In real tests with metallic materials, the engineering stress vs. engineering strain curve shows monotonic increase in stress with the increase in strain, starting with the slope of the Young's modulus. When approaching plastic deformation, the stress still increases, but with decreasing slope, until a maximum stress value is reached. Then, the measured stress starts to decrease. At this point of maximal stress, localization occurs, where a certain area continues plastic deformation without an increase in tension force.

Stability Investigation
In the explicit FEA, the numerical stability plays an important role, since the numerical method is unstable if the timestep is above a critical timestep [23,24]. Another instability results from the material behavior. This material instability is investigated further since the tested separator material shows a large strain in the TD direction.
In real tests with metallic materials, the engineering stress vs. engineering strain curve shows monotonic increase in stress with the increase in strain, starting with the slope of the Young's modulus. When approaching plastic deformation, the stress still increases, but with decreasing slope, until a maximum stress value is reached. Then, the measured stress starts to decrease. At this point of maximal stress, localization occurs, where a certain area continues plastic deformation without an increase in tension force. This important point is so called the necking point. For metals, the material instability is reached, and the material starts to damage and soften until final separation and failure will occur.
In the MD, the specimens reach their terminal strain shortly after the maximum stress is reached. With a well calibrated failure criterium, the finite elements, where this point is reached, will be deleted from further calculation.
In the TD of the tested material, the strain increases significantly after reaching the maximum stress value, without fracture. Therefore, the stability will be investigated more thoroughly. Inspired by the Drucker stability condition [25], which is widely used in hyper-elastic material models [26], the slope of the stress vs. strain curve will be analyzed. Using the true strain: where: ε eng = engineering strain (−); ε true = true strain (−); L 0 = initial length (mm); L 1 = final length (mm); L = measured Length (mm); ∆L = relative displacement (mm) of the parallel length and the true stress: where: σ true = true stress (MPa); F = measured force (N); S 0 = initial cross section (mm 2 ); S = actual cross section (mm 2 ) of the parallel length.
The assumption of constant volume is made due to the lack of measurement of volume change, because of the small thickness of the separator material. The extension of this work, including consideration of the volume change during the test, is the topic of ongoing work and will be the subject of future investigation. Assuming a stable behavior of a material with a positive slope in the engineering stress-strain relationship and isochoric behavior in plastic region, the limit is reached where the slope becomes zero: This indicates the necking point. Proceeding with the slope in the true stress vs. true strain domain, used in the elasto-plastic material models: dσ true dε true = dσ eng · 1 + ε eng + σ eng ·dε eng 1 (1+εeng) ·dε eng = 1 + ε eng 2 ·dσ eng dε eng + 1 + ε eng ·σ eng ·dε eng dε eng , we obtain as a condition for the necking point: dσ true dε true = σ true (8) Figure 11 shows the true stress vs. true plastic strain and the first derivative in the same chart. According to Equation (8), the material instability occurs when the 1st derivative of the true stress vs. true strain curve equals the true stress value, which is the first crossing point. This corresponds to the maximum stress in the engineering stress vs. engineering strain diagram.
Energies 2021, 14, x FOR PEER REVIEW 12 of 18 Figure 11. Extrapolated, fitted hardening curve in TD (green curve) and its 1st derivative, indicating the material stability range (red curve).
From this point, the first derivative is below the true stress vs. true plastic strain curve. This indicates an unstable material behavior. Proceeding further with the elongation up to a strain value of approximately 0.6, the derivative crosses again the true stress vs. true strain curve and continues above this curve. This indicates a stabilization of the material due to potential orientation of the molecular chains in the loading direction.
A similar effect is observed in real tests in TD direction. Figure 12 sows two examples in TD, showing the localization through the rest of the specimen (a) or a localization in a certain area with a second localization occurring while proceeding with the tension test (b).

Material Modeling
In a recent work, a nonlinear visco-elastic-orthotropic and visco-plastic-orthotropic material model was developed, based on the investigation of the PE material used as separator in Li-Ion batteries. In recent publications and commercial FEA products and solvers, there exists several material models to model the behavior of polymers [27][28][29]. From this point, the first derivative is below the true stress vs. true plastic strain curve. This indicates an unstable material behavior. Proceeding further with the elongation up to a strain value of approximately 0.6, the derivative crosses again the true stress vs. true strain curve and continues above this curve. This indicates a stabilization of the material due to potential orientation of the molecular chains in the loading direction.
A similar effect is observed in real tests in TD direction. Figure 12 sows two examples in TD, showing the localization through the rest of the specimen (a) or a localization in a certain area with a second localization occurring while proceeding with the tension test (b).
Energies 2021, 14, x FOR PEER REVIEW 12 of 18 Figure 11. Extrapolated, fitted hardening curve in TD (green curve) and its 1st derivative, indicating the material stability range (red curve).
From this point, the first derivative is below the true stress vs. true plastic strain curve. This indicates an unstable material behavior. Proceeding further with the elongation up to a strain value of approximately 0.6, the derivative crosses again the true stress vs. true strain curve and continues above this curve. This indicates a stabilization of the material due to potential orientation of the molecular chains in the loading direction.
A similar effect is observed in real tests in TD direction. Figure 12 sows two examples in TD, showing the localization through the rest of the specimen (a) or a localization in a certain area with a second localization occurring while proceeding with the tension test (b).

Material Modeling
In a recent work, a nonlinear visco-elastic-orthotropic and visco-plastic-orthotropic material model was developed, based on the investigation of the PE material used as separator in Li-Ion batteries. In recent publications and commercial FEA products and solvers, there exists several material models to model the behavior of polymers [27][28][29].

Material Modeling
In a recent work, a nonlinear visco-elastic-orthotropic and visco-plastic-orthotropic material model was developed, based on the investigation of the PE material used as separator in Li-Ion batteries. In recent publications and commercial FEA products and solvers, there exists several material models to model the behavior of polymers [27][28][29].
Within the nonlinear and crashworthiness FEA models, used by the automotive and other industries, the by far most used material model is based on elasto-visco-plasticity, using the Von Mises plasticity formulation [30]. This is due to the fact that the user can easily obtain the necessary values to fit the material model by using results from various performed tensile tests. Another very important advantage of this model is based on the efficiency and performance in huge FE models, which lead to its wide usage in industrial environments [9]. Therefore, the Von Mises plasticity model was chosen as a basis and enhancements were implemented to handle the visco-elastic and orthotropic behavior of the PE separator.
For the nonlinear interpolation, between the three orthotropic directions, the cosineinterpolation function is used, as described by Bulla [16].
With the assumption of no sudden changes in the orthotropy of neighboring elements, the same method is applied to the Von Mises plastic potential: where: f = Flow potential (MPa); J 2 = deviatoric stress tensor (MPa); k = yield stress obtained from the uniaxial tension test (MPa).
Using a similar approach for viscosity and orthotropy, as used in the elastic region, the yield stress k is dependent on the filtered strain rate and orthotropy direction measured in the real test. With the assumption of no sudden changes in the orthotropy, the associated flow rule is used within this new material model.

FE-Model
In present work the new developed material model is used and validated by the tested results. The orthotropic visco-elastic properties, obtained from test results (see Table 1) are used, to model the orthotropic-visco-elasticity. During the loading phase, a precise modeling of the elastic properties is important to accurately predict the stress response of the entire model since all elements are participating in this mechanical region. Figure 13 shows the three visco-plastic hardening curves for the MD direction obtained by the test carried out during this study and the DD direction result curves taken from Sahraei [17]. The engineering stress vs. engineering strain curves in MD direction are fitted to the Swift and Voce equations [31] and extrapolated to 100 % strain, using an interpolation function as described in Bulla [16]. Figure 14 shows the FE results of the material behavior in MD direction. It can be seen that after the first localization occurs, the localized area remains at the same group of elements. These elements are further extending whereas the neighboring elements remain at the strain level, they had before the localization started. For these material curves (shown in Figure 13, left side) the slope of the true stress vs. true plastic strain curve is decreasing, which means, the material will not become stable back again. where: = 1st invariant on the stress tensor (MPa). Using a similar approach for viscosity and orthotropy, as used in the elastic region, the yield stress k is dependent on the filtered strain rate and orthotropy direction measured in the real test. With the assumption of no sudden changes in the orthotropy, the associated flow rule is used within this new material model.

FE-Model
In present work the new developed material model is used and validated by the tested results. The orthotropic visco-elastic properties, obtained from test results (see Table 1) are used, to model the orthotropic-visco-elasticity. During the loading phase, a precise modeling of the elastic properties is important to accurately predict the stress response of the entire model since all elements are participating in this mechanical region. Figure 13 shows the three visco-plastic hardening curves for the MD direction obtained by the test carried out during this study and the DD direction result curves taken from Sahraei [17]. The engineering stress vs. engineering strain curves in MD direction are fitted to the Swift and Voce equations [31] and extrapolated to 100 % strain, using an interpolation function as described in Bulla [16]. Figure 14 shows the FE results of the material behavior in MD direction. It can be seen that after the first localization occurs, the localized area remains at the same group of elements. These elements are further extending whereas the neighboring elements remain at the strain level, they had before the localization started. For these material curves (shown in Figure. 13, left side) the slope of the true stress vs. true plastic strain curve is decreasing, which means, the material will not become stable back again.  Figure 15 shows the loading direction in the TD direction. As it can be seen in the hardening curve (Figure 13 right side), the slope decreases until ~11% plastic strain and then increases again, which results in stabilization of the material. For materials with such hardening curves, the material becomes unstable as soon as the 1st derivative decreases and undergoes the hardening true-stress-strain curve. Then, the first localization can occur ( Figure 15b) and become stable again. During further loading ( Figure 15c) the material stiffness increases, and the first derivative raises above the true-stress-strain hardening curve. The material becomes stable again and localized strain occurs in another area (Figure 15d) during further loading of the structure. This numerical behavior, which can be simulated only with a sufficient complex material model, is also observed in real tests ( Figure 12). Figure 16 shows the corresponding force vs. displacement simulation result curves. The influence of the strain rate dependent Young's moduli (see Table 1) is clearly visible (Figure 16, right side).   Figure 15 shows the loading direction in the TD direction. As it can be seen in the hardening curve ( Figure 13 right side), the slope decreases until~11% plastic strain and then increases again, which results in stabilization of the material. For materials with such hardening curves, the material becomes unstable as soon as the 1st derivative decreases and undergoes the hardening true-stress-strain curve. Then, the first localization can occur (Figure 15b) and become stable again. During further loading ( Figure 15c) the material stiffness increases, and the first derivative raises above the true-stress-strain hardening curve. The material becomes stable again and localized strain occurs in another area (Figure 15d) during further loading of the structure. This numerical behavior, which can be simulated only with a sufficient complex material model, is also observed in real tests ( Figure 12). Figure 16 shows the corresponding force vs. displacement simulation result curves. The influence of the strain rate dependent Young's moduli (see Table 1) is clearly visible (Figure 16, right side).  For materials with such hardening curves, the material becomes unstable as soon as the 1st derivative decreases and undergoes the hardening true-stress-strain curve. Then, the first localization can occur (Figure 15b) and become stable again. During further loading ( Figure 15c) the material stiffness increases, and the first derivative raises above the true-stress-strain hardening curve. The material becomes stable again and localized strain occurs in another area (Figure 15d) during further loading of the structure. This numerical behavior, which can be simulated only with a sufficient complex material model, is also observed in real tests ( Figure 12). Figure 16 shows the corresponding force vs. displacement simulation result curves. The influence of the strain rate dependent Young's moduli (see Table 1) is clearly visible (Figure 16, right side).
The simulations were conducted with the explicit FEA solver Radioss (version 2021) by Altair Engineering Inc. (1820 E. Big Beaver Rd., Troy MI 48083, USA), on a Windows64 computer using four Intel i7-6820 CPUs at 2.7 GHz with 32GB RAM. The material model was developed in FORTRAN and linked to the FEA solver as a dynamic linked library (dll) compiled using the freeware MinGW-W64 project with a gcc 7.2.0 compiler [32]. The developed validation models were meshed using 4-node shell elements with a mesh size of 1 × 1 mm and a defined thickness of 0.025 mm. The simulations were conducted with the explicit FEA solver Radioss (version 2021) by Altair Engineering Inc. (1820 E. Big Beaver Rd., Troy MI 48083, USA), on a Windows64 computer using four Intel i7-6820 CPUs at 2.7 GHz with 32GB RAM. The material model was developed in FORTRAN and linked to the FEA solver as a dynamic linked library (dll) compiled using the freeware MinGW-W64 project with a gcc 7.2.0 compiler [32]. The developed validation models were meshed using 4-node shell elements with a mesh size of 1 × 1 mm and a defined thickness of 0.025 mm.

Discussion
The objective of the presented work was the development of a finite element material model, for nonlinear, structural CAE simulation and predicting the mechanical behavior of a PE separator material under dynamic and high mechanical load, used for the prediction of stiffness and failure in Li-ion batteries. Nevertheless, our focus so far is set on crashworthiness analysis where batteries undergo very complicated state of deformation where bending plays an important role. Therefore, we set our focus on tensile loading until fracture for the separators.
Based on the real test results, performed with different velocities and different orientations of the specimen, a novel material model was developed. It considers the visco-elastic orthotropic, visco-plastic orthotropic behavior under high mechanical loadings. The developed FE material model accounts for the strain-rate effect in elastic loading phase and different material orientations, as well as the different elasto-plasticorthotropic behavior in the plastic material domain, which was measured in the real tests. Within the material model, the Von Mises potential was used, which is widely used in the industry and suitable for modeling elasto-plastic materials. For material behavior characterization the only necessary information was the results from real tests in the three

Discussion
The objective of the presented work was the development of a finite element material model, for nonlinear, structural CAE simulation and predicting the mechanical behavior of a PE separator material under dynamic and high mechanical load, used for the prediction of stiffness and failure in Li-ion batteries. Nevertheless, our focus so far is set on crashworthiness analysis where batteries undergo very complicated state of deformation where bending plays an important role. Therefore, we set our focus on tensile loading until fracture for the separators.
Based on the real test results, performed with different velocities and different orientations of the specimen, a novel material model was developed. It considers the visco-elastic orthotropic, visco-plastic orthotropic behavior under high mechanical loadings. The developed FE material model accounts for the strain-rate effect in elastic loading phase and different material orientations, as well as the different elasto-plastic-orthotropic behavior in the plastic material domain, which was measured in the real tests. Within the material model, the Von Mises potential was used, which is widely used in the industry and suitable for modeling elasto-plastic materials. For material behavior characterization the only necessary information was the results from real tests in the three directions (MD, DD, TD), using a simple tension test under different loading velocities. Microscopic investigation of the battery separator after high mechanical loading along the MD and TD directions are published by Zhang [33]. Once the numerical material reaches instability, a further consideration of damage propagation and failure must be considered.
Compared to existing orthotropic material models, as developed by Hill or Barlat [34,35], the modeling is much easier due to the fact, that the orthogonal directions are obtained directly by using the real tests as input for this model. It is also possible to apply the Hill or Barlat yield surface if further studies show this need.

Conclusions and Outlook
Within this work, our novel material model and the developed FE model succeeded in predicting the test for all three directions (MD, DD, TD) with different velocities in the linear elastic as well as the plastic domain within a range of 5%. This lies within the scatter range of the tested specimen. For predicting the damage during this dynamic and high mechanical load, a suited failure model should be used in addition to the material model, as proposed by Bulla [16].
This should be useful in the design process of batteries and applications that contain batteries and will serve as an important new computational tool for assessing the safety of lithium-ion batteries against high mechanical loading and crashworthiness. The material model will be implemented in the commercial FE solver Altair Radioss and will thus be available for practical application in e-mobility simulation in near future. For porous polymers such as separators, pore orientation and connectivity affect their effective mechanical properties and performance in a battery, too [36]. However, this is the topic of further investigation and the usage of an extended model.