Optimized Multilevel Time Discretization

Suppose u(t) is the solution of Eq. (3.1). In a multilevel time discretization method, the time axis is divided into a uniform grid with time step At. It will be assumed that the values of u and du/dt are known at time level n, n – 1, n – 2, n – 3. (Note: In CAA, du/dt is given by the governing equations of motion.) To advance to the next time level, the following four-level finite difference approximation is used:

(3.9)

The last term on the right side of Eq. (3.9) may be regarded as a weighted average of the time derivatives at the last 4 mesh points. There are four constants, namely, b0, b1, b2, and b3, that are to be selected. To ensure that the scheme is consistent, it is suggested to choose three of the four coefficients, say bj (j = 1, 2, 3), so that when the terms in Eq. (3.9) are expanded in a Taylor series in At they match to the order (At)2. This leaves one free parameter, b0. The relation between the other coefficients and b0 are as follows:

ro – i s’ лп

b1 =-3b0 + 12 ’ b2 = 3b0 3, b3 =-b0 + 12. (3.10)

The remaining coefficient b0 will now be determined by requiring the Laplace trans­form of the finite difference scheme (3.9) to be a good approximation of that of the partial derivative.

The Laplace transform and its inverse of a function f(t) are related by (see Appendix A) the following:

(Ю

f(m) = 2П j f (t)elmtdt, f (t) = j f(a)e~imtda. (3.11)

0 г

The inverse contour Г is a line in the upper half m plane parallel to the real m­axis above all poles and singularities. Before Laplace transform can be applied to Eq. (3.9), it is necessary to generalize the equation to one with a continuous variable. This can easily be done. The result is as follows:

3

^ du(t – jAt) ,„„„4

u(t + At) = u(t) + At bj—————– . (3.12)

d t

Eq. (3.12) reduces to Eq. (3.9) by setting t = nAt. It will be assumed, for the time being, that u(t) satisfies the following trivial initial conditions (the case with general initial conditions will be discussed later); i. e.,

In Appendix A, the following shifting theorem for Laplace transform is proven:

Transform of f (t + A) = e-imAf; ~ = Laplace transform (3.14)

On applying the shifting theorem to Eq. (3.12), it is easy to find

Thus,

.i(e~mAt -1) _ ( f tdU

– i——————– u = transform of — . (3.15)

3 ,. a V d t f

AtJ2 bjeijmAt

j=0

The Laplace transform of the time derivative of u is – i’«U. Thus, by comparing the two sides of Eq. (3.15), the quantity

i(e-iaAt – 1)

3

At ^2 b. eijmAt
j=0 j

is the effective angular frequency of the time marching scheme (3.9). The weighted error Ex incurred by using ш to approximate ш will be defined as

z

E1 = j [a[Re(«At – rnAt)]2 + (1 – a)[Im(«At – rnAt)]2}d(«At), (3.17)

-z

where Re() and Im() are the real and imaginary part of (). a is the weight and z is the frequency range. One expects ш to be a good approximation of ш. On substituting ш from Eqs. (3.16) and (3.10) into Eq. (3.17), Ex becomes a function of b0 alone. b0 will be chosen so that Ex is a minimum. Therefore, b0 is given by the root of

A range of possible values of a and z has been used in Eq. (3.18) to determine b0 and, hence, b1, b2, and b3. After considering the effect on the range of useful frequency and numerical damping rate, the values a = 0.36 and z = 0.5 were selected. For this value of a and z, the values of the coefficients bj (j = 0,1,2,3) are

b0 = 2.3025580888383, b1 = -2.4910075998482 b2 = 1.5743409331815, b3 = -0.3858914221716.