4. Model simulations of mechanical properties of selected samples

The analysis and measurements of mechanical properties of PU foam samples are generally limited only to certain information, and therefore, it cannot tell us the immediate distribution of deformation and tension in the material structure. This is because the options for measurements are limited; options to place sensors and some properties cannot be well measured (for example, distribution of main tension and deformation of the cellular structure). In this case, a compilation

of appropriate model simulation according to the numerical method becomes a viable option. Programming a model simulation in FEM setting is the most significant to our purpose, but certain ways also offer other numerical methods, for example, method of discrete elements (MDE), method of boundary elements (BEM), or method of finite volume (FVM). The FEM method was used exclusively in this work. Mechanical compression of the selected PU foam samples creates many heterogeneous attributes in its inner structure, which change with the size of deformation, as was stated. While modeling such structures, Refs. [8, 20] agree that it is necessary to simplify or ignore certain characteristic attributes, and both state that a great problem of modeling non-linear attributes is mostly the description of the main tensions in short time differentiations Δt ¼ tiþ<sup>1</sup> � ti. The solution of peripheral problem of great deformations created by compression of the sample further lies not only in the correct input of peripheral conditions and material properties, but also in the construction of proposed net of finite elements. FEM programs are currently very well worked out and allow for the conversion of continuous problem solution to the final solution, where it is possible to suggest appropriate geometrically simple partial sub-regions (finite elements) for approximate solution in the preprocessor.

m�<sup>3</sup>

96 Aspects of Polyurethanes

in force.

and they can be summarized as follows:

dampening of PU foam.

of the finite element method.

,the packing density Ψ (1) increases while influencing the average cell size of the material structure (Table 1). Cellular structure of the PU foam sample influences strain-stress curve during compression. Typical course is characterized by the initial rigidity depending on the strain rate, given stable area so-called plateau and the final steep exponential increase

Mechanical properties are significantly influenced by the character of the cellular structure,

• Deformation ε is the function of not only tension σ, but also time t, and in the examined

• Material deformation is dampened by inner viscose resistances, i.e., dampening of the

• The faster the deformation happens, the more intense the dampening effect of the material viscosity materializes, but also the air viscosity, which cannot be squeezed out of the cellular structure immediately and, therefore, manifests as an initial significant increase in rigidity (the slower the compression, the lower the initial rigidity ! the initial rigidity is

• Recovery manifests (recovery—the ability to immediately recover after deformation)

• Relaxation manifests, i.e., decrease in tension in prestressed condition <sup>σ</sup>ðtÞj<sup>z</sup>�konst <sup>&</sup>lt; <sup>σ</sup>ð0<sup>Þ</sup>

• Creep manifests, i.e., the growth in deformation under constant pressure <sup>ε</sup>ðtÞj<sup>σ</sup>�konst <sup>&</sup>gt; <sup>ε</sup>ð0<sup>Þ</sup> • Mathematic-physical description of PU foam mechanical properties can be described by a constitutive relations and by rheologic models, for example, modified n-parametric Tucket model according to which it is possible to express the coefficient of rigidity and

• For a qualitative analysis of the PU foam during compression, or for wider knowledge of mechanical properties that cannot be sufficiently measured nor mathematically described (distribution of main tensions and transformation in individual directions, contact pressures), it is suitable to construct model simulation of mechanical properties in the setting

The analysis and measurements of mechanical properties of PU foam samples are generally limited only to certain information, and therefore, it cannot tell us the immediate distribution of deformation and tension in the material structure. This is because the options for measurements are limited; options to place sensors and some properties cannot be well measured (for example, distribution of main tension and deformation of the cellular structure). In this case, a compilation

4. Model simulations of mechanical properties of selected samples

area (up to 95 � 3% transformation), it can be considered reversible.

cellular structure, ηtPU, therefore cannot be realized immediately.

nearing the behavior of the so-called plateau).

given viscoelastic properties determined by hysteresis.

Let ℵ⊂R<sup>3</sup> is a continuous area of three-dimensional space. Their boundaries will be Γ, where Γ is the so-called Lipschitz boundary, and let an approximation of the selected basis functions are derived over each finite element with dimension l, because any continuous function may be represented by linear combinations of algebraic polynomials converging to a continuous solution i.e. lim<sup>l</sup>!<sup>0</sup> ! ξ ≈ 1. Thus, FEM can be understood as a special type of a variation method using mathematical description of solved problem. Current substantial commercial FEM software (for example, Ansys, Abaqus, Permas, LS-Dyna, Marc, and PAM CRASH) allow to assemble and then solve the loading of nonlinear materials not only with viscoelastic properties using mathematical relationships based on continuum mechanics and rheological model (for example, Kelvin model, Maxwell model, and so on). Also they allow with sufficient accuracy the studying and modeling of qualitatively more complex problems such as contacts of two or more bodies (i.e., the interaction between the material and probe or human body).

### 4.1. Selection of suitable software for the assembly of FEM model

In this study for all model simulations, software PAM CRASH was selected. It is the FEM software from ESI-Group company (http://www.esi-group.com/) that is used for the study of nonlinear isotropic and anisotropic materials and contact problems in quasi-static and dynamic states. Similar to FEM software as LS-Dyna, Abaqus Explicit, and ANSYS Explicit. The basic principle of explicit methods is second Newton's law, which may be expressed in a matrix form by Eq. (31).

$$M \cdot \ddot{\mu} = F^E - F^I{}\_{\prime} \tag{31}$$

where M is a matrix of mass, u€ is the acceleration matrix of node vectors, FE is the matrix of the vector of external forces acting on the node, and FI is a matrix of vectors of internal forces (volume forces).

The matrix of acceleration vectors, where the acceleration expresses second derivation of the searched (unknowns) displacements (Eq. (32)), can be obtained according to Eq. (31). Then, vector matrix of internal and external forces can be expressed by Eqs. (33) and (34).

$$
\ddot{\mu} = M^{-1} \cdot F^{E} - F^{l} \tag{32}
$$

$$F^{l} = \sum\_{\iota=1}^{N\_{\iota}} \left( \int\_{\varOmega} B^{T} \cdot \sigma\_{n} d\varOmega + F^{Hurg} + F^{kont} \right), \tag{33}$$

$$F^{E} = \sum\_{\varepsilon=1}^{N\_{\varepsilon}} \left( \int\_{\upsilon\_{\varepsilon}} \rho \cdot \kappa\_{\nu\circ}^{\varepsilon} d\upsilon\_{\varepsilon} + \int\_{s\_{\varepsilon}} \chi\_{i} \cdot \mathfrak{d}\_{\circ\circ}^{\varepsilon} d\varsigma\_{\varepsilon} \right) \tag{34}$$

where B is matrix of basic functions of the strain, Fkont is the vector of the contact forces, FHurg is a vector of hourglassing damping forces, σ<sup>n</sup> is the matrix of acting stress in member, ρ is density of the member, κ<sup>i</sup> is vector of volume forces, and χ<sup>i</sup> is vector of surface forces.

Founded matrix of displacement vectors u (Eq. (36)) can be expressed by an integration of the acceleration u€ or velocity of displacements u\_ (Eq. (35)) in accordance with following formulas:

$$
\dot{\mu} = \dot{u}\_{t + \Delta t/2} = \dot{u}\_{t - \Delta t/2} + \ddot{u}\_t \cdot \frac{\Delta t\_t + \Delta t\_{t + \Delta t}}{2} \tag{35}
$$

$$
\mu = \mu\_{t+\Delta t} = \mu\_t + \dot{\mathfrak{u}}\_{t+\Delta t/2} \cdot \Delta t\_{t+\Delta t} \tag{36}
$$

where ut is a vector of instantaneous velocity and ut�Δ<sup>t</sup> and utþΔ<sup>t</sup> are vectors of previous or subsequent displacements.

The software has sophisticated algorithms for complex nonlinear contacts [21], where ongoing simulation model is divided into a selected sequence of m-intervals (where m ≤ t and mmin ¼ 1). For each time step, the displacement vector ut is calculated. This vector describes that in following time step, the origin geometry A0 is changed to current geometry AtþΔ<sup>t</sup> in consequence of the change of displacement vector utþΔ<sup>t</sup> in Eq. (37).

$$A\_{t+\Delta t} = A0 + \mu\_{t+\Delta t} \tag{37}$$

In further steps, instantaneous Cauchy stress σ<sup>t</sup>þΔ<sup>t</sup> can be expressed using constitutive relations according to Eq. (38) that the algorithm processor expresses as strain change of elements dε ¼ ∂u=∂Xi (i ¼1,.., 3). Subsequently, a new vector of internal forces for individual nodes is calculated. Variable t þ Δt is overwritten with t, and the calculation proceeds to the next step.

$$
\sigma\_{t+\Delta t} = f(\sigma\_{t\prime} \, d\varepsilon),
\tag{38}
$$

The resulting simulation time step Δt is described with Eq. (39) and relates to the calculation time, which is proportional to the size of the smallest element lmin, the square root of the material density ρ, and inversely proportional to the square root of elastic modulus E.

The advantage of the explicit method is an order of magnitude faster than explicit method, because the implicit method the time step becomes a quadratic function [21].

$$
\Delta t \le \Delta t^{krit} = l\_{\text{min}} \cdot \sqrt{\frac{\rho}{E}} \tag{39}
$$

where Δt krit expresses the minimal (critical) time step for the simulation.

The matrix of acceleration vectors, where the acceleration expresses second derivation of the searched (unknowns) displacements (Eq. (32)), can be obtained according to Eq. (31). Then,

<sup>B</sup><sup>T</sup> � <sup>σ</sup>nd<sup>Ω</sup> <sup>þ</sup> FHurg <sup>þ</sup> Fkont

ð

<sup>χ</sup><sup>i</sup> � ϑε <sup>∞</sup>ds<sup>ε</sup>

Δtt þ ΔttþΔ<sup>t</sup>

u ¼ utþΔ<sup>t</sup> ¼ ut þ u\_ <sup>t</sup>þΔt=<sup>2</sup> � ΔttþΔ<sup>t</sup> (36)

AtþΔ<sup>t</sup> ¼ A0 þ utþΔ<sup>t</sup> (37)

σ<sup>t</sup>þΔ<sup>t</sup> ¼ fðσt, dεÞ, (38)

sε

<sup>i</sup><sup>∞</sup>dv<sup>ε</sup> þ

where B is matrix of basic functions of the strain, Fkont is the vector of the contact forces, FHurg is a vector of hourglassing damping forces, σ<sup>n</sup> is the matrix of acting stress in member, ρ is

Founded matrix of displacement vectors u (Eq. (36)) can be expressed by an integration of the acceleration u€ or velocity of displacements u\_ (Eq. (35)) in accordance with following formulas:

where ut is a vector of instantaneous velocity and ut�Δ<sup>t</sup> and utþΔ<sup>t</sup> are vectors of previous or

The software has sophisticated algorithms for complex nonlinear contacts [21], where ongoing simulation model is divided into a selected sequence of m-intervals (where m ≤ t and mmin ¼ 1). For each time step, the displacement vector ut is calculated. This vector describes that in following time step, the origin geometry A0 is changed to current geometry AtþΔ<sup>t</sup> in

In further steps, instantaneous Cauchy stress σ<sup>t</sup>þΔ<sup>t</sup> can be expressed using constitutive relations according to Eq. (38) that the algorithm processor expresses as strain change of elements dε ¼ ∂u=∂Xi (i ¼1,.., 3). Subsequently, a new vector of internal forces for individual nodes is calculated. Variable t þ Δt is overwritten with t, and the calculation proceeds to the

The resulting simulation time step Δt is described with Eq. (39) and relates to the calculation time, which is proportional to the size of the smallest element lmin, the square root of the

material density ρ, and inversely proportional to the square root of elastic modulus E.

<sup>u</sup>€ <sup>¼</sup> <sup>M</sup>�<sup>1</sup> � FE � FI (32)

A, (33)

A (34)

<sup>2</sup> (35)

1

1

vector matrix of internal and external forces can be expressed by Eqs. (33) and (34).

ð

0 @

Ω

0 @

ð

<sup>ρ</sup> � <sup>κ</sup><sup>ε</sup>

density of the member, κ<sup>i</sup> is vector of volume forces, and χ<sup>i</sup> is vector of surface forces.

u\_ ¼ u\_ <sup>t</sup>þΔt=<sup>2</sup> ¼ u\_ <sup>t</sup>�Δt=<sup>2</sup> þ u€<sup>t</sup> �

vε

FI <sup>¼</sup> <sup>X</sup> N<sup>ε</sup>

ε¼1

<sup>F</sup><sup>E</sup> <sup>¼</sup> <sup>X</sup> N<sup>ε</sup>

consequence of the change of displacement vector utþΔ<sup>t</sup> in Eq. (37).

subsequent displacements.

98 Aspects of Polyurethanes

next step.

ε¼1

Subsequently, the processor for viscoelastic structure expresses Cauchy stress σij by the tensor of a nominal stress σnom ij , which is inversely proportional to vectors extension λ<sup>i</sup> (Eq. (40)) as described [22].

$$
\sigma\_{ij} = \frac{\sigma\_{ij}^{\text{nom}}}{\lambda\_i \cdot \lambda\_k} \tag{40}
$$

where λ<sup>i</sup> expresses vectors of extension to the principal directions and λ<sup>k</sup> is a permutation index.

### 4.2. FEM simulation of mechanical properties of selected samples of PU foam

Simulations were performed in FE software for sample of the PU foam sample. Simulations were performed for a complete assessment of the selected material because the experimental methods cannot give an explanation of behavior and change of the shape, especially under dynamic loading. Model simulations were performed in the following steps:


Applied types and sizes of elements that affect the final time step Δt of the model (Eq. (39)) are shown in Table 3.

Figure 20. FEM model of sample dynamically compressed without the initial deformation.

Figure 21. Modification of the input parameters in the FEM model dataset (left); visualization of the simulation model in PAM CRASH (right).


Table 3. FEM model dynamically compressed sample.
