A comparison of constitutive models for describing the flow of uncured styrene-butadiene rubber

Uncured styrene-butadiene rubber (SBR) can be modelled as a viscoelastic material with at least two different relaxation mechanisms. In this paper we compare multi-mode constitutive models combining two viscoelastic modes (linear and/or nonlinear) in three possible ways. Our particular choice of the two modes was inspired by models originally developed to describe the response of asphalt binders. We select the model that best fits the experimental data obtained from a modified stress relaxation experiment in the torsional configuration of the plate–plate rheometer. The optimisation of the five model parameters for each model is achieved by minimising the weighted least-squares distance between experimental observations and the computer model output using a tree-structured Parzen estimator algorithm to find an initial guess, followed by further optimisation using the Nelder–Mead simplex algorithm. The results show that the model combining the linear mode and the nonlinear mode is the most suitable variant to describe the observed behaviour of SBR in the given regime. The predictive capabilities of the three models are further examined in changed experimental and numerical configurations. Full data and code to produce the figures in this article are included as supplementary material.


Introduction
Styrene-butadiene rubber (SBR) is the most widely used synthetic rubber and forms a key constituent of many industrial products. It is important to understand its mechanical behaviour in order to accurately predict its response in various settings, both in manufacturing during extrusion, see Rauwendaal [1], and in situ in a finished automotive tyre, see Nakajima [2].
The prediction of the viscoelastic behaviour of rubbers such as SBR is a research area spanning back decades with many contributions. A recent overview with a particular focus on filled SBR compounds was given by Carleo et al. [3]. Montes et al. [4] studied the behaviour of carbon black compounds in various shear flow histories, e.g. stress relaxation and transient shear flow. Ovalle Rodas et al. [5] examined SBR compounds under cyclic loading in a wide range of temperatures. In the context of extrusion, Jugo Viloria et al. [6] characterised the et al. [11,12] that are in turn based on the earlier theoretical work of Rajagopal and Srinivasa [13]. We remark that our work is based on a phenomenological approach to developing viscoelastic models; cf. the popular micromechanical approach, see e.g. Freund et al. [14]. The first model was derived by Málek et al. [11] to describe asphalt binders. The second model can be found in [12], where it is discussed as an equivalent of the classical Burgers model, cf. Burgers [15]. Both models were developed to describe incompressible materials, meaning that all admissible deformations of a body constituted by the given material have to be isochoric. A key differentiating feature between the two models is that the former maintains the instantaneous purely elastic and purely dissipative parts of the total deformation to be isochoric, while the latter admits the instantaneous parts to be compressible (yet the total deformation remains isochoric). The third model is obtained as a combination of the two. This newly proposed phenomenological model describes the material as being composed of different substances, where some of them respond as an incompressible material and some of them as a compressible material.
All of the models we consider have the capability to capture two relaxation mechanisms, although in principle this could be extended further at the expense of introducing extra parameters that need identifying. All variants considered in this study have five parameters, one related to the additional viscous dissipation along with two pairs of parameters (dynamic viscosity and elastic shear modulus) associated with the two relaxation mechanisms. We show that these parameters can be identified using torque and normal force data from a plate-plate rheometer experiment. In the plate-plate configuration each model leads to a coupled system of ordinary differential equation (ODE) that can be solved numerically to compute estimates of the torque and the normal force for a given parameter set, see [12]. For each model the optimal set of parameters is found by minimising the loss between the data and the ODE model in a weighted least-squares sense. An advantage of the chosen approach is that the model is developed directly in a three-dimensional setting allowing parameter identification using experimental data, and/or subsequent direct use of the model and its parameters in numerical simulations of flows in complex geometries discretised using an appropriate numerical method.
The main contributions of this paper are as follows: • We propose a constitutive model with two relaxation mechanisms that is in turn shown to be superior to two existing models of the same type in predicting the behaviour of SBR in a plate-plate rheometry experiment. • In contrast with Málek et al. [11], we use a least-squares loss function to define the optimal or best parameter set for each model. Furthermore, we define a procedure for model selection based on choosing the model with minimum loss function at the optimised parameters. • We demonstrate how a tree-structured Parzen estimator (TPE) algorithm, see Bergstra et al. [16], can be used to find a good initial guess for the parameters. The initial guess is then used in a subsequent Nelder-Mead simplex minimisation algorithm, see Nelder and Mead [17], to find the optimal parameter set. • We use a Monte-Carlo analysis to show that a relatively low variance in the radius of the sample leads to a large variance in the output torque. This finding confirms comments from experimental papers, e.g. Hellström et al. [18], that it is critical to control the sample radius for repeatability.
An outline of this paper is as follows; in Section 2 we give an overview of the general governing and constitutive equations of a viscoelastic fluid. In Section 3 we derive a system of ODEs describing the flow of a viscoelastic fluid in the plate-plate rheometry configuration. In Section 4 we describe the experimental setup. In Section 5 we explain the procedure for finding the optimal parameter set for each model and selecting the best model. We show the results of our numerical study in Section 6, before concluding with a discussion in Section 7.

Governing equations
The kinematics of the isothermal flow of a fluid-like material in a given space-time domain × (0, end ] is mathematically described by the following system of partial differential equations (PDEs), namely This system needs to be supplemented by an appropriate set of initial conditions at time = 0 and a set of time-dependent boundary conditions defined on × (0, end ] which we will discuss later. In the above equations, is the velocity, is the density, T is the symmetric Cauchy stress tensor and is the vector of gravitational acceleration. The gravitational force is assumed to be the only body force acting on the fluid. Introducing the cylindrical coordinate system { , , } with the coordinates ( , , ), such that is aligned with but pointing in the opposite direction, we may write = −∇( ).
In this paper we will limit ourselves to homogeneous, incompressible fluids for which the Cauchy stress tensor T may be written as Here, I is the unit tensor, S is the extra stress tensor and̃is the pressure. More precisely,̃is the Lagrange multiplier that enforces the constraint of incompressibility, see Rajagopal [19]. As such, it must be treated as a primitive variable in the description of the flow. Let denote the pressure including the hydrostatic contribution in the sense of the definition This allows us to rewrite the balance of linear momentum (1a) in the form Furthermore, due to the assumptions of homogeneity and incompressibility, the balance of mass (1b) reduces to div = 0. (4b)

Constitutive equations
In order to supply the dynamic response of the material into the mathematical model, we need to characterise the extra stress tensor S by means of a constitutive equation which relates it to the strain history. Before we discuss some specific examples of such relations, let us make a few comments on the notation used in the remainder of the document.

Remarks on notation
The symbols L and D are reserved for the velocity gradient and its symmetric part, that is, The quantity D is also known as the rate-of-deformation tensor. Its norm, |D| = √ D ∶ D, is used to introduce the shear ratėthrough the formula 2 The upper convected time derivative of an arbitrary second-order tensor A is given by while its deviatoric (traceless) part is defined as

Models studied in this work
We are interested in the specific class of viscoelastic rate-type fluid models that can be derived in a thermodynamically consistent manner using the framework proposed by Rajagopal and Srinivasa [13]. This thermodynamic framework has been further extended and successfully used to develop models for describing the mechanical response of materials as complicated as asphalt binders; see Krishnan and Rajagopal [20] and related papers by Narayan et al. [21,22] or Málek et al. [11,23]. Let us briefly discuss the main ideas of the chosen approach in the following paragraph; consult Málek et al. [12] and the references therein for the detailed explanation.
In general, the observed fluid-like material can be associated with an abstract body which may occupy regions of three-dimensional Euclidean space. The approach by Rajagopal and Srinivasa [13] builds on the assumption of the existence of a natural configuration associated with the placement of a body that evolves as the body deforms, see Fig. 1. This concept allows one to split the total deformation into that associated with the purely elastic response and the dissipative response. Such a decomposition becomes useful when prescribing the way in which the body stores and dissipates energy. This piece of information is supplied into the model by means of two constitutive assumptions specifying the Helmholtz free energy and the rate of entropy production . Depending on the precise form of these two scalar functions, one can obtain different forms of the Cauchy stress tensor T including its evolution equation.
If the material is known to possess more relaxation mechanisms, typically due to its complex microstructure, then the fundamental idea applied in the derivation of a suitable model is to associate these mechanisms with different underlying natural configurations. The present numerical study compares several models which can be written in the canonical form Here, each B = F F ⊤ ( ∈ {1, 2}) denotes the left Cauchy-Green tensor associated with the elastic response between the th natural configuration and the current configuration of the deformed body, while  represents a tensor function which takes the respective B as its single argument. 3 Finally, and are model parameters with the meaning of dynamic viscosity and elastic shear modulus respectively. The term with 0 appearing in (9a) describes the additional viscous dissipation. 4 Each of the two relaxation mechanisms, which can be captured by the present class of models, is characterised by the relaxation time according to the formula Table 1 specifies the three models of our interest. Each model is labelled by  with a pair of subindices indicating either Linearity or Table 1 Variants of viscoelastic models for SBR given by constitutive equations in the form (9).

Table 2
Expressions representing the shear viscosity and the normal stress differences 1 , 2 derived from models in Table 1. Models  NN and  LN capture shear thinning, model  LL does not.
Nonlinearity of the functions  1 and  2 respectively. Table 2 summarises the basic rheometric functions, obtainable in the simple shear flow and/or capillary flow configurations, for the present models. The first model ( LL ) is formally equivalent to the classical Burgers model, see [12]. The second model ( NN ) was previously used to describe the response of asphalt binders, see [11]. Both models describe incompressible fluids in the sense that the total deformation remains isochoric. The difference between them is given by the fact that the instantaneous purely elastic and purely dissipative parts of the total deformation -associated with individual natural configurations -are not necessarily isochoric in the first case ( LL ), but the opposite holds true in the second case ( NN ). The third model ( LN ) is a combination of both.
To best of our knowledge, the last model has not been explicitly considered within the context of modelling the response of materials such as SBR. That said, let us emphasise that all variants covered in the present work are special cases of general multi-mode viscoelastic models with their individual modes corresponding to simpler, well established models (see Remark 2.1). There are multiple ways to derive such multi-mode models. The basic ideas for a unified derivation of the variants from Table 1 within the framework discussed above are presented in Section 2.3 for the specific case  LN . Remark 2.1. Constitutive equations in the canonical form (9) represent a general multi-mode viscoelastic model limited to the case with two modes and a separate viscous contribution. The linear mode L corresponds to the well-known upper convected Maxwell (UCM) model or, more precisely, to the classical Oldroyd-B model taking into account the separate viscous contribution; see Oldroyd [24]. On the other hand, the nonlinear mode N can be recognised as the simplest Leonov model; see Leonov [25]. Indeed, if we formally disable the linear mode (and the separate viscous contribution) in the combined model  LN , we will end up with the single-mode model where we have omitted the superfluous subscripts. 5 5 We should write B ( ) instead of B to emphasise the fact that this quantity characterises the deformation experienced by the ''spring'' part of the virtual spring-dashpot system shown in Fig. 1. 1 ( ) and 2 ( ) , the total deformation F is split into dissipative parts and purely elastic parts corresponding to the mappings

Derivation of the model  LN
Recall that the Eqs. (1a) and (1b) describe the balance of mass and the balance of linear momentum for the given physical system. The balance of angular momentum is in our case equivalent to the requirement that the Cauchy stress tensor T is symmetric. Let us rewrite this basic set of balance equations in the forṁ wherėd ef = + ⋅ ∇ is the shorthand notation for the material time derivative, which is used exclusively within this section. Following the framework outlined in [12, Section 2.2], we consider the additional equation which can be obtained -under circumstances where only isothermal processes are allowed -by combining the balance of energy and the formulation of the second law of thermodynamics.
As the material is supposed to be incompressible, it can undergo only isochoric motions, meaning that Tr D = div = 0. With this constraint the balance of mass reduces tȯ = 0.
Next, we assume that which is the formula that describes the way in which the material stores the energy. It says that the elastic response from the first natural configuration corresponds to that of compressible neo-Hookean solid, while the response from the second one is the same but incompressible. Inserting (15) into (13), using also (14) and doing some further algebraic manipulations 6 (see [27,Section 4.4.1]), we end up with 6 In particular, we make use of the identitiesṪr where C = F ⊤ F ( ∈ 1, 2) denotes the right Cauchy-Green tensor associated with the elastic part of the total deformation (with respect to the th natural configuration) and D is defined by means of its The dissipative mechanism associated with the viscous response from the second natural configuration is assumed to be that for an incompressible fluid, meaning that Tr D 2 = 0.
Using the incompressibility constraints Tr D = Tr D 2 = 0, and the decompositions (2) and (8), we rewrite the expression (16) to take the form As a next step we postulate that These constitutive relations ensure that Each of the last two relations can be rewritten in terms of the corresponding B = F F ⊤ using the identity see [27,Eq. (163)]. We proceed as follows. First, multiplying (18b) from the left by F 1 and from the right by F −1 1 we obtain Second, multiplying (18c) from the left by F 2 and from the right by F ⊤ 2 , and recalling Tr C 2 = Tr B 2 , we obtain Finally, we use (19) to replace the tensor product on the right hand side in the above relations by the upper convected time derivative of the corresponding left Cauchy-Green tensor. After some additional algebraic manipulations, we rewrite the constitutive relations (18)   take their final form with the term that can be incorporated into the pressure variable similarly as the hydrostatic contribution in the definition (3).

Plate-plate rheometry
A rotational shear flow is examined in the plate-plate geometry for the constitutive equations discussed previously in Section 2.2. In this configuration, a sample of a fluid-like material is placed between two cylindrical plates, where the lower one is fixed and the upper one can rotate around its axis, see Fig. 2. The diameter of the sample, = 2 , is supposed to coincide with the diameter of the plates. The fixed distance between the plates is denoted by . The material is assumed to stick to the plates, and the flow induced by the rotation is considered to be laminar and axisymmetric.
We shall focus on the experimental setting with the controlled shear rate, that is, with the controlled rotational speed. We suppose that the material flows only in the azimuthal direction , satisfying the no-slip boundary condition at the plates, and the no-traction condition applies on the lateral surface which remains vertical during the deformation. More precisely, we put where denotes the angular velocity that is controlled in time, see (43). The above ansatz satisfies (4b) and further implies Then, according to (6), we havė The corresponding strain can be obtained by integrating the previous relation in time. In this way, one obtains where denotes the amount of rotation of the upper plate in the sense of the relation = d . Supposing that we are dealing with the creeping flow, so that the inertial effects can be neglected, the balance of linear momentum reads These equations are further coupled with constitutive equations in the form (9). In particular, the components of B ( ∈ {1, 2}) must satisfy the system of evolution equations (28a)-(28f) given in Box I. Let us emphasise that none of these equations explicitly depends on the variable . We shall assume that the material is initially at rest and relaxed, meaning that where O denotes the zero tensor. Since the initial conditions do not depend on , it follows from (28) that the same must hold for the components of B (for > 0) regardless of the choice of  . As a consequence, the components of S are also independent of , see (9a) and (24). Moreover, the choice of initial conditions in (29) ensures that | = | = 0, see (28e) and (28f), and thus also = = 0. Based on these observations, we see that the left hand side of (27c) is identically equal to zero, so that is independent of , (27b) is identically satisfied and (27a) can be rewritten as 7 We are interested in computing the torque and the normal force on the upper plate ⊂ . These quantities are defined by where we have explicitly denoted the dependence of the computed torque and normal force on time, but also on the choice of model The expression for the computation of the normal force can be rewritten in terms of normal stress differences using Eq. (30), see [11]. The alternative formula reads The integrals in (31) and (33) are computed using a three-point Newton-Cotes quadrature scheme (Simpson's rule). At each quadrature point the system of ODEs (28) is solved using a fourth-order adaptive Runge-Kutta method. We implement the solver using the routines available in the Python library SciPy, see Jones et al. [28].

Sample preparation
The sample under investigation is SBR containing 27% of styrene with functionalised chain-ends that are designed to enhance the interaction with silica fillers. The SBR has a molecular weight of 310 000 g mol −1 measured with gel permeation chromatography (GPC) relative to polystyrene standard. The samples under investigation do not contain any fillers (e.g. silica, carbon black) or curatives. However, some other ingredients are incorporated to improve the processability, see Gansen et al. [29] for more details. After mixing and milling, the rubber has the form of a large rubber sheet of an approximate thickness of 4 mm with a highly corrugated surface ( Fig. 3(a)). From this sheet, two discs with a diameter of 20 mm are stamped out. These two discs are positioned one above the other and placed between the two plates of the Haake Mars plate-plate rheometer. A sealed oven encloses the sample. The oven is heated up to 120 • C and the sample is slowly compressed to a thickness of around 5 mm. This procedure typically takes between one and two hours. After cooling down the sample to room temperature, it is carefully removed and due to the compression the excess rubber is removed by stamping the sample out again ( Fig. 3(b)).

Experimental procedure
A modified stress relaxation test is conducted to produce experimental data. The standard stress relaxation test requires the application of a step increase in strain, which is afterwards held constant for a given time. The stress that builds up inside the material, as a consequence of its sudden deformation, will relax over time due to various molecular mechanisms. In the torsional configuration described in Section 3, the application of the step strain corresponds to the instantaneous change of the deflection angle , see (26). It is of course physically impossible to twist the sample instantly. The reason is twofold. First, the equipment ramps to the required steady state in a finite time. Second, even if we admit an ''instantaneous'' deformation to be a fast motion occurring on a time scale that is much shorter than the time of observation, it is not feasible to conduct the experiment under the simplifying conditions discussed in Section 3. Therefore, we tried to use a different testing scheme where the ramping of the deflection angle is controlled to gradually increase to the specified value at = 0 (0 < 0 < end ).
To study the stress relaxation of asphalts, Narayan et al. [30] suggested to use the template with the controlled shear rate in terms of the angular velocity of the form see (25). Here, 0 is a specified constant value. The Haake Mars plateplate rheometer allows one to set up the required shear ratė0 = 0 ∕ .

Experimental data
The experiment was performed at a fixed temperature of 120 • C. Figs. 4-5 show the experimental input/output data for four measurements {A1, A2, B1, B2}. Different samples were used for measurements labelled by different letters. Each of these measurements was repeated twice with the given sample. The thickness of the sample was slightly different in each case as a result of the complicated preparation stage, see Section 4.1. In particular, the mean thickness was = 4.82 mm with the standard deviation of approximately 0.02 mm. The shear ratė 0 = 0.03 s −1 was specified for ∈ (0, 0 ] with 0 = 100 s. After that, the upper plate was held at the same position and the relaxation of the sample was measured over the additional time of 300 s. The normal force and the torque were recorded simultaneously. Looking carefully at Fig. 4, we see that the ramping of the deflection angle is not linear during the time interval (0, 0 ], meaning that the angular velocity is not constant as we have required. We believe that this behaviour is caused by the complex material properties already discussed above, which make it technically difficult to set up and maintain the specified shear rate. Nevertheless, the ramping was done in a repeated path in most cases. It thus makes sense to take such data, fit it by a customised function and use it as an input for the numerical simulation; see relations (41)-(43) in Section 6.
The output data are expected to overlap to a certain extent, taking into account the measurement error as well as the experimental error including variations in the input data for individual measurements (see Fig. 6). The accuracy of the equipment for measuring the torque and the normal force was estimated to be 0.05 mN m and 0.03 N, respectively, following methods suggested by Narayan et al. [30]. In Fig. 5, we see that the measured normal force is almost perfectly consistent in the previous sense. However, a kind of inconsistency is observed in the measured torque for the two different samples. We believe that the observed behaviour is related with the fact that it has not been possible to prepare the two samples with exactly the same radius (with the precision of 0.1 mm). Our conjecture is supported by the numerical experiment in Section 6.2, where we study a sensitivity of the numerical result on the radius of the sample. Fig. 3. SBR sample before preparation as a rubber sheet with dimensions of ∼ 100 mm × 120 mm, and after preparation as a disc with the diameter of ∼ 20 mm.

The optimisation procedure
Our primary objective is to fit the parameters of the three viscoelastic models proposed in Table 1, using the experimental data produced by the plate-plate rheometer described in Section 4.
The observed experimental data consists of sets of discrete torque and normal force measurements, that is, ) and obs = ( 1 , 2 , … , ) for = 1, 2, … , . These discrete observations are taken at times obs = ( 1 , 2 , … , ). We further define an observation operator  obs ∶ → R that takes a function in some function space and evaluates it at obs . In particular, for ∈ we define With this notation in hand we can define the following cost functional where ∈ R is the generalised set of constitutive parameters for a specific constitutive model , and ‖⋅‖ 2 denotes the standard 2 norm.
( ; , ) and ( ; , ) are the torque and normal force functions computed by numerical solution of (31) and (33) for a specified and model . Let us emphasise that the numerical solution depends on several other parameters, including , and , which may vary from case to case. The coefficients and , defined by ensure that approximately equal weighting is given to fitting the torque and normal force measurements. However, by considering the multiplication factor 2 in the first term on the right hand side of (36), we give some additional preference to fitting the torque measurements.
Remark 5.1. A slightly different cost functional was used by Málek et al. [11,Eq. (66)] to fit the data for asphalt binders. In particular, the authors have considered with the weighting coefficients Here, ‖⋅‖ 1 and ‖⋅‖ ∞ denote the standard 1 and ∞ norms, respectively.
The best set of constitutive parameters for each model  is then defined as the argument ∈ R that minimises (36), namely * label = arg min The output of the above process are three sets of optimal parameters , the minimisation problem (38) is overconstrained by the experimental data. Hence, no regularisation term such as Tikhonov (e.g. [31]) or Bayesian (see e.g. [32,33]) is needed. However, the problem is in general non-convex and multiple local minima may exist.
Finally, we define the optimal model  * as the model with the lowest value of the functional (, ) evaluated at the associated optimal parameter * , that is, Remark 5.3. Because the three models have roughly comparable complexity, involving five parameters and the solution of a similar system of ODEs, we do not penalise model complexity in (39).
The minimisation (38) of the functional (36) is performed in two stages. Because of the presence of multiple local minima of , we first use the TPE algorithm as a global optimiser to efficiently search across a large subspace of the whole parameter space. The TPE algorithm requires the specification of a joint prior probability density function on the parameters. In a Bayesian context this prior reflects the analyst's belief about the parameters prior to the optimisation. We choose the variables to be independent and with uniform prior spanning two orders of magnitude 1 , 2 ∼  (10 4 , 10 6 ) and 1 , 2 ∼  (10 2 , 10 4 ), where  ( , ) is a uniform random distribution with lower and upper bounds and . This corresponds to a weakly informative prior. We use the same prior for all models. We note that our prior does not give preference to a particular ordering of the two relaxation timescales in the proposed  LN model.
The output of the TPE algorithm is then denoted (0) and is subsequently used as an initial guess for a restarted Nelder-Mead algorithm, see [17,34]. The purpose of this second step is to improve the estimate of the local minima in the neighbourhood of (0) . In addition, as suggested by Málek et al. [11], in order to eliminate any further chances of finding local minima, we use the Nelder-Mead algorithm in a restarted fashion in the following sense: 1. we find a solution, that is, a seemingly optimal set of model parameters ( ) using Nelder-Mead, 2. we round the values of ( ) and use them as the new initial guess to find ( +1) .
We have found the above two-stage method (TPE followed by Nelder-Mead) to be highly robust in finding a 'reasonable' minimum * . The TPE algorithm does not require the user to input a specific initial guess and is reliable in finding a (0) where the Nelder-Mead algorithm can then work effectively. We use the implementation of TPE found in the Python library Hyperopt, see Bergstra et al. [16], and the Nelder-Mead algorithm from the Python library SciPy.

Results
We use the following formula with four constant parameters { , , , } to fit 8 the experimental deflection angle in the first time segment (see Fig. 4), namely The change of the angle in the whole time domain is then captured by the piecewise defined function with 0 ≈ 100 s. The corresponding change of the angular velocity, which is the input datum for simulations, is then given in the form of the discontinuous function where the formula for the first time segment is obtained as the time derivative of f it .
Remark 6.1. The reason why we fit the deflection angle data is due to the fact that the rheometer has a sensor for recording the angular displacement of the upper plate. These data are thus expected to be more precise than the angular velocity data, especially at low rotational speeds. This is clearly visible when we explore the data in Fig. 4 for > 100 s. (41) is designed in such a way that it provides us with the possibility to capture two dissimilar experimentally observed scenarios, see Section 6.3.

Comparison of models
Let us start with the identification of parameters for the three models proposed in Table 1 using the experimental data from Section 4 and following the optimisation procedure suggested in Section 5.
We take the average of the four measurements as the base for finding the initial guess of parameters (0) for each of the models. In particular, we fit the mean deflection angle by (42) to get the angular   Table 3.  Table 3. With the initial guess obtained by fitting the mean data, we run the ''restarted'' Nelder-Mead algorithm with the aim to simultaneously  Table 3. fit the data coming from two different measurements, 9 namely {A2, B2}. Our main objective remains the same: to minimise the cost functional (36), but this time using the two sets of observed data ( = 2). We fit the deflection angle data by (42), for each measurement separately, to get a couple of angular velocities in the form (43), see Fig. 6. These functions are subsequently considered as inputs for the numerical computation of ( ; , ) and ( ; , ) corresponding to the individual measurements, see (31) and (33). Table 3 lists the optimised parameter values together with their initial guess found by the TPE algorithm. The corresponding values of the cost functional (36) are contained in the rightmost column. Based on these values, we are tempted to choose the optimal model  * =  LN in accordance with (39). Before doing so, let us explore the functions ( ;  label , * label ) and ( ;  label , * label ) for label ∈ {LL, NN, LN} which are shown in Figs. 7-9.
First and foremost, all three models capture the jump response in the torque that is experimentally observed at the critical time 0 as a consequence of the step change of the angular velocity. The model  LL fits the experimental data almost perfectly with only a slight overshoot of the torque data in the vicinity of the critical time 0 . There is no such overshoot in case of the model  NN which has the shear thinning property in contrast to  LL , see Table 2. On the other hand, we clearly see that the model  NN undershoots the torque data in the second time segment during the relaxation stage. Finally, the model  LN combines the properties of both previous models and fits the experimental data the best. Remark 6.3. Using the alternative cost functional introduced in Remark 5.1, we end up with qualitatively identical results. 9 We intentionally choose measurements characterising the variance in the data, in order to mitigate the risk of fitting data with higher experimental error.

Sensitivity analysis
In this section, we study the influence of the radius of the sample on the expected output using the Monte-Carlo approach, see e.g. Caflisch [35]. The experimental data shown in Fig. 10 represent the average data from Fig. 5. The mean value of the recorded deflection angle from Fig. 4 was used to get the input for the simulation. The results support our previous statement regarding the observed inconsistency in the measured torque, see Fig. 5 and the discussion at the end of Section 4.3. Indeed, we see that a relatively small change of the radius of the sample, which is still assumed to coincide with the radius of the plates, leads to a variance in the simulated torque that is both qualitatively and quantitatively comparable with the variance in the experimental data.

Tests with modified configurations
The aim of this section is to illustrate the capability of the models to make predictions in configurations other than the one used for the parameter fitting, see Fig. 6. These so-called testing configurations are captured in Fig. 11.
The testing scheme is the same as in the previous case, but different times for the ramping of the deflection angle were used in case of measurements C1 and D1. On the other hand, in the remaining case we used exactly the same setting as in the fitting configuration. However, the ramping was done in a completely different -but also repeatablepath in this particular case. Indeed, looking at Fig. 11, we see that the upper plate rotated fast at the beginning of the measurement E1 (for ≪ 0 ) and slowed down as the time increased. This is in contrast with all other cases, where the rotation was gradually accelerated. Remark 6.4. The testing configurations are in the same rotational rheometer geometry as the one used in the parameter fitting. To test the full predictive power of the model it would be prudent for future work to test in a different geometry, e.g. capillary rheometer.

Table 3
Output of TPE algorithm subsequently used as the initial guess in the restarted Nelder-Mead algorithm, and the final optimised parameter values for models listed in Table 1 Table 3.  Table 3.  Table 3.
Comparison of the simulation results with the experimental data for individual models is shown in Figs. 12-14. Recall that the computation of ( ; , ) and ( ; , ) is based on the numerical solution of the system of ODEs (28) that was derived under certain simplifying assumptions. The acceleration observed at the start-up of measurement E1 makes some of these assumptions irrelevant. Moreover, the choice of the initial conditions (29) is not suitable in this specific case. In order to make them acceptable, we allow a fast ramping of the angular velocity to a specified value within the first few time steps. Due to these inaccuracies one cannot expect an ideal fit of the experimental data for this particular measurement in the time segment (0, 0 ]. Let us conclude with the observation that the model  LN provides us with a satisfactory prediction of the response of the SBR in different torsional flow configurations.

Conclusions
By performing the experiment described in Section 4, we have observed that the response of the uncured SBR is qualitatively comparable with the response of asphalt binders, especially in the context of the work by Narayan et al. [30]. Such an observation led us to the idea to describe the complicated nonlinear viscoelastic response of SBR by one of the models recently discussed by Málek et al. [11,12]. We compared two existing models, here denoted by  LL and  NN , with the newly proposed variant  LN (the derivation is sketched in Section 2.3), see Table 1.
The first two models showed a relatively good match between the simulated data and the experimental data, with only a few imperfections especially in the torque data. Possible reasons which make the two models imperfect in the given configuration may be related to the following facts. First,  LL cannot describe shear thinning behaviour which was experimentally proven for SBR, see e.g. Mourniac et al. [36]. Second,  NN describes the material as being ''fully incompressible'' which may be inappropriate in case of SBR. On the other hand,  NN captures shear thinning and the derivation of  LL is based on the assumption that the material stores the elastic energy as a compressible neo-Hookean solid. The new model combines physical properties of both in an appropriate way. In particular, it can describe shear thinning behaviour and a part of the elastic energy is stored as in the case of a compressible neo-Hookean solid, see (15). Probably this makes  LN superior to the other two models in predicting the behaviour of SBR in the chosen fitting and testing configurations.
For models investigated in this work, we have observed that step changes (discontinuous jumps) in the input angular velocity lead to the corresponding step changes in the simulated torque, see Figs. 6-9 and 11-14. It seems to plausibly reflect the experimentally observed behaviour. Rigorous analysis of the response of nonlinear materials to step inputs is feasible within the framework of Colombeau algebra, see Řehoř et al. [37] and Průša et al. [38].
Another line for future investigation would be to investigate the validity of the new model and the fitted parameters in a completely different geometrical configuration, e.g. a capillary rheometer. Furthermore, only unfilled SBR at a fixed temperature was studied in this paper. The influence of different fillers in SBR compounds as well as the influence of the temperature on model parameters should be subject to further intensive research.