Computer methods

in applied
mechanics and
englaoerlng
ELSEYIER Comput. Methods Appl. Mech. Engrg. 142 (1997) 203-214

A time marching integral equation method for unsteady state

problems
Jichao Su”, Behrouz Tabarrok
Department of Mechanical Engineering, University of Victoria, Victoria, B.C., V8W 3P6, Canada

Received 15 October 1993; revised 17 May 1996

Abstract

A time marching integral equation method is developed for solving three classes of differential equations--elliptic, hyperbolic
and parabolic with the elliptic equation as the steady state limiting case of hyperbolic or parabolic equations. The method
described is semi-analytical with its numerical aspect in the time discretization of the differential equation and the space
discretization of the integral equation and its analytical aspect in the derivation of the integral equation via the Green’s function
method. The analyses of the time marching integral equation for three limiting cases indicate that both the steady integral
equation and the explicit algorithm are special cases of the time marching integral equation. Finite element ideas have been used
to discretize the derived integral equations. Some physical problems are analysed to illustrate the method and demonstrate its
reliability and efficiency.

1. Introduction

Most mathematical models of practical problems in engineering science take the form of differential
equations accompanied by appropriate boundary conditions and/or initial conditions. Generally,
analytical solutions of these problems are practically impossible to obtain due to complicated boundary
and/or initial conditions, irregular geometry or nonlinearity of the problems, and we have to rely on
numerical techniques to obtain approximate solutions.
The most widely used numerical techniques are the finite difference method (FDM) and the finite
element method (FEM). The FDM approximates the derivatives in the differential equations by some
type of truncated Taylor expansion in terms of the values at a number of discrete mesh points. This
results in a set of algebraic equations to which the boundary and/ or initial conditions are applied in
order to obtain an approximate solution of the problem. The FEM approximates the unknown
functions over each subdomain called an element in terms of polynomial interpolation functions. This
also results in a set of algebraic equations which are to be solved with appropriate boundary and/or
initial conditions when the weak formulation, derived from the differential equations using some
variational principles or a weighted residual technique, is discretized.
The boundary element method (BEM), also called the boundary integral equation method, is
another numerical technique which has been receiving increasing attention among scientists and
engineers, and has become an important alternative to the finite difference method and the finite
element method for the solution of a variety of engineering problems.

* Corresponding author, currently at Institute for Aerospace Research, National Research Council of Canada KlA 0R6.

00457825/97/$17.00 0 1997 Elsevier Science S.A. All rights reserved

PII SOO45-7825(96)01143-7
204 J. Su, B. Tabarrok I Comput. Methods Appl. Mech. Engrg. 142 (1997) 203-214

The modern theory of boundary integral equation can be traced back to 1903 with Fredholm [l], who
established the existence of solutions on the basis of his limiting discretization procedure. The
emergence of computer in 1950s stimulated the development of numerical methods including the
integral equation method. Much work has been done by Hess and Smith [2] for potential flow
problems, Jaswon and Ponter [3] for a torsion problem, Jaswon [4] and Rizzo [5] for a electrostatic
capacitance problem, Wu et al. [6], Tasaka and Onishi [7] and Dargush and Ganerjee [8] for
incompressible viscous flows, and many others. Great contributions have been made to bring this
method to the level of engineering applications by Brebbia [9,10], Zienkiewicz [ll] and others. Today,
the method under the title of Boundary Element Method is being used to deal with engineering
problems in many fields.
The BEM in most cases is applied for linear problems. When a nonlinear problem is encountered, a
field integral of unknown quantities is included in the integral equation whose discretization generally
yields a set of algebraic equations requiring more computational cost to solve. In this case the BEM is
more suitable to be called a field integral equation method.
While for an unsteady problem, most researchers in the field deal with the integral equation derived
by using a time-dependent Green’s function (i.e. a fundamental solution), Tosaka and Onishi [7] (to
authors’ knowledge) first introduced an approach of implementing time differencing before deriving an
integral equation. The similar technique was developed independently by Su [12,13] to obtain time
marching integral equations for the solutions of two- and three-dimensional unsteady transonic flows
around wings. The main procedures of the method are as follows. Firstly, time discretization of the
differential equation is made by replacing the time derivatives by finite difference. Secondly, the
time-discretized differential equation is transformed into an equivalent integral equation by applying
the Green’s function method. Finally, the integral equation is discretized in space and solved in each
time step using a numerical integration technique. In this paper further work along this direction has
been done to improve the method. Three limiting cases are discussed so as to better understand the
features of the method. To test the reliability and accuracy of the method, the calculations are carried
out for three examples including an initial and boundary value problem of the nonlinear Burger’s
equation.

2. Wave equation

An initial and boundary value problem of the one-dimensional wave equation is described by

I$=$+f(x,t) (O<x<l,t>O)
(1)
L
u(0, t) = u(1, t) = 0
U(X,0) = 0, r&, 0) = g(x)

Now let us expand the time derivative a2u/dt2 in the wave equation as

=- l (u ,,+’ - 2u, + u,_~) + O(At”) (2)

At2
where At denotes a time step. The space derivative is expanded in time as

+ (1 - 20)O(Aht) + O(At”) (3)

By substituting Eqs. (2) and (3) in the wave equation and omitting the higher-order terms, we obtain
 .I. Su, B. Tabarrok I Comput. Methods Appl. Mech. Engrg. 142 (1997) 203-214 205

d*u,+l
0 ~--p,+l
1
=--&U~_l-2u.)-(l-~)~-f(X,tn) (4)
It is easy to see that when 8 = 0, an explicit algorithm is obtained with a first-order accuracy in time;
when 8 = 1 a full implicit algorithm is obtained with a first-order accuracy in time; when 8 = l/2 an
implicit algorithm (Crank-Nicolson method) is obtained with a second-order accuracy in time.
Thus, by discretizing problem (1) in time, we obtain a boundary value problem of the ordinary
differential equation for u, +1(x)

= F(x)
(5)
%+1(O) = u,+,(l) =0
with

d*u,_,
(u,_, - 224,) - At’( 1 - 0) - du2 - At* f(x, t,)
I
uo(x) = 0
u,(x) = At - ut(x, 0) = At - g(x)

where k, = 1 /(A&).
To transform Eq. (4) into an integral form, we need to obtain the Green’s function which satisfies

(5-k;
) G=S(x-x,) (6)
with the boundary conditions

G(O,x,)=G(l,x,)=O (7a)
where 8(x - x0) denotes the Dirac delta function. Solving Eq. (6) subject to the boundary condition
(7a), we obtain
(e-*%x - I)[1 - e*h’~o-“]
G(x, x0) = &1 ek~‘“-“O’H(xo- x) 1 _ e-*k,

I 2i Lkl(XO-X)H(X _ xo) (e-2k1*o

- 1)11- e2k1(x-1)l (7)
1 1 _ e-*k,

where H( ) is the Heaviside function. Evidently, G(x, x0) = G(x,, x). By multiplying Eq. (4) with
G(x, x,,) and integrating over the interval (0,l) (integration by part is required), the solution of
problem (5) is obtained as

U ..,(x)=~G(x,x,)[u._,-2u,-At’(l-s)+-

with G(x, x0) = ktG( x,x0). The initial condition for Eq. (8) is
0
At2 f(x,, t,,)
1
% (8)

uo(x) = 0
ur(x) =At.g(x)

3. Analyses for three limiting cases

To gain insight into the features of integral expression (8), we analyse three limiting cases:
(i) k,+O(At-+~,f!l ZO);
(ii) k, --, m(At-+ 0,8 # 0);
206 J. Su, B. Tabarrok I Comput. Methods Appl. Mech. Engrg. 142 (1997) 203-214

(iii) k,+w(O+O, At#O).

3.1(i). kl+O(At+m,e #O)

It is convenient to define an operator for the parameter k, approaching zero. Let such an operator be
denoted by A,. Then

A, = lib&t

ekla
_1
We have

A,
[ 11k =ff

and when A, operates on the Green’s function we obtain

A,[G(x,x,)] +(x0 -x)(-2x)[-2(x0 - 1)]/2 +$Z(X -x0)(-2x,)[-2(x - 1)]/2

= H(xO -x)x(x0 - 1) + x()(x - l)H(x - x0)

= (x - xJH(x - x(J + x(x0 - 1)

Set

x0) = (x - x,&7(x - x0) + x(x,, - 1)

G,(x, (9)

This expression is the Green’s function for the following (steady state) boundary value problem

$ =-f(x) (10)
{ u(0) = U(1) = 0

In this limiting case, integral expression (8) reduces to

0. u,+,(x)+ (1 - eb,-,(x) = 44 = - i,’G,(x,x,)f(xo) % (11)

which may be easily derived from Eq. (10) directly. Here, the source term f(x,) now is independent of
time. Therefore, as k, +O, the solutions for unsteady problems reduce to the solutions for steady
problems, as expected. From a physical point of view, parameter k, should be a characteristic
parameter reflecting the features of the unsteady problems (for example, a frequency parameter for
harmonic oscillatory problems), since for a given unsteady motion the time step At (thus k,) should be
bounded in order to obtain a reasonable numerical result. For oscillatory motion, the number of time
steps in every cycle is obviously controlled by the accuracy required. As frequency increases and period
decreases, At should decrease and k, increases. Hence, k, = 0 yields solutions for steady problems,
indicating the continuity of the solutions from unsteady states to steady states as k,*O.

344. k,+a (At-O, 8 ZO)

In this case we define an operator A, as follows

A, =lii’-t
 J. Su, B. Tabarrok I Comput. Methods Appl. Mech. Engrg. 142 (1997) 203-214 207

Thus, for any x E (0,l)

A&,q,) = L[k:G(x,x0)1
=-A,
[
kl h(X-Xo)~(XO _
ye x) + k, e%‘“O-“‘~(x
2
_ xo)
1
Although C.?(x,x0) is not defined at x0 = x (since H(x, - x) is not defined at x0 = x), we may redefine the
above expression such that it holds at X, = x.
Let us investigate the nature of G(x, x0) in the case of k, +m. From the above expression, we have
that

-G,(x, x0) = A, 2
kl e-klix-xOi
[ 1
0 x,#x
=
1 m x0=x

For any interval (x - E, x + E) including x,, = x where E is any positive quantity (certainly, x - E, x + E E
(0, l)), we have

X+E X+E k

I G,(x, x0) dx, = A,

I -f e-kl’x-xo’ d.x,

-
x-e x--E

m
_1 e-1ri dt =
=
I
--m2 I0
e-‘dt = 1

where t= k,(x, -x). It is therefore evident that -G,(x,xo) is the one-dimensional S-function, i.e.

kl
-G’m(~,~o)=A, ~e-kl’x-xo’
[ 1
=8(x,-x) (12)

By the above property of -@x, x0), we may derive a form for Eq. (8) in the limiting case of k, + 03.
For any point x E (0, l), taking the limit of k, + CQof Eq. (8)

u .+lo.=A~~o1d(x,xo)[(U.,-2U.)-*t2(l-B)~-*t~f(XU,t~)]dlo
0

= o1[-6(x, - x)](z&-r - 2uJ dX,

I
=24x) - &-l(X)

Now, it is clear that whether or not the original differential wave equation is homogeneous, the time
marching integral expression reduces to

%+1(X)
=2%(x)-%-l(x) (13)

There are two cases for Eq. (13):

(a) If U,(X) =u.-r(x), then u,+~ (x) = u,(x). This simply reflects the continuous dependence of
physical quantity U(X,t) on the time variable t since k, + 0 is equivalent to At-, 0;
(b) If u,(x) # u,-,(x), then u,+~ (x) # u,(x). This is a singular solution corresponding to the case of
infinite frequency which does not exist physically.
208 .I. Su, B. Tabarrok I Comput. Methods Appl. Mech. Engrg. 142 (1997) 203-214

3.3&i). k, + ~0 (0 --, 0, At Z 0)

As before, taking the limit of k,+m of Eq. (S), we obtain

U n+l(~)=n,j-o’
G(x,x,) u._,-2u.-At’@-B)+ At2 f(xo, r,)
1
h,

=01[-6(xo-x)]
0

(u._,-2u.)-At2+- At2 f(xo, L) ho

I [ 0 I

= 2u,(x) - Us_, + At2 2d2z-’ + At2 f(x, t,)

This is the case when 8 = 0 in Eq. (4). So an explicit algorithm is reached in this limiting case.

4. Heat conduction equation

The same idea may be employed to solve the heat conduction initial and boundary value problem as
follows

1
y$+f(x,t) (O<X<l, t>O)
(14)
u(0, t) = Co(t), u(1, t) = l&(t)
4% 0) = d-4

By discretizing the equation in time, we obtain

e
d2u,+11 1
u n - (1 -e&-f(x, tn+o) (15)
T-g%+1 = -bt

This leads to a boundary value problem of the ordinary differential equation for u,,+~(x)

= F(x)
(16)

with

d2u
-u,-At(l-e)-&Atf(x,tn+o)
I
uo(4= g(x)
where k, = l/X&-%. It is clear that we obtain the same form of boundary value problem for both the
wave equation initial boundary value problem and the heat conduction initial boundary value problem
once the problems are discretized in time. Hence, the Green’s function for problem (16) is as given in
Eq. (7) but with k, = l/G. In the same way, we obtain an integral form of the solution of problem
(16)

(17)
 J. Su, B. Tabarrok I Comput. Methods Appl. Mech. Engrg. 14.2 (1997) 203-214 209

with G(x, x,,) = k:G(x, x,,). The initial condition for Eq. (17) is

44 = g(x)

From the above analyses, we note that the dependence of k, on At is determined by the type of the
governin e uation. For the wave equation k, = l/(AtV%), while for the heat conduction equation
k, = l/ Y-0.At.

5. Advection-diffusion equation

An initial and boundary value problem of the time-dependent advection-diffusion equation is

expressed as follows

J
2

~+Tg=p$ (--m<X<~,?>O)

(18)

l u(-a, t) = 0,
4-5 0) = g(x)
u(y t) = 0

where the parameters 7 and CLare constants. Here, the purpose of the choice for an infinite span is to
show one of the advantages of the time marching integral equation method compared with other
numerical methods.
In contrast to the above sections, we apply a leap-frog algorithm to the advection term and a
generalized 8 method to the diffusion term in problem (18). This gives

du, d’u,_,
(%+1 -zL,)=rz-CC &2 (19)

The boundary value problem of Eq. (19) for u,+~ (x) is written in the standard form as

- q-l) = F(x)
(20)
u,+,(-) = u,+,(c=9
=0
with

uo(x)= g(x)
where k, = 1 /j/-t.
The Green’s function in this case is

1
G(x,x,) = - Te h(~-Xo)~(xo _ x) + & ewO-“‘H(x _ xo)]

1 -k*(x-x01
z--e

24 (21)
The integral form of the solution for problem (20) is

(22)
210 J. Su, B. Tabarrok I Comput. Methods Appl. Mech. Engrg. 142 (1997) 203-214

with

uo(x)= m
G(x,x0) = k;G(x, x0)

If u(x, t) has any order derivatives with respect to x, we may obtain, by integration by part, the
following expansion for u, +1(x)

(23)

with a condition that the above series converges,

6. Numerical discretization and result discussion

Eq. (8) for wave problems, Eq. (17) for heat conduction problems and Eq. (22) for advection-
diffusion problems are discretized by using linear interpolation approximation for u(x, I) on each
element. When the first- or second-order derivatives are encountered, integrations by part are
employed. Then, a two-point Gaussian numerical quadrature is used for calculating the resulting
integrals and thus the solution.
From Figs. 1 and 2, we observe that as time step At decreases, k, increases, the function G(x, x0)
becomes more and more concentrated around x0 = x. From the integral equation we derived, it is clear
that G(x, x0) accounts for the percentage of contribution of integrand at x0 to the solution u,+~ at x.
Thus, such contribution becomes more concentrated near x0 =x as At decreases.
In Section 3, we have shown that as k, +O, Green’s function G(x, x0) approaches steady state
Green’s function G,(x, x0). In Figs. 3 and 4, the variations of -G(x, x,,) with respect to x0 at x = 0.2
and x = 0.5 are shown for various values of k,. It is evident that when k, < 1, G,(x, x0) (steady state
counterpart of G(x,x,)) already is a good approximation for G(x,x,). That is, when a low frequency
(slowly varying) unsteady state problem is to be solved, a quasi-steady approach may be employed, i.e.
in each time step, a steady state approximation is made and solved.
As the first numerical example [12,13], we choose problem (1) with g(x) = sin ~FX.The exact solution
is

1
2.4(x,t) = _rr sin 7rt sin ~FX

43 -G

6 6

0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.1 0.9 I 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.6 0.9 1

Fig. 1. The variation of -C&r, x,,) vs. x0 with various values of k, (x = 0.5).

Fig. 2. The variation of -&x,x0) vs. x0 with various values of k, (x = 0.2).

 J. Su, B. Tabarrok I Compui. Methods Appl. Mech. Engrg. 142 (1997) 203-214 211

0.1 0.2 0.3 0.4 as 0.6 0.7 0.8 0.9

Fig. 3. The variation of -G(x, x0) vs. x0 with various values of k, (x = 0.5).
Fig. 4. The variation of -G(x, x0) vs. x,, with various values of k, (x = 0.2).

-0.1 -

-0.2 -

-0.3 -

Fig. 5. The comparison of numerical solutions and exact solution of the wave equation. Solid line--exact solution; dashed
line-numerical solution, 0 = 0.5; dashdot line-numerical solution, f3= 1 .O.

Fig. 5 gives the comparison of numerical solutions and analytical solution at x = 0.5. The calculation is
carried out with uniform mesh of 50 elements and At = 0.05. From Fig. 5 it is noted that even after five
cycles the present method still provides excellent agreement with the exact solution when 0 = 0.5 (which
leads to an algorithm with second-order accuracy in time).
We calculate, for the second test case [14], the heat conduction initial and boundary value problem
(14) with u,(t) = 1, ul(t) = 0 and g(x) = 0. The exact solution is

m sinnrrx _-n~r~,
u(x,t)=(l---X)--i: ye
n-1

Comparisons of the present numerical solution and exact solution are presented in Figs. 6 and 7. A
good agreement is obtained for a uniform mesh of 20 elements. The time step At = 0.002.
One important goal of most numerical methods is to solve nonlinear problems. Here, we select the
well-known Burgers’ equation, which serves as a one-dimensional model for viscous compressible flows
P51

g++p$ (O<x<l,t>O)
212 J. Su, B. Tabarrok I Comput. Methods Appl. Mech. Engrg. 142 (1997) 203-214

Fig. 6. The comparison of numerical solution and exact solution of the heat conduction equation. Solid line-exact solution;
dashed line-numerical solution. 1: at x = 0.8; 2: at x = 0.6; 3: at x = 0.4; 4: at x = 0.2.

0 0.1 02 0.9 OA 05 0.6 0.7 Od 0.9 I

X

Fig. 7. The comparison of numerical solution and exact solution of the heat conduction equation. Solid line-exact solution;
dashed line-numerical solution.

with the boundary and initial conditions selected to agree with the exact solution

0.1 e-” + 0.5 eeb + e-’

U(X,t) = -n + e-b + e-c
e

where a(~, t) = 0.05(x - 0.5 + 4.95t)/p, b(x, t) = 0.25(x - 0.5 + 0.75t) /P, C(X,t) = 0.5(x - 0.375t) /P
with p = 3 x 10w3. We may simply set f(x, t) = --u au/ ax in the governing equation of problem (14) or
T = u in the governing equation of problem (18) with suitable modifications of the boundary and initial
conditions. The results are given in Fig. 8. A good agreement of the present numerical result with the
exact result is obtained by using a nonuniform mesh of 40 elements. The time step in this case is 0.005.
 J. Su, B. Tabarrok I Comput. Methods Appl. Mech. Engrg. 142 (1997) 203-214 213

U 1.1
141

0.9 -

0.8 -

0.7 -

0.6 -

0.5 -

0.4 -

0.3 -

0.2 -

0.1 .
0 0.1 0.2 0.3 0.4 05 0.6 0.7 0.8 0.9 1

Fig. 8. The comparison of numerical solution and exact solution of the Burger’s equation. Solid line-exact solution; plus and
circle-numerical solution.

7. Conclusions

The most obvious conclusion to be reached for the time marching integral equation method is that a
same standard form of the semi-discretized differential equation is obtained for both the wave
(hyperbolic) and heat conduction (parabolic) equations once the primitive governing differential
equations are discretized in time. The nature of the primitive differential equation is reflected by the
dependence of the parameter k, on time step At. The analysis of the limiting cases indicates that both
the steady integral equation and the explicit algorithm are special cases of the time marching integral
equations. It is noted that there are no inherent difficulties for extending the method to two- and
three-dimensional time-dependent nonlinear problems. In fact, a time marching integral equation
method was developed by Su [12,13] for two- and three-dimensional unsteady transonic flows.
Recently, the method was extended to two-dimensional incompressible thermally coupled viscous flows
[16]. The present work gives general discussions and further development of the -method for one-
dimensional problems.

References

ill I. Fredholm, Sur une classe d’equations fonctionelles, Acta Math. 27 (1903) 365-390.
PI J.L. Hess and A.M.O. Smith, Calculation of potential flow about arbitrary bodies, in: D. Kuchemann, ed., Progress in
Aeronautical Sciences, Vol. 8 (Pergamon Press, London, 1967).
[31 M.A. Jaswon and A.R.S. Ponter, An integral equation solution of the torsion problem, Proc. Roy. Sot. A273 (1963)
237-246.
[41 M.A. Jaswon, Integral equation methods in potential theory, I, Proc. Roy. Sot. A275 (1963) 23-32.
PI F.J. Rizzo, An integral equation approach to boundary value problems of classical elastostatics, Quart. App. Math. 25
(1967) 83-95.
161J.C. Wu and J.F. Thompson, Numerical solutions of time-dependent incompressible Navier-Stokes equations using an
integro-differential formulation, Comput. Fluids 1 (1973) 197-215.
[71 N. Tosaka and K. Onishi, Boundary integral equation formulations for unsteady incompressible viscous fluid flow by
time-differencing, Engrg. Anal. 3 (1986) 101-104.
181 G.F. Dargush and L.K. Ganerjee, A boundary element method for steady incompressible thermoviscous flow, Int. J.
Numer. Methods Engrg. 31 (1991) 1605-1626.
191 CA. Brebbia, The Boundary Element Method for Engineers (Pentech Press, London, 1978).
WI CA. Brebbia and S. Walker, Boundary Element Techniques in Engineering (Newnes-Butterworths, London, 1980).
ill1 O.K. Zienkiewiu, Marriage a la mode-finite element and boundary integral method, in: International Symposium on
Innovative Numerical Analysis in Applied Engineering Science, Versailles, CETIM Publications, Paris, 1977.
214 J. Su, B. Tabarrok I Comput. Methods Appl. Mech. Engrg. 142 (1997) 203-214

(121 J. Su, Time marching integral equation method for unsteady transonic flows around airfoils, Int. J. Numer. Methods Engrg.
34 (1992) 455-467.
[13] J. Su, Unsteady transonic flow integral equation methods, Ph.D. Dissertation, Beijing University of Aeronautics and
Astronautics, Sept., 1989.
[14] P.M. Gresho and R.L. Lee, Don’t suppress the wiggles--they’re telling you something, Comput. Fluids 9 (1981) 223-253.
(1.51 G.F. Carey and J.T. Oden, Finite Elements: Computational Aspects, Vol. 3 (Prentice Hall, Englewood Cliffs, NJ, 1984).
[16] J. Su, B. Tabarrok and S. Dost, A time marching boundary element method for two-dimensional incompressible thermally
coupled flows, Comput. Methods Appl. Mech. Engrg. 128 (1995) 11-23.