Masataka Tanaka
email: dtanaka@gipwc.shinshu-u.ac.jp

Wen Chen
email: chenw@homer.shinshu-u.ac.jp

Dual reciprocity BEM applied to transient elastodynamic problems with dierential quadrature method in time

Keywords: Elastodynamics, Dual reciprocity BEM, Dierential quadrature method, Time-stepping integration, Lyapunov matrix equation

This paper is concerned with numerical analysis of the transient elastodynamic problems by the dual reciprocity BEM (DRBEM) in space combined with the dierential quadrature method (DQM) in time. Emphasis is placed on a comparative study of various timemarching schemes. Three numerical examples considered are for the longitudinal forced vibration of plates and the transverse free vibration of membranes. To authorsÕ best knowledge, this is the ®rst attempt to apply the DQM to the second-order time derivative in the DRBEM elastodynamic formulations. A recent approach using boundary conditions in the DQM was here extended to handle the initial conditions of elastodynamic problems. The resulting algebraic formulation is a Lyapunov matrix equation, which can be very eciently solved by the Bartels±Stewart algorithm. It is revealed that the DQM is an unconditionally stable algorithm and gives much better accuracy than the standard ®nite dierence schemes such as the Wilson h, Newmark and Houbolt methods, for the same time step size for the cases considered.

Introduction

Recently, the introduction of the dual reciprocity boundary element method (DRBEM) [1±3] has been bringing a major breakthrough in the BEM analysis of time-dependent problems. The DRBEM has been widely recognized as a technique which can retain the boundary-only merit in the BEM analysis of general nonlinear and linear time-dependent problems [START_REF] Partridge | The Dual Reciprocity Boundary Element Method[END_REF]. Nowadays, the application of the DRBEM to various transient problems has been a subject of growing interest. The boundary integral equations for the DRBEM are dependent only on geometrical data and free of interior cells. The resulting DRBEM formulation of initial-boundary value problems is therefore expressed in the standard form of ordinary differential equations of initial value problems, which can easily be solved by the mature time integrators.

In conjunction with the DRBEM discretization in space, the predominant numerical procedures currently used for time integration in analysis of elastodynamic problems are ®nite dierence approximations such as the Newmark, Houbolt and Wilson methods, which is similar to the situation in the FEM. In compound sources of errors due to approximate representations of spatial and time derivatives, it is seen that the accuracy and computing eciency of the resulting solutions depend greatly on the proper choice of time-marching schemes. Some pioneer works on the DRBEM analysis of elastodynamics are due to Nardini and Brebbia [START_REF] Nardini | A new approach to free vibration analysis using boundary elements[END_REF] and Loeer and Mansur [START_REF] Loeer | Analysis of time integration schemes for boundary element applications to transient wave propagation problems[END_REF]. They have concluded that the Houbolt method is preferred in the DRBEM since the arti®cial damping inherent with this method can eectively depresses in¯uence of higher modes in the response. Kontoni and Beskos [5] also followed this idea to choose the Houbolt method in the DRBEM analysis of dynamic elastoplastic problems. However, this conclusion is not in agreement with the common sense in the FEM dynamics analysis that the Newmark method is normally preferred to other methods [START_REF] Kim | From the trapezoidal rule to higher-order accurate and unconditionally stable time-integration method for structural dynamics[END_REF]. There exist some de®nite shortcomings in the Houbolt method; it is recognized that high numerical damping in the Houbolt method often impairs the accuracy of the solutions due to undesirable amplitude attenuation if large time step is used, especially for structural dynamic problems in which the response is often dominated by the low-frequency components of the system [START_REF] Subbaraj | A survey of direct time-integration methods in computational structural dynamics ± II. Implicit methods[END_REF]. Furthermore, Hilber and Hughes [START_REF] Hilber | Collocation, dissipation and ÔovershootÕ for time integration schemes in structural dynamics[END_REF] have pointed out that the self-starting is one of essentially desirable properties of a competitive numerical integrator for most elastodynamic problems, while the Houbolt method necessarily requires a distinct starting procedure and fails to satisfy this condition. The time integrators of such type not only increase the programming labor but also may cause some complexity of computing. Numerical experiment which preferred to the Houbolt method in [START_REF] Partridge | The Dual Reciprocity Boundary Element Method[END_REF][START_REF] Loeer | Analysis of time integration schemes for boundary element applications to transient wave propagation problems[END_REF] was an elastic wave propagation problem in which the contribution of high frequency structural modes to the response is important. However, in the cases of such type, the explicit algorithms are very ecient and rather frequently used in practice [START_REF] Subbaraj | A survey of direct time-integration methods in computational structural dynamics ± II. Implicit methods[END_REF]. In addition, the time step size used in [START_REF] Partridge | The Dual Reciprocity Boundary Element Method[END_REF][START_REF] Loeer | Analysis of time integration schemes for boundary element applications to transient wave propagation problems[END_REF] was very small, and therefore, those experiments cannot adequately provide an explicit comparison among a variety of integrators. It is clear from the preceding discussions that there is still a strong need to further investigate various numerical integrators in the DRBEM elastodynamic analysis.

On the other hand, Singh and Kalra [START_REF] Singh | Time integration in the dual reciprocity boundary element analysis of transient diusion[END_REF] presented a comprehensive comparative study on the time integrators in the context of the DRBEM formulation of transient diusion problems involving only the ®rstorder time derivative. They concluded that the one step least squares algorithm was the most accurate and ecient numerical integrators among all the ones assessed. However, for the problems involving the Dirichlet boundary condition, all time integration methods used in [START_REF] Singh | Time integration in the dual reciprocity boundary element analysis of transient diusion[END_REF] encountered a sharp drop in accuracy and eciency of computation. Very recently, the present authors [START_REF] Tanaka | Analysis of transient diusion problems by a combined use of dual reciprocity BEM and dierential quadrature method[END_REF] employed the dierential quadrature method (DQM) to approximate the ®rst-order time derivative in the DRBEM analysis of the same transient diusion problems investigated in [START_REF] Singh | Time integration in the dual reciprocity boundary element analysis of transient diusion[END_REF]. It was found that the obtained solutions of the Dirichlet problems are very accurate with comparatively much less computing eort than the other existing methods. This work has also revealed that the ineciency of the DRBEM analysis of the Dirichlet problems reported in [START_REF] Singh | Time integration in the dual reciprocity boundary element analysis of transient diusion[END_REF] are due to the numerical integrator in time rather than the DRBEM discretization in space.

The success in applying the DQM for the diusion problems has been encouraging the present authors to further extend the method along with the DRBEM to the elastodynamic problems with the second-order time derivatives. This paper also places emphasis on a comprehensive comparative study of various standard time-marching schemes in the context of the DRBEM formulation of elastodynamic problems. The in-plane forced vibration of plates and the transverse free vibration of membranes are taken as numerical examples because their analytical solutions are easily obtainable.

The DQM [START_REF] Bert | Dierential quadrature method in computational mechanics: a review[END_REF][START_REF] Chen | Dierential quadrature method and its applications in engineering ± applying special matrix product to nonlinear computations[END_REF] may not be well-known in computational community due to its recent origin. The method can be regarded as the ``direct approach'' of the traditional collocation method in the sense that the governing equations are analogized in terms of physical variables instead of usually ®ctitious expansion coecients. Some details on this method are presented. It has been well known that the DQM as well as the collocation method shows the exponential rate of convergence. The weakness of this method is its ¯exibility for irregular geometry. Since the time domain has no such diculty, the strengths of higher-order accuracy of the DQM can be fully exploited in approximating the time derivatives. In addition, the DQM also holds the unconditionally stable feature for accuracy of order more than two. This is due to the fact that the method is not a traditional time step scheme and circumvents the rigorous limitation of accuracy for unconditionally stable algorithms due to Dahlquist theorem [START_REF] Dahlquist | Special stability problem for linear multistep methods[END_REF]. On the other hand, the DQM advances progressively in time domain element-wisely from the initial state, and thus keeps the simplicity and ¯exibility of the standard time step methods. The resulting set of algebraic equations is in fact the Lyapunov matrix equation and can be very eciently solved by the well-established Bartels±Stewart algorithm [START_REF] Bartels | A solution of the equation AX XB C[END_REF]. The total computing eort required in the DQM is comparable to the common implicit step methods. Following the ideas of recent work by Wang and Bert [START_REF] Wang | A new approach in applying dierential quadrature to static and free vibrational analyses of beams and plates[END_REF] for boundary value problems, this paper also develops a simple technique to exactly satisfy all the initial conditions in the DQM analogue of the dynamic systems with the second-order time derivatives.

Vibration of plates and membranes

The governing dierential equation for the longitudinal vibration of plates and the transverse vibration of membranes can be expressed as

r 2 uxY t 1 c 2 o 2 uxY t ot 2 Y x P XY 1
subject to the initial conditions:

uxY 0 u 0 xY 2a uxY 0 v 0 xY 2b
and the displacement and traction boundary conditions:

uxY t uxY tY x & C u Y 3
T xY t T xY tY x & C T Y 4
where T ouaon, in which n is the unit outward normal. It is assumed that the domain X P R 2 is bounded by a piece-wise smooth boundary C C u C T . In this study, the three typical vibration problems of plates and membranes with analytical solutions available in the literature are considered to evaluate the performances of the DQM, Newmark, Houbolt and Wilson methods in conjunction with the DRBEM spatial discretizations. The spatial variable domains of the test problems are square. It is noted that the analytical solutions of the longitudinal vibration of square plates considered are in agreement with those of one-dimensional rod cases [START_REF] Partridge | The Dual Reciprocity Boundary Element Method[END_REF][START_REF] Mansur | Further development on the solution of the transient scalar wave equation[END_REF].

Vibration of plate subjected to periodic in-plane force

The initial conditions of this case are assumed as uxY 0 0Y 5a uxY 0 0Y 5b and the boundary conditions are speci®ed as

uxY t 0Y x 1 1Y 6a
T xY t 0Y x 1 0Y x 2 0X1X 6b
When a periodic in-plane force

p P cos 3p 8a E q s t 7
is applied at the left side of plate as shown in Fig. 1, we can derive the analytical solution [START_REF] Timoshenko | Vibration Problems in Engineering[END_REF][START_REF] Kondou | Theory of Vibration[END_REF]:

u 8 p 2 Pa EA I i1 À1 iÀ1 2i À 1 2 À 9a16 3 sin 2i À 1 2a px 1 cos 3p 8a E q s t 2 À cos 2i À 1p 2a E q s t 3 Y 8
where a denotes the length of the plate edges; A is area of cross-section; E is the Young modulus and p the density of the plate material. The traction can be obtained by

T EA ou ox X 9

Free vibration of membranes

It is assumed that a square membrane of a  b is released from rest in an initial position and velocity, i.e.,

uxY 0 ax 1 À À x 2 1 Á bx 2 À À x 2 2 Á Y 10a uxY 0 0X 10b
The boundary conditions are given by uxY t 0 at x 1 ÀaY a and x 2 ÀbY bX 11

By using formulas given in [START_REF] Sneddon | Elements of Partial Dierential Equations[END_REF], the analytical solution for this case is given by

uxY t I mYn1 16 2 a 2 b 2 p 6 1 À À1 m m 3 1 À À1 n n 3 sin mpx 1 2a sin npx 2 2b cos k mn ctY 12
where c T q s Y k 2 mn p 2 m 2 4a 2 n 2 4b 2 Y
T and q are tension and density, respectively.

Vibration of plate subjected to Heaviside-type impact load

In this case, the initial and boundary conditions are the same as Eqs. (5a)±(6b) of example 1. An impact load of Heaviside-type p P Y t P 0 13 is enforced at the left side of plate as shown in Fig. 1. The analytical solution is given by [START_REF] Kondou | Theory of Vibration[END_REF]

u 8 p 2 Pa EA I i1 À1 iÀ1 2i À 1 2 sin 2i À 1px 1 2a 1 À cos 2i À 1p 2a E q s t 3 X 14
The traction can be evaluated by Eq. ( 9).

DRBEM discretization in space

The governing equation ( 1) can be weighted by the fundamental solution u à of Laplace operator as follows:

X 1 c 2 o 2 u ot 2 À r 2 u u à dX 0X 15
Applying GreenÕs second identity to Eq. ( 15) yields

d i u i C T à u Àu à T dC À X 1 c 2 o 2 u ot 2 u à dXY 16
where subscript i denotes the source point, T à ou à aon, and d i X dfY x dXX The dual reciprocity method transforms the domain integral on the right-hand side of Eq. ( 16) by means of a set of coordinate function f j x:

uxY t % N L j1 f j x a j tY 17
where the superimposed dot represents the time derivative, a j are unknown functions of time, and N and L are the numbers of the boundary and selected internal nodes, respectively. The coordinate functions used in this paper were presented by Wrobel and Brebbia [START_REF] Wrobel | The dual reciprocity boundary element formulation for nonlinear diusion problems[END_REF]. These functions are also linked with w j x through r 2 w j f j X 18

Therefore, we have

X o 2 u ot 2 u à dX N L j1 X u à r 2 w j dXX 19
Eq. ( 15) can ®nally be reduced to

d i u i C T à u Àu à T dC N L j1 d i w j i C T à w j À À u à g j Á dC ! a c 2 Y 20
where g j ow j aon. Note that w j and f j are known functions. The resulting DRBEM formulation for the present transient elastodynamic problems is given by

M u Hu À GT 0Y 21
where M GE À HWF À1 is the mass matrix; H and G denote the whole matrices of boundary element with kernels T à and u à , respectively. F, W and E are comprised of the coordinate function column vectors f j , w j and g j . The discretization procedure in detail can be found in [START_REF] Partridge | The Dual Reciprocity Boundary Element Method[END_REF].

If displacement boundary conditions are involved in the problems, Eq. ( 21) is a dierential algebraic system. By using an approach of matrix partition [START_REF] Nardini | A new approach to free vibration analysis using boundary elements[END_REF][START_REF] Wrobel | The dual reciprocity boundary element formulation for nonlinear diusion problems[END_REF]

u u 1 u 2 & ' and T T 1 T 2 & ' Y 22
where u and T are divided into two parts corresponding to C u and C T parts of the boundary, we can in general have

M u 2 Hu 2 Mu 1 Hu 1 GT 2 X 23
In particular, we have

M u 2 Hu 2 GT 2 24
for forced vibrations when only external traction is applied, which is the case for examples 1 and 3 in this study. Note that all excitation sources are placed in the right-hand side of Eq. [START_REF] Strize | Harmonic dierential method and applications to structural components[END_REF]. All the coecient matrices of these equations are dependent only on the geometric data of the problem. The remaining solution procedure is the same as the treatment of the standard initial value problems. The desired traction can be easily calculated after the solutions of the above dierential system are accomplished. The linear element (DC 0X1) was employed in the present DRBEM discretization as shown Figs. 1 and2. For examples 1 and 3, one internal point is placed in the center of domain, as shown in Fig. 1. For example 2, 33 internal points are used due to the fact that a homogeneous boundary condition u 0 is applied at all the boundary (Fig. 2).

DQM approximation in time

The DQM is in fact a variant of the standard collocation methods. The advantages of the DQM over the latter were well established in [START_REF] Bert | Dierential quadrature method in computational mechanics: a review[END_REF][START_REF] Quan | New insights in solving distributed system equations by the quadrature methods ± II[END_REF]. First, the practical physical values are directly computed in the DQM, whereas the collocation methods use the indirect expansion (spectral) variables. This greatly simpli®es implementations and manifests the DQM in easy-to-choose starting solutions of nonlinear iterations. It is noted that the ®ctitious expansion variables in the collocation methods usually have not physical meanings and are therefore dicult to choose initial iterative solutions for nonlinear problems [START_REF] Burka | Solution of sti ordinary dierential equations by decomposition and orthogonal collocation[END_REF]. Second, the DQM provides more ¯exibility to choose grid points to enhance the rate of convergence [START_REF] Quan | New insights in solving distributed system equations by the quadrature methods ± II[END_REF]. For more details of the method see [START_REF] Bert | Dierential quadrature method in computational mechanics: a review[END_REF] and references cited therein.

The DQM analogue of the ®rst-and second-order derivatives of function f t can be expressed as

df t dt ti N j1 A ij f t j Y 25a
d 2 f t dt 2 ti N j1 B ij f t j Y i 1Y 2Y F F F Y N Y 25b
where t j Õs are the discrete points in the temporal variable domain. f t j is the function values at these points, A ij and B ij are the related DQM weighting coecients for the ®rst-and second-order derivatives, respectively. In the present study, the Chebyshev±Gauss±Lobatto collocation points are used in each time element of the DQM, namely,

t i s 2 1 À cos i À 1 N À 1 p ! Y i 1Y 2Y F F F Y N Y 26
where s denotes the length of DQM time element and N is the number of grid points.

The DRBEM formulations ( 23) and ( 24) can be expressed as the standard form

u Ku f Y 27
where K M À1 H is stiness matrix. By using simple algebraic transformation

z u À u 0 À v 0 t v 0 t 0 Y 28
where t 0 is the initial instance of each DQM time element, Eq. ( 27) is restated as

z Kz FY 29
where F f À Ku 0 À Kv 0 t Kv 0 t 0 , and the respective two initial conditions (2a) and (2b) are set zero, namely,

zj tt 0 0Y 30a zj tt 0 0X 30b
The above transformation is a key step to apply initial conditions exactly in the DQM approximation of the second-order time derivative. The DQM analogue of the ®rst-order time derivative can be stated as

Afzg f zgY 31a

according to the initial condition (30a), and Af zg f zg 31b

according to the initial condition (30b), where A is yielded by removing the ®rst column of the original DQM weighting coecient matrix A in Eq. (25a). Substituting Eq. (31a) into Eq. (31b), we have

AAfzg Bfzg f zgY 32

where the modi®ed DQM coecient matrix B is N À 1 Â N À 1 dimension. It is stressed that the initial conditions speci®ed in Eqs. (30a) and (30b) have been built into the modi®ed coecient matrix B. The above scheme is an analogy with the recently developed technique in applying boundary conditions for the DQM solution of high-order boundary value problems presented by Wang and Bert [START_REF] Wang | A new approach in applying dierential quadrature to static and free vibrational analyses of beams and plates[END_REF]. In terms of approximate formula (32), Eq. ( 29) can be analogized as

ZB T KZ FtY 33
where B T is a transpose of modi®ed DQM weighting coecient matrix B with the inclusion of two initial conditions given in Eqs. (30a) and (30b). It is noted that Z in Eq. ( 33) is a rectangular matrix rather than a vector. Therefore, Eq. ( 33) is a Lyapunov matrix equation. The present DQM procedure discretizes the time variable in the element-wise way, which is somehow dierent from the standard step-wise solution of single step or multi-step integration. In other words, the present methodology is to advance element by element and the multiple grid points are employed in each time element to keep high accuracy of solution. It is worth stressing that the DQM advances progressively in time and thus can maintain the simplicity and ¯exibility of the common time step methods.

Solver of Lyapunov matrix equations

It is easy to transform the resulting algebraic formulation of the Lyapunov equation (33) into a standard form of simultaneous algebraic equations applicable to being solved by the LU decomposition method. However, such procedure will fail to fully utilize the special structure inherent with the Lyapunov equation so that the computing eort is not necessarily high. In the present study, the so-called Bartels±Stewart algorithm [START_REF] Bartels | A solution of the equation AX XB C[END_REF] was utilized to solve the Lyapunov algebraic matrix equation (34). The performances of this method are very ecient, stable and accurate. The solution procedures include the following four steps:

Step 1. Reduce K and B T of Eq. (34) into certain simple form via the similarity transformations

G P À1 KP and R V À1 B T VX
Step 2. Q P À1 FV for the solution of Q.

Step 3. Solve the transformed equation GY YR Q for Y.

Step 4. Z PYV À1 . The time-consuming calculation of OM 3 scalar multiplication is required only in Step 1 for Eq. ( 29) of M dimension, while all implicit step methods also demand OM 3 operations using LU decomposition. Moreover, operation in Step 1 need be done only once. On the other hand, Steps 2±4 of the Bartels±Stewart algorithm need be executed repeatedly in each time element, which requires OM 2 multiplication. The standard step methods also demand analogous computational effort in each time step. Therefore, the present DQM scheme is comparable in computing effort to that of the normal step implicit methods.

A-stability

It is centrally important whether or not an algorithm is stable in the integration solution of ordinary dierential equations. It is known that the collocation method is A-stable [START_REF] Burka | Solution of sti ordinary dierential equations by decomposition and orthogonal collocation[END_REF] which is, in the terminology of computational structural dynamics, unconditionally stable. Therefore, the DQM is also unconditionally stable due to the actual equivalence to the collocation method.

Error estimation and accuracy

The accuracy of the algorithms is of vital importance in computing eciency and closely related to the truncation error. The error estimator of the DQM approximation of the second-order derivative of function f t is given by Chen [START_REF] Chen | Dierential quadrature method and its applications in engineering ± applying special matrix product to nonlinear computations[END_REF]

jR i j 6 K Á err i Á Dt N À1 Y i 1Y 2Y F F F Y N Y 34
where K maxfjf N njY jn i f N 1 njgY err i denotes the error constants dependent on grid spacing and can be obtained easily, N the number of grid points in the DQM time element and Dt is time step size.

According to formula (34), the accuracy of the DQM is ODt N À1 , for example, ten-order of accuracy when N 11. Due to Dahlquist [START_REF] Dahlquist | Special stability problem for linear multistep methods[END_REF], there is no third-order accurate unconditionally stable linear multistep method, and the maximum order of accuracy in the step methods up to two in order to preserve the Astability. The DQM is not a traditional multi-step algorithm and therefore circumvents this rigorous limitation of solution accuracy. The DQM can produce accurate solutions by using larger time step, while still attaining the desirable A-stability merits.

Results and discussion

In this section, the numerical results of three examples in Section 2 are provided and discussed based a performance comparison of the DQM, Newmark, Houbolt and Wilson methods in conjunction with the DRBEM spatial discretization. Parameters a 0X25 and d 0X5 are taken in the Newmark method and h 1X4 in the Wilson method as in [START_REF] Bathe | Numerical Methods in Finite Element Analysis[END_REF]. The coordinates of the displacement, traction and time are dimensionless as uE/Pa and TA/P for plate, and ct where c 1aa Eaq p for plates and c T aq p for membrane. In this study,

Dt c s N Y 35
where s is the length of the DQM time element, N denotes the number of grid points in the time element and is taken as 11. Dt represents the average two-point spacing, which is analogous to the time step size in the usual time step methods. Figs. 3 and4 display the solutions of time±displacement curves at point A of Fig. 1 for the longitudinal vibration of square plate subjected to a periodic in-plane force. It is found from Fig. 3 that all methods provide the exact results using suciently small time step Dt 0X1. It is also noted that various ®nite dierence methods confront some small attenuation of amplitude and overshoot for long-term response, while DQ method always gives the very accurate solutions. On the other hand, this also reveals high degree of accuracy of the DRBEM spatial discretization. When a larger time step D 0X5 is employed, Fig. 4 illustrates quite distinct performances of various dierent time-marching schemes. It is observed that the DQM produces strikingly much more accurate solutions than all other time schemes. This demonstrates the superb converging rate and accuracy of the DQM. The Wilson method is found obvious overshoot. In fact, the solutions in long-term response by the various ®nite dierence schemes using the large time step Dt 0X5 is not acceptable in engineering accuracy for this case.

The time±traction response curves at ®xed point B of Fig. 1 for the same example 1 are shown in Figs. 5 and6. As in the displacement situations, the response by all these methods closely traces the analytical solution curve for the small time step size Dt 0X1 as shown in Fig. 5. However, as the time step becomes larger Dt 0X5, the solutions of the traditional ®nite dierence schemes including the Newmark, Houbolt and Wilson methods encounter a sharp drop of accuracy for the long-term response as shown in Fig. 6, while the DQM results always remain in very good agreement with analytical solutions. The Wilson solution is found a strong oscillation at the initial response and evident overshoot in long-term response. Also, all ®nite dierence schemes produce a manifest phase shift.

Figs. 7 and 8 illustrate the analytical and numerical displacement±time curves at central point C of free vibration of a square membrane. As is expected, all methods produce the exact solutions using the small Fig. 3. Displacement curve at point A of a square plate subjected to a periodic in-plane force Dt 0X1. time step Dt 0X1 shown in Fig. 7. Among them, the Houbolt solutions have visible phase shift. Signi®cant dierences in numerical accuracy, amplitude attenuation, and phase shift are easily observed from Fig. 8 when using larger time step Dt 0X5. Except for the DQM, all other methods have a great loss of accuracy and evident phase shift. The amplitude attenuation is clearly noted in the Houbolt and Wilson methods, especially in the former due to the excessive arti®cial damping. To provide more insights into the distinct performances of these methods under coarse time step, Fig. 9 also plotted the normal derivative curves at the boundary point D in Fig. 2. using Dt 0X5. Again, the DQM shows the high accuracy, while the Newmark, Houbolt and Wilson methods have a great drop in the solution accuracy and apparent phase shift. In particular, undesirable high numerical damping of the Houbolt method is also observed from Fig. 8.

In conclusion of the above two examples, there are essential dierences in the accuracy, amplitude attenuation and phase shift behaviors between the DQM and the standard dierence schemes if large time step is chosen for computing economy. This is easily explained by the fact that the DQM is characterized as high accuracy and rate of convergence we referred to din Section 4. It is commonly known that for periodic system analysis, the major factors aecting accuracy of a given method should be truncation error, can yield the exact solutions of the time±displacement curve. However, as was seen in the previous examples 1 and 2, Fig. 11 shows that the Newmark, Houbolt and Wilson methods confront an obvious fall in solution accuracy when using large Dt 0X5. Also, an evident amplitude attenuation and phase shift is observed in the Houbolt method. The DQM remains the exact solutions for Dt 0X5 and behaves very well over any other methods.

High-order modes have a strong eect on the traction response behavior of example 3. This fact is easily observed from Figs. 12±14. It is seen from Fig. 12 that the solutions of all methods have oscillations when very small step size Dt 0X05 is employed. Among them, the Houbolt method gave the best solutions with the smallest ¯uctuations. The high arti®cial damping of the Houbolt method becomes bene®cial in this case. This is one of the major reasons in [START_REF] Partridge | The Dual Reciprocity Boundary Element Method[END_REF][START_REF] Loeer | Analysis of time integration schemes for boundary element applications to transient wave propagation problems[END_REF] to conclude that the Houbolt method is preferred over the other ®nite dierence schemes. However, as mentioned previously, it was generally recommended that the explicit method is more ecient for this type of problems in which the high and intermediate frequency components have important aect on the response [START_REF] Subbaraj | A survey of direct time-integration methods in computational structural dynamics ± II. Implicit methods[END_REF]. The explicit methods are much cheaper and simpler than the implicit algorithms. The weakness in the explicit methods lies in that the numerical stability requires employing the very small time steps. In such case, the explicit methods are much more advantageous in applying the small time step size than the implicit Houbolt method. We con®ne our attentions in this paper within the implicit methods. For more details on explicit methods see [START_REF] Dokainish | A survey of direct time-integration methods in computational structural dynamics ± I. Explicit methods[END_REF].

To investigate the aect of step size on the oscillations of the traction solutions, Figs. 13 and 14 illustrate respectively the traction curves using time step Dt 0X1 and Dt 0X4. It is found from Fig. 13 that the Houbolt method loses some solution accuracy, while the oscillations of the DQM and Newmark method become somewhat weak when using time step Dt 0X1. The Wilson method emerges comparatively the best time scheme with Dt 0X1. Furthermore, Fig. 14 depicts the traction curves of the DQM and the Houbolt method. We can ®nd that the oscillation of the DQM solutions using Dt 0X4 is much less than using Dt 0X05 in the reference of Fig. 12. Mansur and Brebbia [START_REF] Mansur | Further development on the solution of the transient scalar wave equation[END_REF] also pointed out the evident relationship between time step size and traction solutions noise of this case, in which a dierent BEM technique rather than the DRBEM was used to handle domain integral. It is also noted here that the accuracy of the DQM using Dt 0X4 is almost the same as that of the Houbolt method using Dt 0X05, while the Houbolt method confronts a great loss of accuracy employing Dt 0X4. By comparing behaviors of various methods in Figs. 12±14, it is concluded that for the Houbolt method, the time step should necessarily be small enough to produce adequately accurate solutions. In contrast, the DQM method can smooth the solutions and alleviate the oscillation eect of high-order modes by using the larger time step. From the viewpoint of computational economy using coarse time step, the DQM appears the method of choice for this case.

The foregoing discussions indicate that, for the cases where response is primarily dominated by low and intermediate frequency modes, the DQ method exhibits an impressive advantage in the solution accuracy. For the systems on which high modes have important aect such as the traction analysis of example 3, the DQM is preferred if the large time step is employed and the Houbolt method is the best for the very small time step. However, an explicit algorithm should be considered at ®rst for such type of dynamic systems. It is conceivable that the DQM with other type of basis functions such as harmonic basis functions [START_REF] Strize | Harmonic dierential method and applications to structural components[END_REF] may yield better characteristics traction analysis of example 3. This is left for further study. Whether the numerical oscillation in computing traction of example 3 is partly due to the spatial discretization also requires further investigation.

Conclusions

Some important characteristics of numerical integrators for elastodynamics problems have been examined in detail through numerical experiments. It is found that although the Houbolt method seems predominant currently in the solution of the DRBEM formulation of elastodynamics problems [1,3±5], excessive numerical damping in the method can have a very detrimental eect on the accuracy of solutions if large time step is employed, and strengthens the case of applying an alternative method. Based on the present study, we conclude that the Newmark method should be in general preferred in the context of the DRBEM formulations of elastodynamics problems in comparison with the Wilson and Houbolt methods.

In this paper, we have applied the DQM to approximate temporal derivative in conjunction with the DRBEM spatial discretization for three typical elastodynamics problems. It is validated through theoretical analysis and numerical experiments that the DQM holds the desirable attributes of unconditional stability yet has higher-order of accuracy than the standard ®nite dierence integrators such as the Newmark, Houbolt and Wilson methods. An eective approach applying the initial conditions was developed in this paper by analogy with a recent work by Wang and Bert [START_REF] Wang | A new approach in applying dierential quadrature to static and free vibrational analyses of beams and plates[END_REF] in the DQM solution of high-order boundary value problems. The use of the Bartels±Stewart algorithm greatly reduces the computational eort of the DQM to comparable level of the normal implicit ®nite dierence methods. However, the special procedure of the Bartels±Stewart algorithm increases the complexity of the programming. The robustness and superior accuracy of the DQM over the traditional ®nite dierence schemes are clearly observed by comparing the numerical results of three examples. The DQM appears to be a promising technique in practical engineering computations. Further evaluation of this method should be bene®cial.

How to satisfy all initial conditions exactly is a key step to successfully implement the DQM to analogize the second-order derivatives in time. There exist two competitive approaches employing multiple boundary conditions in the DQM solution of high-order boundary problems [START_REF] Bert | Dierential quadrature method in computational mechanics: a review[END_REF][START_REF] Chen | Dierential quadrature method and its applications in engineering ± applying special matrix product to nonlinear computations[END_REF][START_REF] Wang | A new approach in applying dierential quadrature to static and free vibrational analyses of beams and plates[END_REF][START_REF] Shu | A generalized approach for implementing general boundary conditions in the GDQ free vibration analysis of plates[END_REF]. This paper follows the basic idea in [START_REF] Wang | A new approach in applying dierential quadrature to static and free vibrational analyses of beams and plates[END_REF] to incorporate initial conditions into the modi®ed DQM coecient matrix. It is still possible to develop a dierent approach by analogy with the strategy in [START_REF] Chen | Dierential quadrature method and its applications in engineering ± applying special matrix product to nonlinear computations[END_REF][START_REF] Shu | A generalized approach for implementing general boundary conditions in the GDQ free vibration analysis of plates[END_REF]. In addition, the choice of grid spacing in each DQM time element should have some eect on the solution accuracy as in the DQM solution of boundary value problems. The study of the above these problems are now under way.

Fig. 1 .

 1 Fig. 1. Linear elements and interior points of a square plate.

Fig. 2 .

 2 Fig. 2. Linear elements and interior points of a square membrane.

Fig. 4 .

 4 Fig. 4. Displacement curve at point A of a square plate subjected to a periodic in-plane force Dt 0X5.

Fig. 5 .

 5 Fig. 5. Traction curve at point B of a square plate subjected to a periodic in-plane force Dt 0X1.

Fig. 6 .

 6 Fig. 6. . Traction curve at point B of a square plate subjected to a periodic in-plane force Dt 0X5.

Fig. 7 .

 7 Fig. 7. Displacement curves at center point C of free vibration of a square membrane with initial displacement Dt 0X1.

Fig. 8 .

 8 Fig. 8. Displacement curves at center point C of free vibration of a square membrane with initial displacement Dt 0X5.

Fig. 9 .

 9 Fig. 9. Normal derivative curve at edge point D of free vibration of a square membrane with initial displacement Dt 0X5.

Fig. 10 .

 10 Fig. 10. Displacement curve at point A of a square plate subjected to a Heaviside-type impact Dt 0X1.

Fig. 11 .

 11 Fig. 11. Displacement curve at point A of a square plate subjected to a Heaviside-type impact Dt 0X5.

Fig. 12 .

 12 Fig. 12. Traction curve at point B of a square plate subjected to a Heaviside-type impact Dt 0X05.

Fig. 13 .

 13 Fig. 13. Traction curve at point B of a square plate subjected to a Heaviside-type impact Dt 0X1.

Fig. 14 .

 14 Fig. 14. Traction curve at point B of a square plate subjected to a Heaviside-type impact.

Acknowledgements

This work was carried out as a part of the research program supported by the Japan Society for the Promotion of Science. Additional ®nancial support was provided as Grant-in-Aid for JSPS fellows by the Ministry of Education, Science, Sports and Culture, Japan.


P Wapperom
email: wapperom@math.vt.edu

A Leygue

R Keunings

Numerical simulation of large amplitude oscillatory shear of a high-density polyethylene melt using the MSF model

Keywords: LAOS, Integral MSF model, Deformation field method, Linear polymer melts

Numerical simulation of large amplitude oscillatory shear of a high-density polyethylene melt using the MSF model

P. Wapperom, Adrien Leygue, Roland Keunings

Introduction

Molecular models have become increasingly popular to describe the complex rheological behavior of entangled polymeric liquids. Practically all recently developed molecular models are extensions of the Doi-Edwards theory. Although the basic Doi-Edwards (DE) model [START_REF] Doi | The Theory of Polymer Dynamics[END_REF] can successfully predict the damping function and the plateau modulus of linear viscosity, it has some important deficiencies like excessive shear thinning in fast shearing flows. Recent extensions of the Doi-Edwards theory have alleviated these shortcomings. We consider one of such extensions, the molecular stress function (MSF) model which has recently been proposed by Wagner et al. [START_REF] Wagner | The molecular stress function model for polydisperse polymer melts with dissipative convective constraint release[END_REF]. The non-linear response is captured through the introduction of an extra evolution equation for the molecular stress function. The model has two additional material parameters, one to describe extensional and one to describe shear flows. With only two non-linear parame-ters, the MSF model is able to accurately predict the non-linear response in start up of shear and extension for a commercial linear high-density and branched low-density polyethylene [START_REF] Wagner | The molecular stress function model for polydisperse polymer melts with dissipative convective constraint release[END_REF].

The linear spectrum of polymer melts is normally determined in small amplitude oscillatory shear flow. At low strains, the material response is approximately linear. At larger strains, however, non-linear effects do play an important role. Early large amplitude oscillatory shear (LAOS) experiments date already from almost 50 years ago [START_REF] Philippoff | Vibrational measurements with large amplitudes[END_REF]. A recent overview of LAOS has been given by Giacomin and Dealy [START_REF] Giacomin | Using large-amplitude oscillatory shear[END_REF]. Wilhelm et al. [START_REF] Wilhelm | Fourier-transform rheology[END_REF][START_REF] Wilhelm | High sensitivity Fouriertransform rheology[END_REF] introduced high-sensitivity Fourier-transform rheology and analysed large amplitude step shear oscillations of polymer melts [START_REF] Van Dusschoten | Two-dimensional Fourier transform rheology[END_REF]. Very recently, Debbaut and Burhin [START_REF] Debbaut | Large amplitude oscillatory shear and Fourier transform rheology for a high-density polyethylene: experiments and numerical simulation[END_REF] performed LAOS experiments on a commercial high-density polyethylene melt up to high strains of 10. Their simulations with a Giesekus model showed good agreement for moderate non-linear regimes and larger deviations for the most non-linear regimes that were experimentally achievable.

In this paper, we study the rheometrical response of the MSF model in large amplitude oscillatory shear flow. In particular, we investigate whether such models derived from molecular theory are able to accurately predict experimental data of a commercial linear polymer melt. In order to perform the simulations in an efficient and accurate manner, we modify the deformation field method of Hulsen et al. [START_REF] Hulsen | A new approach to the deformation fields method for solving complex flows using integral constitutive equations[END_REF] (Section 4). After validation of the numerical technique in Section 6, we show that the MSF model, having only one relevant material parameter in LAOS, can accurately predict experimental results for a linear highdensity polyethylene melt [START_REF] Debbaut | Large amplitude oscillatory shear and Fourier transform rheology for a high-density polyethylene: experiments and numerical simulation[END_REF]. At medium strains, the performance is comparable to that of a six-mode Giesekus fluid having six material parameters to describe the non-linear regime. At these strains the basic Doi-Edwards model already underpredicts the experimental data. At the higher strains, both the MSF and Giesekus model overpredict the experimental data, although the MSF model renders slightly to significantly better predictions depending on the value of the non-linear material parameter. Completely unexpectedly, it is the Doi-Edwards model that shows excellent agreement with the experimental data at the high strains. Then, the DE model, having no non-linear material parameters, outperforms both non-linear models.

Governing equations

For the MSF model the stress is related to the deformation history by

T (t) = 5 t -∞ m(t -t )f 2 t (t)Q [B t (t)] dt , (1)
where B t (t) is the Finger tensor which measures the deformation of a fluid particle at the current time t with respect to a reference time t . The memory function m assumes the classical multimode Maxwell form

m(t -t ) = i G i λ i t -∞
e -(t-t )/λ i , [START_REF] Wagner | The molecular stress function model for polydisperse polymer melts with dissipative convective constraint release[END_REF] where G i are the moduli and λ i are the relaxation times of the fluid.

The tensor Q is the strain measure for which the independent alignment approximation was used in the original derivation in [START_REF] Wagner | The molecular stress function model for polydisperse polymer melts with dissipative convective constraint release[END_REF]. Instead, we use the Currie approximation to the Doi-Edwards deformation tensor [START_REF] Currie | Constitutive equations for polymer melts predicted by the Doi-Edwards and Curtiss-Bird kinetic theory models[END_REF]. The orientation tensor Q is then directly related to the Finger strain B t and the Cauchy strain B -1 t by

Q = 1 J -1 B t - 1 (J -1)(I 2 + 3.25) 1/2 B -1 t . (3)
Here, J = I 1 + 2(I 2 + 3.25) 1/2 and I 1 and I 2 are the first and second invariants of B t , respectively. For large amplitude oscillatory shear flow this is a very good approximation to the Doi-Edwards tensor using the independent alignment approximation, as we show in Section 7.

The MSF model is completed by an evolution equation for the Finger tensor and the molecular stress function f t . The Finger tensor is governed by

DB t Dt = κ • B t + B t • κ T , ( 4
)
where κ is the transpose of the velocity gradient. For the MSF model, the evolution equation for the square of the molecular stress function takes different forms for linear and branched polymers. For linear polymers, we have

Df 2 t Dt = f 2 t κ : Q - 1 f 2 t -1 CR (5)
while for branched polymers the right-hand side is slightly modified [START_REF] Wagner | The molecular stress function model for polydisperse polymer melts with dissipative convective constraint release[END_REF][START_REF] Wagner | Quantitative assesssment of strain hardening of low-density polyethylene melts by the molecular stress function model[END_REF]. The dissipative constraint release CR is expressed as

CR = 1 2 (f 2 t -1) 2 a 1 A 2 1 : Q + a 2 |A 2 : Q -A 2 1 : Q| , ( 6
)
where A 2 1 and A 2 are second-order Rivlin-Erickson tensors which are related to the rate-of-deformation tensor d = (κ + κ T )/2 and rate-of-rotation tensor w = (κκ T )/2 by

A 2 1 = 4d 2 , A 2 = DA 1 Dt + A 2 1 + 2w • d + 2d • w T . ( 7
)
The parameters a 1 and a 2 in Eq. ( 6) are the only two non-linear parameters in the MSF model. The only parameter that is relevant in steady non-rotational flows is a 1 , since A 2 : Q -A 2 1 : Q vanishes in that case. The value of a 1 can be determined by fitting extensional flow data. Next, a 2 can be obtained from shear viscosity and first normal stress data.

The MSF theory is an extension of the basic Doi-Edwards tube theory [START_REF] Doi | The Theory of Polymer Dynamics[END_REF]. The difference between the models is the inclusion of the molecular stress function for the MSF model. To evaluate the improvement of this model in large amplitude oscillatory shear, we use the Doi-Edwards model for comparison. The governing equations are then Eqs. ( 1)-( 4) with f 2 t = 1 in the expression for the polymer stress. The DE model has no non-linear parameters and is fully characterized by the linear spectrum.

Large amplitude oscillatory shear flow

In large amplitude oscillatory shear flow a fluid is subject to a periodic shear deformation γ with amplitude γ 0 and frequency ν,

γ(t) = γ 0 sin(2πνt). ( 8
)
The deformation is applied for t > 0 and up to t = 0 the fluid is assumed at rest. The corresponding periodic shear rate γ equals

γ(t) = 2πνγ 0 cos(2πνt). (9)
With this shear rate, the polymeric stress is obtained as a function of time from the constitutive equations ( 1)-( 4) and the evolution equation for the molecular stress function for linear polymers, Eqs. ( 5)- [START_REF] Van Dusschoten | Two-dimensional Fourier transform rheology[END_REF]. Since the stress response attains a steady periodic state, a Fourier transform facilitates a detailed quantitative analysis. The shear stress is decomposed into an infinite sum of trigonometric functions,

T xy (t) = A 0 2 + ∞ n=1 A n sin(2nπνt) + ∞ n=1
B n cos(2nπνt), [START_REF] Currie | Constitutive equations for polymer melts predicted by the Doi-Edwards and Curtiss-Bird kinetic theory models[END_REF] where A n and B n are the Fourier coefficients. For a real signal, these coefficients are given by

A n = 2ν c+(1/ν) c T xy sin(2nπνt) dt, B n = 2ν c+(1/ν) c
T xy cos(2nπνt) dt [START_REF] Wagner | Quantitative assesssment of strain hardening of low-density polyethylene melts by the molecular stress function model[END_REF] for an arbitrary period ranging from c to c + (1/ν). The Fourier coefficient A 0 vanishes in view of the two-fold symmetry of the periodic shear stress signal. Standard numerical integration has been used to evaluate the integrals in Eq. [START_REF] Wagner | Quantitative assesssment of strain hardening of low-density polyethylene melts by the molecular stress function model[END_REF].

Recently, Debbaut and Burhin [START_REF] Debbaut | Large amplitude oscillatory shear and Fourier transform rheology for a high-density polyethylene: experiments and numerical simulation[END_REF] have described a new viscometric device to characterize polymer melts in large amplitude oscillatory shear flows. The experimental equipment consisted of an oscillatory device with a closed chamber to allow for higher frequencies. At a temperature of 170 • C, LAOS experiments have been performed at various frequencies and amplitudes for the commercially available high-density polyethylene melt Finathene ® 3802 YCF. Experimental data are available for the frequencies ν = 0.1, 0.3, 1, and 3 Hz and for the amplitudes γ 0 = 0.5, 1, 2.5, 5, and 10. For the higher frequencies, not all amplitudes can be achieved experimentally. For ν = 1 Hz, the highest achievable amplitude is γ 0 = 5, while for ν = 3 Hz this is γ 0 = 1. For further reference, we note that the maximum shear rates that occur at the highest amplitude for each frequency are γmax ≈ 6.28 s -1 for ν = 0.1 Hz, γmax ≈ 18.85 s -1 for ν = 0.3 and 3 Hz, and γmax ≈ 31.42 s -1 for ν = 1 Hz.

Numerical method for simulation in LAOS

The evolution equations for the Finger tensor B t are easily integrated analytically for large amplitude oscillatory shear, i.e. for the shear rate specified by Eq. [START_REF] Hulsen | A new approach to the deformation fields method for solving complex flows using integral constitutive equations[END_REF]. We obtain for the nonconstant components of the Finger tensor

B xy t (t) = γ 0 sin(2πνt) t ≤ 0 -γ 0 sin(2πνt ) + γ 0 sin(2πνt) t > 0 B xx t (t) = 1 + (B xy t ) 2 (t) . ( 12
)
The strain measure Q can thus be obtained analytically and the only evolution equation that remains to be solved numerically in LAOS is the evolution equation for f 2 t , Eq. ( 5) for linear polymers. To solve the evolution equation for f 2 t , we follow the idea of Hulsen et al. [START_REF] Hulsen | A new approach to the deformation fields method for solving complex flows using integral constitutive equations[END_REF], who used the age τ = tt as independent variable instead of t to solve the evolution equation for the Finger tensor B(t, τ) = B t (t). The introduction of τ as an independent variable modifies the time derivative in the evolution equation, but leaves the right-hand side unaltered. The proper equation for f 2 (t, τ) = f 2 t (t) involves a derivative with respect to t and τ resulting from the material derivative while leaving the right-hand side of an evolution equation unaltered. For the molecular stress function describing linear polymers, Eq. ( 5), we obtain

Df 2 Dt + ∂f 2 ∂τ = f 2 κ : Q - 1 f 2 -1 CR (13)
subject to the boundary condition f 2 (t, 0) = 1 and the initial condition f 2 (0, τ) = 1, since the fluid is assumed to be at rest for t < 0.

In [START_REF] Hulsen | A new approach to the deformation fields method for solving complex flows using integral constitutive equations[END_REF], the discontinuous Galerkin (DG) method is used for the discretization of the differential equation for the Finger tensor in the τ direction. For transient flow problems involving large and small time scales, the DG discretization in τ is not very appealing. Large relaxation times imply a long time before a steady periodic state has been reached. For the simulations in Section 7 we used 200 s. The small relaxation times on the other hand require a fine discretization for small τ. The DG method requires a time step of t < τ/6 for reasons of stability. For our simulations in Section 7, we use 1520 τ intervals which have a minimum length of 9.69 × 10 -6 s. For the discontinuous Galerkin method to remain stable, the maximum time step is approximately 1.6 × 10 -6 . To avoid very long computation times, we proceed differently.

Discretization of the age τ is performed in a similar way as in [START_REF] Hulsen | A new approach to the deformation fields method for solving complex flows using integral constitutive equations[END_REF]. We replace the semi-infinite age interval τ ∈ [0, ∞) by a finite interval [0, τ c ] with τ c the cut-off age. This value has to be large compared to the largest relaxation time λ max of the fluid. For the LAOS simulations in Section 7, we have used τ c = 20λ max , which is a rather conservative value. The interval [0, τ c ] is divided into N subintervals [τ j , τ j+1 ] for each j = 0, . . . , N -1. These subintervals are of increasing size to take advantage of the fast decaying memory function for large values of τ/λ i . For a one mode upper-convected Maxwell integral model, Hulsen et al. [START_REF] Hulsen | A new approach to the deformation fields method for solving complex flows using integral constitutive equations[END_REF] used a stretched mesh and determined an optimal stretching factor. We found that this was not the optimal τ discretization for our multimode MSF model. To determine the mesh for the age discretization, we define a cut-off time τ c,i = 20λ i for every relaxation time. For the interval [0, τ c,1 ], we compute the τ discretization exactly as in [START_REF] Hulsen | A new approach to the deformation fields method for solving complex flows using integral constitutive equations[END_REF]. For other intervals [τ c,i , τ c,i+1 ], we use less subintervals since the region up to τ c,i has already been discretized based on the smaller relaxation times.

In view of the evaluation of the stress integral Eq. ( 1), solutions to Eq. ( 13) are computed in the two-point Gauss points on each subinterval. The 2N Gauss points τ G k for each k = 1, . . . , 2N define the mesh on which we compute the molecular stress function. To include the end points of the τ domain we define τ G 0 = 0 and τ G 2N+1 = τ c . At each Gauss point τ G k for k = 1, . . . , 2N + 1, we need to solve Eq. ( 13). Note that for LAOS, the Finger tensor and thus the deformation tensor is known as a function of time and τ via Eq. [START_REF] Bird | Dynamics of Polymeric Liquids[END_REF].

Integration of Eq. ( 13) from time t i to t i+1 and over a τ interval between two Gauss points [τ G k , τ G k+1 ] is performed using a trapezoidal rule for the time and τ direction. This results, for each k = 0, . . . , 2N, in a difference equation for the molecular stress function at the new time level i + 1,

f 2 i+1,k+1 = f 2 i,k + t -τ t + τ (f 2 i+1,k -f 2 i,k+1 ) + 1 2 t τ t + τ × (r i+1,k+1 + r i,k+1 + r i+1,k + r i,k ), ( 14
)
where t = t i+1t i , τ = τ G k+1τ G k , and r denotes the righthand side of Eq. ( 13). All quantities f 2 i, * and r i, * are at the previous time level and are known at the start of a new time step. The quantities f 2 i+1,0 and r i+1,0 are known from the boundary condition at τ = 0, f 2 i+1,0 = 1. For the first τ interval [0, τ 1 ], only the quantities f 2 i+1,1 and r i+1,1 , which is a non-linear function of f 2 , are unknown. To handle the non-linearity of the right-hand side, we use a predictor-corrector scheme, for which we use f 2 i,j+1 as a predictor. Once f 2 i+1,1 is known, we can apply the same procedure to obtain f 2 i+1,2 and so on, till we have computed

f 2 i+1,2N+1
at the cut-off age of the largest relaxation time.

The molecular stress function for a linear polymer melt has to remain in the range (0, f 2 max ). When we use Eq. ( 14) in LAOS simulations, we found that f 2 can easily exceed the lower and upper bound under strongly non-linear flow conditions, particularly at large ages τ where the difference between two consecutive ages is large. This makes the method unstable. For example, once f 2 becomes negative for some value of τ it remains negative and rapidly increases in magnitude. The instability disappeared when we applied a transformation that ensures that the molecular stress function always remains between its lower and upper bound. For all our computations, we used the transformation

h = ln f 2 f 2 max -f 2 which maps (0, f 2 max ) to (-∞, ∞).
Instead of solving the differential equation for f 2 , we solve the corresponding differential equation for h. Since -∞ < h < ∞, we cannot violate any constraints. After obtaining the value of h at a new time level, the value of f 2 is recovered using the inverse mapping

f 2 = e h f 2 max 1 + e h
which indeed ensures that 0 < f 2 < f 2 max . In terms of time t and age τ, the integral for the stress, Eq. ( 1), becomes

T (t) = 5 ∞ 0 m(τ)f 2 (t, τ)Q[B(t, τ)] dτ. ( 15
)
Once the molecular stress function is computed at the new time level, the polymer stress T (t i+1 ) can be computed by integrating over τ. On the interval [0, τ c ], the integral is approximated by a finite sum and on [τ c , ∞) we assume f 2 (t, τ) = f 2 (t, τ c ) and B(t, τ) = B(t, τ c ), so that the integral can be integrated exactly.

The resulting stress at the new time level is of the form

T (t i+1 ) = 5 2N k=1 w k m(τ G k )f 2 (t i+1 , τ G k )Q[B(t i+1 , τ G k )] + 5M(τ c )f 2 (t i+1 , τ c )Q[B(t i+1 , τ c )], (16)
where w k , k = 1, . . . , 2N are the weights corresponding to a two-point Gauss quadrature rule and

M(τ c ) = i G i exp(-τ c /λ i ).
The numerical simulation of the Doi-Edwards model in LAOS is more straightforward since no evolution equation needs to be solved for f 2 . For every time step, only the stress needs to be computed using Eq. ( 16) with f 2 = 1 and the analytic solution for the Finger tensor Eq. [START_REF] Bird | Dynamics of Polymeric Liquids[END_REF]. For this we use exactly the same age discretization as for the MSF model.

Melt properties

The linear spectrum of the Finathene melt that we use in all our viscoelastic simulations is the six-mode spectrum identified in [START_REF] Debbaut | Large amplitude oscillatory shear and Fourier transform rheology for a high-density polyethylene: experiments and numerical simulation[END_REF] where it was used for LAOS simulations with the Giesekus model.

In addition to the linear spectrum, the MSF model only contains two non-linear parameters for the molecular stress function, a 1 and a 2 . The parameter a 1 describes the melt rheology in nonrotational flows and its value can be obtained from fitting elongational data. Next, the value of a 2 , which is relevant for rotational flows, can be determined from the shear viscosity and first normal stress difference. For the Finathene melt, however, only shear data are available. This is not a major limitation since the impact of the parameter a 1 in large amplitude oscillatory shear flows is very small as we show in Section 7. In this section, we use the value a 1 = 0.02 which was identified in [START_REF] Wagner | The molecular stress function model for polydisperse polymer melts with dissipative convective constraint release[END_REF] for another high-density polyethylene melt.

The second parameter a 2 is identified using steady shear viscosity data. These shear data were obtained from small amplitude oscillatory shear measurements and the Cox-Merz rule which is valid for the Finathene fluid [START_REF] Debbaut | Large amplitude oscillatory shear and Fourier transform rheology for a high-density polyethylene: experiments and numerical simulation[END_REF]. For another highdensity polyethylene melt, the value of a 2 = 2.3 was identified in [START_REF] Wagner | The molecular stress function model for polydisperse polymer melts with dissipative convective constraint release[END_REF]. Using this parameter value, steady shear viscosity predictions also agree well with experimental data of the Finathene fluid as can be observed from Fig. 1. At shear rates of order unity, the experimental data are slightly overpredicted while for γ > 10 s -1 a slight underprediction is apparent. Increasing the value of a 2 leads to better agreement with experimental data at shear rates around γ ≈ 1 s -1 . The underprediction of the data for γ > 10 s -1 , however, persists and even slightly increases when a 2 is increased. We conclude from Fig. 1 that the steady shear viscosity predictions are not very sensitive to changes in a 2 and that a rather wide range of parameter values a 2 fits the data equally well. For this reason, we will consider both a 2 = 2.3, 4, and 8 in the LAOS simulations in Section 7.

Results of the MSF model in large amplitude oscillatory shear will be compared with the basic Doi-Edwards model, the Giesekus model, and the inelastic Carreau-Yasuda model [START_REF] Bird | Dynamics of Polymeric Liquids[END_REF]. For both viscoelastic models, we use the same linear spectrum as for the MSF model. For the Giesekus model, we use in addition the non-linear parameter values identified in [START_REF] Debbaut | Large amplitude oscillatory shear and Fourier transform rheology for a high-density polyethylene: experiments and numerical simulation[END_REF]. The viscosity of the Carreau-Yasuda model is given by

η = η 0 (1 + [λI 2 ] a ) (n-1)/a , ( 17
)
where I 2 is the second invariant of the rate-of-strain tensor d. The four adjustable parameters in the model are the zero-shear viscosity η 0 , a time constant λ, the power-law index n, and a numerical parameter a. The parameters are obtained from fitting the shear viscosity. We found a good fit using η 0 = 8.66 × 10 4 Pa s, λ = 18 s, n = 0.5, and a = 0.85 The steady shear viscosity predictions of the Doi-Edwards and Carreau-Yasuda models are displayed in Fig. 2. Surpris-ingly, the predictions of the Doi-Edwards model, which has no adjustable parameters, agree well with the data up to moderate shear rates. At higher shear rates, γ > 10 s -1 , the Doi-Edwards model underpredicts the experimental data. Up to shear rates that can be reached in the LAOS experiments ( γ < 32 s -1 ), however, differences are relatively small. In fact, the predictions of the DE model are only slightly lower than those of the MSF model with a 2 = 8 as can be observed by comparing Figs. 1 and2. The Carreau-Yasuda fit overpredicts the viscosity at large shear rates. However, for the range of shear rates in the LAOS experiments, γ ≤ 32 s -1 , the purely viscous model shows good agreement. For further reference, we have also included in Fig. 2 the steady shear viscosity predictions of the Giesekus model. The Giesekus model shows excellent agreement up to the largest experimental shear rate.

To analyse shear predictions for a transient flow like LAOS, it is also important to know the response in transient rheometrical flows. Fig. 3 displays the viscosity in start up of shear and start up of uniaxial elongation at various deformation rates. The differences between the results of the MSF models with different values of the non-linear parameter a 2 are small. Both a 2 = 2.3, 4, and 8 show small time overshoots at moderate shear rates. The magnitudes of these overshoots and the steady viscosity slightly decrease with increasing value of a 2 . For the largest value of a 2 = 8, the predictions lie only just above those of the Doi-Edwards model. For further reference, we also include in Fig. 3a the model predictions of the Giesekus model. There are two differences with the molecular models. First, the larger time overshoots. This becomes apparent at relatively low shear rates of γ = 1 and becomes more pronounced when the shear rate is increased. Second, the Giesekus model has a higher steady shear viscosity at high shear rates. This is in better agreement with the data in Fig. 1. For the LAOS experiments, however, such high rates cannot be achieved experimentally.

Fig. 3b shows that the differences between the models is much more pronounced in start up of uniaxial extension, at significantly large values of the extension rate ˙ . The MSF model shows considerable strain hardening which increases when a 1 is increased. This behaviour is absent for the Doi-Edwards model.

Validation of the numerical results

For the validation of the numerical technique and various numerical parameters, we use the MSF model for linear polymers with parameter values a 1 = 0.02 and a 2 = 2.3. Other parameter settings, not shown in this section, gave identical results.

We first validate our numerical technique described in Section 4, by comparing with the deformation field method using the discontinuous Galerkin method to discretise the τ direction in the equation governing the molecular stress function. This is exactly the same approach as discussed in [START_REF] Hulsen | A new approach to the deformation fields method for solving complex flows using integral constitutive equations[END_REF] for the Finger tensor. For both techniques we take the same τ discretization using 1520 τ subintervals. Fig. 4a shows the results of both techniques in LAOS with ν = 1 Hz and γ 0 = 5 which has the highest achievable shear rate in the experiments. At the scale of the plot the methods are indistinguishable. The main difference lies in the time step employed. For the new method we have used a time step of t = 2.5 × 10 -4 s while for the DG method a time step of t = τ/6 is necessary for reasons of stability [START_REF] Hulsen | A new approach to the deformation fields method for solving complex flows using integral constitutive equations[END_REF]. For a τ discretization using 1520 subintervals, the minimum length of a subinterval equals τ min = 9.6 × 10 -6 s. This requires a time step smaller than t = 1.6 × 10 -6 s. In combination with the long time required to reach a periodic steady state, such small time steps lead to long simulation times for the DG method. The new technique, however, remains stable for much larger time steps leading to substantially smaller computation times.

Second, it needs to be checked whether the discretization in τ using N = 1520 subintervals is sufficient to capture the steady periodic stress response. In order to verify this, we compare with a τ discretization that has twice as many subintervals, N = 3040, which implies a smallest τ subinterval which is halved, τ min = 4.8 × 10 -6 s. For both computations, we used a time step of t = 2.5 × 10 -4 s. The steady periodic shear stress for both τ discretizations is displayed in Fig. 4b. On the scale of the figure there are no differences between the two τ discretizations. Henceforth, we use the smallest number of subintervals, N = 1520, for all simulations with the MSF model.

A time step of t = 2.5 × 10 -4 s is sufficiently small to capture the steady periodic regime as can be observed from Fig. 5a. When the time step is halved to t = 1.25 × 10 -4 s, both time steps produce identical results on the scale of the figure. For both simulations, we used N = 1520 subintervals for the τ dis- cretization. Henceforth, we use t = 2.5 × 10 -4 s for all LAOS simulations.

In Section 7, LAOS simulations are performed for a time of t s = 200 s. Fig. 5b compares shear stresses using t s = 200 and 1000 s for the highest experimentally achievable amplitudes at a low (ν = 0.1 Hz) and a high frequency (ν = 1 Hz). For both calculations, we used t = 2.5 × 10 -4 s and N = 1520. We conclude from Fig. 5b that t s = 200 s is sufficiently large to reach the steady periodic regime. This value of t s is a rather conservative choice. The stress maxima only differ 0.1% from the steady periodic value after 3 and 15 cycles for ν = 0.1 and 1 Hz, respectively.

Large amplitude oscillatory shear results

The MSF theory uses the strain measure based on the independent alignment approximation. In view of the large number of time steps that have to be performed, it is advantageous to use the computationally much more efficient Currie approximation instead. This is a good approximation in steady and start up of shear for the Doi-Edwards model [START_REF] Currie | Constitutive equations for polymer melts predicted by the Doi-Edwards and Curtiss-Bird kinetic theory models[END_REF]. It is, however, a priori not clear how well both strain measures correspond in large amplitude oscillatory shear flows of an MSF fluid. For this we compared the steady periodic shear stress at a low (ν = 0.1 Hz) and high frequency (ν = 1 Hz) for the two highest amplitudes γ 0 that are experimentally achievable for each frequency. As can be seen from Fig. 6, the flow curves of the Currie and independent alignment approximation agree very well. Only near the minima and maxima we observe that the Currie approximation slightly overpredicts the independent alignment approximation. This is, however, negligible compared to the differences between the various model predictions and the experimental data as we discuss shortly. The approximation is equally good for the other two frequencies, ν = 0.3 and 3 Hz, which are not shown in Fig. 6. We conclude from Fig. 6 that in large amplitude oscillatory shear simulations of an MSF fluid, the Currie approximation is a very good approximation to the independent alignment approximation of the Doi-Edwards strain mea- sure. Henceforth, we use the Currie approximation in all LAOS simulations.

To accurately determine the parameter a 1 of the MSF model for a polymeric fluid, extensional data are required. For the Finathene melt, however, only shear data are available. Therefore, we first compare the periodic steady shear response for two realistic values of a 1 . For the first value, we take a 1 = 0.02 (f 2 max = 51) as identified in [START_REF] Wagner | The molecular stress function model for polydisperse polymer melts with dissipative convective constraint release[END_REF] for another high-density polyethylene melt. The linear polymers considered in [START_REF] Wagner | Determination of elongational viscosity of polymer melts by RME and Rheotens experiments[END_REF] have a lower f 2 max . For the second value, we take f 2 max = 11 which corresponds to a 1 = 0.1. Fig. 7 shows the impact of a 1 on the periodic steady shear response for ν = 1 Hz at the two highest experimentally achievable amplitudes γ 0 = 2.5 and 5. At the scale of the figure, both parameter values of a 1 result in identical model predictions. We conclude that the only relevant MSF model parameter for large amplitude oscillatory shear flow of linear polymers is a 2 and henceforth we use a 1 = 0.02 for all simulations.

The LAOS experiments have been performed at various frequencies and amplitudes. At low values of ν and γ 0 , the flow does not deviate much from low amplitude oscillatory shear. It is therefore not surprising that all viscoelastic models predict similar results that do not differ much from the experimental data. As an example, we display for γ 0 = 1 the steady periodic response for ν = 0.1 and 1 Hz in Fig. 8. All viscoelastic models considered predict the correct phase shift of the periodic shear stress. The MSF model, for all parameter values of a 2 considered, also correctly predict the amplitude. The Doi-Edwards model, however, clearly underpredicts the amplitude in the mildly non-linear regime. To investigate the impact of viscoelasticity, we also display results of an inelastic model possessing only shear-thinning behaviour. This model is only able to predict the correct amplitude and fails to predict the phase shift. At ν = 1 Hz and γ 0 = 1, which corresponds to a maximum shear rate of γ ≈ 6.3 s -1 , the inelastic model also significantly overpredicts the amplitude of the periodic shear stress signal.

At larger values of γ 0 , differences between the various parameter settings in the MSF model and the Doi-Edwards model become more apparent. In Fig. 9, we display for each frequency the steady periodic shear stress for the two largest experimentally achievable values of γ 0 . At the lowest strains, the Doi-Edwards model consistently underpredicts the experimental data, while the MSF predictions still show good agreement for all values of a 2 considered, particularly a 2 = 4. For the larger γ 0 (and thus larger shear rates) the Doi-Edwards predictions are again fairly close to the experimental data. Surprisingly, the results are in better agreement with the experiments than most of the predictions of the MSF model which has two extra parameters to capture the non-linear behavior. The best MSF parameter for the largest strains seems to be a 2 = 8. The parameter values a 2 = 4 and particularly a 2 = 2.3 overpredict the shear response. This is most apparent at ν = 1 Hz and γ 0 = 5 which has the highest maximum shear rate. We note in passing that the Giesekus model also overpredicts the experimental signal at more non-linear flow conditions [START_REF] Debbaut | Large amplitude oscillatory shear and Fourier transform rheology for a high-density polyethylene: experiments and numerical simulation[END_REF]. At ν = 1 Hz and γ 0 = 5, the predicted shear stress amplitude is, for example, 117 kPa for the Giesekus model, while the largest amplitude that occurs for the MSF model is 113.5 kPa for a 2 = 2.3. At this amplitude and frequency, the amplitude of the experimental signal 106 kPa is best predicted by a 2 = 8 which gives an amplitude of 106.8 kPa. The Doi-Edwards model slightly underpredicts the amplitude, 101.8 kPa. Similar trends are visible at the lower frequencies ν = 0.1 and 0.3 Hz. At ν = 3 Hz, no data are available at high enough γ 0 to confirm the counter intuitive result that the Doi-Edwards model gives better predictions at large strains (large shear rates). It is remarkable that at corresponding values of the shear rate, the steady shear viscosity data in Fig. 1 are better predicted by the MSF and Giesekus model than by the Doi-Edwards model which shows a stronger underprediction of the shear viscosity at high shear rates. We also observe that the LAOS simulations are more sensitive to changes in the parameter a 2 than steady shear flow at corresponding shear rates. This indicates that large oscillatory shear might be a better flow to determine the non-linear parameters that are important for shearing flows.

The Lissajous plot for a frequency ν = 1 Hz and all experimentally available strain amplitudes are depicted in Fig. 10. All viscoelastic models predict the deviation from an ellipsoidal shape that is characteristic for the non-linear response at large strains. Quantitative differences are observed at the largest strain γ 0 = 5 for the Giesekus model and the MSF model, particularly with a 2 = 2.3 and to a lesser extent with a 2 = 4. At this strain, predictions of the MSF model with a 2 = 8 and Doi-Edwards model are in good agreement with the experimental data. At lower strains, however, some discrepancies are noticeable for the Doi-Edwards model. For all viscoelastic models, deviations are in accordance with the over and underprediction of the amplitude of the periodic shear stress signal as a function of time observed in Fig. 9.

For the first normal stress difference N 1 , no experimental data are available. For completeness, we have included the model predictions at experimentally achievable strains for ν = 1 Hz in Fig. 11. The periodic N 1 response has a period that is half of the corresponding shear stress signal and has a non-zero average. All models predict the same phase shift for all strains and nearly the same minimum value of N 1 . The difference between the models lie in the maxima of the amplitudes. Similar to the periodic shear stress, the Doi-Edwards model predicts the lowest and the MSF model with smallest value of a 2 the Fig. 9. Comparison of MSF model using values of a 2 as in the legends and Doi-Edwards model with periodic shear stress data for Finathene melt [START_REF] Debbaut | Large amplitude oscillatory shear and Fourier transform rheology for a high-density polyethylene: experiments and numerical simulation[END_REF]. Left column: one period; right column: zoom around minimum. largest amplitude. Different magnitudes of the amplitudes become already apparent at relatively low strains of γ = 1 where the shear stress predictions are still very similar, particularly for the various values of a 2 of the MSF model (Fig. 8). The periodic normal stress, however, can clearly be distinguished and differ-ences between the models continue to grow when the strain is increased. It would be interesting to compare these results with experimental data to establish whether the Doi-Edwards model also better predicts the first normal stress difference at high strains. Differences between the model predictions are better quantified in the frequency domain than in the time domain. Under more non-linear flow conditions, the higher odd harmonics become more significant. The even harmonics should vanish except for some numerical noise. For the following computations, the order of magnitude of the even harmonics was at least 6 orders lower than that of the largest odd harmonics. We focus on the two highest experimentally achievable amplitudes at a low frequency of ν = 0.1 Hz and a high frequency of ν = 1 Hz.

We first consider the two highest experimentally achievable strains at the lowest frequency ν = 0.1 Hz. Table 1 shows the Fourier coefficients A n and B n for the MSF model using various values of a 2 . At this low frequency the response in phase with the shear rate is dominating and becomes more dominant when the strain is increased. This is correctly predicted by the MSF model for all parameter values considered. Furthermore, the MSF model predicts the signs of all Fourier coefficients correctly. There are only small quantitative differences. At the low strain γ 0 = 5, the low harmonics agree best for a 2 = 4, particularly the dominating amplitude B 1 which corresponds to the lost work per cycle [START_REF] Giacomin | Using large-amplitude oscillatory shear[END_REF]. At the high strain a 2 = 8 gives the best agreement for the dominating amplitude B 1 . This explains the better agreement of a 2 = 4 at γ 0 = 5 and of a 2 = 8 at γ 0 = 10 in Fig. 9.

The amplitudes of the first four non-zero harmonics for the Doi-Edwards, Carreau-Yasuda, and Giesekus model are tabulated in Table 2. The amplitudes of the Doi-Edwards model closely follow the MSF results with a 2 = 8, with the exception of B 1 which is somewhat lower. At γ 0 = 5, the underprediction of the amplitude of the experimental shear stress signal in Fig. 9 is caused by an underprediction of the first harmonics A 1 and B 1 . The higher harmonics of the Doi-Edwards model are, however, in good agreement with the experimental data, while the Giesekus model predicts the signs of some higher harmonics incorrectly. It is remarkable that at the larger strain of γ 0 = 10, the differences with the experimental A 1 and B 1 are considerably smaller, resulting in a slightly better agreement with the experimental data than the best MSF model at this γ 0 . Furthermore, the higher harmonics of the Doi-Edwards model are in excellent agreement with the experimental data while the Giesekus model shows again much larger deviations and predicts some signs incorrectly.

At the frequency ν = 0.1 Hz and at the strain values in Table 2, the Carreau-Yasuda model is still able to predict correctly the coefficient B 1 , i.e. the lost work per cycle. The amplitudes B 3 , B 5 , and B 7 , however, are off by a factor 2-10 and these might have the wrong sign. Furthermore, the inelastic model only predicts non-zero values for the coefficients B n and is therefore incapable to predict the phase shift in the shear stress response.

At the frequency ν = 1 Hz, the response in phase with the strain is much larger than for the low frequency of ν = 0.1 Hz. This is correctly predicted by the MSF model. The deviations between the experiments and the MSF results are caused by an overprediction of B 1 , as can be observed from Table 3. For the parameter value a 2 = 2.3, the value of B 1 is overpredicted by at least 5% and deviations become larger when the strain increases. For larger values of a 2 , the deviations from the experimentally obtained B 1 become smaller. Particularly, a 2 = 8 agrees well with the experimental data. For some other amplitudes like A 3 and B 3 , however, we then find larger differ- ences with the experiments. These amplitudes are off by at least 25%.

Table 4 shows the corresponding results for the Doi-Edwards, Giesekus, and Carreau-Yasuda model. In this more non-linear regime, the purely viscous model is not capable to predict any of the amplitudes correctly. Also, B 1 is now overpredicted considerably. The Doi-Edwards model underpredicts both A 1 and B 1 at the lower strain γ 0 = 2.5. For this strain, the amplitudes of the low harmonics are better predicted by the Giesekus and MSF model. The higher harmonics, which arise due to non-linear effects, are again much better predicted by the Doi-Edwards than by the Giesekus model. For the more non-linear regime at the higher strain, the Doi-Edwards model does not only show better agreement with experimental data for the higher harmonics. Also, the dominating term B 1 is much better predicted than the Giesekus model does. This confirms the counter intuitive result found for the lower frequency of ν = 0.1 Hz that the Doi-Edwards model, which has no nonlinear parameters and underpredicts the steady shear viscosities at high rates, is more accurate in the more non-linear regimes. Fig. 12 shows the model predictions of the MSF and Doi-Edwards model outside the experimental window. To compare the model predictions with a macroscopic stress model, results for the Giesekus model have been included as well. All models predict the same trends at large strains, for example, the same sign is predicted for all Fourier coefficients B n and A n at large strains. The dominating amplitude is B 1 , while the magnitude of A 1 decreases at large strains and reaches a plateau. For the higher harmonics, all amplitudes B n in phase with the shear seem to reach a higher magnitude at large strains than the corresponding amplitudes A n that are in phase with the strain. The Fourier coefficient B n seems to continue to increase in magnitude for a much longer range of strains than the corresponding A n , which reaches a plateau value. For the two lowest harmonics, the plateau is reached at medium strains of 10. For the higher harmonics much higher strains are necessary. Quantitatively, however, there are differences between the models. Only for A 1 , all models predict the same values. For all other Fourier coefficients, the Doi-Edwards model predicts the lowest magnitudes. At large strains, the amplitudes are considerably lower than those of the MSF model with a 2 = 8. For the MSF model, the magnitudes of the odd harmonics increase when the parameter a 2 is decreased. The magnitudes predicted by the MSF model with a 2 = 2.3 are, however, still considerably lower than those predicted by the Giesekus model. Particularly for the dominant amplitude B 1 the molecular models predict significantly lower magnitudes at large strains. This is consistent with the smaller time overshoots for these models in the transient shear viscosity in Fig. 3a. Whether the trends at large strains are correct and whether the Doi-Edwards predictions are still more accurate than those of the non-linear models remains to be established.

Concluding remarks

We have evaluated the rheometrical response of the integral molecular stress function model in large amplitude oscillatory shear. For reasons of efficiency, we modified the deformation field method to avoid the restrictive time step constraint resulting from the discontinuous Galerkin method used in [START_REF] Hulsen | A new approach to the deformation fields method for solving complex flows using integral constitutive equations[END_REF]. The new technique is both accurate and efficient for the LAOS computations we performed.

In LAOS, the MSF model only has one relevant material parameter to describe the non-linear material response. By fitting this parameter to steady shear data of a linear high-density polyethylene melt, the MSF model is able to predict qualitatively and quantitatively the response in large amplitude oscillatory shear. Up to medium strains, results agree very well with experimental data and are comparable with those of a six-mode Giesekus fluid, having six parameters to describe the non-linear material response. At medium strains, the non-linear response becomes important. As expected, the introduction of the molecular stress function results in better predictions than the basic Doi-Edwards theory which has no material parameters to describe the non-linear response. Although correctly predicting the phase shift, the DE model underpredicts the amplitude of the experimental signal. This is caused by an underprediction of the amplitude of the most dominant odd harmonics.

At the highest experimentally achievable strains, both the MSF and Giesekus model overpredict the amplitude of the periodic shear stress, although this is more significant for the Giesekus model. At these strains, however, it is the Doi-Edwards model that shows excellent agreement with the experimental data. In the time domain it only slightly underpredicts the amplitude of the periodic shear stress. In the frequency domain the correspondence of both the dominant harmonic and the higher odd harmonics is striking. We recall that the DE model does not have any non-linear parameters and underpredicts the steady shear stress at high shear rates. This is a remarkable and unexpected result indeed. Whether this trend continues at higher strains remains to be established. Since the Doi-Edwards model is well known to underpredict experimental data in start up of shear and steady shear flows, this also raises the question whether the experimental data are somewhat inaccurate at high strains. Further experimental results are clearly needed to confirm our counter intuitive findings at high strains.

Fig. 1 .

 1 Fig. 1. Steady shear viscosity of the Finathene melt: experimental data of [8] and model predictions of the MSF model using various values of a 2 (a 1 = 0.02).

Fig. 2 .

 2 Fig. 2. Steady shear viscosity of the Finathene melt: experimental data of [8] and model predictions of the Doi-Edwards, Carreau-Yasuda, and Giesekus model.

Fig. 3 .

 3 Fig. 3. Transient viscosity of MSF model for various parameter settings. The Doi-Edwards and Giesekus model are included for comparison. Deformation rates (s -1 ) are indicated in the figures. Start up of: (a) shear for a 1 = 0.02 and values of a 2 indicated in the legend and (b) uniaxial elongation for values of a 1 indicated in the legend (a 2 irrelevant).

Fig. 4 .

 4 Fig. 4. Validation of the numerical technique. Comparison of: (a) new method with the DG method of [9] and (b) various τ discretizations.

Fig. 5 .

 5 Fig. 5. Validation of sufficiently small time step (a) and sufficiently large t s (b).

Fig. 7 .

 7 Fig. 7. Steady periodic shear stress at ν = 1 Hz for the MSF model using a 2 = 2.3 and values of a 1 indicated in the legend.

Fig. 6 .

 6 Fig. 6. Comparison of the steady periodic shear stress prediction of the MSF model using the Currie and independent alignment approximation at amplitudes γ 0 indicated in the figures: (a) ν = 0.1 Hz and (b) ν = 1 Hz.

Fig. 8 .

 8 Fig. 8. Periodic shear stress at moderate shear rates predicted by the viscous and viscoelastic models together with the Finathene data. Frequencies and strains indicated in the figures.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Lissajous plots for ν = 1 Hz for Finathene data of [8] and model predictions. Strain amplitudes γ 0 are 0.5, 1, 2.5, and 5 (from innermost to outermost loops).

Fig. 12 .

 12 Fig. 12. (a-d) Amplitudes A n and B n of the odd harmonics (kPa) outside the experimental window for MSF model using values of a 2 as in the legends. The Doi-Edwards and Giesekus model are included for comparison.

Table 1 Amplitudes

 1 A n and B n of the odd harmonics (kPa) for ν = 0.1 Hz (Comparison of experimental data and MSF results using various values of a 2 .)

	γ 0	n	A n (kPa)

Table 2 Amplitudes

 2 A n and B n of the odd harmonics (kPa) for ν = 0.1 Hz (Comparison of experimental data, Doi-Edwards, inelastic Carreau-Yasuda model, and Giesekus model. Data for the Giesekus model are from[START_REF] Debbaut | Large amplitude oscillatory shear and Fourier transform rheology for a high-density polyethylene: experiments and numerical simulation[END_REF].)

	γ 0	n	A n (kPa)				B n (kPa)
			Experimental	DE	CY	G	Experimental	DE	CY	G
	5	1	1 0 .9	8 .06	-	8.19	38.8	3 6 .6	3 9 .3	3 9 .2
		3	-2.84	-3.02	-	-2.90	-0.34	-0.799	-5.29	-0.163
		5	0 .259	0.297	-	-0.124	-0.534	-0.499	2.30	-0.704
		7	0 .027	0.103	-	0.151	0.116	0.106	-1.32	-0.062
	10	1	8.44	6.70	-	6.51	55.3	5 4 .1	5 6 .2	5 9 .4
		3	-4.15	-4.34	-	-4.34	-3.43	-3.71	-7.76	-3.27
		5	1 .38	1.42	-	1.20	-0.219	-0.284	3.43	-1.11
		7	-0.295	-0.152	-	0.388	0.37	0.411	-2.01	0.419

Table 3 Amplitudes

 3 A n and B n of the odd harmonics (kPa) for ν = 1 Hz (Comparison of experimental data and MSF results using various values of a 2 .)

	γ 0	n	A n (kPa)

Table 4 Amplitudes

 4 A n and B n of the odd harmonics (kPa) for ν = 1 Hz (Comparison of experimental data, Doi-Edwards, inelastic Carreau-Yasuda model, and Giesekus model. Data for the Giesekus model are from[START_REF] Debbaut | Large amplitude oscillatory shear and Fourier transform rheology for a high-density polyethylene: experiments and numerical simulation[END_REF].)

	γ 0	n	A n (kPa)				B n (kPa)
			Experimental	DE	CY	G	Experimental	DE	CY	G
	2.5	1	3 8 .0	3 4 .1	-	3 8 .2	7 4 .6	7 1 .6	8 9 .5	7 7 .8
		3	-2.92	-4.52	-	-3.45	2.40	1.14	-12.6	2 .31
		5	-0.04	-0.274	-	-0.474	-0.512	-0.538	5.65	-0.227
		7	0 .098	0.063	-	-0.008	0.296	-0.025	-3.35	-0.059
	5	1	3 3 .1	3 0 .7	-	3 5 .1	108	106	127	119
		3	-9.66	-11.2	-	-10.2	-1.23	-2.46	-18.0	-0.211
		5	0 .817	0.996	-	-0.492	-1.62	-1.85	8.11	-2.42
		7	0 .075	0.391	-	0.492	0.412	0.385	4.83	-0.242

Acknowledgements

We wish to thank Henri Burhin and Christian Bailly for making the Finathene data available to us. We also thank Manfred Wagner for fruitful discussions. This work is supported by the programme Action de Recherche Concertée, Communauté Franc ¸aise de Belgique.