The equations Laplacian operator omega = 0, (1.1a) and omega = Laplacian operator Chi, (1.1b) describe, in suitable units, 2-D Stokes flow of an incompressible fluid occupying a domain D in which omega is the vorticity and Chi is the stream function. The flow is uniquely determined by specifying the velocity on the boundary B of D, a condition which leads to specifying the stream function Chi and its normal derivative Chi sub n on B. A mathematically similar problem arises in describing the equilibrium of a flat plate in structural mechanics where a related 1-D problem by finite difference or finite element methods is to introduce effective methods for imposing the boundary conditions through which (1.1a) is coupled to (1.1b). These models thus provide a simple starting point for examining the general treatment of boundary conditions for more general time dependent Navier-Stokes incompressible flows. For the purpose of discussion it is assumed that D is a square domain. A standard finite difference method to solve (1.1) is to introduce a uniform grid and then use standard five point finite difference operators to express each equation in (1.1). At any point on the boundary B a value of Chi is specified by the boundary conditions but a value of omega at the same boundary mesh point will also be required to complete the computation. Methods are discussed which overcome the difficulty in solving these problems

Rose, M. E.

English

NASA Technical Reports Server

A Finite Difference Treatment of Stokes-Type Flows:
(Preliminary Report)
M.E. Rose
Department of Mechanical Engineering
North Carolina A&T State University
(NASA Contract No. NAG- 1-812)
1.Introduction
The equations
V2co = O, (1. la)
to - -V2X (1. Ib)
describe, in suitable units, two-dimensional Stokes flow of an incompressible fluid
occupying a domain D in which co is the vorticity and X is the stream function. The flow is
uniquely determined by specifying the velocity on the boundary B of D, a condition which
leads to specifying the stream function X and its normal derivative Xn on B. A
mathematically similar problem arises in describing the equilibrium of a flat plate in
structural mechanics where a related one-dimensional problem describes the equilibrium of
a clamped beam. A key to treating these simple problems by finite difference or finite
element methods is to introduce effective methods for imposing the boundary conditions
through which (1. la) is coupled to (1. lb). These models thus provide a simple starting
point for examining the general treatment of boundary conditions for more general time-
dependent Navier-Stokes incompressible flows.
K[-P'_d,T (_ior,;:_ _ .-_Jrolin,, A_riculfura] qnd
i,_chnit._l Jf._t;, Univ. ) ii 7 Cqt.L ?0 r3
N90-27970
r. 31_a
Sat, Sep 16, 1989
https://ntrs.nasa.gov/search.jsp?R=19900018654 2020-03-19T21:04:42+00:00Z
For the purpose of our discussion we may assume D is a square domain. A
standard fmite difference method to solve (1.1) is to introduce a uniform grid and then
employ standard five-point finite difference operators to express each equation in (1.1). At
any point on the boundary B a value of× is specified by the boundary conditions but a
value of 0o at the same boundary mesh point will also be required to complete the
computation. Peyret and Taylor [ 1] review the use of extrapolation methods to achieve this
using Taylor series arguments. However, this technique can be expected to be of limited
value when time-volume methods are used to treat curved boundaries or when geometrical
singularities arise when using curvilinear coordinates. As we shall also see, its use can be
expected to result in a loss of accuracy even in simple cases. The method discussed in this
report can be expected to overcome these difticultites.
2. A Compact Difference Scheme
We first describe a compact finite difference scheme for solving _" "= g which has
been described in Rose[2]. It is a specialization to one-dimension of a more general finite
volume scheme for solving div u = g, grad _ = u on general domains.
The scheme expresses a relationship between certain primaryvariables _i, u-(xi),
u+(xi) which are associated with the endpoints of the non-overlaooin_ intervals (x i, xi+ 1),
i = 0,1 ,...,M-1, and dualvadables _ which are associated with the interval centerpoints
xj, j = 1/2,3/2,...,M-1/2, which we regard as forming a dual (orstaggared) grid. (In the
following, the index i will be associated with the primary grid, while j will indicate the dual
grid.) Indicating the cell endpoints by the symbol x and the cell midpoints by the symbol
o, the points form the pattern
XOXOXOXOXOX .... XOX
Sat, Sep 16, 1989
totalling 2M+ 1 points.
For uniform grids the scheme to solve 4)"= g can be described as follows: write the
equation as the first order system u' = g, 4)" -- u and express the first equation by the
difference equations
[ U-(Xj+l/2) - U+(Xj_I/2)] lax = gj j = 1/2, 3/2,..., M-1/2
(2.1)
in each interval of length Ax = 2h. We interpret the second equation as relating u-and u ÷
with forward and backward differences involving the variables 4) at primary and secondary
points of the grid:
u-(xi) = (_-__l/2)/h, i = 1,2,..., M (2.2)
u+(xi ) = (_+1/2- _)/h, i = 0,1,..., M-1
Introduce the central average and difference operators
ta _- ( 4)j+ 1/2 + _-1/2 )/2' A _j--- ( 4)j+1/2 - _-1/2)"
and impose the condition that the variables u-(x i) and u+(xi) be continuous at interior
endpoints of the primary grid, i.e.,
ui -- u-(xi) = u+(xi), i = 1,2, .... M-1. (2.3)
Using the definitions (2.2) in (2.1) we find, with 1<-- 1/Ax 2,
gj = 41< (la4)j - 4)j), j --- 1/2, 3/2,..., M- 1/2
whike use of the continuity conditions (2.3) leads to
0 =/a M -_ i = 1,2,..., M-1
This tri-diagonal system is to be solved for the variables 4) with specified boundary
conditions of either Dirichlet or Neumann type.
(2.4)
(2.5)
A cyclic (odd-even) reduction technique leads to the system of equations
_tg i = K A2_ i, i = 1,2,..., M-1 (2.6)
3
Sat, Sep 16, 1989
for the primary variables and, separately, to
gj = 1<A24pj, j = 3/2,,5/2.., M-3/2 (2.7)
for the dual variables. Under Dirichlet-type boundary conditions, the solution of the first
set of equations can be solved directly for the primary variables; the values qbl/2 , ¢_vI-1/2
for the dual variables can be found in terms of these primary values by solving each of the
equations
gj = r A2_j, j = I/2, M-1/2 (2.8)
and these, in turn, can then be used to provide Dirichlet data to solve (2.7). This reduction
is less useful when Neumann-type data are imposed; in this case it is best to solve (2.4)-
(2.5) directly.
The following table compares numerical results obtaained by solving _" "= g by a
standard finite difference scheme and by the compact scheme just described.
Error Norm Comparison of Standard and Compact Schemes
for the Equation ¢"=g
¢=x(1 -x):
(The precision of results for this example is questionable because of machine limitations)
endpoint error norms
# intervals solution
standard compact
derivative solution derivative
6 2.77556e-17 .166667 5.55112e-17 5.55112e-17
12 5.55112e-17 8.33333e-2 8.32667e-17 3.33067e-16
24 1.11022e-16 4.16667e-2 2.49800e-16 9.99201e-16
48 2.63678e-16 2.08333e-2 1.83187e-15 5.77316e-15
4
Sat, Sep 16, 1989
¢ = X2(1-X)2:
endpoint error norms
standard compact
# intervals solution derivative solution derivative
6 6.17284e-3 6.48148e-2 1.23457e-2 1.85185e-2
12 1.7361 le-3 4.16667e-2 3.47222e-3 9.25926e-3
24 4.34028e-4 2.40162e-2 8.68056e-4 2.89352e-3
48 1.08507e-4 1.62399e-2 2.17014e-4 7.95718e-4
3. The Clamped Beam Problem
The deflection _(x) of a uniform, straight beam under a load -fix) per unit length on an
interval [1_, 1+] is, in dimensionless form, governed by the simple fourth order differential
equation
d_.... = -f (3.1)
If the slope u, bending moment v, and shear force w are given by
u = _', v = -u" , w = v" (3.2)
then typical well-posed boundary conditions allow one to prescribe pairs of values among
(qb,u,v,w) at each endpoint of the interval. We may also write (3.1) as a coupled system of
second-order equations for _ and v
(a) v'" = f,
(b) dp"'= -v.
(3.3)
For a problem in which the boundary conditions involve the pair of values (_, v) at each
endpoint, for example, a simple Green's function technique allows (a) to be solved for v(x)
in terms of f(x) together with the boundary conditions v(l+), while a similar construction
gives _(x) in terms of the boundary conditions qb(l+) and v(x). However, for clamped
Sat, Sep 16, 1989
boundary conditions, qb= u = 0 at both endpoints and this Green's function construction
fails.
A standard finite difference technique for handling clamped boundary conditions
involves the use of a Taylor series expansion at each endpoint in order to express v(1 +) in
terms off(l+) and u(l+), both values of which are prescribed, as well as one or more
values ofd_ at interior mesh points of the interval. We shall describe this technique first,
using a rather direct finite difference argument which is suggested by a method described
in Peyret and Taylor [ 1] in connection with a treatment of the Navier-Stokes equations in
vorticity-stream function variables. We then discuss and illustrate another, closely related,
technique which arises from an application of a compact finite difference method (Rose[2])
to this problem.
4. Some Standard Finite Difference Approaches.
We adopt the standard finite-difference notations x i -- i Ax, h -- Ax/2, u(x i) = u i.
Divide [1_, 1÷] into M non-ovedappin_ intervals Ij -= {xl Xj_l/2 < x < xj+ 1/2}with
centerpoints xj, j= 1/2, 3/2,..., M- 1/2. Also, recall the central average and difference
operators
la qbj---( ¢j+1/2 + _)j-1/2 )/2, A (_j-- ( _+1/2- (_j-1/2).
introduced earlier.
(4.1)
With t<= l/Ax2,a standard finite difference approach to solving (4.3) is to consider
the coupled difference equations
(a) fi= _A2vi
(b) -v i = K A2_
i = 1, 2, ..., M-1 (4.2)
6
Sat, Sep 16, 1989
eachof which may be separately solved with Dirichlet-type data. For the damped beam,
the boundary conditions _1 +) = 0, u(l +) = 0 translate into
(a) _ = _ = 0, (4.3)
(b) u0 = U M = 0.
Were v i known, (4.2b) could easily be solved for _ under either pair of these boundary
conditions by a standard tri-diagonal solver (in case (b) we can add thc condition ([_0= 0)
In order to solve (4.2a) we will specify values ofv 0 and v M Consistency requires that
thesc values be related to values of_ and u at points on or ncar thc boundary points and we
may write
v0 -- B0(qb,u), (4.4)
V M = BM(qb,u )
where B is a suitable boundary operator which incorporates the prescribed boundary
conditions for.dp and u.
One method to obtain a boundary operator B is to use a Taylor series approximation
as follows: Write
dpl =qb0 + Ax u 0 - L(Ax) 2 v0 +... (4.5)
2
_'l-1 = _ - &X U M - L(Ax) 2 v M +...,
2
so that by imposing the homogeneous clamped beam boundary conditions on _ and u we
find
2 _1 = - (Ax)2 v0 +""
2 qbM_ 1 = - (Ax)2 VM + ....
Higher order extrapolations are discussed in Peyret and Taylor[ 1].
(4.6)
7
Sat, Sep 16,1989
5. A Compact Scheme for the Clamped Beam.
Our objective will be to describe a difference scheme based upon these ideas which
solves (1) under the clamped beam boundary conditions stated in (4.3), i.e.,
t_0 = _vt = 0, u0 = uM = 0. The boundary conditions forv will not require the additional
use of a Taylor series other than that which is implicit in the proposed difference equations.
We will consider the coupled system
-!tlv i = r A2_,
=,,a2vj,
The first of these corresponds to (4.5) for the primary variables associated with a compact
scheme for solving dp'"= -v and the second corresponds to (4.6) for the dual variables
associated with a compact scheme for solving v'" = f.
Using the definitions ofu given by (4.2) the values ofv at the endpoints of the dual
i = 1,2,..., M-1 (5.1)
j = 3/2,,5/2.., M-3/2.
grid are found to satisfy
- Ax vj = (uj+ 1/2 - uj-1/2) j = 1/2, M- 1/2 (5.2)
and imposing the Neumann-type boundary conditions for u gives
u 1 = - Ax Vl/2, UM_ 1 = Ax VM_l/2.
Using the definition ofu leads to
21,: (qb1 - dPi/2)= - Vl/2 (5.3)
2K (¢M-1/2 - ¢_I-1 )= VM-1/2"
One difference between this and the treatment of boundary conditions for the
standard finite difference method described earlier lies in the position on the mesh of at
which the boundary values assigned to v are imposed. In the standard problem v i was
considered a primary variable and a Taylor expansion was required in order to fumish an
additional equation which allowed boundary values of v at the endpoints of the domain to
Sat, Sep 16, 1989
be determined bye. In the compact scheme, on the other hand, the vj are dual variables
and the boundary conditions, which determine their values at the boundary of the dual
mesh, are a consequence of applying the difference equations lay i -- - K A2_ when
i = 1/2, M-1/2. The cyclic reduction technique allowed the dual variables to be solved
directly with such data.
This distinction between variables defined on a primary and dual grid is a common
feature in treating the Navier-Stokes equations in primitive variables where the pressure
term is commonly associated with the dual grid. This connection is worth further
exploration.
Numerical Examples for the Clamped Beam.
The following table compares numerical results for the clamped beam problem
using the standard finite difference method with extrapolation techniques to set the
boundary conditions and the compact scheme described above.
Error Norm Comparison of Standard and Compact Schemes
for a Clamped Beam
test solution : qb=x2(l-x)2:
Standard scheme using 1st order extrapolated BC;
# intervals f u w
12
24
48
96
1.07692 .213675 1.92901e-2
.324713 7.24377e-2 1.17071e-2
.124513 5.87587e-2 7.65804e-3
5.41195e-2 3.55117e-2 4.37644e-3
Sat, Sep 16, 1989
Standard scheme using 2nd order extrapolated BC;
# intervals f u w
12 2.80765e-2 .14098 .863799
24 1.77672e-2 9.04946e-2 .313276
48 9.78894e-3 6.49779e-2 .1223
96 5.07409e-3 3.71921e-2 5.17336e-2
Compact scheme ;
# intervals f u w
12 2.82377e-2 3.8407e-5 6.94711e-3
24 6.9977e-3 2.25417e-6 1.73624e-3
48 1.74305e-3 1.58858e-7 4.34032e-4
96 4.35881e-4" 1.39912e-7" 1.08463e-4"
*precision doubtful because of machine limitations.
Note that, as predicted, the compact scheme fumishes second-order accuracy for
the variables.
6. A Time-dependent Stokes-type Problem
The time-dependent equations
oat = V2oa,
to = -V2x
can serve as a model for studying the effect of handling boundary conditions and differs
10
Sat, Sep 16, 1989
fromtheNavier-Stokes equations in that only the convective terms have been omitted. A
method for adapting ADI techniques to solve this problem by a compact scheme is outlined
in [2].
A code to test the accuracy of solving this problem by a compact scheme has been
developed by J. M. Klimkowski and is reported upon separately.
References
[ 1] Peyret, R. and Taylor, T.D., Computational Methods for Fluid Flow, Springer-
Verlag, 1983.
[2] Rose, Milton E., Compact Finite Volume Methods for the Diffusion Equation,
submitted to J. Comp. Phys.
11
Sat, Sep 16, 1989


A finite difference treatment of Stokes-type flows: Preliminary report

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server