Citation
Fluid motion due to a bending wave

Material Information

Title:
Fluid motion due to a bending wave
Creator:
Raskob, Jr., Anthony William
Publication Date:
Language:
English
Physical Description:
xi, 73 leaves : ; 28 cm

Subjects

Subjects / Keywords:
Hydrodynamics ( lcsh )
Fluid mechanics ( lcsh )
Navier-Stokes equations ( lcsh )
Wave-motion, Theory of ( lcsh )
Fluid mechanics ( fast )
Hydrodynamics ( fast )
Navier-Stokes equations ( fast )
Wave-motion, Theory of ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaves 72-73).
General Note:
Department of Mechanical Engineering
Statement of Responsibility:
by Anthony William Raskob, Jr.

Record Information

Source Institution:
University of Colorado Denver
Holding Location:
Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
41471095 ( OCLC )
ocm41471095
Classification:
LD1190.E55 1998m .R37 ( lcc )

Full Text
Fluid Motion Due to a Bending Wave
by
Anthony William Raskob, Jr.
B.S., University of Colorado, 1982
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Mechanical Engineering
1998


1998 by Anthony William Raskob, Jr.
All rights reserved.


This thesis for the Master of Science
degree by
Anthony William Raskob, Jr.
has been approved
Ken Ortega
/cV
Date


Raskob, Anthony William Jr. (M.S., Mechanical Engineering)
Fluid Motion Due to a Bending Wave
Thesis directed by Professor Samuel Welch
ABSTRACT
The fluid motion caused by a bending wave on the boundary of a viscous, incom-
pressible fluid is presented and analyzed. The motion is input in the form of a stand-
ing wave applied to the boundary of a semi-infinite fluid. A scaling law analysis
of the Navier-Stokes and continuity equations is presented which gives the con-
ditions under which a linearized form of the Navier-Stokes equations may be ex-
pected to yield useful approximate solutions. A linearized solution is obtained and
presented. A numerical solution to the full Navier-Stokes equations is obtained us-
ing the commercial solver FLUENT running an implementation of the SIMPLE al-
gorithm. Comparison of the numerical solution to the linearized solution reveals
significant differences in the near-boundary fluid motions. These differences are
attributed to the non-linear fluid behaviour not accounted for in the linearized equa-
tions.
This abstract accurately represents the content of the candidals
mend its publication.
Signed
IV


CONTENTS
Chapter
1. Introduction 1
1.1 Problem Statement................................................ 1
1.2 Previous Research by Others...................................... 3
1.2.1 Inviscid fluid................................................... 3
1.2.2 Viscous fluid 5
2. Linear Theory Solution 7
2.1 Scaling Analysis 7
2.1.1 Non-dimensionalization using bending wavelength.................. 8
2.1.2 Non-dimensionalization using viscous wavelength................. 11
2.2 Linearized Theory Results 11
2.2.1 Solution to linear momentum equations........................... 12
2.2.2 Fluid behaviour according to linear solution 14
2.3 Consideration of Non-linear Effects............................. 21
v


3. Numerical Solution
23
3.1 SIMPLE Algorithm................................................ 23
3.2 FLUENT Implementation of Bending Wave Problem................. 23
3.2.1 Selection of mesh parameters.................................... 24
3.2.2 Use of symmetry; mesh boundary conditions...................... 26
3.2.3 Moving mesh forcing function 26
3.3 FLUENT Results.................................................. 28
3.3.1 Grid refinement................................................. 28
3.3.2 Time step refinement............................................ 30
3.3.3 Starting transient 30
3.4 Comparison of Linear Theory and FLUENT Results................ 32
3.4.1 Velocity profile results 32
3.4.2 Vector field results 38
4. Conclusions 42
4.1 Directions for further research 42
Appendix
A Development of Linearized Solution using Mathematica 43
B. NEWMESH Subroutine Listing 69
vi


FIGURES
1.1 Semi-infinite fluid bounded by undulating plate......................... 1
1.2 Sectional view of undulating boundary 2
1.3 Near-field coupling of sub-acoustic plate bending waves 3
1.4 Bending wave coupling to fluid.......................................... 4
2.1 Expression for real portion of u solution to linearized equation set 15
2.2 Pressure contours calculated using linearized theory. 16
2.3 Wall-normal (v) velocity contours calculated using linearized theory. 17
2.4 Wall-normal (u) velocity profile calculated using linear theory (near
wall region)..................................................... 17
2.5 Wall-normal (v) velocity profile calculated using linear theory (over-
all)........................................................... 18
2.6 Spanwise (u) velocity profile calculated using linear theory (near
wall region)..................................................... 19
2.7 Spanwise (u) velocity profile calculated using linear theory (overall). 19
2.8 Spanwise (u) velocity contours calculated using linear theory (near
wall region)..................................................... 20
vii


2.9 Spanwise (u) velocity contours calculated using linear theory (over-
all)................................................................... 20
2.10 Velocity vector field calculated using linear theory (overall). 21
2.11 Velocity vector field calculated using linear theory (near wall region). 22
3.1 Meshing scheme employed in fluent.............................. 25
3.2 Boundary conditions used with fluent mesh. 26
3.3 Geometry updated after solution of momentum and pressure equa-
tions.................................................................. 27
3.4 FLUENT solution time histories for various degrees of mesh reso-
lution (u velocity component, x = Ab/2)................................ 29
3.5 FLUENT solution time histories for various degrees of mesh reso-
lution (v velocity component, x Ab/4)................................ 29
3.6 FLUENT solution time histories for various degrees of mesh reso-
lution (v velocity component, x = Ab/2)................................ 30
3.7 FLUENT solution time histories for various degrees of time step re-
finement............................................................... 31
3.8 FLUENT solution time histories for various degrees of time step re-
finement (detail)...................................................... 31
3.9 Difference in first and second cycle u velocities (y = 7.2AS). ... 32
viii


3.10 Comparison between first and second cycle u velocities (y = 0.2AS). 33
3.11 Difference in first and second cycle u velocities (y = 0.2AS). 33
3.12 Comparison of linear theory solution and fluent result, u velocity
component at 180. 34
3.13 Comparison of linear theory solution and FLUENT result, u velocity
component at 315. 34
3.14 Comparison of linear theory solution and FLUENT result, v velocity
component at 180. 35
3.15 Comparison of linear theory solution and FLUENT result, v velocity
component at 315. 35
3.16 Comparison of linear theory solution and fluent result (dashed),
u velocity component near 270, x = Ab/2................ 36
3.17 Comparison of linear theory solution and FLUENT result, near-wall
(y = 0.2Ab) u velocity component versus time............ 36
3.18 Comparison of linear theory solution and FLUENT result, u velocity
component near 270, x = Xh/4........................... 37
3.19 Comparison of linear theory solution and FLUENT result, v velocity
component near 270. x = Ab/2............................. 37
3.20 Schematic illustration of convection-driven near-wall motion. ... 38
IX


3.21 Linear theory solution velocity field at several points through wall
motion cycle....................................................... 39
3.22 FLUENT results region showing near-wall effect. 39
3.23 FLUENT calculation near-wall velocity field at several points through
wall motion cycle (0-75)......................................... 40
3.24 FLUENT calculation near-wall velocity field at several points through
wall motion cycle (90-l 65)...................................... 41
x


TABLES
2.1 Typical dimensional values for bending wave experiments........... 10
3.1 Levels of mesh resolution investigated. 25
3.2 Gas law properties of FLUENT. 26
xi


Chapter 1
Introduction
This work considers the problem of fluid motions induced by wave motion on
the surface or boundary of a viscous, incompressible fluid. Specifically, it considers
a fluid in contact with a vibrating solid boundary. Because the solution of this prob-
lem pertains to fluid motions caused by bending waves on immersed plates, for ease
of reference I refer to the boundary vibration as a bending wave. Thus, I consider
a semi-infinite fluid, bounded on one side by a boundary which is undulating in a
manner similar to that of a vibrating plate (Figure 1.1).
Figure 1.1: Semi-infinite fluid bounded by undulating plate
1.1 Problem Statement
The bending wave motion takes the form of a single mode wave with charac-
teristic wavelength A and radian frequency of vibration u;. In general a sinusoidal
wave will be assumed, although other shapes could be considered. Because only a
single mode of vibration is considered, and because of the semi-infinite geometry
1


| y deformed shape
Figure 1.2: Sectional view of undulating boundary
of the problem, coordinate axes can be chosen without loss of generality as shown
in Figure 1.1, with the x axis in the plane of the boundary and oriented perpendic-
ular to the crests and troughs of the wave pattern, the y axis perpendicular to the
plane of the boundary, and the ~ axis in the plane of the boundary and aligned with
the waves crests and troughs. A section of the x,y plane at the boundary is shown
in Figure 1.2. With this choice of coordinate system, the boundary conditions for
the x direction velocity component (it) and y direction velocity component (v) can
be expressed:
u(.r, y,t)\y=d = 0,
i>(.r,t/,i )!>,=,/ = v0g(x) cosut,
with
d(x,t) = d0g{x)smu>t
being the out of plane deflection of the boundary, d0 representing the maximum out
of plane deflection, and g(x) is the mode shape function. Note that v0 = d0uj. The
first boundary condition arises from the no-slip condition, while in the second is
contained the assumption that the boundary is non-porous. I am considering a two
dimensional formulation only, therefore iv = 0 for simplicity.
Of particular interest is the motion very close to the wall. Experiments have
shown that the undulating motion described above can effect significant changes in
the structure of the near wall boundary layer. It is in part the purpose of this work to
aid in understanding this effect. Since the lengthscales important in viscous flows
may span a large range and may include very small lengths, it will be necessary to
accurately model the near wall motions. This in turn necessitates accurate consider-
ation of these viscous effects which will govern the fluid behaviour in regions close
to a solid boundary.
2


1.2 Previous Research by Others
This problem touches on several areas of interest, and various aspects of it have
been previously considered. Working in the area of sound radiation and control,
Cremer [3], Brillioun [1], and others [8, 11] have developed formulations for the
far and near field radiation patterns of undulating boundaries. As will be shown in
the discussion of Section 1.2.1, these fail, however, to address the effects of fluid
viscosity. Likewise, Segal [16] has extensive treatment of fluid surface waves, but
employs an inviscid formulation, while Lighthill [7] considers both surface and in-
ternal waves, but with limited treatment of viscous effects.
Lamb [5] on the other hand does consider internal and surface waves, including
the effects of viscosity, as discussed in Section 1.2.2. As will be shown, his treat-
ment of free surface waves needs to be modified to account for the solid boundary
I am considering.
1.2.1 Inviscid fluid
is shorter than the corresponding Figure L3: Near-field coupling of sub-
acoustic wavelength, the vibration acoustlc Plate bending waves
is sub-acoustic, and creates near-
field pressure oscillations in the fluid. These oscillations are partly in a direction
perpendicular to the local surface normal of the plate or shell. Thus, for properly
oriented bending waves in a shell, spanwise near-field oscillations result. In addi-
tion, there is motion of fluid in the wall-normal direction, so that the coupled bend-
ing energy produces near-field motions having an elliptical profile (Figure 1.3).
For an infinite plate, an expression for the normal velocity due to the existence
of a bending wave (Figure 1.4) can be given by (making use of Eulers identity e'e =
cos 6 + i sin 6, and use of the real portion of each expression for extracting physical
Cremer et al. [3] and others [8,
11] have derived expressions for
the pressure field created in a fluid
by the action of sub-acoustic vibra-
tion modes in a structure in con-
tact with the fluid. Vibration in
the structure generates vibrations
of the same frequency in the fluid,
and if, for that frequency, the wave-
length of the structural vibration
near-field oscillations induced in fluid
surface normal (exaggerated)
3


X
Figure 1.4: Bending wave coupling to fluid
properties is implied):
rjy(j) = I>0e tkBl.
Cremer [3] gives for the sound pressure radiated into an ambient medium or fluid
in contact with the plate an expression of the form:
p(.v. y) = p0e~kBTe~lkyy. (1.1)
He also gives an expression for the wave equation in the fluid:
V2p + k2p = 0. (1.2)
Equation 1.1 must satisfy the wave equation, and at the same time, the normal com-
ponent of the fluid velocity at the plate surface must be equal to the corresponding
normal velocity of the surface (for a non-porous surface with no voids in the fluid).
Substituting Equation 1.1 into Equation 1.2 and enforcing the plate-fluid bound-
ary condition yields information about the sound pressure above the plate and the
wavenumber ky, namely:
p(x,y)= , Vol>c e~il(1.3)
yj1 *i/*2
with
K = k2 kl
One can clarify the physical meaning of Equation 1.3 by introducing an angle 6
at which the sound is radiated. This angle must be such that the trace1 of the acoustic
wave along the plate has the same wavelength as the plate bending wave [3], or
M.e., the projection of the acoustic wave onto the structure
4


sin 6 = A/AB (i.e., sin 9 = k^/k). When AB > A, Equation 1.3 describes a set of
plane acoustic waves (Figure 1.4) radiating from the plate. However, when AB < A,
there is no real angle at which sound is radiated, and Equation 1.3 instead describes
a pressure field which decreases with increasing distance from the plate. In this
case, the particle motions can be found to have x and y components:
= !kBVo
yJk'Q k2
vy = v0e-lkBl'e-'/kl-k2y.
This may be interpreted as indicating that the fluid avoids acoustic compression (re-
quired for the generation of a sound wave) by moving out of the way laterally, re-
sulting in a back and forth sloshing of the fluid [3].
The above represents a potential function solution to the wave equation, and
as such, does not enforce the no-slip condition at the wall. Lamb [5] considers a
combined potential and stream function solution which more correctly describes
the near-wall spanwise motion, as well as accounting for fluid viscosity, but for a
certain case of boundary conditions representative of water waves at an air inter-
face.
1.2.2 Viscous fluid
Lamb begins with the momentum equations linearized in this form:
du 1 dP
dt
dv
di
p dx
1dP_
p dy
+ i/V n,
g + uV2v
in which g is the gravitational constant and v is viscosity. By setting
u
v
p
p
d(j>
dx
d d'tk
~ dy'
dy + dx
dcj)
dt
- gy-
5


where and 0 are potential and stream functions, respectively, the equations of
momentum can be cast in the form of the auxiliary conditions:
V26 = 0,
^ (1.4)
dt
Assuming a solution periodic in x with a prescribed wavelength A = 2k/k and a
time factor ent, Lamb gives solutions of (1.4) as:
0 = {Atky + Be-kv)e{ikT+nt),
0 = (Cemy + De-my)e(lkx+nt).
provided m2 = k2 + n/u.
In his particular case, Lamb deals with water waves at a water/air interface, and
as such his boundary conditions, which include the effects of surface tension, are
inconsistent with the problem I am considering. However, in Section 2.2 I again
take up Lambs analysis and apply the boundary conditions of my problem to obtain
a pertinent solution.
Before leaving Lambs solution, I would like to make reference to the quantity
n/u which appears in the equation for the wavenumber m used in the expression for
0. This quantity appears in the solution to Stokes first and second problems [17]
as the square of the wavenumber representing the spatial decay of penetration of
viscous motion into the fluid due to spanwise boundary motion. As such I refer to
this wavenumber as ks = \Jnjv, that is the Stokes wavenumber. Evidently the
solution proposed by Lamb possesses similar viscosity-related lengthscale charac-
teristics, and it will be shown that the bending wave problem does as well.
6


Chapter 2
Linear Theory Solution
To develop my solution, I begin with a scaling analysis of the governing equa-
tions to investigate the relative importance of each term (Section 2.1). From this, I
extend the work of Lamb to account for my particular set of boundary conditions,
and solve a linearized set of the governing equations (Section 2.2). Finally, 1 briefly
discuss the limitations of my solution with reference to its exclusion of the nonlin-
ear terms (Section 2.3).
2.1 Scaling Analysis
In his derivation, Lamb drops the nonlinear terms of the Navier-Stokes equa-
tions for mathematical convenience. I performed a scaling analysis to more closely
examine the effect of dropping these terms on the solution for the flow I am con-
sidering.
Of primary importance in non-dimensionalizing the problem is the choice of
lengthscale. The work of Cremer and Lamb both show the role of the bending
wavelength At = 2-k/k^ in defining lengthscale in both x and y. In addition, Lamb
shows the role played by the viscous lengthscale As = 2n/ks. In my analysis, I
consider these two governing lengthscales in turn first the bending wavelength
to determine the essential momentum equation terms for the fluid region y > At,
and second the viscous lengthscale to do likewise for the near-wall region y < As.
To make use of the scaling results it will be necessary (eventually) to substitute
dimensional quantities so that the relative significance of the terms in the momen-
tum equations may be assessed. For this I will make use of conditions which have
previously been used in experiments, as will be detailed below.
7


2.1.1 Non-dimensionalization using bending wavelength
I begin with a standard expression for the Navier-Stokes equations in two di-
mensions:
du dn du
dv dv dv
m+uai + vdi,
i dP
-^-+9 + v
p dx
1 dP
-+g + V
P oy
d2u d2u\
dx2 dy2) '
d2v d2v\
dx2 dy2 J
I normalized all quantities to the following scales (the starred terms represent di-
mensionless quantities):
U = Uoll
V = V0V*
y = hym
x = Abx
t = Tt*
P = -P\
T
with the choice of pressure normalization discussed subsequently. This results in:
du* uqt xdu*
dt* Xb dx*
dv* vqt , dv*
dt* Xb U dx*
UqT Mdu
+ -V
A b
dy*
dv
, V0T .
+ ~T~V TT~~
At. dy
v dP
u0Xb dx*
v dP
v0Xb dy
r l (d2u* d2u*
+ 9 u 0 + ReAb \<9.r*2 + dy*2
T 1 / d2v d2v*'
+ g Vo + ReAb V<9.r2 + dy*2 /
(2.1)
The velocities u0 and v0 are generated by the motion of the wall in y as it oscillates
with period r. Thus, for example, u0 = d0/r, where the amplitude of the wall
motion is represented by d0. By continuity, u0 ~ t>0, so that the Equations 2.1 can
be written:
du*
W*
+ Ell
dv
du ,du*
d7*+£v dr
dv ndv*
dx* dy*
VT dP T2 1 / d2ll' d2u\
Xbd0 dx* + 9 d0 + Re.\b \d.r*2 ^ dy*2)
vt dP t2 1 / d2v* d2i>*\
At^o dy + 9d0 + Re.\b \9x*2 + dy*2) 7
(2.2)
8


where:
do
A,-
From the above it can be readily seen that if the normal displacement of the
boundary dQ is small with respect to the bending wavelength (i.e. £ the convective terms are negligible in comparison to the inertial and viscous terms
(provided Re,\t is not too large). The body force terms will be negligible provided
t2/d0 < (). The pressure terms are negligible if ut < Afec/0. This ability to ne-
glect the nonlinear terms if d0 is sufficiently small was presumed by Lamb in his
development.
At this point I would like to recognize that there are two other choices for nondi-
mensionalizing the pressure terms namely to use the stagnation pressure brought
about by the fluid velocity, or a pressure derived from the fluid inertia terms. In
other words, use
P Piuo + vo)P*
in the case of stagnation-dominated pressure or use
p PuXb p
T
in the case of inertially driven pressure. The choice between this form of pressure
normalization and the viscous form used in obtaining Equations 2.2 is made based
on the relative magnitude of the inertial, stagnation and viscous pressures. That is,
based on which is larger; piioXb/r, p(ul + Cq) or p/r. Assuming u0 ~ u0> and that
the peak velocity is due to the amplitude of the wall motion divided by the period
(u0 = do/r), these expressions can be cast in the forms shown below.
inertial stagnation viscous
d0\h T 2^ T V
0.012 m2/s 80 x 106 m2/s 15.9 x 10"6 m2/s
To make this determination it will be necessary to consider typical values for
v, d0 and t. Table 2.1 presents a nominal set of values for these and the other im-
portant variables, as determined from conditions found to be of interest in the pre-
viously mentioned wind tunnel experiments. From these, it can be seen that in fact
9


Table 2.1: Typical dimensional values for bending wave experiments.
Variable Value
d0 1 O X TS
0.06m
T O.OOlsec
V 15.9 x 10_6m2/s
1.3 x 104m
( A s j A b 0.002
£ 0.003
1.5
Redo 2.5
ReA 225, 000
ReA> 1.1
r-/d0 0.005s2/m
v is smaller than 2d%/r, but not greatly smaller, and that the inertially-derived pres-
sure greatly dominates both of these. This would seem to imply that one should use
the inertial pressure for normalization. If this is done, the multipliers for the pres-
sure terms in Equations 2.2 become unity and the pressure terms are found to be of
the same order as the inertial terms. Because of this, I retain the pressure terms in
my linear analysis to better account for the case in which inertial forces drive the
pressure.
Placing the typical values of Table 2.1 into the Equations 2.2 reveals the follow-
ing for the region y > Af>:
1) convective terms are negligible,
2) viscous terms are negligible, and
3) body force terms are negligible
in comparison with the inertial and pressure terms. This result tends to reinforce
the central supposition of both Cremer and Lamb regarding the importance of the
inertial and pressure terms for their respective problems, but tends to open to ques-
tioning their retention of the viscous terms.
10


2.1.2 Non-dimensionalization using viscous wavelength
In the case As is used as the .r and y lengthscale, the Equations 2.2 become:
du
Ot
+ A .U
+ c-A,l
du
dy
(hP
dt-
+ £\,u ~
dr'
d.r
+ S\sV
dr
dy
VT dP T2 1
Asd0 dx + J d0 + Re.\a
jn-_dP^ r2 1
Asd0 dy" + J (l0 Re,\a
d2u d2u\
dx2 dy2 J
d2v d2t?\
d^ + d^2)1
(2.3)
where:
Re.\s
c As
VT
d0
K'
The chief difference between Equations 2.2 and 2.3 is embodied in the differ-
ence in magnitude between Xb and Xs: since both Re,\s and cAa approach unity, the
viscous and convective terms can no longer be neglected in comparison to the in-
ertia and pressure terms. This shows that for the region y < Xs only the body force
terms are negligible. This change in the relative emphasis of the terms in the mo-
mentum equations implies that a linearized solution which may model the flow at
moderate distances from the boundary may not model it as well in the region very
close to the boundary.
2.2 Linearized Theory Results
The scaling analysis showed that for y > h a linearized set of equations cap-
tures the important terms of the momentum equations. Applying this, and making
use of the work of Lamb, I solved a set of linearized equations for my particular
boundary conditions. Once I found and verified a solution, I made several plots
to help elucidate the chief features. In developing this solution I made extensive
use of the symbolic mathematics program Mathematica [18] to find and check so-
lutions, manipulate and simplify expressions, perform analysis, and make plots of
results. A Mathematica notebook which contains this development is included with
this thesis as Appendix A.
11


2.2.1 Solution to linear momentum equations
The linearized equations determined from the scaling analysis are:
du
Â¥
dv
dt
1 dP
p dx
i dP
- + |/V2v.
P dy
The pressure terms were estimated to be significant both near and far from the wall,
while the viscous terms were estimated to be significant only near the wall. I re-
tain both terms and thereby consider both near- and far-wall regions together. In
this way I avoid the necessity of developing a matching solution criterion at the
boundary between two dissimilar governing equation regions. Following Lambs
technique, I make use of a combined stream (0) and potential (0) function repre-
sentation of the flowfield:
do <90
<90 <90
dy ^ Ox
dp
dt
The use of these functions casts the momentum equations into the auxiliary form:
V20 = 0, (2.4)
w = "V*- (2-5)
I assume a solution periodic in x with a prescribed wavelength At = 2k/kb and a
time factor eiw< which takes the form:
0 = (Aekby + Be-kby)e(lkbX+wt\
0 = (Cemy + De-my)e(lkb£+ut).
with m2 = yjkf + ih'i (using the terminology k] = u/u). This last condition
comes from satisfying Equation 2.5. In order that the velocities remain finite as
y > oo, it is clear that A = C = 0. This gives for the velocities u, v and pressure
12


P/p:
P{x,y,t)/p = ibue-ykbet{kbX+ujt).
Of course it remains to take the real portion of these equations to find the physical
expressions for the velocities, but that must wait until after the constants b, d are
found.
Now I apply the boundary conditions
the first of which is an expression of the no-slip condition, and the second is a state-
ment of the fact that the wall is non-porous. This yields expressions for the con-
stants:
u(x,0,t) = 0,
v(x,0,t) = v0e'{kbT+ujt)
I'o l'o
kb h \jH + ik2s
d
iv o
Substitution of these into the expressions for u, v and Pjp yield:
(2.6)
(2.7)
13


^e~y^be(kbX + ajt)
(2.8)
These expressions are rather complicated and it is difficult to envision the nature of
the velocity field. To begin with I solve for the real portion of these expressions. As
with all of the above development, the details of this process are contained in Ap-
pendix A, which is the Mathematica [18] notebook I used as an aid. The result for
the u velocity is given in Figure 2.1 as an example (in which t = \S/Xb kb/ks).
As can be seen from this example however, the expressions are still complicated
and it is extremely difficult to clearly grasp any spatial form or pattern they may be
describing.
2.2.2 Fluid behaviour according to linear solution
In order to provide a clearer picture of the underlying nature of this solution, I
took advantage of the large difference in the magnitudes of the two wavenumbers
involved (k and kb) to considerably simplify the solution. As can be seen from Ta-
ble 2.1 the ratio e = Xs/Xb = kb/ks is on the order of 10-3, and Figure 2.1 shows
that this ratio appears in the velocity equations to the second as well as fourth pow-
ers. Therefore it may be instructive to examine the form of the solution as ks be-
comes infinitely large in comparison to kb, that is as e > 0. In considering this
asymptotic condition, it is useful and important to recognize that both ks and kb are
in fact independent: ks is determined by the frequency of the input vibration (as-
suming viscosity is given), while kb is determined by the bending wavelength. It is
true that for a given solid structure there may be well defined relationships between
vibration frequency and bending wavelength; for the fluid however there is no sim-
ilar a priori relationship. Thus there is a basis for using the asymptotic condition
to shed some light on the underlying fluid structure implied by the solution. Taking
the limit e > 0 in the expressions for u, v and P/p yields:
u(.i\yj)
(2.9)
e"<*<-+*') {kb ks)
v{x. y, t) =
(2.10)
u sin(A:(, .r) ks v0 cos(u>t)
(2.11)
14


u(x,y,t)

Jl+t^yco5(ar-Ila;ti21)
/ . , /arctan(e2)\ /arctan(e2)\ . A .
I cos(a' kb) sin I -------I + cos I ----------- 1 sin(a- kb) I kb
sykb
,os{Zki + (1 + £f sin k.) sin +
COS
arctan(e2)
2 / \ 2
2' ' 1 . (arctan(e2)
^ j sin{or kb + (\ + e4) 4 y sin
(1+£y Ls(^ tykb sinj.r hj + (l + rA'j 1 y sin A fcl
kS)
^ + sin(arCta2n(t2)yj
A'f,+

W
fe2 2 (1 + £"); cos falctan(<2)j h k. + v^TTC k2
Figure 2.1: Expression for real portion of u solution to linearized equation set


2
Position (x)
Figure 2.2: Pressure contours calculated using linearized theory.
These expressions are much simpler and easier to understand. Below I consider
each expression in turn, working from the simplest (the pressure) to the more in-
volved (the u velocity). In each of the flow field figures below I show the field at
the point in time corresponding to maximum velocity (or pressure) for clarity.
Pressure The expression for pressure/density represents a field varying sinu-
soidally in x (wavelength kb) and decaying exponentially in y (decay constant also
kb). A contour plot of pressure is presented in Figure 2.2.
Velocity v Like the pressure, this expression varies sinusoidally in x, and is in
quadrature with the pressure field (i.e. cosine x dependance versus pressures sine
dependence, see Figure 2.3). In y, the function is a maximum at the wall and decays
exponentially from there. The exponential with ks as the decay constant is signifi-
cant only very near the wall, and rapidly yields significance to the other exponen-
tial, which governs the subsequent overall behavior. This is shown in Figures 2.4
and 2.5. Figure 2.4 shows the region very close to the wall, while Figure 2.5 shows
a region extending farther from the wall.
Velocity u This expression varies sinusoidally in x and is in phase with the pres-
sure field. In y, unlike the v velocity expression which is at its maximum at the wall,
u is zero at the wall (in keeping with the no-slip boundary condition), but exponen-
tially builds (with spatial growth factor ks) as distance from the wall is increased.
16


Figure 2.3: Wall-normal (v) velocity contours calculated using linearized theory.
Figure 2.4: Wall-normal (v) velocity profile calculated using linear theory (near
wall region).
17


10/ks 1/kb 2/kb 3/kb
distance from wall
Figure 2.5: Wall-normal (v) velocity profile calculated using linear theory (overall).
At a distance not far from the wall it reaches a maximum value and then begins
an exponential decay (with decay constant kb). The distance at which it reaches its
maximum is easily determined:
_ i = / h- ks r1
Vc kc l In it*-in it J
I characterize this distance as the reciprocal of a wavenumber (kc) which repre-
sents the wavenumber of the crossover point between viscous and bending wave-
length dominated effects. The maximum u velocity is subsequently found as:
h ^ s
ks t'o
At first this may appear to suggest that under the asymptotic approximation
which 1 am using that the peak n is zero. Taking the limit kb/ks > 0 as before,
it can be seen that in fact in the limit the maximum approaches v0. It is of note that
this is the same peak velocity in u predicted by the inviscid solution of Cremer.
The u velocity behaviour is shown in Figures 2.6 and 2.7, the first of which
shows the region very close to the wall, and the second shows a region extending
farther from the wall.
Another view of the u velocity is shown in the contour plots shown in Fig-
ures 2.8 and 2.9. These highlight the periodicity of the field in x and the decay in
y-
18


Figure 2.6: Spanwise (u) velocity profile calculated using linear theory (near wall
region).
Figure 2.7: Spanwise (u) velocity profile calculated using linear theory (overall).
19


Figure 2.8: Spanwise (u) velocity contours calculated using linear theory (near wall
region).
o
Position (x)

Figure 2.9: Spanwise (u) velocity contours calculated using linear theory (overall).
20


_2_
kb
3
S 2 kb
flj
$
e
o
'o
i 4 / r
i 4 4 +
i 4 < r
\ 1 f ?/>-*>* v
^ ^ > V w
^ ^ 1 f ^ >
J_
i*
<0 ------
w 2 kt,
b
* / / ^ k
J t / ^ ^ ^
t t / s .*
\ t //^
"V'...-.
_-** *
^ ... v. y *
^ v. \ \
\
-w, //-^x \
' ' 1 ! f / "
^ \ ^ , f f *
^ \ \ , ,' / ^
f I
*
I f
1
4A
2 Ab
Position (x)
Figure 2.10: Velocity vector field calculated using linear theory (overall).
Velocity vector field Finally, Figures 2.10 and 2.11 show the velocity vector field
as predicted by linear theory. Note how the 90 spatial phase relationship between
the u and v velocity components results in regions (near x = 0 and x = Ab/2)
with essentially only wall-normal velocity, and other regions (near x = Ab/4 and
x = 3Ab/4) with essentially only wall-parallel velocity. At the wall of course the
wall-parallel velocity component is zero due to the no-slip condition. The build-up
in the u component as one moves away from the wall is apparent in Figure 2.11.
2.3 Consideration of Non-linear Effects
Schlichting [15] presents a solution for a periodic flow which considers non-
linear effects using a technique referred to as the method of successive approxima-
tions. In this approach, the velocities are assumed to be the sum of a zeroeth order
or linear term and an additional term which represents a perturbation to the linear
solution:
u(x,y,t) = u0(x,y,t) + ufx,y,t) (2.12)
In this approach, the linear term is the solution to the linear set of equations. This
term is replaced into the complete momentum equations and the resulting perturbed
solution becomes u\.
As pointed out by Schlichting, the first-order perturbation term will involve
terms such as e2'"'; the frequency doubling arising from the nonlinear terms.
It is possible to formulate using this method a set of equations which can be
21


ks
= 1
i kT
E
i A
I
1
! / /^_
t / /^_
f / /^_
: / /^_
* /
t /^_
! t t ,
"A ) ,
^\w/^
) / /^
''A; / /^
'\\\i//^
1 ;
,^\N \ t
^\\ \ t
,^\\ \ t
\ \ *
,.v\ \ t
_^\\ \ 1
* t I
! i
1
4Ab
1
2Xb
Position (x)
3
4Ab
Figure 2.11: Velocity vector field calculated using linear theory (near wall region).
solved to yield the first perturbation terms u\, vi. However, these equations would
still be required to meet the boundary conditions of the linear problem (since bound-
ary conditions do not distinguish between the method of solution, driven as they are
by physical constraints). The boundary conditions, by definition, will not contain
any frequency doubled terms. This presents a problem, namely that without the re-
quired source terms how can the perturbation solution capture this behaviour?
This remains an area for further investigation.
22


Chapter 3
Numerical Solution
In order to more fully understand the flowfield created by the studied motions,
and to provide some perspective on the utility of the linearized equations, a solu-
tion to the full Navier-Stokes equations was found using a numerical technique as
implemented in the commercial solver FLUENT. FLUENT employs the SIMPLE al-
gorithm (with an overview of this algorithm given in Section 3.1), and contains a
variety of features which allows it to be used in a wide variety of fluids problems.
Of particular importance for the present problem was the moving mesh feature (dis-
cussed more fully in Section 3.2 along with other solution and implementation de-
tails). The results of the FLUENT analysis are presented in Section 3.3.
3.1 SIMPLE Algorithm
fluent employs a semi-implicit method for pressure-linked equations (SIM-
PLE), the general form of which was pioneered by Patankar and Spalding [12, 13].
The governing Navier-Stokes equations are reduced to algebraic form using fi-
nite analytic method spatial and first-order forward finite-difference temporal dis-
cretization. Spatial discretization is by finite volume integration over the compu-
tational cells into which the domain is divided [14]. The resulting equations are
solved using alternating direction application of the Line-Gauss-Seidel (LGS) tech-
nique with underrelaxation. The semi-implicit solution uses the momentum equa-
tions and previous pressures to solve for an estimate of current velocities, and then
applies a pressure-correction derived from continuity and linearized momentum
equations to obtain current velocities and pressure. FLUENT monitors the imbal-
ance in a generalized conservation equation as a residual. Iteration proceeds until
the value of this residual decreases below a (user) set value.
3.2 FLUENT Implementation of Bending Wave Problem
The major challenge presented in modelling the bending wave problem was
creating a mesh which adequately represented the range of lengthscales involved.
23


Recall that the ratio of the bending wavelength to the viscous wavelength was on
the order of 103. To first order, for a two-dimensional Cartesian space, this im-
plied that in order to resolve viscous forces as well as correctly capture overall mo-
tion some 106 cells would be required. Investigations with the software installed
on the Marangoni node at CU Denver revealed that memory restrictions prevented
calculations with more than about 111.000 cells. Creative mesh generation using
the non-structured form of the code was considered but ruled out. Ultimately a
structured grid was defined which allowed for acceptable resolution near the solid
boundary where viscous effects would be most important, maintained a reasonable
cell aspect ratio, and limited the number of cells to a workable level. This definition
is briefly covered in Section 3.2.1. Of great assistance in enabling these calculations
was the use of symmetry in the mesh, as discussed in Section 3.2.2. Finally, the fi-
delity of the numerical solution was enhanced by taking advantage of the moving
mesh feature of FLUENT. This allowed the code to realistically mimic the defor-
mation of the solid boundary. In order to implement this feature, a user subroutine
was written and a special executable version of FLUENT was compiled and linked
with this subroutine. The use of the moving mesh feature and the development of
the user subroutine is detailed in Section 3.2.3.
3.2.1 Selection of mesh parameters
The built-in structured grid generation programs GeoMesh and pmesh were
used to produce a grid for use with FLUENT. GeoMesh takes dimensional input
and creates patches which are then meshed using pmesh. A code named Leo was
used to monitor cell aspect ratio and skewness of trial meshes.
The basic challenge was to create very y-compact cells near the solid boundary
to provide adequate calculational resolution in the viscous region, and allow for in-
crease in the y-thickness as the distance from the wall grew, while maintaining an
overall workable number of cells with a reasonable aspect ratio. The scheme I se-
lected is presented in Figure 3.1. I elected to use two patches. The patch near the
wall was meshed using a geometrically growing cell Ay, while the other patch em-
ployed a fixed Ay. A.r was fixed for both patches. The number of cells in the x di-
rection dictated the Ay at the wall because I limited cell aspect ratio to 1:5 (Ay: A.r).
Similarly, the number of cells in x dictated the fixed Ay in the patch further from
the wall. The choice of maximum cell aspect ratio was based on recommendations
in the FLUENT manual. In order that there be no dramatic change in cell size at the
boundary of the two patches, I required that Ay at the bottom of the top patch be
roughly equivalent to Ay at the top of the wall-adjacent patch. Thus the Ay growth
factor for the bottom patch was essentially dictated by the choice of A.r, as sum-
24


a y
15 mm
, i 1
+ 1 |
i !

2 1
pater |
T 1 | :
1
j i j I !
1 I i
! ! i
1 ' TT i
pateni
Lfc

computational domain
Ay=5a
geometrically increasing
Ay from wall in patch 1,
constant Ay in patch 2
Ay=a/5 at wall
/ x
(0,0)
Ax=a
15 mm
Figure 3.1: Meshing scheme employed in FLUENT.
Table 3.1: Levels of mesh resolution investigated.
Level Aywaii(^ni) A. total cells
a 60.549 4.95 1240
b 27.051 11.09 4514
c 12.804 23.43 17182
d 6.231 48.15 66998
marized in Figure 3.1.
Various degrees of mesh spatial resolution were investigated in order to study
the effect of cell size. A summary of the resolutions investigated is presented in
Table 3.1. Consideration of various mesh resolutions was facilitated by a built in
set of commands in FLUENT which allow the user to increase or decrease cell size
by a factor of two, in x or y or both x and y. In my case, I needed to simultaneously
change x and y in order to maintain proper cell aspect ratios.
A ratio of ks/kb = 100 was investigated with FLUENT to be representative
of experimental investigations. The radian frequency of input motion was chosen
(again based on experiment) as u; = 5817.771/sec (~ 926 Hz) and a bending wave-
length of 0.015m (kb = 209.441/m) was selected. These values are not only repre-
sentative of experimental values, but they also tend to be consistent with the asymp-
totic approximation used to derive simplified linearized equations.
The gas law option of FLUENT was used to represent air as the fluid. The default
25


Table 3.2: Gas law properties of fluent.
Property value units
working pressure molecular weight viscosity density kinematic viscosity 1.0133 x 105 28.970 1.72 x 105 1.293 13.302 x 10~6 Pa moles kg/m s kg/m3 m2/s
Figure 3.2: Boundary conditions used with FLUENT mesh,
properties for this option were used and are presented in Table 3.2.
3.2.2 Use of symmetry; mesh boundary conditions
The extent of the calculational zone was positioned to take advantage of sym-
metry as shown in Figure 3.2. It was predicted by Lamb and Cremer (as well as the
linearized solution I computed) that the u velocity would be zero at bending wave
motion antinodes, thus a symmetry boundary condition was enforced at the left and
right boundaries of the problem. An outflow boundary condition was enforced at
the top of the mesh. The bottom boundary is represented as the type z-wall which
is the designation for solid walls to be used with the deforming mesh model.
3.2.3 Moving mesh forcing function
The moving mesh feature of FLUENT proved to be quite useful for this prob-
lem. For one, it allowed greater fidelity to the physical problem since for real plate
vibration the fluid-solid interface is moving. In addition, the use of a user-defined
subroutine, required for moving mesh problems, allowed for fairly precise spatial
specification of input motion.
26


start cycle
solve momentum
equations
update velocity
finish cycle
calculate
geometrical
parameters
update mesh
solve mass balance
(pressure correction)
equation
update velocity
and pressure
(adapted from
FLUENTmanual)
""f-
check convergence
Figure 3.3: Geometry updated after solution of momentum and pressure equations.
User-defined moving mesh subroutine The fluent code incorporates a mov-
ing mesh flag which, when set, causes the mesh position to be defined by the user-
defined subroutine NEWMESH. Updating the mesh position is accomplished dur-
ing the solution cycle as shown in the flow chart of Figure 3.3.
Mesh motion description A standing spatial, time harmonic deformation of the
z-wall boundary was used. I modified a subroutine NEWMESH to input this motion
using an equation of the form:
y{x,y,t) = dQs(y)m cos(kbx)sm(Lut), (3.1)
where s(y) is a y-scale factor which smoothly adjusts the deformation from its max-
imum at the lower boundary to zero at the upper boundary:
s(y) =
H-y
H
where H is the overall extent in y of the mesh. The parameter m controlled how
strongly the deformation was concentrated near the lower boundary. Various values
of m were tried and it was found that too large a value resulted in slower conver-
gence. A value which produced reasonable convergence rates over a broad range of
conditions was m = 0.5. Equation 3.1 was written in a discretized form for imple-
mentation into the NEWMESH subroutine. Note that the x positions for the mesh
points were not changed.
27


A printout of the NEWMESH subroutine is presented in Appendix B.
3.3 FLUENT Results
As the FLUENT solution was developed, it was important to verify that the mesh
resolution or other ancillary effects were not driving response. To investigate this
a number of analyses were conducted. Results from a series of grid refinements to
verify that cell size was sufficiently small are presented in Section 3.3.1. Results
from a corresponding time step study are presented in Section 3.3.2. Results from
a study looking at the duration of the startup transient are presented in Section 3.3.3.
In all of these studies, comparison with the linear theory prediction is made.
3.3.1 Grid refinement
Time histories at several locations were made to compare the effect of the var-
ious grid resolutions discussed in Section 3.2.1. Figures 3.4 through 3.6 present
these results. In each of these figures, a small schematic shows the point from which
the presented data are taken. Also shown for comparison purposes is the corre-
sponding data from the linear theory solution.
It is clear from these figures that there is generally good agreement between the
linear theory and FLUENT, even for the coarsest meshes. Figure 3.4 presents the u
velocity as a function of time for the point x = Ab/2,y = 7.2AS. Spatially, this
point is placed in a region where u is relatively large. Figure 3.5 presents the v
velocity for the point x = Ab/4, y = 7.2AS, at which v is relatively significant. In
both these figures the numerical solution can be seen to closely track the linear one.
Refining the mesh clearly helps, but does not seem to dramatically affect results.
A slightly different result emerges when considering regions in which the linear
theory predicts insignificant velocities. For example, Figure 3.6 presents v velocity
results for the point a- = Ab/2, y = 7.2As, which is the same point the u velocity was
presented for in Figure 3.4. At this point in the flowfield the linear theory predicts
v = 0. As can be seen, the FLUENT solution shows a very small v velocity, but not
identically zero. Furthermore, the degree of mesh refinement clearly has significant
impact. The coarsest mesh considered (mesh a) appears to acquire a steady-state
positive drift as the solution proceeds. This drift becomes less apparent as the mesh
becomes more refined. The finest mesh (mesh d) exhibits small oscillations of the
v velocity at this point, but no drift. Because of this, the finest mesh was chosen for
the remainder of the analyses.
28


0 0.0005 0.001 0.0015 0.002
Time (seconds)
Figure 3.4: FLUENT solution time histories for various degrees of mesh resolution
(u velocity component, x = Ab/2).
0 0.0005 0.001 0.0015 0.002
Time (seconds)
Figure 3.5: FLUENT solution time histories for various degrees of mesh resolution
(v velocity component, x = Ab/4).
29


o
>
>
o
o
0.04
0.03
0.02
0.01
0
£ -0.01
-0.02
-0.03 f
^x
/ \
/ V
XX.
\
-.'-55-
' /
/
/
1/y
A
-''X
/ \
/ \
/ X
1
! ^
If ,'-x X
* ///;- \ -x^
VX
\
_\
(0.5^,7.2^)
linear theory j
mesh d (402x139)
-- mesh c (242x71) j
--mesh b (122x37) !
- mesh a (62x20) )
0.0005 0.001 0.0015
Time (seconds)
0.002
Figure 3.6: FLUENT solution time histories for various degrees of mesh resolution
(v velocity component, x = Ab/2).
3.3.2 Time step refinement
Various levels of time resolution were examined as well. The degree of time
resolution was quantified in terms of the number of time steps per cycle of bend-
ing wave motion. Five levels were investigated: 6, 12, 18, 24, and 30 time steps
per cycle. Velocity histories for these cases are plotted versus the linear solution
in Figures 3.7 and 3.8. Both figures reveal that the solution results are essentially
independent of time step size, at least for the step sizes considered here. Note how
all solutions pass through essentially the same points. The only difference is the
resolution in time with which the resulting motion can be characterized. In order
to obtain reasonable time resolution with reasonable computational effort, 24 time
steps per cycle was chosen for the remainder of the analyses (in addition, the mesh
resolution studies were made using 24 time steps per cycle).
3.3.3 Starting transient
In the FLUENT analysis all of the fluid begins at rest. Therefore it is necessary,
prior to comparing results to the linear theory, to allow the FLUENT solution to come
to steady state motion. Figures 3.9 through 3.11 present a comparison of the initial
portion of a FLUENT time history with a portion of the same history after one cycle
of motion. In general, after the first half cycle of motion, the two portions are nearly
identical, implying that the starting transient dies out after about 180 of motion.
Figure 3.9 presents the difference in the u velocity between the first and second cy-
30


0 0.0005 0.001 0.0015 0.002
Time (seconds)
Figure 3.7: FLUENT solution time histories for various degrees of time step refine-
ment.
0.425 0.450 0.475 0.500 0.525 0.550 0.575
Time (milliseconds)
Figure 3.8: FLUENT solution time histories for various degrees of time step refine-
ment (detail).
31


o
"D
>%
O
o
(D
>
0.0015
0.00125
0.001
0.00075
0.0005
0.00025
0
-0.00025
(0.5Xb,7.2Xs)
mesh d (482x139)
mesh c (242x71)
--- mesh b (122x37)
mesh a (62x20)
-V
J '

;i
/:
/;
; /
0.0002 0.0004 0.0006 0.0008 0.001
Time (seconds)
Figure 3.9: Difference in first and second cycle u velocities (y = 7.2AS).
cles of motion for the point x = Ab/2, y = 7.2AS. Only a very slight difference is
apparent. This behaviour is representative of the fluent solution in general, with
the exception of a region very near the wall. An example of this region is shown in
Figure 3.10. In this figure the velocities are plotted (as opposed to the differences
in the velocities). As can be seen there is a significant difference initially, which
disappears by the end of the first cycle. This is shown perhaps more clearly in Fig-
ure 3.11, which is the velocity difference for the data of the previous figure. In the
case of the very near wall velocities it is clear that a full cycle of motion is required
to eliminate the starting transient, and for this reason all subsequent results present
data from the second (or sometimes later) cycles of motion.
3.4 Comparison of Linear Theory and FLUENT Results
In addition to the comparisons presented in the section previous, fluent re-
sults are compared to those of the linear theory making use of velocity profiles (Sec-
tion 3.4.1) and vector field plots (Section 3.4.2).
3.4.1 Velocity profile results
Figures 3.12 and 3.13 show u velocity profiles in y at three x locations along
the bending wave for both the FLUENT calculation as well as the linear theory. Fig-
ure 3.12 captures the velocity profile at a point in time half of the way through a
cycle of wall motion (180), while Figure 3.13 presents results from 315. Both
32


0 0.0002 0.0004 0.0006 0.0008
Time (seconds)
0.001
Figure 3.10: Comparison between first and second cycle u velocities (y = 0.2AS).
or
3
I
d)
O
c
a>
a>
*0
o
o
a>
>
-0.05 |
-0.1
-0.15
-0.2
-0.25
-0.3
-0.35
Jr
in.
:n
a
i
i
> (0.5Xt),0.2X,)
mesh d (482x139)
mesh c (242x71)
- mesh b (122x37)
-- mesh a (62x20)
0.0002 0.0004 0.0006 0.0008 0.001
Time (seconds)
Figure 3.11: Difference in first and second cycle u velocities (y = 0.2AS).
33


0 0.002 0.004 0.006 0.008 0.01 0.012 0.014
Position (meters)
Figure 3.12: Comparison of linear theory solution and FLUENT result, u velocity
component at 180.
show generally good agreement with linear theory, the only notable difference be-
ing in the near-wall region of the profile in Figure 3.13. This will be explored in
more detail subsequently.
Figures 3.14 and 3.15 present similar data for the v velocity component. These
indicate as well that the FLUENT solution is in generally good agreement with linear
theory, again with the possible exception of certain near-wall discrepancies.
As mentioned above, a difference emerges between the fluent calculation and
the linear theory for comparisons at or near the points in time for which the linear
theory predicts no motion in the fluid (such as 270), for regions very close to the
Position (meters)
Figure 3.13: Comparison of linear theory solution and FLUENT result, u velocity
component at 315.
34


Position (meters)
Figure 3.14: Comparison of linear theory solution and FLUENT result, v velocity
component at 180.
Position (meters)
Figure 3.15: Comparison of linear theory solution and FLUENT result, v velocity
component at 315.
35


2.0
*

240 255" 270' 285' 300'
Velocity profile (u/uma,) at cycle phase
Figure 3.16: Comparison of linear theory solution and fluent result (dashed), u
velocity component near 270, x = Ab/2.
Time (seconds)
Figure 3.17: Comparison of linear theory solution and FLUENT result, near-wall
(y = 0.2Ab) u velocity component versus time.
wall. Figure 3.16 presents a series of u velocity profiles (spanning the points in
cycle motion from 240 to 300 in 15 increments) at x = Ab/2. As the linear
theory approaches the quiescent portion of the cycle, the FLUENT calculation shows
a near-wall bulge in the velocity profile in the direction opposite to the overall
motion. After passing through the quiescent portion, this bulge continues as a near
wall protrusion in the velocity profile. Looking at time history data (Figure 3.17),
this bulge can be characterized as a phase-leading zero crossing of the u velocity
profile in the very near wall region of the flow.
Figure 3.18 presents a similar series of u velocity profiles for the location x =
Ab/4 while Figure 3.19 presents representative profiles for the v velocity compo-
nent.
36


2.0
1.5 -
I 1.0 -
Q- 0.5
>.
As

240 255 270 285' 300'
Velocity profile (u/um) at cycle phase
Figure 3.18: Comparison of linear theory solution and FLUENT result, u velocity
component near 270, x = Ab/4.
Velocity profile (v/vma,) at cycle phase
Figure 3.19: Comparison of linear theory solution and FLUENT result, v velocity
component near 270, x = Ab/2.
37


fluid at wall moving
with peak velocity
as boundary passes
through neutral
position
wall motion slows,
moving fluid is decelerated,
creating significant near-wall
velocity gradient in wall-normal
direction (selected region
shown for clarity)
large gradient in y gives
rise to large gradient in
x (by continuity); resulting
convective coupling
creates near-wall flow
in direction opposite to
prevailing u velocity component
Figure 3.20: Schematic illustration of convection-driven near-wall motion.
This near-wall difference between the linear theory and the FLUENT results,
coupled with the findings of the scaling analysis (namely that the convective terms
are not negligible for the region V < As) suggests that this near-wall behaviour
arises from convective forces. In fact it is reasonable to expect that as the mov-
ing fluid is decelerated by the wall (in the wall normal direction), there would be
coupling by convection of fluid motion in the wall-parallel direction. This is shown
schematically in Figure 3.20, in which the presumption that the convective terms
are significant follows from the results of the scaling analysis. As the moving fluid
encounters the slowing motion of the wall (referring to the figure) a strong wall-
normal velocity gradient is created. This dv/dy gradient couples to du/dx by con-
tinuity, creating significant convective forces in both the x and y momentum equa-
tions.
3.4.2 Vector field results
Figure 3.21 presents the linear theory result for the velocity vector field through
one cycle of wall motion. This summarizes the motion of which various aspects
have previously been presented in some detail, with the position of the wall bound-
ary superimposed. It also shows clearly the portions during the wall motion cycle
for which the fluid is at rest (these portions correspond to the points of maximum
38


Figure 3.21: Linear theory solution velocity field at several points through wall mo-
tion cycle.
deflection of the wall boundary, when it too has come to rest and is about to reverse
its direction of motion).
As noted above, the linear theory and FLUENT calculations agree very well for
regions y > As. Focusing on the region of the computational domain near the wall
(identfied in Figure 3.22), a series of FLUENT results which highlight the convective
motion near the quiescent phase angle of 90 are presented in Figures 3.23 and 3.24.
These plots show results for the near-wall region at x = Ab/2 at various points
through the wall motion cycle. They clearly show the nature of the near-wall phase
leading motion in u which is not predicted by the linear theory.
Figure 3.22: FLUENT results region showing near-wall effect.
39


*
0
15
Ul'tUUli
H.l bST Clod hie f COVHGl RATIOS
FLl'ENT Ond File f CONFKJl RATIOS a
Veloaty Vecure (M/S)
Ijw 9 419E-01 ljae.4 M2F.-0 Time
45

R.l ENT Ond Kite f COS FIGURATION > me*03
Velocity Vectcri (MA)
Ijnex 6.971B-01 Ijwa 1-20IE-02 Time.2.2%
. CkiOS l*M
| FlMitt
: FtMlK.
]4S
IW
) DM!
u
ia
75
^ / j 7 n n n j
U-iU M t i i I i
FLUENTGnd File /t'ONFIGl RATIOS .
| Vcluoty Vedcei (MA)
| lmi 2^29fc-0l U .1J90B-02 Tit
metbOI */
e 2.JME-0J
Oct 0$ 1991
Fheeat 4 41
FWe* lac.
Figure 3.23: FLUENT calculation near-wall velocity field at several points through
wall motion cycle (0-75).
40


w*
t-
i. V' ^ Ij £ f (- £ f j- p f
FULST C*W hie f a>SFKU RATIOS me*0.< /
Vclocay Vector* (M/Si
|jMi>2.2J5fc-OI lai I CiK-03 Time 2-4XX-.-03
{ Oct 05 IWI
l Flwa 4 41
Fteeatlac
T05^

KU'LNT(Snd hie / COSFKil'RATIOS me*h03 > Octal mi
Veloaty VMutiM/S) ' Hmi4-U
Lmi ) 973fc-0l Um 4 S52E-02 time 147M -03 Fte* Ik.

T35CT
P r f t f?
f Pf Ff f
R.rtNT God file /COSHGI'RATIOS > ee*03 /
Velucay Vector* (M/Si
Iah 5 745F.-01 liiuiSJUh 02 Time = 2.5201.03
! Oa05 IW* |
Hue* 4 4* j
I hW Ik i
R J fcNTOnd File / CONFK31 RATIOS me*03 ,
Velocity Vech*t IMA/
Inu 7.37OF-01 l.nui
S.V42F-02 Time
LWI-03
ua as i94
I1W4.II
Fluea lac.
150
FUHSTGod Hie f (T)SHCI RATIOS *
Veloaty Vvn
Lm** 16391.-01 Ijmi a ) 902B-Q2 Tim
T'
Oct 03 IW*
Hues 4.44
FUteat Ik.
| FU F.STGnd HU / COSF1GIRATIOS raefaOA /
' Velocity VeOtxnWSr
Itau *9 VMI -OI lmi = t_Mlh "2 Titm J.655F-03
FktCM 4.41
Fkcm Ik.
Figure 3.24: fluent calculation near-wall velocity field at several points through
wall motion cycle (90-165).
41


Chapter 4
Conclusions
Fluid motions resulting from a bending wave motion of a solid wall were inves-
tigated using both theoretical and analytical techniques. An analysis of the motion
predicted by linearized Navier-Stokes equations was conducted, using a solution
technique inspired by Lamb. Theoretical results were based on a scaling analysis
which showed that convective terms are negligible provided the region considered
is sufficiently removed from the near-wall region. Theory showed that the bending
wave motion created a fluid velocity field which exponentially decayed with dis-
tance from the wall with a decay constant equal to the bending wave wavenumber
(A'b). Because of the no-slip condition at the wall, the wall-parallel velocity compo-
nent at the wall is zero. Theory showed however that this component exponentially
increases with distance from the wall with constant ks = up to a peak value
which is governed by both kh and ks, as well as the peak wall-normal boundary
motion velocity.
An analytical study was performed using the computational fluids dynamics
code FLUENT to compare and contrast with the linearized theory. Results from this
study showed very good agreement with the theory overall with one exception -
analytical results indicate a phase-leading near-wall motion of the fluid in the time-
frame adjacent to times at which the theory predicts no fluid motion at all. This
discrepancy is attributed to convection, consistent with the findings of the scaling
analysis. At all other times and in regions greater than y ~ 2ir/ks, there was re-
markable agreement.
4.1 Directions for further research
There is some potential that the near-wall fluid motion may be fruitfully ex-
plored using either a perturbation solution or method of successive approximations
to develop a better understanding of the convective fluid motion manifested in this
region.
42


Appendix A: Induced Fluid Motions
of Beam Bending Waves Solution to Linearized Navier-Stokes using Stream
and Potential Functions
Mathematics notebook developed by Anthony W. Raskob 3-98 through 9-98 in support
of Master's Thesis.
copyright 1998 A.W. Raskob, Jr.
This notebook describes the solution of a linear form of the Navier-Stokes equations to model the
near- and far-field fluid motions induced by a bending wave at one boundary. Effects of fluid
viscosity are included. The approach follows Lamb's (Sir Horace, Art. 349 of Hydrodynamics)
solution to a water wave problem.
Linearized Navier-Stokes Equations
Preliminaries
First some preliminaries. We will use two-dimensional cartesian space, with x and y as the
coordinate axes. V2 is the Laplacian operator:
Ubprotect [Power] ; v2[q_] := d(x,2) q+a(y,2) q; Protect [Power] ;
Well also need the Plot Field package which will allow us to make vector field plots.
<< Graphics'PlotField'
The nature of the Mathematica notebook presents each section as a group when printed. Because
of this grouping some page space seems empty. This is merely due to this unfortunate aspect of
Mathematica when printing.
In contrast, the symbolic manipulation capability of Mathematica contributed greatly to the
reduction of algebraic errors, and the ability to consider a wide range of different solution
techniques to the problem considered.
Throughout this notebook I make use of several of the standard Mathematica packages. These are
usually available with the standard installation of the program. Complete descriptions of these
43


packages and their capabilities are contained in the Mathematica book by Steven Wolfram as well
as in the online help libraries which come with the code. In the next portion of this section I will
give an overview of the form of Navier-Stokes equations which I used to obtain my linearized
solution.
Form of Navier-Stokes
I am using a linearized form of the Navier Stokes equations, the particular form of which was
arrived at using dimensional analysis. These equations are similar to those used by Lamb, the
major difference being that I am neglecting body forces in y, which turns out to be appropriate for
my conditions. Here are the momentum equations:
du 1 dp 2 dv
57" = -----5 + V V U and ~r~
dt p ox at
P
dp
dy
- ---r + v vz v.
Lamb formulates a solution to this set of equations (along with the continuity equation) in terms of
stream and potential functions:
u[x_, y_, t_] = -dx$[x, y, t] -dyi)r[x, y, t] ;
v[x_, y_, t_] = -dy[x, y, t] +axi/r[x, y, t];
P[x_, Y_ t_] = 3t [x, y, t];
which will meet the N.S. and momentum equations provided V2

suggests the following forms for the potential and stream functions, respectively:
y, t_] := b Exp 11 kx] Exp[-ky] Exp [I w t] ;
i/r[x_, y, t_] := dExp[I kx] Exp[-ny] Exp [I a> t] ;
The form of these equations is dictated by the periodic form of the imposed input in both space
and time (accounting for the Exp[I k x] and Exp [I oj t] terms), and the decay in y with increasing
distance from the moving boundary suggested by behaviour observed in similar periodic
excitation cases (see for example a discussion of Stokes 2nd problem in Lamb).
Now I check to see that these expressions meet the auxilliary conditions:
Sinplify[v2 [<£ [x, y, t] ] ]
0
Siiqplify[v v2 [<)r[x, y, t] ] -dt *[x, y, t] ]
dEI (kx.lny.tu) (-^2 v + n2 v J u)
44


As can be see, the expression for

to solve for n such that the auxilliary condition is met:
Solve[% == 0, n]
r r v/k2 v + I u , r Vk2 v + I ui ,
{{.-----^}. 1--^1)
Using the latter for n, the potential and stream functions become:
4>[x_, y t_] := bExp[Ikx] Exp[-ky] Exp[Iu t] ;
*[x_, y_, t_] := dExp[I kx] Exp[-VkA2 tlu/v y] Exp[It] ;
Checking the auxilliary conditions, we see that these are indeed solutions
Sinplifyfv2 [$ [x, y, t] ] ]
0
Sim?lify[v v2[*[x, y, t] ] at ^[x, y, t] ]
0
These then are the expressions for u,v, and p/p:
u[x, y, t]
-IbEIltx'ky*It"k + dEIkx,It'Wlt2*if /k2 +
v[x, y, t]
bEIkx-ky.Itu k+ j d glkx.I tu-y ^
P[x# Y. t]
i bEIkxkyu ui
Do these expressions conserve momentum?
Substituting these expressions in the N.S. equations:
Sim>lify[dtu[x, y, t] + a*p[x, y, t] -w2[u[x, y, t] ] ]
0
45


Simplify[dtv[x, y, t] + dyp[x, y, t] -w2[v[x, y, t] ] ]
0
This shows that they do.
How about continuity?
Sinplify[dx u[x, y, t] + ay v[x, y, t] ]
0
These expressions meet the governing equations for the problem as I have proposed it. Now I will
move on to finding the value of the constants b,d.
Application of Boundary Conditions
The boundary conditions will be used to determine the constants b,d. At y=0, u=0 due to the
no-slip condition, and also at y=0 v=v0 Exp[/ k jc] Exp[/ a>t].l will make use of the symbol £v2to
denote what I call the Stokes wavenumber, after Stokes 2nd problem. k/=y.
constants=
Substituting the above into the expressions for u, v, and pIp gives the following (which for ease of
future manipulations I call usol, vsol, and psol):
usol[x_, y t_] =
u[x, y, t] /. { constants! [1, 1] ], constants! [1# 2] ], u / v -> k.2}
Solve [{u[x, 0, t] == 0, v[x, 0, t] ==v0 Exp [I kx+ lot]}, {b, d}] /.
{id / v -> k.2}
,-y.y^.Ik! 2 + k2 Vo
I E1 Xx-ky*I t u k Xl
Vo
k Vk2 + I k2s
k Vk2 + I k|
vsol [x_f y_, t_] =
v[x, y, t] / { constants! [1? 1] ] / constants! [1, 2] ], m / v -> k.2}
46


psol[x_, y_, t_] =
p[x, y, t] / { constants! [1, 1] ], constants! [1# 2] ], u / v -> k,3}
k Vk2 + I k|
In order to visualize this velocity field, it is necessary to find the real parts of the above
expressions.
Results
Extraction of real portion of expressions for u, v, p/
Firstly, it is convenient to recognize that because I am considering a standing wave solution, the
common time factor Exp[I to t] in these expressions can be dropped, leaving expressions for the
spatial patterning of the velocity field, the magnitude of which would vary in time sinusoidally:
umag[x_, y_] = Simplify[usol[x, y, t] /Exp[Iwt]]
Vk2*Iki (Eky EyVk2-IkI ) Vk2 + I k2 v0
-k + Vk2 + I k§
vmag[x y_] = Simplify[vsol [x, y, t] / Exp [Iwt]]
EIkx |_E-y A2.ikl k + E-ky ^/kTTYkf | Vo
-k + Vk2 + I k2
I*nag[x_# y_] =psol[x, y, t] /Exp[It]
The next step involves replacing the complex radical Vk2 + I k| with an expression in which the
real an imaginary parts can be easily separated.
Cos [ y ArcTan [ ] ] (k4 + k^ )1/4 + I Sin [ y ArcTan [ ] ] (k4 + k^
1/4
47


I will make use of the following shorthand and set f=( j-)2such that Vk2 + I k2 would become
Cos [ -5- ArcTan [ f ] ] ks ( f2 + 1)1/4 + I Sin [ -5- ArcTan [ f ] ] ks ( f2 + 1)1/4
Load Relm package
The following loads a package which automates the solution of real and imaginary parts and will
facilitate considerably the next steps:
<< Algebra'Relm'
This next lines attach rules to the indicated variables which tells Mathematica that they are
real-valued:
lm[v0] A= 0; Im[k] A = 0; Im[k.] A= 0; Im[x] A= 0;
Im[y] A= 0; Im[f] A= 0; Unprotect [Power]; Im[ (f2 + 1)1/4] A= 0;
Protect[Power]; Im[catf] A= 0; Im[satf] A= 0;
Now I take each expression in turn, beginning with that for u:
To make the expressions easier to read, I will let calf = Cos [ ArcTan [ f ] ] and
satf=Sin [ y ArcTan [ f ] ]. umag represents the spatial magnitude function of the u velocity:
umag[x_, y_] =
umag[x, y] / {Vk2 + I kj -> k {f2 + 1)1/4 (catf + I satf )}
/j glkx-ky-(lf:l1/4 (catf + 1 satf) y ks

(catf.lsatf) yksj (1 + f2}^ (catf + J satf) ks V() ) j
-k+(l + f2) (catf+Isatf)ks
Now I manipulate the expression above, taking the complex portions out of the exponential using
Euler's identity:
48


umag[x_, y_] = umag[x, y] /.
{Exp[g_ (catf + I satf) ] -> Exp[gcatf] (Cos [g satf] + I Sin[g satf]),
Exp [I h_ + j_ + g_ (catf + I satf) ] ->
(Cos[h] +ISin[h]) Exp[j] Exp [g catf] (Cos[gsatf] + I Sin[g satf ])}
(i Eky-catf l.£2,1,4yks + f 2 ) 1/4
(catf + I satf) (Cos[kx] + I Sin[kx] ) (Eky Ecatf a-f2)1/4 y ks
(Cos [ (1 + f2) *'4 satf y ks ] + I Sin [ (1 + f2)1,4 satf y ks ] ) )
(Cos [ (1 + f2)1/4 satf yks] I Sin [ (1 + f2) 1 4 satf y ks ] ) kE v0) /
(-k + (1 + f2)1/4 (catf + Isatf) ks)
Next I multiply the numerator by the conjugate of the denominator:
unum[x_, y_] = Numerator[umag[x, y] /. {y-> -yy, -y->yy}]
Conjugate [Denominator [umag [x, y] /. (y -> -yy, -y -> yy} ] ] ;
unura[x_, y_] = unum[x, y] / {-yy -> y, yy -> -y};
unum[x_, y_] = Siinplify[unum[x, y] ]
I
g-y (k*cat£ (1 + £2 )
gcat£ (1 + f2 ) 1/4
(Cos[(1+f2)
(-k + (1+f2)
,1,4ks) (1 + f2 )1/4 (catf + Isatf) (Cos [k x] +ISin[kx]) (Eky-
yks (Cos [ (1 + f2)1/4 satf yks] + I Sin [ (1 + f2)1/4 satf yks] ) )
1/4 satf yks] I Sin [ (1 + f2) 1/4 satf yks]) ks
1/4 (catf I satf) ks) v0
Next I multiply the denominator by its conjugate. With this, all the complex portions of the
expression are linear factors in / = V -1
udan[x_, y_] =
Sinplify[Expand[DenQfninator[umag[x, y] /. (y-> -yy, -y -> yy) ]
Conjugate[Dencnlnator[uinag[x, y] /. {y-> -yy, -y -> yy} ] ] ] ] /.
{catfA2 + satf A2 -> 1, catf -> Cos[-^- ArcTan[f] ] }
k2 2 (1 + f2)1/4 k Cos [ ArcT^nt£] ] ks + Vl + f2 k2
Next I expand the numerator and set the complex portions to zero. Then I replace the satf and catf
variables with the actual functions. This yields an expression for the numerator of the real portion
of the u velocity solution:
ureal [x y_] = Expand [unum[x, y] ] ; ureal [x y_] =
ureal[x, y] /. {I-> ii, -I > -ii>; ureal[x_, y_] =ureal[x, y] /.
{ 11 -> 0, catf -> Cos [ ArcTan[f ] ], satf -> Sin[ ArcTan[f ] j }
49


Finally, I put the numerator over the denominator:
ureal [x_, y_] = ureal [x, y] /uden[x, y]
jEky-y(k*,i*f>,1'co.[i=iSntL.]lt.) (1 + f2ji/kCos[ ArcTanlf^j
Cos [ (1 + f2)1/4 y Sin[ ArcT^n[£] ] kg] sin[kx] ks v0 -
Ell.f;)1/ (1 + f2)1/4 k Cos [ ^rcT^n [ £ 1 ]
Cos [ (1 + f2)1/4 y Sin [ ArcT^n^^ ] ks] Sin [k x] ks v0 +
Eky-y (k*(l.£2)l/4 cos [] ks) (1 + f2 ) ^ k CoS [ ]c X ]
Cos[(l + f2)1/4ySin[-ArcT^[£1 ] ks] sin[^-T^L£] ] k= Vo _
Ea,£2)lMyCos[A£cT^lks.y(kM1,£2)1/lcos[Arc3|nUL]ks) (1 + f2) 1/4 kCOS [kx]
cos [ (1 + f2)1/4 y Sin [ ArcT;n[£] ] ks] Sin[ ^Z|El£L ] ks Vo -
gky-y (k+ (1+f2 )1/4 Cos( ArcTanIf1 1 ks) (1 + f2)1/4 kCos[kx]
COS[ArcTan[f] ] Sint (l.f>,'ySi.,[ArcT"[£1 ]*,]*,
gk y-y (k+(l + f2)1/4 Cos [ ArcTan f f1 ' 1 k=) (1 + f2)1/4 k Sin[k x]
_ r ArcTan[f ] Sln ( 2 ] Sin[ (1 + f2)14 y Sin [ ArcT^nt £1 ]ks]ki
g(l+£2)1/4 y Cos [ ] ks_y Cos [ A^-T|nLfl ] ks)
(1 + f2)1/4 k Cos [ ^£^Z!|£LL£L ]
Sin[k x] Sin[ (1 + f2)1/4 y Sin[ ArcT^n^£^ ] ks] ks v0 -
£(l + f2)1/4 y Cos [ krcTan£^J_ ] ks _y (k+ (1+ f 2 > 1/4 Cos [ ^££l|£lilL ] ks ) ( 1 + f 2 ) 2'4 k COS [ k X]
Sin[ArcT^[£] ] Sin[ (1 + f2)4 y Sin[ ] ks]2 ks Vo .
Eky-y (k.|l,£2)1/4Cos[ArcT^nlfl ] k.) C0= [ ArCT^n [ £ ] ]"
Cos [ (1 + f2)1/4 y Sin [ ArcT^n^£^ ] ks] Sin [k x] k2 v0 +
e(1*£2)1/4 y Cos [ ArcT^ntf J ] k3-y (k- (l,f2 )1/4 Cos [ krc^n[_f J_ ] kg j Cp= [ ArCT^1 [ £ ] ]2
Cos [ (1 + f2)1/4 y Sin[ ArcT^n[£] ] kj] sin[kx] k2 v0 -
Eky-y(kMl*f2,1/4Cos[^iH]k3) COS [
(1 + f2)1/4 y Sin[ ArcT^n[£ 1 ]ks]Sin[kx] sin[^C-T;n[£1 ] k2 v0 +
E Vl + f2 Cos [ (1 + f2)1/4 y Sin [ ArcT^n[f 1 ] k= ]
50


r, r ArcTan [ f ] 1 2 2
Sin [k x] Sin ----- 1
Eky_y
Cos
k* (1 *£2) Cos
ArcTan[f1
Eky-y (k. (l.t2)1/a Cos [ AreT;nlt) ] ks) + f2 COS[kx]
- J k2 V0 +
ks) a/1 + f2 Cos[kx]
Sin[(1 + f2 ) y Sin
ArcTan[f]
Sin
ArcTan[f]
Sin [ (1 + f2 ) y Sin
ArcTan[f]
ks k2 v0 +
ks k2 v0
E(l.f2 )1/J yeos[ ArcTfnIfl
ks-y (k.ll.f2)1'4 Cosl-K£i$aiiL] ks
ArcTanffl n
Sin[k x] Sin (1 + f2) y Sin
r.--TT ^ r ArcTan [ f ]
VI + f2 Cos [-----
2
:s1 kj v0 +
Ea.£2)l/SCos[kIcT|s^]ks.y(k.atE2)lMCos[kIHTai^lks) sin[kx]
r ArcTan[ f ] i 2 . r ,, r ArcTan[f] i i2 2
Sin -------. Sin (1 + f2) y Sin ---------, ks k2 v0
k2 2 (1 + f2
k Cos
ArcTan[f1
] ks + V1 + f2 kj
This is a very long expression. Rather than spend effort in simplifying it directly, I will make use
of the fact that in the case of interest ks k. As will be shown, this will make considerable
simplifications possible. But first, I go on to perform the same manipulations on the expressions
for v and p/p that I did above for u.
Taking real portion of v:
To make the expressions easier to read, I will let catf = Cos [ ArcTan [ f ] ] and
satf=Sin [ \ ArcTan [ f ] ]. vmag represents the spatial magnitude function of the v velocity:
vmag[x_, y_] =
vmag[x, y] /. {V*2 + I kj -> k. (f2 + 1)1/4 (catf + I satf ) }
EIkx (-E'(1,f211 4 (=atf.isatf> k + E'ky (1 + f2)1/4 (catf + I satf) ks) v0)
(-k + (1 + f2)1/4 (catf + I satf) ks)
Now 1 manipulate the expression above, taking the complex portions out of the exponential using
Euler's identity:
51


vmag[x_, y_] =vmag[x, y] /. {Eb?)[Ih_] ->Cos[h] +ISin[h],
Exp[g_ (catf + I satf) ] -> Exp[gcatf] (Cos[gsatf] + I Sin[g satf ]),
Exp [I h_ + j_ + g_ (catf + I satf) ] ->
(Cos[h] +ISin[h]) Exp[j] Exp[gcatf] (Cos [g satf] + I Sin[g satf ]) }
( (Cos [kx] + I Sin [k x] ) (-E~cacf ^ y ks k
(Cos[(l + f2)1,4satfyks] -ISin[(l + f2)1,4satfyks]j +
E'ky (1 + f2) 1/4 (catf + I satf) ks) v0) /
( k + (1 + f2)1,4 (catf + I satf) ks)
Next I multiply the numerator by the conjugate of the denominator:
vnum[x_, y_] = Numerator[vmag[x, y] /. (y-> -yy, -y->yy}]
Conjugate [Denominator [vmag[x, y] / {y -> -yy, -y -> yy} ] ];
vnum[x_, y_] = vnum[x, y] / {-yy -> y, yy -> -y} ;
vnum[x_, y_] = Simplify [vnum[x, y] ]
(Cosjkx] + I Sin[kx] )
( -k + (l + f2)1/4 (catf I satf) ks) (-Ecat£ -l.f2;1/1 y
k (Cos[(l + f2)14satfyks] -ISin[(l + f2)14satfyks]) +
E~ky (1 + f2)174 (catf + I satf) ks)
v0
Next I multiply the denominator by its conjugate. With this, all the complex portions of the
expression are linear factors in / = a/^T .
vden[x y_] =
Sinplify[Expand[Denominator[vmag[x, y] /. {y-> -yy, -y->yy}]
Conjugate[Dencminator[vmag[x, y] / {y -> -yy, -y -> yy} ] ] ] ] /
{catf A2 + satf A2 -> 1, catf -> Cos[ ArcTan[f] ] ]
k2 2
k Cos
ArcTan[f]
2
ks + V1 + f2 k2
Next I expand the numerator and set the complex portions to zero. Then I replace the satf and catf
variables with the actual functions. This yields an expression for the numerator of the real portion
of the u velocity solution:
vreal[x_, y_] = Expand[vnum[x, y] ]; vreal[x_, y_] =
vreal[x, y] /. {I-> ii, -I -> -ii}; vreal[x_, y_] =vreal[x, y] /.
{ ii -> 0, catf -> Cos[ ArcTan[f ] ], satf -> Sin[ ArcTan[f] ] }
52


Finally, I put the numerator over the denominator:
vreal[x_, y_] = vreal [x, y] / vden[x, y]
jE-l.>1'*yC0.[^1U-]k. k2 Cos[kx]
Cos [ (1 + f2)1/4 y Sin [ ArcT^n 1 £ 1 ] k,] Vo + E- u.f2 -1/1 yC[ ^-111] ks
k2 Sin [k x] Sin [ (1 + f2)1/4 y Sin [ ArCT^D ^ £^ ] ks ] v0 -
E'ky (1 + f2)1/4 k Cos [k x] Cos [ -r--^^L£-I ] ks v0 -
E- ll.f2)1,4 y Cos[ ArcT;"lf) ] ks (1 + f2)l'4 k Cos [k x]
cos [ ArCTfn^ ] cos [ (1 + f2)1/4 y Sin [ ] ks ] ks v0 +
Eky (1 + f2)1/4 k Sin[k x] Sin [ ArcT*n[£] j ks Vq _
E- Cos [ (1 + f2)1/4 y Sin [ ArcT^n[£] ] ks] sin[kx] Sin[ ArcT^ntf 1 ] ks Vo _
E- (i.£2)1/4 y cos[nzpiL] ,s (x + f2 j i/ k Cos [ ArgTan_[_fj_ j
Sin [k x] Sin [ (1 + f2)1/4 y Sin [ ArcT^n tf i j ks] ks Vo +
E- (1*£2, 1/4 y Cos[ AJISZJILLLL ] ks (1 + f2)l/ kCos[kx]
Sin [ ArcT^nt£] j sin[ (i + f2)1/4 y sin [ ArcT^n[f 1 ] ks] ks Vg +
E"ky V1 + f2 Cos [k x] Cos [ A.rcT^n[£] ]2 k2 Vq +
Eky Vl + f2 Cos [k x] Sin [ ArcT^n[£] ]2 k| Vq| j
(k2 2 (l + f2)1/4kCos[ArCT;n[£1 ] ks + VITf^ k2)
Now I move on to the pressure term.
Taking real portion of pIp:
To make the expressions easier to read, I will let catf = Cos [ -A ArcTan [ f ] ] and
satf=Sin [ -A ArcTan [ f ] ]. vmag represents the spatial magnitude function of the v velocity:
y_] =
pmag[x, y] / {Vka + I kj -> k. (f2 + l)1/4 (catf + I satf )}


pmag [ x y_] = Together[pmag[x, y] ]
E i k x k y (i + f2)1/4 (-Icatf + satf) w ks v0
k (k catf (1 + f2)1/4 ks I (1 + f2)1'4 satf ks)
Now I manipulate the expression above, taking the complex portions out of the exponential using
Euler's identity:
F*nag[x_, y_] =
E*nag[x, y] / {Exp[Ih_ + g__] -> Exp[g] (Cos[h] + I Sin[h]),
Exp[g_ (catf + I satf) ] -> Exp[g catf] (Cos [g satf] + I Sin[g satf]),
Exp [I h_ + j_ + g_ (catf + I satf) ] ->
(Cos[h] +ISin[h]) Exp[j] Exp[gcatf] (Cos [g satf] + I Sin[g satf]) }
(E~ky (1 + f2)1/4 (-1 catf + satf) w (Cos [ k x] + I Sin (k x] ) ks v0) /
(k (k catf (1 + f2)1/4 ks I (1 + f2)1'4 satf ks) )
Next I multiply the numerator by the conjugate of the denominator:
pnum[x_, y_] = Numerator [pmag [x, y] /. (y-> -yy, -y->yy}]
Conjugate[Dencminator[pmag[x, y] / (y -> -yy, -y -> yy} ] ];
pnum[x_, y_] = pnum[x, y] / {-yy -> y, yy -> -y} ;
pnum[x_, y_] = Sinplify[pnum[x, y] ]
E"ky (1 + f2 )1/4 k ( -1 catf + satf) u (Cos [kx] + I Sin [k x] ) ks
(k (1 + f2)1'4 (catf- I satf) ks) v0
Next I multiply the denominator by its conjugate. With this, all the complex portions of the
expression are linear factors in / = V^-T.
pden[x_, y_] =
Sirrplify[ Expand [Denominator [pmag [x, y] /. {y-> -yy, -y->yy}]
Conjugate[Denominator[pmag[x, y] / {y -> -yy, -y -> yy} ] ] ] ] / .
variables with the actual functions. This yields an expression for the numerator of the real portion
of the u velocity solution:
Next I expand the numerator and set the complex portions to zero. Then I replace the satf and catf
54


preal [x_, y_] = Expand[pnum[x, y] ]; preal[x_, y_] =
preal[x, y] /. {I -> ii, -I -> -ii}; preal[x_, y_] = preal[x, y] /.
{ ii -> 0, catf -> Cos[i ArcTan[f] ], satf -> Sin[ ArcTan[f] ] }
E"k y (1 +
Eky (1
,2,i/ f2 k2 ai Cos -----------
,2, 1/4. 2 r, r ArcTan [f
f2 k2 (j Cos k x Sin ------
Sin [k x] ks Vo +
ks v0 -
E~ky Vl + f2 k u Cos [ ArCT^n ^ ^^ ] Sin [k x] k2 v0 -
E ky a/ 1 + f2 kojSin[kx] Sin [ ArcT^n t ^ j k2 v0
Finally, I pul the numerator over the denominator:
preal[x_, y_] = preal[x, y] /pden[x, y]
E~ky (1 + f2) k2 ai Cos [---- J Sin[k x] ks v0
E ky (1 + f2) k2 a: Cos [ k x] Sin [--- J ks
E'ky V1 + f2 k a; Cos [ ArCT^n ^£ 1 ] Sin [k x] k2 v0
E"ky Vl + f2 k a) Sin [k x) Sinf ArcT^n[£] |2 k2 Vq
v0 -
k2 (k2 2 (1 + f2)1/4 k Cos
ArcTan[f1
ks + Vl + f2 k2
Use of asymptotic condition to simplify expressions and reveal underlying
flowfield structure
Recapitulation: / have found a solution to a linearized set of momentum equations using a
method employing stream and potential functions. / have found expressions for u,v and p/p. I
have solved for the real part of these expressions and have shown them to be quite complicated.
Now I will take advantage of certain conditions in the flowfield I am studying to considerably
simplify these expressions, and in so doing make apparent the underlying flowfield structure
which they represent.
55


Application of asymptotic condition
Recall that ks = yjai/v This shows that the viscous lengthscale depends upon the frequency of
wall motion. The frequency and wavelength of wall motion are independent variables. For the
bending actuator conditions I am studying, the frequency is such that the viscous wavenumber is
much greater than the bending wavenumber, that is ks k. In the expression for velocity and
pressure I have made use of a variable / = (^-) Thus as ks becomes very much larger than k, f
tends toward zero. In the following lines I take advantage of Mathematica's Limit function to find
what the expressions for velocity and pressure tend to in the asymptotic case:
usimp[x_, y_] = Limit[ureal[x, y], f -> 0]
E~y ik*ksi (Eky Eyks) Sin[kx] ks v0
k ks
vsimp[x_, y_] = Limit[vreal[x, y], f -> 0]
Cos [k x] (E~y ks k E'ky ks) v0
k ks
psiin>[x_, y_] = Limit [preal[x, y], f -> 0]
E'ky ai Sin [ k x] ks v0
k2 k ks
As can be seen, in the case where the viscous wavelength is much smaller than the bending
wavelength the form of these expressions is considerable simplified. It is interesting to note the
denominator of these simplified expressions. It seems to warn that these expressions are not valid
if ks ~ k, which is to be expected given the asymptotic condition from which they derive.
statement of simplified solutions for ease of development
E-y usinp[x_, y_] =
vsinp[x_, y_] =
psiirp[x_, y_] =
k k,
Cos[kx] (E-yk k-Ekyk.) v0
k-k.
E'ky oi Sin[kx] k, v0
k2 -kk. '
Interpretation of results
I will deal with the three expressions in order of increasing complexity: p, then v and finally u.
56


Pressure: the expression for pressure/density represents a Held varying sinusoidally in x
(wavelength k) and decaying exponentially in y (decay constant also k). A contour plot of
pressure is shown below:
ContourPlotfAbs[psinp[x, y] /. {ai-> (kJ -kk,) /k}] /
{v0 ->1, k->20, k,-> 20000}, {x, 0, 2 7r / 20}, {y, 0,
Contours ->{0.85 0.7, 0.55 0.4, 0.25, 0.1, 0.05,
FramaTicks-> {{{0, StyleForm[ "0", FontSize-> 12]},
[27t/ (4 20), StyleForm[ Ab", FontSize-> 12]},
{27t/ (2 20), StyleFormf "-^-Ab", FontSize-> 12]},
{2 7r 3 / (4 20), StyleForm[" Ab", FontSize -> 12] }
{2tt/ (20), StyleFormf "Ab", FontSize-> 12]}},
{{0, StyleForm["0", FontSize-> 12]},
{l / (4 20), StyleFormf"
{l / (2 20) StyleFormf"
{3 / (4 20) StyleFormf"
1
4kT
1
2 kt,
3
4^
n
n
n
FontSize -> 12]},
FontSize -> 12] } ,
FontSize -> 12]},
{1/ (20), StyleFormf" ", FontSize -> 12] },
ki,
{3/ (2 20), StyleFormf"----", FontSize-> 12]},
2 kjj
2
{2/20, StyleFormf" ", FontSize-> 12] }}, None,
kb
ColorFunction -> mycf, AspectRatio -> 1 / GoldenRatio,
PlotPoints -> 40]
2/20},
0.02 },
None},
57


Velocity v: this expression also varies sinusoidally in x, and is in quadrature with the pressure
field (i.e. cosine x dependance versus pressure's sine dependence). In y, the function is a
maximum at the wall and decays exponentially from there. The exponential with ks as the decay
constant is significant only very near the wall, and rapidly yields significance to the other
exponential, which governs the subsequent overall behavior. This is shown in the following two
plots, the first of which shows the region very close to the wall, and the second shows a region
extending farther from the wall.
58


velocity (v/v0)
Plot[vsini>[0, y] /. {v0 -> 1, k -> 20, k, -> 20000},
{y, 0, 10 / 20000}, Frame -> True, Axes -> False,
FrameLabel-> {StyleForm["distance from wall", FontSlze-> 16],
StyleForm["velocity (v/v0)", FontSize-> 16] },
FrameTicks-> {{{1/ 20000, StyleForm[ "l/k", FontSize-> 12]},
{5/20000, StyleForm[" 5 / k, ", FontSize-> 12] },
{10/20000, StyleForm["10/k.", FontSize-> 12] }} ,
Automatic, None, None},
Text Style-> { FontFami ly > "Helvetica", FontSize-> 12}]
0.99
0.99
0.99
0.99
1/ks
5/ks
distance from wall
10/ks
- Graphics -
59


Plot[vsinp[0, y] /. {v0 -> 1, k -> 20, k, -> 20000},
{y, 0, 3/20), Frame-> True,
FrameLabel -> {StyleFormf distance from wall", FontSize -> 16],
StyleFormf velocity (v/v0)n, FontSize-> 16] },
FrameTicks-> {{{10/20000, StyleFormf "10/k,", FontSize-> 12]},
{1/20, StyleForm["1/kb", FontSize-> 12] },
{2/20, StyleFormf 2/kb", FontSize-> 12]},
{3/20, StyleFormf 3/kt", FontSize-> 12] }},
Automatic, None, None},
TextStyle-> {FontFamily-> "Helvetica", FontSize-> 12}]
distance from wall
- Graphics -
Velocity u: this expressions varies sinusoidally in x, and is in phase with the pressure field. In y,
unlike the v velocity expression, which is at its maximum at the wall, u is zero at the wall (in
keeping with the boundary condition), but exponentially builds (with spatial growth factor^) as
distance from the wall is increased. At a distance not far from the wall it reaches a maximum
value and then begins an exponential decay (with decay constant kb). The distance at which it
reaches it maximum is easily determined:
peaky = Solve[dy using? [x, y] == 0, y]
{{y
Lg[^]
k ks
}}
I characterize this distance as the reciprocal of a wavenumber kc, or critical wavenumber:
f. __ _____________
l' Logl^-LoglA:,)
which represents the wavenumber of the crossover point between viscous and bending wavelength
dominated effects. The maximum u velocity is subsequently found as:
60


Simplify[usinf>[7r / (2 k), y] /. peaky[[1, 1] ] ]
;s K v0
k ks
At first this may appear to suggest that under the asymptotic condition under which I am working
that the peak u is zero. Making the substitution f=k//cvas before, it can be seen that in fact in the
limit as ks &,the maximum approaches vo. It is of note that this is the same peak velocity in u
predicted by the inviscid solution of Cremer.
% / k. -> k/ f
f (k -f)
Limit [H, f -> 0]
The u velocity behaviour is shown in the following two plots, the first of which shows the region
very close to the wall, and the second shows a region extending farther from the wall.
61


velocity (u/v0)
Plot [usinp[7r / (2 20) y] / {v0 -> 1, k-> 20, k. -> 20000} ,
{y, 0, 20 / 20000}, Frame -> True, Axes -> False,
FrameLabel -> {StyleFoxm[ "distance from wall", FontSize -> 16],
StyleForm["velocity (u/v0) ", FontSize -> 16]},
FrameTicks-> {{{1/ 20000, StyleForm["1/k,", FontSize-> 12]},
{Log[20 / 20000] / (20 20000), StyleForm[ "1/kc", FontSize-> 12] } ,
{3/ 20000, StyleForm[ "3/k,", FontSize-> 12] },
{10/ 20000, StyleFonn["10/k,", FontSize-> 12] } ,
{20/20000, StyleFoxm[ "20/k, ", FontSize -> 12] }} ,
Automatic, None, None},
TextStyle-> {FontFamily-> "Helvetica", FontSize-> 12}]
62


Plot [usimp[7r / (2 20) y] / {v0 -> 1, k-> 20, k, -> 20000},
{y, 0, 2/20}, Frame-> True,
FrameLabel -> {StyleForm["distance from wall", FontSlze -> 16],
StyleFonnf"velocity (u/v0) ", FontSize -> 16]},
FrameTicks -> {{{10/20000, StyleForm[ "10/k,", FontSlze-> 12]},
{1/20, StyleFonn["l/kb", FontSize-> 12]},
{2/20, StyleForm[ "2/k,,", FontSize-> 12]},
{3/20, StyleForm["3/kb", FontSize-> 12]}},
Automatic, None, None},
TextStyle -> {FontFamily-> "Helvetica", FontSize -> 12}]
distance from wall
- Graphics -
Another view of the u velocity is shown in the following contour plot. This highlights the
periodicity of the field in x and the decay in y.
mycf[c_] := GrayLevel[l (1 c) / 1.8]
63


ContourPlot[Abs[usinp[x, y] ] /. {v0 -> 1, k-> 20, k, -> 20000},
{x, 0, 2 7r / 20} {y, 0, 2/20},
Contours >{0.85 0.7, 0.55 0.4, 0.25, 0.1, 0.05, 0.02 },
FramaTicks -> {{{0, StyleForm["0", FontSlze-> 12]},
{2tt/ (420), StyleForm["Ab ", FontSlze-> 12]},
{2 tt/ (2 20), StyleForm[" -^-Ab", FontSize-> 12] },
{2 7r 3 / (4 20), StyleForm["-^-Ab", FontSize-> 12]},
{2 7r / ( 20), StyleForm["Ab", FontSize-> 12]}},
{{0, StyleForm[ "0", FontSize -> 12]},
{1/ (4 20), StyleFormf"-^-", FontSize-> 12]},
{1/ (2 20), StyleForm["---", FontSize-> 12]},
2 kt,
[3/ (4 20), StyleForm["", FontSize-> 12]},
[1/ (20), StyleFonn[" ", FontSize-> 12]},
kb
[3/ (2 20), StyleForm["--", FontSize-> 12] } ,
2 kb
2
{2/20, StyleForm[" ", FontSize -> 12] } }, None, None},
kb
ColorFunction -> mycf, As pact Ratio -> 1 / GoldenRatio,
PlotPoints -> 40]
64


65


ContourPlot [ Abs[usimp[x, y] ] /. {v0 -> 1, k-> 20, k. -> 20000},
{x, 0, 2 7r / 20} {y, 0, 10/20000},
Contours -> {0.85 0.7 0.55 0.4 0.25, 0.1, 0.05, 0.02 },
FrameTicks -> {{{0, StyleFormf"0, FontSize-> 12]},
{2x/ (4 20), StyleFormf "-^-Ab", FontSize-> 12]},
{2x/ (2 20), StyleForm[" Ab", FontSize-> 12]},
{2x3/ (420), StyleForm[" Ab", FontSize-> 12] } ,
{2 x / (20), StyleForm[ "Ab", FontSize-> 12]}},
{{0, StyleForm["0", FontSize-> 12]},
{1/ (20000), StyleFormf" ", FontSize-> 12]},
{3 / (20000), StyleFormf" ", FontSize-> 12]},
{Log[20 / 20000] / (20 20000), StyleFormf" ", FontSize -> 12] },
{10/ (20000), StyleFormf"-^-", FontSize-> 12]},
20
{20/ (20000), StyleFormf"---", FontSize -> 12]}}, None, None},
ColorFunction -> mycf, AspectRatio -> 1 / GoldenRatio,
PlotPoints -> 40]
66


- ContourGraphics-
Finally, I present some velocity vector field plots
pO = PlotVectorField[{using?[x, y] /. {v0 -> 1, k-> 20, k. -> 20000},
vsinp[x, y] / {v0 -> 1, k -> 20, k, -> 20000} } {x, 0, 2 7r / 20} ,
{y, 0, 2/20}, PlotPoints - {16, 9}, AspectRatio-> 1 /GoldenRatio,
ScaleFactor -> (#8 .0&), ScaleFunctlon-> (# .04&), Frame -> True,
FrameTicks-> {{{0, StyleForm["0", FontSlze-> 12]},
{2 7t/ (4 20), StyleFormf" Ab", FontSlze-> 12]},
{2 tr/ (2 20), StyleFonn["yAbn, FontSlze-> 12]},
(2tt3/ (420), StyleForm[" Ab ", FontSlze-> 12] } ,
{2 7r / (20), StyleForm[ nAb", FontSize-> 12] } } ,
{{0, StyleForm[ "0", FontSlze -> 12] } ,
{1 / (4 20), StyleForm["
{l / (2 20), StyleForm["
{3 / (4 20) StyleForm["
4 kt
1
2 kt
3
4kT
{1/ (20), StyleFormf" ",
kb
FontSlze -> 12
FontSize -> 12
FontSize -> 12
FontSize -> 12]} ,
]}
]}
)}'
{3 / (2 20), StyleForm["
(2/20, StyleForm[" "
kb
----", FontSize > 121 },
2 kb J1
, FontSize -> 12] }}, None
None}]
67


_2_
kb
3
2 kb
3kb
4 kb
2k^
r i 4 4 < \ f r f > " > V i
! i 4 / < * \ i r f * -* > V i 1
\ * 4 / < * \ 1 F t -* > V i
i 4 / < ^ \ F F t + > V i i
L * 4 \ \ F F t -* > V i k
( k 4 ' < ^ ^ \ i * t > *4 > V y k
it f y ^ v \ t
U i f / ^ ^ \ \ 1 i * A ^ "V. \ \ t
k t., A \ , 1 / X F i t f i i
1 F 1 i i i F F
0 i, 1 , 3 Ah
TAb yAb 4Ab
- Graphics -
68


Appendix B SUBROUTINE NEWMSH
C sees ID @(#)newmsh.F 1.4 03/21/97
c-----------------------------------------------------------------------c
c c
C NAME : NEWMSH C
C C
C PROGRAM : FLUENT C
C VERSION : V4.3x C
C C
C ARGS DESCRIPTION C
C C
C C
C C
C PURPOSE: GENERATES GRID AT EACH TIME STEP AND COMPUTES GRID C
C VELOCITY. USEFUL FOR MOVING MESH PROBLEMS. C
C C
C COMMENTS: THE SUBROUTINE NEWMSH WHEN CALLED SHOULD DEFINE A NEW C
C GRID WHICH HAS THE SAME NUMBER OF CELLS IN THE I, J AND C
C K DIRECTIONS AS THE GRID IN THE ORIGINAL CASE FILE. C
C THIS AMOUNTS TO DEFINING THE X-Y-Z COORDINATES OF THE C
C GRID POINTS (IN METERS), i.e. THE COORDINATES OF THE C
C INTERSECTION OF THE GRID LINES WHICH YOU SEE WHILE C
C DISPLAYING THE GRID. IN FLUENT, THESE ARE STORED IN THE C
C VARIBALES XC, YC, and ZC. NEWMSH ROUTINE WILL BE CALLED C
C AT EVERY TIME STEP OF A TRANSIENT CALCULATION WHEN THE C
C DEFORMING-MESH OPTION IS ON. C
C TARGET VARIABLES: C
C XC, a real array which stores the x-coords in meters. C
C YC, a real array which stores the y-coords in meters. C
C ZC, a real array which stores the z-coords in meters. C
C C
C SUPPORTING VARIABLES: C
C TCSTEP, a real variable, stores current time-step size C
C NI,NJ and NK, integer variables, which store the number C
C of grid points as used internally by C
C FLUENT, e.g. NI is 2 plus the number of C
C internal computational cells. C
C C
C ALGORITHM SPEC: C
C Compute NI*NJ*NK values in the XC, YC and ZC arrays at C
C time step n+1, given those values at time step n, and C
C the time step size, TCSTEP. The grid coordinates must C
C be defined in meters. C
C C
C C
C CALLED BY: PICALC, TNEXT C
C C
69


c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
CREATED : 02/08/94
AUTHOR : Narayan Behera ,FLUENT, INC., LEBANON, NH
REVISION HISTORY :
NO. BY DATE MODIFICATIONS/COMMENTS
001 AWR 07/07/98 sinusoidal boundary condition
002 AWR 09/04/98 make #time steps/cyc USRPAR3
#include
#include
iinclude
iinclude
#include
#include
#include
IMPLICIT.INC"
SIZE.INC"
DMESH.INC"
GRID.INC"
FLPGEO.INC"
TDTEMP.INC
USPARS.INC"
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C LOCAL VARIABLES
C
INTEGER NIJK,
REAL YTAU
REAL XDIM
REAL K
REAL OMEGA
REAL YAMP
REAL M
REAL TSPERCYC
#include "LINDEX.INC
#include "LINDEX.DEF
I, J, L
! y velocity
! the width of the grid in x
! bending wavenumber
! bending radian frequency
! bending amplitude, USPAR1 is peak velocity
grid y displacement decay exponent-USPAR2
! # of time steps per half cycle of motion-
! note USPAR3 is the number of steps/cycle
C
c---------------------------------------------------------------------------c
c
NIJK = NI NJ NK
C...set full width of structured grid to equal 1/2 of bending wavelength
XDIM=ABS(XC(LINDEX(NI,2,1))-XC(LINDEX(2,2,1)))
K=3.141592654/XDIM
C...set radian frequency to 36 times the timestep (TCSTEP constant)if
C USPAR3=0, otherwise use USPAR3 value
TSPERCYC=USPAR3/2.0
IF(USPAR3.LE.0.0) TSPERCYC=18
OMEGA=3.141592654/(TSPERCYC*TCSTEP)
70


C...set amplitude as determined by user-defined peak velocity
YAMP=USPAR1/OMEGA
M=USPAR2
WRITE(*,*) K,OMEGA,YAMP,M,TIME1
DO J=2,NJ,1
DO 1=2,NI,1
IF(J.LT.(NJ-2)) THEN
YTAU=YAMP*OMEGA*COS(K*XC(LINDEX(I,J,1)))*COS(OMEGA*TIMEl)/
(NJ/(NJ-J+2))**M
ELSE
YTAU=0.0
END IF
YC(LINDEX(I,J,1)) = YC(LINDEX(I,J,1)) + YTAU TCSTEP
END DO
END DO
C...monitor node positions
WRITE(*,*) NINT(NI/2.)+1,XC(LINDEX(NINT(NI/2.)+1,3,1)),
YC(LINDEX(NINT(NI/2.)+1,3,1)),
YC(LINDEX(2,2,1))
RETURN
END
71


Bibliography
[1] Brillioun, M.J., Problems de rayonnement en acoustique du batiment,
Acustica 2, pp. 65 76, 1952
[2] Clark, R.L., and Fuller, C.R., Active structural acoustic control with adaptive
structures including wavenumber considerations, Recent Advances in Ac-
tive Control of Sound and Vibration, Rogers, C.A. and Fuller, C.R., Eds.,
Virginia Polytechnic Institute, April, 1991
[3] Cremer, L., Heckl, M. and Ungar, E.E., Structure-Borne Sound; Structural
Vibrations and Sound Radiation at Audio Frequencies, Springer-Verlag,
Berlin, 2nd edition, 1988
[4] Fox, R.W., and McDonald, A.T., Introduction to Fluid Mechanics, 2ed., Wi-
ley and Sons, NY, 1978
[5] Lamb, H., Hydrodynamics, Dover, New York, 1945
[6] Leissa, A.W., The free vibration of rectangular plates, J. Sound Vib. 31(3),
pp.257-293, 1973
[7] Lighthill, J., Waves in Fluids, Cambridge University Press, Cambridge, 1978
[8] Ljunggren, S., Forced vibrations of infinite plates, J. Sound Vib. 121 (2),
221-236, 1988
[9] Meirovitch, L. and Thangjitham, S., Control of Sound Radiating from an Or-
thotropic Plate, J. Vib. and Acoustics, Vol. 114 October 1992
[10] Meirovitch, L., Analytical Methods in Vibrations, Macmillan Co., London,
1967
[11] Norton, M.P., Fundamentals of Noise and Vibration Analysis for Engi-
neers, Cambridge University Press, 1989
[12] Patanker, S.V., and Spalding, D.B., A calculation procedure for heat, mass
and momentum transfer in three-dimensional parabolic flows, Int. J. Heat
Mass Transfer 15, pp. 1787-1806, 1972
72


[13] Patanker, S.V., Numerical Heat TVansfer and Fluid Flow, Hemisphere, New
York, 1980
[14] Peric, M., Kessler, R., and Scheurer, G., Comparison of finite-volume nu-
merical methods with staggered and collocated grids, Comput. Fluids 16(4),
pp. 389-403, 1988
[15] Schlichting, H., Boundary-layer Theory, 7th Edition, McGraw-Hill, Inc.,
New York, 1987
[ 16] Segal, Lee A., Mathematics Applied to Continuum Mechanics, Macmillan
Publishing Co., New York, 1977
[17] for a discussion of Stokes first and second problems, see (for example), Pan-
ton, R.L., Incompressible Flow, 2ed. pp 263-272, John Wiley and Sons, New
York, 1996
[18] Wolfram, S., Mathematica, Cambridge University Press, 1996
73