- Email: [email protected]
Finite element approach to steady-state vibrations in a fluid of finite depth
COMPUTER METHODS NORTH-HOLLAND
IN APPLIED
MECHANICS
AND
ENGINEERING
46 (1984) 259-276
FINITE ELEMENT APPROACH TO STEADY-STATE VIBRATIONS IN A FLUID OF FINITE DEPTH* Kazimierz Polish Academy
SZMIDT
of Sciences, Institute of Hydroengineering,
Revised
PL-80-953,
Received 12 April 1983 manuscript received 9 September
Gdatisk, Poland
1983
The paper deals with the finite element analysis of waves in a semi-infinite layer of the considerations is to examine the problem of irregular spacing of points in the layer on results of numerical solutions. It is shown that introduction of a nonhomogeneous in the layer introduces a nonhomogeneity into the field. The discrete system may different from the continuous field.
fluid. The aim of and its influence spacing of points have properties
0. Introduction Many problems in hydroelasticity can be treated as steady-state wave propagation problems in a layer of fluid. It is possible to find analytical solutions only for special cases. For many practical problems it is necessary to apply numerical methods. In these methods the continuum is replaced by a discrete space of chosen points. In a general case the discrete model may have properties different from the continuous field. Such differences are important when the discretization leads to solutions which have no physical meaning but result from discretization. The properties of the discrete solutions of the steady-state hydrodynamical problems for a semi-infinite layer of fluid were discussed by Wilde and Szmidt in [3,4]. In [4] the discrete solution was constructed on the basis of finite differences. The accuracy of the description by finite differences was discussed and compared to the analytic solution of the problem. Special attention was focused on solutions in a case of nonhomogeneous spacing of points in the layer of fluid. The paper [3] deals with the solution of the problem by the finite element method. The discrete solution was constructed for a case of regular spacing of points in the layer. Conditions for the spacing of an assumed net necessary to obtain a reasonable approximation of the analytic solution of the problem were given. The present paper deals with the same problems and in a sense it is a continuation of [3]. The solution of the hydrodynamical problem is constructed with the help of the finite element method. Special attention is paid to the problem of a nonhomogeneous spacing of points in the layer of fluid and its influence on results of the numerical solution. As in [3,4], three types of boundary conditions on the upper surface of the layer are considered. The conditions correspond to cases of no flow through the boundary, zero pressure and a surface-gravitational *This research 00457825/84/$3.00
was sponsored
by the Polish Academy
@ 1984, Elsevier
Science
Publishers
of Sciences
within program
B.V. (North-Holland)
05.12.
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
260
wave. The algebraic equations formulation.
of the problem
are derived with the help of a variational
1. Formulation of the problem Let us consider a semi-infinite layer of fluid of depth h, see Fig. l(a). It is assumed that the fluid is compressible, nonviscous and that disturbances are so small that potential theory can be applied. The movement of the fluid is induced by an assumed velocity field at x = 0. The differential equation of the problem has the following form:
(1.1) where @ is the potential of the velocity field, c the velocity of sound; V2 denotes the Laplace operator and dots denote derivatives with respect to time. The potential must satisfy the differential equation (1.1) and corresponding boundary conditions. The boundary conditions for z = 0 and for x = 0 (see Fig. l(a)) are -aQ,
a2
r=O=Oy
g=
x
0
= K,
(1.2)
where V, denotes the given velocity at the boundary x = 0. Three conditions on the upper surface of the layer are considered:
a@
0,
(p-p.&)
r=h
cases of boundary
=o,
(1.3)
where p is a density of the fluid and g the gravitational acceleration. The first condition in (1.3) means that there is no flow through the boundary z = h. In the second equation of (1.3) the pressure is zero on the surface. The third condition in (1.3) corresponds to the surfacegravitational wave. To the conditions (1.2) and (1.3) the Sommerfeld condition must be added that the waves die out or propagate to the right at infinity.
Fig. 1. A layer of fluid and assumed discrete points.
Fig. 2. Finite elements.
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
For steady-state formula @(x, Substitution
2,
it is possible
to introduce
a spatial
potential
cp(x, z) according
t) = cp(x, 2) - elm’.
261
to the
(1.4)
of (1.4) into (1.1) yields V2q + k2, -cp=
0,
k, = o/c.
(1.5)
The classical analytical solution of the equation is obtained by separation of variables [4]. Let us consider the solution of the problem by the finite element method. The assumed net of the discrete model is shown in Fig. l(b). The spacing of the vertical lines is a while the horizontal lines have a spacing b, in the lower part and b2 in the upper part of the layer. Without loss of generality let us assume that bI < bZ. For the differential equation (1.5) and the boundary conditions considered, the appropriate variational formulation is [3]
where f(z) is defined by V,= w -f(z). For the case of no flow through the upper boundary and zero pressure on the boundary the last integral in (1.6) must be omitted. The terms which appear in (1.6) have physical meaning. The surface integral corresponds to the kinetic and strain energy of the fluid. The second integral in (1.6) corresponds to the power introduced by the assumed velocity field at x = 0. The line integral on the free surface is connected with the potential energy of the gravitational waves. The integrals are evaluated over rectangular and linear finite elements, respectively. The finite elements are shown in Fig. 2. For the elements the following shape functions are assumed:
c(x,z)=-&* [w(x-a)+-b)-
(P2.X.(Z-b)+(P3’X’Z-(P4.(X-a).2],
(1.7) ~(z)=~.[-~l.(Z-b)+(P2.z].
~(x)=~.[-~I.(x-a)+~Z.X],
Substitution of the first in (1.7) into the surface integral (1.6) and integration over one rectangular element yields a polynomial of second degree with respect to the values of CJIat the corners of the element. The result may be written as follows: T, =
i(S)t *[K] *(6))
where (8)’ = (cp,, (p2, (Pi, (p4). The elements:
(1.8) symmetric
square
(4 x 4) matrix
[K1] has
the
following
262
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
Substitution of the second of (1.7) into the last integral of (1.6) and integration gives a formula analogous to (1.8) but in this case we have two nodal points and a square (2 X 2) matrix with the elements io2. a
k12= k21=
co2-a $. g
(1.10)
*
If f(z) may be approximated by a linear shape function, then from the second integral of (1.6) and the third of (1.7) we have T2
=
(S)t
* [K2]
’
(f) ,
(1.11)
where (S)t = ((pl, (Pi)! (f)’ = (fl,f2) and the elements of the matrix [K2] are k*,
=
k22
=
4b
7
k12
=
k21
=
&b.
(1.12)
Taking the results (1.8)-(1.12) into account and differentiating the final linear set of equations of the problem,
WI ’ w =(P) >
with respect to qi,j, one obtains
(1.13)
(P~.~, . . . , c2, . . .>, (P> is a column matrix connected with the given velocity at x = 0, and [K] is a square matrix of coefficients in the equations. For a typical point in the region only the elements in the neighbourhood enter into the equation for the considered point. For a point (i, j) (i denotes a horizontal line corresponding to a row of the matrix while j denotes a vertical line corresponding to a column of the matrix) in the lower part of the layer, the equation is
where (S)t= hl,
k 13*qii-1,,-1+
2k14* CP~-I~+ k13. Vi-l,j+l
+ 2k12 *+Iz,j-l+ 4kll . upi,j + 2kl2 +Cp
i,J+l
+
k13 *(Pi+l,j-l + 2k14 *
(pi+l,J
+
k13
*
(Pi+l,j+l
=
0
.
(1.14)
A remark is needed. The terms kiJ (i, j = 1,. . . ,4) in the last equation are described by (1.9)
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
263
but in this case one should substitute b1 instead of b in (1.9). Similarly, in the equations written for points in the upper part of the layer one has to substitute bz for b in (1.9). To distinguish the terms kij for the elements of the upper part of the layer we use capital letters Kij for these terms. A formal approach to the problem leads to an infinite set of equations (1.14) written for all nodal points of the assumed net. The solution has to satisfy the Sommerfeld conditions which mean that no progressive wave comes from infinity and for propagating and standing waves the values do not increase when going to infinity. Thus let us look for a solution of the form: Cpi.j+s =
Cp 1.1
.exp(-s.r*a),
(1.15)
where s is a natural number. The last equation represents a dependance between values of +CJ at the points (Xi, Z, = const.) and (Xj+s, Zi = const.). For real and positive values of r the solution (1.15) represents a standing wave which decays when going to infinity. If r is imaginary a progressive wave may be obtained. Substitution of (1.15) into (1.14) gives (k14+ hk13) * P
i-l,j
+
2
(kll + hk12)
’
’ Cp r,j +
(kl4
+
* (pi+l,j
hk13)
=
0
9
(1.16)
where A = cash ra .
(1.17)
Equation (1.16) refers to the unknown values of pi, in points belonging to one vertical line x, = const. Hence we may omit the index j and instead of the infinite set of equations (1.14) we have the finite set (1.16). The equations are written for all i including the bottom and the upper boundary. The result is a set of equations which can be written in the following form:
WI + A*[Bl) *(cp)= 0 .
(1.18)
The matrices [A]and [B] are shown in Tables 1 and 2. The matrices are real and symmetric. Table 1
Table 2
P
*, .
Bl =
[Al =
I 1
K,< 1 2K,,
[
.
K:, = K,,
p’o
K,: = 2Ko K,: = K,, - ;.“+
2
K,,
9:
1 K,d 1 K;
.
2K,, K,4
K,4 K::
= K,z
264
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
The eigenvalues A are real numbers if at least one of the matrices is positive definite [l, p. 3411. In the following we will show that for the cases of no flow through the upper boundary and zero pressure on the boundary the matrix [A]is positive definite. The matrix [A]will be positive definite if the determinants CA, = det[A,], 1= 1,2, . . . , n, where [A,]means a matrix obtained by crossing the last (n - 1) rows and columns in the matrix [A],are all greater than zero. For the case of a gravitational wave on the upper boundary we will establish a necessary and sufficient condition that [A] be positive definite. Let us transform the matrix [A] according to the formula
[Al = PI *[Cl
(1.19)
’
Let us consider first the case kr4>0. matrix [D] has the following form:
In the case K1, > 0 and the diagonal of the diagonal n-k
k [D]diag=
(LTLFXCb
(1.20)
K,,,...,)7
where k denotes the point where the spacing changes (see Fig. l(b)). The matrix [C] is shown in Table 3. One can see that in this case the (n - 1) first elements of the principal diagonal of the matrix [C] are all larger than 2. Hence we may make the substitutions Rl=p=coshx,
j&=!+
14
cash y , 14
p=fE, 14
(1.21) cu=coshx+&coshy,
x,y>O,
Because the matrix [D] is positive definite it is enough to confine our attention to the matrix [C] only. To evaluate all the determinants CDr = det(Cij\, i, j = 1,2,. . . , 1 of the matrix we apply the well-known theorem which states that the value of a determinant is left unchanged, if the elements of a row are altered by adding to them any constant multiple of the Table 3
p=
2
ct=R,+P.Rn
R*=Ra
.
.
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
265
corresponding elements in any other row [2, p. 4061. According to the theorem the first row of the matrix [C] (Table 3) is multiplied by R;’ and subtracted from the second row, then the changed second row is multiplied by C;,l and subtracted from the third row and so on. The procedure leads to a triangular matrix, the elements of which below the diagonal are all zeros. Due to (1.21) the elements of the diagonal of the resulting matrix can be described as follows:
a, =
cash rx cosh(r - 1)x ’
r= 1,2 ,...,
k,
(1.22)
cosh(s - k)y + G * sinh(s - k) y b, = cosh(s-k-l)y+G.sinh(s-k-1)y’
s=k+l,...,n-1,
where G is a constant. According to the last relations the diagonal of the triangular matrix [C”] may be written in the following form: [C*ldlag = [cash X, cos:h’,“, . . . , w,
ak,
(cash y + G *sinh y) ,
cosh(n-k-l)y+G.sinh(n-k-1)y cash 2y + G - sinh 2y cash y + G - sinh y ’ * * * ‘cosh(n-k-2)y+G*sinh(n-k-2)y’
1’
b n
(1.23)
where (1.24)
ak = sinh x - tanh(k - 1)x + /3 - cash y .
The last element b,, in (1.23) depends on a specified boundary condition on the upper surface of the layer. In the case of no flow through the upper boundary the element is b,, = sinh y -
tanh(n-k-l)y+G l+G.tanh(n-k-1)y’
(1.25)
For the case of zero pressure on the upper boundary the corresponding equal to N = n - 1. In this case the element bN assumes the form bN =
cosh(N-k)y+G.sinh(N-k)y cosh(N-k-l)y+G.sinh(N-k-1)y’
In the case of a gravitational holds:
number of equations is
(1.26)
wave on the upper surface of the layer the following relation
tanh(n-k-l)y+G 1 bn = sinh y ’ 1 + G - tanh(n - k - 1)y --*3K,,
1~2.a g
*
(1.27)
To obtain the determinants CD, = det]Cij], i, j = 1,. . . , 1 one should take into account the first 1 terms in (1.23). It is easy to see that the first k determinants of the matrix [Cj are all greater
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
266
than zero. Let us consider the values of b,. The term may be rewritten as b1 = 2 * cash y - p - $ = cos y + &* [/3 - sinh’y + sinh x * tanh(k - 1)x - cash y] . (1.28) One can see that because x > 0 and y > 0 all the terms appearing on the right-hand side of the last relation are greater than zero. Comparison of the last result with (1.23) shows that the following relation holds: G *sinh y = k. It is easy to Substitution flow through definite. For (1.23) leads follows:
+ sinh x - tanh(k - 1)x * cash y] .
(1.29)
see that G > 0 and all the determinants CDk+l,. . . , CD,_, are greater than zero. of (1.25) and (1.26) into (1.23) leads to the conclusion that for both cases of no the upper boundary and zero pressure on the boundary, the matrix [A] is positive the case of a gravitational wave on the upper boundary, substitution of (1.27) into to a condition that [A] be positive definite. The condition may be written as
bz. co2
-<
[p . sinh’y
tY
[
1+
2 ’ _ f . (b,.
0
kc)21 . tanh y
. tanh(n - k - ‘)y + G 1+ G * tanh(n - k - 1)y ’
(1.30)
To complete the considerations let us discuss the remaining cases. For the case k14= 0, K14 > 0 the first (k - 1) elements of the diagonal of the matrix [D] are equal to one while the last (n - k + 1) elements are equal to K,,. Changes occur in the substitution (1.21). In this case cr = cash y + kll/Kl, and Rr = k,l. The elements C1.2, Ck,k-l and Ci,i-1, Ci,f+l (i = 2,3, . . . , k - 1) of the matrix [C] are equal to zero. To obtain a triangular matrix corresponding to the matrix [Cj one should apply only the second relation in (1.22) for the last (n - k) rows of the matrix. Because kll > 0 and K1, > 0, the first k determinants of the matrix are greater than zero. The other determinants of the matrix may be evaluated in the way we have shown for the previous case. Hence, the inequality (1.30) is still valid but in this case G *sinh y =
y . cash y + sinh’y , o y+cosh y ’
k ’ =e*
(1.31)
The next case to be considered is k14< 0 and K14> 0. The elements Dj (j = 1,2, . . . , k) of the diagonal of the matrix [D] are equal to -k14 while the others remain unchanged. In this case p = -K14/k14 >O and R1 = -kll/k14. One can show that R1 > 1. The elements C1,2, C,,_l and Ci,i-1, Ci,i+l (i = 293, * e * 9 k - 1) of the matrix [q change sign. The corresponding triangular matrix is obtained by a similar transformation as in the first case but now multiplied rows are added to the next ones. The further procedure is as in the first case and leads to the same results. The subsequent case is k14< 0 and K,, = 0. The elements Dj (j = k + 1, . . . , n) of the diagonal of the matrix [D] are equal to one while the remaining elements are as in the previous case. With reference to the previous case, the elements Cn+l, Ci,i-1, Ci,i+l (i = k + 1, ...)n - 1) of the matrix [C] are equal to zero. In this case /3 = 0. To obtain a triangular
K. Szmidt, Finite element approach to steady-state vibrations in a &id of finite depth
267
matrix corresponding to the matrix [q one should apply only the first relation of (1.22) for the first k rows of the matrix. Instead of the inequality (1.30) we have the simpler formula (1.32) The last case, k14c 0 and K14< 0, can be treated as the first case. With reference to (1.20) the elements of the diagonal matrix [D] have opposite signs. The elements situated above and below the principal diagonal of the matrix [C] (see Table 3) change sign. In this case R, > 1 and R2 > 1 so the substitutions (1.22) are valid. The procedure is similar to that in the first case and leads to the same conclusions. With this case all possibilities are considered. For a special case of regular spacing of points (bl = b2 = b) in the layer, one can show that G = tanh kx and instead of the relations (1.30) and (1.32) we have b - co2
< [ 1 + (iy
-
$(b - /cc)21 * tanh
-
x - tanh(n - 1)x.
(1.33)
g
For many cases it is possible to simplify the condition (1.30). If (n - k - 1)y is a large number, say (n - k - 1)y 2 10, we may assume that tanh(n - k - 1)y = 1, and instead of the relation (1.30) we may use the approximated form be
y
2
<
[I + (2)’
-
f(b, - kc)21 . tanh
y.
(1.34)
The value of (b, * kc)’ is very small compared to the other terms in the last relation and may be ignored in practical calculations. In the remainder of this paper we are dealing with cases for which the matrix [A] is positive definite. Thus, the eigenvalues A are real numbers. According to the notation given in (1.17) , if A 3 1, r is real (standing wave), if IAl < 1, r is pure imaginary (progressive wave) and if A s -1 then r=i.Areacosh(-A)+i*(2n-l)T, where n is a natural number. In the last case the solution (1.15) corresponds V
i,j+s
=
i=v--1,
(1.35)
to a standing wave of the form
(- 1)” * cpl,j * exp[- s * Re(ru)] .
(1.36)
Such a solution does not correspond to the one for a continuum. For further considerations it is convenient to rewrite (1.18) in the following form:
[A*].(cp)=
0.
(1.37)
The matrix [A*] is simply obtained by adding the matrices [A] and A . [B], see (1.18). Let us assume that k14+ A *k13Z 0 and K14+ A *K13f 0. If the first assumption is not satisfied then
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
268
the first (k - 1) components of the eigenvector (cp) will be equal to zero. It means that a nontrivial solution may exist only for a layer hl d z < h (Fig. l(b)). If the second assumption is not satisfied then the last (n - k) components of the vector will be equal to zero. In such a case a nontrivial solution may exist only for a layer 0 < z < h1 (Fig. l(b)). These special cases are of no interest and are excluded from the considerations. According to the above assumptions the first k equations (1.37) are multiplied by -l/(k14+ A . ku), the remaining equations by -1/(K14 + A * K13), respectively. To simplify the notation the following substitutions are made: sl=
kll+A-ku -(b
s
s**=
s*=sl+as*,
Kll+A’Kl*
=
*
+ A - kn) ’
-&+AX&
Klz,+A.Kn (y = -(k14 + A . k13) ’
G,+A.fS -(Ku+ A .&)’
(1.38)
From the equations it follows that
kll+&h,_
A
Kll + S2.
Ku
(1.39)
=-k12+S1.k/-K12+S2.K1s*
According to the last notations,
PI *(44= 0
instead of (1.37) we may consider (1.40)
7
where the matrix [S] is shown in Table 4. The matrix is not symmetric. It should be stressed that S, and S2 are not independent. Elimination of A from the two first equations (1.38) leads to the relations s =2(1-&)+S*(1+2&) 1 (2+&)+S**(l-E) Table 4
’
s =2(1-&)-&‘(2+& * s,.(1-8)-(1+2e))’
b1 *<1 . (1.41) &= 06
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
The parameter
269
(Y can be expressed as follows:
a=&;.[l+2
3-G (2+ &)+ S**(l-
*&-&(1-E)]=
The values of Si and S, are real numbers. It is convenient S2 =
s1 = cash K ,
(1.42)
&)’
to make the following substitutions: (1.43)
cash 7 .
The values K and 7 are complex numbers whose real and imaginary parts are not arbitrary. The numbers may assume the following values: K =
ul+i.uZ,
where u1 = 0 or u2 = j
q =
u1+i*v2,
where ul=Oor
* T, j = 0, 1,2, . . . ,
~~=j*7~,j=O,1,2
Pw
,....
It is enough to limit our attention to the positive values of u1 and ul. One can see that if S1 3 1, if u2 = (2 .j - l)n, j = 1,2, . . . , S1 d - 1 and if u1 = 0, then z42=2+n, j=O,l,..., -1 G S1 < 1. The same cases for 7 lead to the same results. The relations (1.43) and (1.44) enable us to represent an arbitrary real number S1 (or S,) in the form given by (1.43). When (p2 is expressed in terms of (pl from the first equation in (1.40) and then (p3 from the second and so on, it results in [4] cpr=
(~1.
cosh(r -
Cps= (Pk * cosh(s
r =
l)K,
- k )q
+A
*sinh(s-
where A is a constant. Substitution of (1.45) into (1.40) for point
1,2, . . . , k
s = k +
k)q,
k (zk = hl,
,
(1.45)
1,. . . , ~1,
see Fig. l(b)) gives (1.46)
Let us consider first the case of no flow through the upper boundary. Substitution of (1.45) into the nth equation in (1.40) yields after simple manipulations the following relation: sinh 77*[ql *cosh(k - 1)~ . sinh(n -
k)q
+ A . cosh(n
-
k)q]
= 0 .
(1.47)
There exists a nontrivial solution of the set of equations (1.46) and (1.47) if the determinant equal to zero. Thus it follows that cosh(k - 1)~ * sinh 7 * cosh(n - k)v
.
[sinh
K *
is
tanh(k - 1)~
+ a * sinh 7 . tanh(n -
k)v]
= 0.
(1.48)
The last equation is the necessary condition for nontrivial solutions of (1.46) and (1.47). It is not easy to calculate in a general case all the roots of the equation, but it is possible to discuss
270
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finitedepth
properties of the solutions (1.45) corresponding to the condition (1.48). Let us discuss all possible cases connected with values of S1 and SZ. Case 1. SZ> 1. From the relations (1.41) and (1.42) it follows that the values of S1 belong to the open interval; S1 E (1, (1+ 2~)/(1- c)) and (YE (0, V/E). In this case K = u1 and n = ol. Now in (1.48) all the terms are positive, thus the equation has no solution. Case 2. SZ = 1. For this case the second of (1.41) leads to the result S1 = 1. Thus the equation (1.48) has the roots K = q = 0. The nontrivial solution corresponding to the roots is cp,= (pl = const. ,
r = 1,2, . . . , n .
(1.49)
Case 3. -lG&
(1.50)
Nontrivial solutions have the following form: r= 1,2 ,...,
cpI= (PI*cos(r - l)Uz ) $&= (Pk’ co+
- k ) v2 + A * sin(s
- k)v2,
s= k +
k,
(1.51)
1, . . . , n ,
Case 4. SZ< -1.
In this case (1.41) result in two solutions: S, > (1 + 2&)/(1- .Y)and S1 < (1 - 4&)/(1+ 2~). For the first case (1.42) yields (YCO. One can show that there are no solutions of (1.48) for this case. For the second case instead of the inequality S1 < (1- 4&)/(1+ 2~) we may consider two possibilities: S1 < -1 and S1 E (-1, (1 - 4&)/(1+ 28)). If S1 -< - 1 the determinant has no solution. Thus a nontrivial solution is possible only for S1 E (-1, (1 - 4&)/(1+ 28)). The case corresponds to S1 = cos u2 and S2 = -cash ul. The determinant is as follows: cos(k - 1)~~. sinh u1 - cosh(n - k)vI - [sin u2 - tan(k - l)u2
+ (Y- sinh u1 * tanh(n - k)q] = 0. The case corresponds
(1.52)
to the following solutions:
cp, = (PI*cos(r - l)uz )
r=1,2
CD,= (-l)“-k[V k.cosh(s-k)2)1+A.sinh(s-k)u1],
s=k+l,...,n.
,...,
Such solutions do not correspond to the ones for a continuum. Let us now consider the case of zero pressure on the upper boundary. into the nth equation in (1.40) yields the relation
k,
(1.53)
Substitution
of (1.45)
271
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
(~1. cosh(k - 1)K - cosh(n - k)q + A - sinh(n - k)q = 0. The determinant
(1.54)
of (1.46) and (1.54) assumes the form
cosh(k - 1)~ -cosh(n - k)q - [sinh
K *
tanh(k - 1)~ - tanh(n - k)q
+&*sinhq]=O.
(1.55)
As in the previous case for S, > 1, the equation has no solution. For the case Sz = 1 (n = 0) the last equation is satisfied but (1.54) leads to the trivial solution. The discussion of the cases - 1s S2 < 1 and S, < - 1 is similar to the previous case and does not give new basic ideas. The last case is one of a gravitational wave on the free surface of the layer. The last equation of (1.40) is transformed into the following one: ql *cosh(k - 1)~ * [sinh 7 - sinh(n - k)v +AThe determinant
[
sinhn*cosh(n-k)q-$7.
-f 7.
(2 + cash 7) - cosh(n - k)g]
(2+coshv)*sinh(n-k)v]=O.
(1.56)
of (1.46) and (1.56) is
cosh(k - { sinh
1)K
K -
* cosh(n
- k)q -
tanh(k - 1)~ * [ sinh 77-$ y
sinhq.tanh(n-k)T-$v.
* (2 +
cash r)) - tanh(n - k)?]
(2+coshq)
11
=O.
(1.57)
Let us consider first the case S2 > 1. In the case K = ul, r] = u1 and S1 = cash ul, S2 = cash ul. All the functions appearing in the last equations are positive. Thus the solution may exist only if
fb.23
(2+
cash ul) - tanh(n - k)u, < sinh v1 < $9
’
-
2+coshu1 tanh(n - k)vl ’
(1.58)
For the case tanh(n - k)ul = 1 it is easy to calculate the root of (1.57). From the equation follows that sinh v1 1 b 2’02 2+cosh ol=‘g’ Simple manipulations
lead to the solution
it
(1.59)
272
K. Szmidt, Finite element approach to steady-state
vibrations in a fluid of finite depth
The corresponding eigenvalue A may be obtained from the relations (1.39). The nontrivial solution is described by (1.45) in which one has to substitute K = u1 and 77= ul. For the case S, = 1 (7 = 0), (1.57) is satisfied, but (1.56) leads to a trivial solution. In the case -1s SZ< 1 the determinant (1.57) assumes the form cos(k - l)u, . cos(n - k)uz . . sin u2. tan(k - l)u2 * sin v~--$ !C-&. g I
[
+ (Y * sin
V2- sin u2. tan(n -
(2+ cos 02) - tan(rt - (2+cosuz)
k)V2 + f y
I]
-
k)U2 I
=o.
(1.61)
It is not easy to calculate the roots of (1.61) and (1.41). The corresponding expressed by (1.51). The last case S2 < -1 results in the following determinant:
solutions
are
cos(k - l)u, - cosh(n - k)ul .
I
sin u2 - tan(k - l)u, * sinh v1 + 3 g” ’ O2 - (2 - cash u,) * tanh(n [
+ (Y - sinh
The corresponding
v1 - sinh v1 * tanh(n -
k)ul
+ i y.
k)vl]
(2 - cash ul)]} = 0 .
(1.62)
solutions are described by (1.53).
2. Detailed solutions for special cases The analysis presented up to now gives an insight into the solutions corresponding to each of the considered boundary conditions on the upper surface of the layer. In the following we are dealing with a special case for which it is possible to assume that tanh(n - k)q = 1. For this case all the determinants (1.48) (1.55) and (1.57) h ave a common multiplier. By equating the multiplier to zero we have the following equation:
sinh
K -
tanh(k - 1)~ + (Y* sinh q = 0 ,
(24
where q = o1 + iu 2, o2 = 2j - T or u2 = (2j - l)~, j = 1,2, . . . . If there exist solutions of (1.41) and (2.1), then the solutions do not depend on the specified boundary conditions on the upper surface of the layer. In other words, the solutions are insensitive to the boundary conditions. Such solutions do not correspond to the continuum and follow the discretization of the problem. According to the previous considerations the last equation leads to nontrivial solutions only if S,<-1
and
S,E
(2.2)
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
273
For the last cases (2.1) assumes the following form: sin u * tan(k - 1)~ + (Y. sinh ZJ= 0,
(2.3)
where u = u2 and v = vl. The cases (2.2) correspond to the solutions (1.53). The last equations of (2.3) must be supplemented by the relations (1.41). It may be seen that for an unknown value of (Y* sinh v > 0 there are k - 1 roots of the equation in the interval (0, n), The roots are in the intervals 2j- 1 2’ 2 - (k - 1) .~
j=l,2
,...,
k-l.
(2.4)
If we continue the calculation of roots for (k - 1). ZA> n, the important values will lead to the same solutions. Thus it is enough to restrict our attention to the interval (0, 71). It does not mean that there are k - 1 solutions of (1.41) and (2.3). According to the second relation in (2.2) it may happen that
(k-l).arccose>s.n,
(2.5)
where s is a natural number. For such a case the first s roots of (2.3) and (1.41) should be crossed out and the number of eigenvectors (1.53) must be reduced to the number k - l- s. According to the considerations it is possible to define bounds for eigenvalues. To make the discussion easier let us neglect the term k2, in (1.9). This term is very small compared to the other terms in the equations. Simple calculations show that for w = 100 set-‘, k: = 4.46 - lo-‘; thus the simplification seems to be justified for our purposes. Substitution of (1.9) into (1.38) and (2.2) gives _l
s 2
=Z.(l+[)+A.(l-25)<1-4s 2-5+/U(1+5)
=2*(1+l)+h*(1-25) 2-5+A*(l+<)
Let us consider three cases: Case 1.2-l>O. For this case the expressions 4+l --<*c-z. 2-5 Case 2. 2-6~0 and 2-5>0. This case follows the solutions
<-1,
1+2&’ +)*,
b, 2 5= (-)a ’ !$=&a[.
(2.6)
(2.6) lead to the following result:
P-7)
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
274
A<-%
or
(2.8)
A>%.
Case 3. 2- 5<0. The solution for this case assumes the form 4+5 -
(2.9)
l-2
According to (1.39) and the above relations it is possible to get a more precise description of intervals for values of S1 and Sz. To illustrate the considerations, an example is presented below. For the case shown in Fig. l(b) let us assume the spacing of points: a = 1 cm, br = 0.4cm, b2 = 2 cm, it = 29, k = 11. It is seen that 5 = 0.16, 6 = 4.0, E = 0.4 and the case corresponds to the relations (2.8). For the parameters, the relations (2.8) give A < -2.26087 or A b4.0. If A < -2.26087, then from (1.39) it follows: - 1 < S, < 0.5862069 and - 1.9411763 < S2< -1.40. Tacking into account the upper bound for S1 = cos u, one can get (k - 1). u = (k - 1) - arccos 0.5862069 = 3.00625 * n. Hence, the number of roots of (2.3) is reduced to the number k - 1 - s = 7. The corresponding eigenvalues Ai are in the intervals: Al E (-9.6802, -5.4469), . . . , AgE (-2.3106, -2.2823), A7E (-2.2661, -2.2609). The solution obtained by a standard numerical procedure (1.18) results in: A1= -7.4048, . .. . , A6= -2.3088, A, = -2.2660. Comparison of the above results shows that the intervals obtained by analytical considerations are in good agreement with the result of numerical computations. Let us consider the case A > 4.0. In this case l-4& 0.5862069 c S, = cos u < -1 + 2E = 0.777777
and
- 1.40 < Sz < - 1.0 .
The bounds result in 2.1634689 * T < (k - 1) - u < 3.00625 - IT. It is easy to see that in this interval there is only one root of (2.3). For the case the first in (1.39) results in the relation 7.2653695 < A < 676.43221. The obtained interval is very wide. The eigenvalue A found by a
Fig. 3. Eigenvectors
corresponding
to equation
(2.1).
K. Szmidt. Finite element approach to steady-state vibrations in a fluid of finite depth
275
standard numerical procedure is equal to 17.0621. Thus, in conclusion, there are 8 eigenvectors corresponding to the roots of (2.3). Some of the corresponding eigenvectors are plotted in Fig. 3. From the graphs it is seen that the components of the eigenvectors are going to zero when one moves along a vertical line in the z-direction. The influence of the solutions on the final results depends of course on the character of a given problem. In cases when higher harmonics are important the influence may be significant.
3. Example of application To get a better insight into the discussed problem let us consider numerical examples. For comparison let us consider two cases of excitements of the fluid motion: V,(z, t) = w * d * cos ot
o - d - cos wt K(z,
t) =
o
forOsz
(3.1)
for 0 s t G hl , forhl
Our aim is to determine the pressure distribution along the wall OA, see Fig. 4. To simplify the solution the zero pressure condition on the upper surface of the layer is assumed. To describe the influence of an assumed irregular net on the final results, the solutions are constructed for two different nets. Both the two nets have an equal number of points in the vertical direction. The first net is a regular one while the second has a denser spacing of points in the lower part of the layer. For comparison the analytical solution of the problem is also given. The obtained results of numerical calculations are plotted in Fig. 4. For a regular distribution of the assumed velocity field at x = 0 (first of (3.1)) the discrete solutions for both regular and irregular nets are in very good agreement with the analytical results. For the
a)
b)
T*P=O WC 24
f’
hs4orm,
-
Fig. 4. Pressure
ANALYTICAL
distribution
SOLUTIONS
along the wall CA.
276
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
second of (3.1) the numerical solutions are different from the analytical one. The difference is about 15% with reference to the analytical solution. It should be stressed that the solution for the regular spacing of points is better than in the case of the nonhomogeneous spacing of points in the layer. From the graphs it is easy to see how important for obtaining a reasonable solution is a proper spacing of points in the layer. In cases when higher harmonics are important (second of (3.1)) the better way of discrete solution seems to be applying a homogeneous spacing of points instead of a net with a denser spacing in one part of the layer.
4. Conclusions For the three cases of boundary conditions on the upper surface of the layer considered in this paper the analytical solutions obtained by separation of variables in Cartesian coordinates exist and are well known [4]. If the continuum is replaced by a discrete space of chosen points an analogous solution may be obtained. In the paper the discrete solution is constructed with the help of the finite element method in cases of nonhomogeneous spacing of points in the layer. The discrete solution leads to velocity profiles which do not appear in the solution of the differential equation of the problem. Thus the introduction of a nonhomogeneous spacing of points introduces a nonhomogenity into the field. One should be careful in evaluating results of the finite element method because the discrete system may have properties different from the continuous field. The numerical model gives good results when one considers typical distributions of velocities in the layer. In cases when higher harmonics are important the influence of discretization of the problem on final results may be significant if the number of points is small.
Acknowledgment The author wishes to thank Professor P. Wilde for his valuable remarks during preparation of the work.
References [l] R.E. Bishop, G.M. Gladwell and S. Michaelson, Matrix Analysis of Vibrations (Wyd. Nauk. Techn., Warszawa, 1972) (in Polish). [2] E. Kreyszig, Advanced Engineering Mathematics (Wiley, New York, 1962). [3] P. Wilde and K. Szmidt, F.E. analysis of waves in a semi-infinite layer of fluid, Proc. 4th Conf. Finite Element Methods in Water Resources, Hannover, 1982. [4] P. Wilde and K. Szmidt, Numerical analysis of waves in a semi-infinite layer of fluid, Comput. Meths. Appl. Mech. Engrg. 36(l) (1983) l-23.
IN APPLIED
MECHANICS
AND
ENGINEERING
46 (1984) 259-276
FINITE ELEMENT APPROACH TO STEADY-STATE VIBRATIONS IN A FLUID OF FINITE DEPTH* Kazimierz Polish Academy
SZMIDT
of Sciences, Institute of Hydroengineering,
Revised
PL-80-953,
Received 12 April 1983 manuscript received 9 September
Gdatisk, Poland
1983
The paper deals with the finite element analysis of waves in a semi-infinite layer of the considerations is to examine the problem of irregular spacing of points in the layer on results of numerical solutions. It is shown that introduction of a nonhomogeneous in the layer introduces a nonhomogeneity into the field. The discrete system may different from the continuous field.
fluid. The aim of and its influence spacing of points have properties
0. Introduction Many problems in hydroelasticity can be treated as steady-state wave propagation problems in a layer of fluid. It is possible to find analytical solutions only for special cases. For many practical problems it is necessary to apply numerical methods. In these methods the continuum is replaced by a discrete space of chosen points. In a general case the discrete model may have properties different from the continuous field. Such differences are important when the discretization leads to solutions which have no physical meaning but result from discretization. The properties of the discrete solutions of the steady-state hydrodynamical problems for a semi-infinite layer of fluid were discussed by Wilde and Szmidt in [3,4]. In [4] the discrete solution was constructed on the basis of finite differences. The accuracy of the description by finite differences was discussed and compared to the analytic solution of the problem. Special attention was focused on solutions in a case of nonhomogeneous spacing of points in the layer of fluid. The paper [3] deals with the solution of the problem by the finite element method. The discrete solution was constructed for a case of regular spacing of points in the layer. Conditions for the spacing of an assumed net necessary to obtain a reasonable approximation of the analytic solution of the problem were given. The present paper deals with the same problems and in a sense it is a continuation of [3]. The solution of the hydrodynamical problem is constructed with the help of the finite element method. Special attention is paid to the problem of a nonhomogeneous spacing of points in the layer of fluid and its influence on results of the numerical solution. As in [3,4], three types of boundary conditions on the upper surface of the layer are considered. The conditions correspond to cases of no flow through the boundary, zero pressure and a surface-gravitational *This research 00457825/84/$3.00
was sponsored
by the Polish Academy
@ 1984, Elsevier
Science
Publishers
of Sciences
within program
B.V. (North-Holland)
05.12.
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
260
wave. The algebraic equations formulation.
of the problem
are derived with the help of a variational
1. Formulation of the problem Let us consider a semi-infinite layer of fluid of depth h, see Fig. l(a). It is assumed that the fluid is compressible, nonviscous and that disturbances are so small that potential theory can be applied. The movement of the fluid is induced by an assumed velocity field at x = 0. The differential equation of the problem has the following form:
(1.1) where @ is the potential of the velocity field, c the velocity of sound; V2 denotes the Laplace operator and dots denote derivatives with respect to time. The potential must satisfy the differential equation (1.1) and corresponding boundary conditions. The boundary conditions for z = 0 and for x = 0 (see Fig. l(a)) are -aQ,
a2
r=O=Oy
g=
x
0
= K,
(1.2)
where V, denotes the given velocity at the boundary x = 0. Three conditions on the upper surface of the layer are considered:
a@
0,
(p-p.&)
r=h
cases of boundary
=o,
(1.3)
where p is a density of the fluid and g the gravitational acceleration. The first condition in (1.3) means that there is no flow through the boundary z = h. In the second equation of (1.3) the pressure is zero on the surface. The third condition in (1.3) corresponds to the surfacegravitational wave. To the conditions (1.2) and (1.3) the Sommerfeld condition must be added that the waves die out or propagate to the right at infinity.
Fig. 1. A layer of fluid and assumed discrete points.
Fig. 2. Finite elements.
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
For steady-state formula @(x, Substitution
2,
it is possible
to introduce
a spatial
potential
cp(x, z) according
t) = cp(x, 2) - elm’.
261
to the
(1.4)
of (1.4) into (1.1) yields V2q + k2, -cp=
0,
k, = o/c.
(1.5)
The classical analytical solution of the equation is obtained by separation of variables [4]. Let us consider the solution of the problem by the finite element method. The assumed net of the discrete model is shown in Fig. l(b). The spacing of the vertical lines is a while the horizontal lines have a spacing b, in the lower part and b2 in the upper part of the layer. Without loss of generality let us assume that bI < bZ. For the differential equation (1.5) and the boundary conditions considered, the appropriate variational formulation is [3]
where f(z) is defined by V,= w -f(z). For the case of no flow through the upper boundary and zero pressure on the boundary the last integral in (1.6) must be omitted. The terms which appear in (1.6) have physical meaning. The surface integral corresponds to the kinetic and strain energy of the fluid. The second integral in (1.6) corresponds to the power introduced by the assumed velocity field at x = 0. The line integral on the free surface is connected with the potential energy of the gravitational waves. The integrals are evaluated over rectangular and linear finite elements, respectively. The finite elements are shown in Fig. 2. For the elements the following shape functions are assumed:
c(x,z)=-&* [w(x-a)+-b)-
(P2.X.(Z-b)+(P3’X’Z-(P4.(X-a).2],
(1.7) ~(z)=~.[-~l.(Z-b)+(P2.z].
~(x)=~.[-~I.(x-a)+~Z.X],
Substitution of the first in (1.7) into the surface integral (1.6) and integration over one rectangular element yields a polynomial of second degree with respect to the values of CJIat the corners of the element. The result may be written as follows: T, =
i(S)t *[K] *(6))
where (8)’ = (cp,, (p2, (Pi, (p4). The elements:
(1.8) symmetric
square
(4 x 4) matrix
[K1] has
the
following
262
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
Substitution of the second of (1.7) into the last integral of (1.6) and integration gives a formula analogous to (1.8) but in this case we have two nodal points and a square (2 X 2) matrix with the elements io2. a
k12= k21=
co2-a $. g
(1.10)
*
If f(z) may be approximated by a linear shape function, then from the second integral of (1.6) and the third of (1.7) we have T2
=
(S)t
* [K2]
’
(f) ,
(1.11)
where (S)t = ((pl, (Pi)! (f)’ = (fl,f2) and the elements of the matrix [K2] are k*,
=
k22
=
4b
7
k12
=
k21
=
&b.
(1.12)
Taking the results (1.8)-(1.12) into account and differentiating the final linear set of equations of the problem,
WI ’ w =(P) >
with respect to qi,j, one obtains
(1.13)
(P~.~, . . . , c2, . . .>, (P> is a column matrix connected with the given velocity at x = 0, and [K] is a square matrix of coefficients in the equations. For a typical point in the region only the elements in the neighbourhood enter into the equation for the considered point. For a point (i, j) (i denotes a horizontal line corresponding to a row of the matrix while j denotes a vertical line corresponding to a column of the matrix) in the lower part of the layer, the equation is
where (S)t= hl,
k 13*qii-1,,-1+
2k14* CP~-I~+ k13. Vi-l,j+l
+ 2k12 *+Iz,j-l+ 4kll . upi,j + 2kl2 +Cp
i,J+l
+
k13 *(Pi+l,j-l + 2k14 *
(pi+l,J
+
k13
*
(Pi+l,j+l
=
0
.
(1.14)
A remark is needed. The terms kiJ (i, j = 1,. . . ,4) in the last equation are described by (1.9)
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
263
but in this case one should substitute b1 instead of b in (1.9). Similarly, in the equations written for points in the upper part of the layer one has to substitute bz for b in (1.9). To distinguish the terms kij for the elements of the upper part of the layer we use capital letters Kij for these terms. A formal approach to the problem leads to an infinite set of equations (1.14) written for all nodal points of the assumed net. The solution has to satisfy the Sommerfeld conditions which mean that no progressive wave comes from infinity and for propagating and standing waves the values do not increase when going to infinity. Thus let us look for a solution of the form: Cpi.j+s =
Cp 1.1
.exp(-s.r*a),
(1.15)
where s is a natural number. The last equation represents a dependance between values of +CJ at the points (Xi, Z, = const.) and (Xj+s, Zi = const.). For real and positive values of r the solution (1.15) represents a standing wave which decays when going to infinity. If r is imaginary a progressive wave may be obtained. Substitution of (1.15) into (1.14) gives (k14+ hk13) * P
i-l,j
+
2
(kll + hk12)
’
’ Cp r,j +
(kl4
+
* (pi+l,j
hk13)
=
0
9
(1.16)
where A = cash ra .
(1.17)
Equation (1.16) refers to the unknown values of pi, in points belonging to one vertical line x, = const. Hence we may omit the index j and instead of the infinite set of equations (1.14) we have the finite set (1.16). The equations are written for all i including the bottom and the upper boundary. The result is a set of equations which can be written in the following form:
WI + A*[Bl) *(cp)= 0 .
(1.18)
The matrices [A]and [B] are shown in Tables 1 and 2. The matrices are real and symmetric. Table 1
Table 2
P
*, .
Bl =
[Al =
I 1
K,< 1 2K,,
[
.
K:, = K,,
p’o
K,: = 2Ko K,: = K,, - ;.“+
2
K,,
9:
1 K,d 1 K;
.
2K,, K,4
K,4 K::
= K,z
264
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
The eigenvalues A are real numbers if at least one of the matrices is positive definite [l, p. 3411. In the following we will show that for the cases of no flow through the upper boundary and zero pressure on the boundary the matrix [A]is positive definite. The matrix [A]will be positive definite if the determinants CA, = det[A,], 1= 1,2, . . . , n, where [A,]means a matrix obtained by crossing the last (n - 1) rows and columns in the matrix [A],are all greater than zero. For the case of a gravitational wave on the upper boundary we will establish a necessary and sufficient condition that [A] be positive definite. Let us transform the matrix [A] according to the formula
[Al = PI *[Cl
(1.19)
’
Let us consider first the case kr4>0. matrix [D] has the following form:
In the case K1, > 0 and the diagonal of the diagonal n-k
k [D]diag=
(LTLFXCb
(1.20)
K,,,...,)7
where k denotes the point where the spacing changes (see Fig. l(b)). The matrix [C] is shown in Table 3. One can see that in this case the (n - 1) first elements of the principal diagonal of the matrix [C] are all larger than 2. Hence we may make the substitutions Rl=p=coshx,
j&=!+
14
cash y , 14
p=fE, 14
(1.21) cu=coshx+&coshy,
x,y>O,
Because the matrix [D] is positive definite it is enough to confine our attention to the matrix [C] only. To evaluate all the determinants CDr = det(Cij\, i, j = 1,2,. . . , 1 of the matrix we apply the well-known theorem which states that the value of a determinant is left unchanged, if the elements of a row are altered by adding to them any constant multiple of the Table 3
p=
2
ct=R,+P.Rn
R*=Ra
.
.
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
265
corresponding elements in any other row [2, p. 4061. According to the theorem the first row of the matrix [C] (Table 3) is multiplied by R;’ and subtracted from the second row, then the changed second row is multiplied by C;,l and subtracted from the third row and so on. The procedure leads to a triangular matrix, the elements of which below the diagonal are all zeros. Due to (1.21) the elements of the diagonal of the resulting matrix can be described as follows:
a, =
cash rx cosh(r - 1)x ’
r= 1,2 ,...,
k,
(1.22)
cosh(s - k)y + G * sinh(s - k) y b, = cosh(s-k-l)y+G.sinh(s-k-1)y’
s=k+l,...,n-1,
where G is a constant. According to the last relations the diagonal of the triangular matrix [C”] may be written in the following form: [C*ldlag = [cash X, cos:h’,“, . . . , w,
ak,
(cash y + G *sinh y) ,
cosh(n-k-l)y+G.sinh(n-k-1)y cash 2y + G - sinh 2y cash y + G - sinh y ’ * * * ‘cosh(n-k-2)y+G*sinh(n-k-2)y’
1’
b n
(1.23)
where (1.24)
ak = sinh x - tanh(k - 1)x + /3 - cash y .
The last element b,, in (1.23) depends on a specified boundary condition on the upper surface of the layer. In the case of no flow through the upper boundary the element is b,, = sinh y -
tanh(n-k-l)y+G l+G.tanh(n-k-1)y’
(1.25)
For the case of zero pressure on the upper boundary the corresponding equal to N = n - 1. In this case the element bN assumes the form bN =
cosh(N-k)y+G.sinh(N-k)y cosh(N-k-l)y+G.sinh(N-k-1)y’
In the case of a gravitational holds:
number of equations is
(1.26)
wave on the upper surface of the layer the following relation
tanh(n-k-l)y+G 1 bn = sinh y ’ 1 + G - tanh(n - k - 1)y --*3K,,
1~2.a g
*
(1.27)
To obtain the determinants CD, = det]Cij], i, j = 1,. . . , 1 one should take into account the first 1 terms in (1.23). It is easy to see that the first k determinants of the matrix [Cj are all greater
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
266
than zero. Let us consider the values of b,. The term may be rewritten as b1 = 2 * cash y - p - $ = cos y + &* [/3 - sinh’y + sinh x * tanh(k - 1)x - cash y] . (1.28) One can see that because x > 0 and y > 0 all the terms appearing on the right-hand side of the last relation are greater than zero. Comparison of the last result with (1.23) shows that the following relation holds: G *sinh y = k. It is easy to Substitution flow through definite. For (1.23) leads follows:
+ sinh x - tanh(k - 1)x * cash y] .
(1.29)
see that G > 0 and all the determinants CDk+l,. . . , CD,_, are greater than zero. of (1.25) and (1.26) into (1.23) leads to the conclusion that for both cases of no the upper boundary and zero pressure on the boundary, the matrix [A] is positive the case of a gravitational wave on the upper boundary, substitution of (1.27) into to a condition that [A] be positive definite. The condition may be written as
bz. co2
-<
[p . sinh’y
tY
[
1+
2 ’ _ f . (b,.
0
kc)21 . tanh y
. tanh(n - k - ‘)y + G 1+ G * tanh(n - k - 1)y ’
(1.30)
To complete the considerations let us discuss the remaining cases. For the case k14= 0, K14 > 0 the first (k - 1) elements of the diagonal of the matrix [D] are equal to one while the last (n - k + 1) elements are equal to K,,. Changes occur in the substitution (1.21). In this case cr = cash y + kll/Kl, and Rr = k,l. The elements C1.2, Ck,k-l and Ci,i-1, Ci,f+l (i = 2,3, . . . , k - 1) of the matrix [C] are equal to zero. To obtain a triangular matrix corresponding to the matrix [Cj one should apply only the second relation in (1.22) for the last (n - k) rows of the matrix. Because kll > 0 and K1, > 0, the first k determinants of the matrix are greater than zero. The other determinants of the matrix may be evaluated in the way we have shown for the previous case. Hence, the inequality (1.30) is still valid but in this case G *sinh y =
y . cash y + sinh’y , o y+cosh y ’
k ’ =e*
(1.31)
The next case to be considered is k14< 0 and K14> 0. The elements Dj (j = 1,2, . . . , k) of the diagonal of the matrix [D] are equal to -k14 while the others remain unchanged. In this case p = -K14/k14 >O and R1 = -kll/k14. One can show that R1 > 1. The elements C1,2, C,,_l and Ci,i-1, Ci,i+l (i = 293, * e * 9 k - 1) of the matrix [q change sign. The corresponding triangular matrix is obtained by a similar transformation as in the first case but now multiplied rows are added to the next ones. The further procedure is as in the first case and leads to the same results. The subsequent case is k14< 0 and K,, = 0. The elements Dj (j = k + 1, . . . , n) of the diagonal of the matrix [D] are equal to one while the remaining elements are as in the previous case. With reference to the previous case, the elements Cn+l, Ci,i-1, Ci,i+l (i = k + 1, ...)n - 1) of the matrix [C] are equal to zero. In this case /3 = 0. To obtain a triangular
K. Szmidt, Finite element approach to steady-state vibrations in a &id of finite depth
267
matrix corresponding to the matrix [q one should apply only the first relation of (1.22) for the first k rows of the matrix. Instead of the inequality (1.30) we have the simpler formula (1.32) The last case, k14c 0 and K14< 0, can be treated as the first case. With reference to (1.20) the elements of the diagonal matrix [D] have opposite signs. The elements situated above and below the principal diagonal of the matrix [C] (see Table 3) change sign. In this case R, > 1 and R2 > 1 so the substitutions (1.22) are valid. The procedure is similar to that in the first case and leads to the same conclusions. With this case all possibilities are considered. For a special case of regular spacing of points (bl = b2 = b) in the layer, one can show that G = tanh kx and instead of the relations (1.30) and (1.32) we have b - co2
< [ 1 + (iy
-
$(b - /cc)21 * tanh
-
x - tanh(n - 1)x.
(1.33)
g
For many cases it is possible to simplify the condition (1.30). If (n - k - 1)y is a large number, say (n - k - 1)y 2 10, we may assume that tanh(n - k - 1)y = 1, and instead of the relation (1.30) we may use the approximated form be
y
2
<
[I + (2)’
-
f(b, - kc)21 . tanh
y.
(1.34)
The value of (b, * kc)’ is very small compared to the other terms in the last relation and may be ignored in practical calculations. In the remainder of this paper we are dealing with cases for which the matrix [A] is positive definite. Thus, the eigenvalues A are real numbers. According to the notation given in (1.17) , if A 3 1, r is real (standing wave), if IAl < 1, r is pure imaginary (progressive wave) and if A s -1 then r=i.Areacosh(-A)+i*(2n-l)T, where n is a natural number. In the last case the solution (1.15) corresponds V
i,j+s
=
i=v--1,
(1.35)
to a standing wave of the form
(- 1)” * cpl,j * exp[- s * Re(ru)] .
(1.36)
Such a solution does not correspond to the one for a continuum. For further considerations it is convenient to rewrite (1.18) in the following form:
[A*].(cp)=
0.
(1.37)
The matrix [A*] is simply obtained by adding the matrices [A] and A . [B], see (1.18). Let us assume that k14+ A *k13Z 0 and K14+ A *K13f 0. If the first assumption is not satisfied then
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
268
the first (k - 1) components of the eigenvector (cp) will be equal to zero. It means that a nontrivial solution may exist only for a layer hl d z < h (Fig. l(b)). If the second assumption is not satisfied then the last (n - k) components of the vector will be equal to zero. In such a case a nontrivial solution may exist only for a layer 0 < z < h1 (Fig. l(b)). These special cases are of no interest and are excluded from the considerations. According to the above assumptions the first k equations (1.37) are multiplied by -l/(k14+ A . ku), the remaining equations by -1/(K14 + A * K13), respectively. To simplify the notation the following substitutions are made: sl=
kll+A-ku -(b
s
s**=
s*=sl+as*,
Kll+A’Kl*
=
*
+ A - kn) ’
-&+AX&
Klz,+A.Kn (y = -(k14 + A . k13) ’
G,+A.fS -(Ku+ A .&)’
(1.38)
From the equations it follows that
kll+&h,_
A
Kll + S2.
Ku
(1.39)
=-k12+S1.k/-K12+S2.K1s*
According to the last notations,
PI *(44= 0
instead of (1.37) we may consider (1.40)
7
where the matrix [S] is shown in Table 4. The matrix is not symmetric. It should be stressed that S, and S2 are not independent. Elimination of A from the two first equations (1.38) leads to the relations s =2(1-&)+S*(1+2&) 1 (2+&)+S**(l-E) Table 4
’
s =2(1-&)-&‘(2+& * s,.(1-8)-(1+2e))’
b1 *<1 . (1.41) &= 06
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
The parameter
269
(Y can be expressed as follows:
a=&;.[l+2
3-G (2+ &)+ S**(l-
*&-&(1-E)]=
The values of Si and S, are real numbers. It is convenient S2 =
s1 = cash K ,
(1.42)
&)’
to make the following substitutions: (1.43)
cash 7 .
The values K and 7 are complex numbers whose real and imaginary parts are not arbitrary. The numbers may assume the following values: K =
ul+i.uZ,
where u1 = 0 or u2 = j
q =
u1+i*v2,
where ul=Oor
* T, j = 0, 1,2, . . . ,
~~=j*7~,j=O,1,2
Pw
,....
It is enough to limit our attention to the positive values of u1 and ul. One can see that if S1 3 1, if u2 = (2 .j - l)n, j = 1,2, . . . , S1 d - 1 and if u1 = 0, then z42=2+n, j=O,l,..., -1 G S1 < 1. The same cases for 7 lead to the same results. The relations (1.43) and (1.44) enable us to represent an arbitrary real number S1 (or S,) in the form given by (1.43). When (p2 is expressed in terms of (pl from the first equation in (1.40) and then (p3 from the second and so on, it results in [4] cpr=
(~1.
cosh(r -
Cps= (Pk * cosh(s
r =
l)K,
- k )q
+A
*sinh(s-
where A is a constant. Substitution of (1.45) into (1.40) for point
1,2, . . . , k
s = k +
k)q,
k (zk = hl,
,
(1.45)
1,. . . , ~1,
see Fig. l(b)) gives (1.46)
Let us consider first the case of no flow through the upper boundary. Substitution of (1.45) into the nth equation in (1.40) yields after simple manipulations the following relation: sinh 77*[ql *cosh(k - 1)~ . sinh(n -
k)q
+ A . cosh(n
-
k)q]
= 0 .
(1.47)
There exists a nontrivial solution of the set of equations (1.46) and (1.47) if the determinant equal to zero. Thus it follows that cosh(k - 1)~ * sinh 7 * cosh(n - k)v
.
[sinh
K *
is
tanh(k - 1)~
+ a * sinh 7 . tanh(n -
k)v]
= 0.
(1.48)
The last equation is the necessary condition for nontrivial solutions of (1.46) and (1.47). It is not easy to calculate in a general case all the roots of the equation, but it is possible to discuss
270
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finitedepth
properties of the solutions (1.45) corresponding to the condition (1.48). Let us discuss all possible cases connected with values of S1 and SZ. Case 1. SZ> 1. From the relations (1.41) and (1.42) it follows that the values of S1 belong to the open interval; S1 E (1, (1+ 2~)/(1- c)) and (YE (0, V/E). In this case K = u1 and n = ol. Now in (1.48) all the terms are positive, thus the equation has no solution. Case 2. SZ = 1. For this case the second of (1.41) leads to the result S1 = 1. Thus the equation (1.48) has the roots K = q = 0. The nontrivial solution corresponding to the roots is cp,= (pl = const. ,
r = 1,2, . . . , n .
(1.49)
Case 3. -lG&
(1.50)
Nontrivial solutions have the following form: r= 1,2 ,...,
cpI= (PI*cos(r - l)Uz ) $&= (Pk’ co+
- k ) v2 + A * sin(s
- k)v2,
s= k +
k,
(1.51)
1, . . . , n ,
Case 4. SZ< -1.
In this case (1.41) result in two solutions: S, > (1 + 2&)/(1- .Y)and S1 < (1 - 4&)/(1+ 2~). For the first case (1.42) yields (YCO. One can show that there are no solutions of (1.48) for this case. For the second case instead of the inequality S1 < (1- 4&)/(1+ 2~) we may consider two possibilities: S1 < -1 and S1 E (-1, (1 - 4&)/(1+ 28)). If S1 -< - 1 the determinant has no solution. Thus a nontrivial solution is possible only for S1 E (-1, (1 - 4&)/(1+ 28)). The case corresponds to S1 = cos u2 and S2 = -cash ul. The determinant is as follows: cos(k - 1)~~. sinh u1 - cosh(n - k)vI - [sin u2 - tan(k - l)u2
+ (Y- sinh u1 * tanh(n - k)q] = 0. The case corresponds
(1.52)
to the following solutions:
cp, = (PI*cos(r - l)uz )
r=1,2
CD,= (-l)“-k[V k.cosh(s-k)2)1+A.sinh(s-k)u1],
s=k+l,...,n.
,...,
Such solutions do not correspond to the ones for a continuum. Let us now consider the case of zero pressure on the upper boundary. into the nth equation in (1.40) yields the relation
k,
(1.53)
Substitution
of (1.45)
271
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
(~1. cosh(k - 1)K - cosh(n - k)q + A - sinh(n - k)q = 0. The determinant
(1.54)
of (1.46) and (1.54) assumes the form
cosh(k - 1)~ -cosh(n - k)q - [sinh
K *
tanh(k - 1)~ - tanh(n - k)q
+&*sinhq]=O.
(1.55)
As in the previous case for S, > 1, the equation has no solution. For the case Sz = 1 (n = 0) the last equation is satisfied but (1.54) leads to the trivial solution. The discussion of the cases - 1s S2 < 1 and S, < - 1 is similar to the previous case and does not give new basic ideas. The last case is one of a gravitational wave on the free surface of the layer. The last equation of (1.40) is transformed into the following one: ql *cosh(k - 1)~ * [sinh 7 - sinh(n - k)v +AThe determinant
[
sinhn*cosh(n-k)q-$7.
-f 7.
(2 + cash 7) - cosh(n - k)g]
(2+coshv)*sinh(n-k)v]=O.
(1.56)
of (1.46) and (1.56) is
cosh(k - { sinh
1)K
K -
* cosh(n
- k)q -
tanh(k - 1)~ * [ sinh 77-$ y
sinhq.tanh(n-k)T-$v.
* (2 +
cash r)) - tanh(n - k)?]
(2+coshq)
11
=O.
(1.57)
Let us consider first the case S2 > 1. In the case K = ul, r] = u1 and S1 = cash ul, S2 = cash ul. All the functions appearing in the last equations are positive. Thus the solution may exist only if
fb.23
(2+
cash ul) - tanh(n - k)u, < sinh v1 < $9
’
-
2+coshu1 tanh(n - k)vl ’
(1.58)
For the case tanh(n - k)ul = 1 it is easy to calculate the root of (1.57). From the equation follows that sinh v1 1 b 2’02 2+cosh ol=‘g’ Simple manipulations
lead to the solution
it
(1.59)
272
K. Szmidt, Finite element approach to steady-state
vibrations in a fluid of finite depth
The corresponding eigenvalue A may be obtained from the relations (1.39). The nontrivial solution is described by (1.45) in which one has to substitute K = u1 and 77= ul. For the case S, = 1 (7 = 0), (1.57) is satisfied, but (1.56) leads to a trivial solution. In the case -1s SZ< 1 the determinant (1.57) assumes the form cos(k - l)u, . cos(n - k)uz . . sin u2. tan(k - l)u2 * sin v~--$ !C-&. g I
[
+ (Y * sin
V2- sin u2. tan(n -
(2+ cos 02) - tan(rt - (2+cosuz)
k)V2 + f y
I]
-
k)U2 I
=o.
(1.61)
It is not easy to calculate the roots of (1.61) and (1.41). The corresponding expressed by (1.51). The last case S2 < -1 results in the following determinant:
solutions
are
cos(k - l)u, - cosh(n - k)ul .
I
sin u2 - tan(k - l)u, * sinh v1 + 3 g” ’ O2 - (2 - cash u,) * tanh(n [
+ (Y - sinh
The corresponding
v1 - sinh v1 * tanh(n -
k)ul
+ i y.
k)vl]
(2 - cash ul)]} = 0 .
(1.62)
solutions are described by (1.53).
2. Detailed solutions for special cases The analysis presented up to now gives an insight into the solutions corresponding to each of the considered boundary conditions on the upper surface of the layer. In the following we are dealing with a special case for which it is possible to assume that tanh(n - k)q = 1. For this case all the determinants (1.48) (1.55) and (1.57) h ave a common multiplier. By equating the multiplier to zero we have the following equation:
sinh
K -
tanh(k - 1)~ + (Y* sinh q = 0 ,
(24
where q = o1 + iu 2, o2 = 2j - T or u2 = (2j - l)~, j = 1,2, . . . . If there exist solutions of (1.41) and (2.1), then the solutions do not depend on the specified boundary conditions on the upper surface of the layer. In other words, the solutions are insensitive to the boundary conditions. Such solutions do not correspond to the continuum and follow the discretization of the problem. According to the previous considerations the last equation leads to nontrivial solutions only if S,<-1
and
S,E
(2.2)
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
273
For the last cases (2.1) assumes the following form: sin u * tan(k - 1)~ + (Y. sinh ZJ= 0,
(2.3)
where u = u2 and v = vl. The cases (2.2) correspond to the solutions (1.53). The last equations of (2.3) must be supplemented by the relations (1.41). It may be seen that for an unknown value of (Y* sinh v > 0 there are k - 1 roots of the equation in the interval (0, n), The roots are in the intervals 2j- 1 2’ 2 - (k - 1) .~
j=l,2
,...,
k-l.
(2.4)
If we continue the calculation of roots for (k - 1). ZA> n, the important values will lead to the same solutions. Thus it is enough to restrict our attention to the interval (0, 71). It does not mean that there are k - 1 solutions of (1.41) and (2.3). According to the second relation in (2.2) it may happen that
(k-l).arccose>s.n,
(2.5)
where s is a natural number. For such a case the first s roots of (2.3) and (1.41) should be crossed out and the number of eigenvectors (1.53) must be reduced to the number k - l- s. According to the considerations it is possible to define bounds for eigenvalues. To make the discussion easier let us neglect the term k2, in (1.9). This term is very small compared to the other terms in the equations. Simple calculations show that for w = 100 set-‘, k: = 4.46 - lo-‘; thus the simplification seems to be justified for our purposes. Substitution of (1.9) into (1.38) and (2.2) gives _l
s 2
=Z.(l+[)+A.(l-25)<1-4s 2-5+/U(1+5)
=2*(1+l)+h*(1-25) 2-5+A*(l+<)
Let us consider three cases: Case 1.2-l>O. For this case the expressions 4+l --<*c-z. 2-5 Case 2. 2-6~0 and 2-5>0. This case follows the solutions
<-1,
1+2&’ +)*,
b, 2 5= (-)a ’ !$=&a[.
(2.6)
(2.6) lead to the following result:
P-7)
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
274
A<-%
or
(2.8)
A>%.
Case 3. 2- 5<0. The solution for this case assumes the form 4+5 -
(2.9)
l-2
According to (1.39) and the above relations it is possible to get a more precise description of intervals for values of S1 and Sz. To illustrate the considerations, an example is presented below. For the case shown in Fig. l(b) let us assume the spacing of points: a = 1 cm, br = 0.4cm, b2 = 2 cm, it = 29, k = 11. It is seen that 5 = 0.16, 6 = 4.0, E = 0.4 and the case corresponds to the relations (2.8). For the parameters, the relations (2.8) give A < -2.26087 or A b4.0. If A < -2.26087, then from (1.39) it follows: - 1 < S, < 0.5862069 and - 1.9411763 < S2< -1.40. Tacking into account the upper bound for S1 = cos u, one can get (k - 1). u = (k - 1) - arccos 0.5862069 = 3.00625 * n. Hence, the number of roots of (2.3) is reduced to the number k - 1 - s = 7. The corresponding eigenvalues Ai are in the intervals: Al E (-9.6802, -5.4469), . . . , AgE (-2.3106, -2.2823), A7E (-2.2661, -2.2609). The solution obtained by a standard numerical procedure (1.18) results in: A1= -7.4048, . .. . , A6= -2.3088, A, = -2.2660. Comparison of the above results shows that the intervals obtained by analytical considerations are in good agreement with the result of numerical computations. Let us consider the case A > 4.0. In this case l-4& 0.5862069 c S, = cos u < -1 + 2E = 0.777777
and
- 1.40 < Sz < - 1.0 .
The bounds result in 2.1634689 * T < (k - 1) - u < 3.00625 - IT. It is easy to see that in this interval there is only one root of (2.3). For the case the first in (1.39) results in the relation 7.2653695 < A < 676.43221. The obtained interval is very wide. The eigenvalue A found by a
Fig. 3. Eigenvectors
corresponding
to equation
(2.1).
K. Szmidt. Finite element approach to steady-state vibrations in a fluid of finite depth
275
standard numerical procedure is equal to 17.0621. Thus, in conclusion, there are 8 eigenvectors corresponding to the roots of (2.3). Some of the corresponding eigenvectors are plotted in Fig. 3. From the graphs it is seen that the components of the eigenvectors are going to zero when one moves along a vertical line in the z-direction. The influence of the solutions on the final results depends of course on the character of a given problem. In cases when higher harmonics are important the influence may be significant.
3. Example of application To get a better insight into the discussed problem let us consider numerical examples. For comparison let us consider two cases of excitements of the fluid motion: V,(z, t) = w * d * cos ot
o - d - cos wt K(z,
t) =
o
forOsz
(3.1)
for 0 s t G hl , forhl
Our aim is to determine the pressure distribution along the wall OA, see Fig. 4. To simplify the solution the zero pressure condition on the upper surface of the layer is assumed. To describe the influence of an assumed irregular net on the final results, the solutions are constructed for two different nets. Both the two nets have an equal number of points in the vertical direction. The first net is a regular one while the second has a denser spacing of points in the lower part of the layer. For comparison the analytical solution of the problem is also given. The obtained results of numerical calculations are plotted in Fig. 4. For a regular distribution of the assumed velocity field at x = 0 (first of (3.1)) the discrete solutions for both regular and irregular nets are in very good agreement with the analytical results. For the
a)
b)
T*P=O WC 24
f’
hs4orm,
-
Fig. 4. Pressure
ANALYTICAL
distribution
SOLUTIONS
along the wall CA.
276
K. Szmidt, Finite element approach to steady-state vibrations in a fluid of finite depth
second of (3.1) the numerical solutions are different from the analytical one. The difference is about 15% with reference to the analytical solution. It should be stressed that the solution for the regular spacing of points is better than in the case of the nonhomogeneous spacing of points in the layer. From the graphs it is easy to see how important for obtaining a reasonable solution is a proper spacing of points in the layer. In cases when higher harmonics are important (second of (3.1)) the better way of discrete solution seems to be applying a homogeneous spacing of points instead of a net with a denser spacing in one part of the layer.
4. Conclusions For the three cases of boundary conditions on the upper surface of the layer considered in this paper the analytical solutions obtained by separation of variables in Cartesian coordinates exist and are well known [4]. If the continuum is replaced by a discrete space of chosen points an analogous solution may be obtained. In the paper the discrete solution is constructed with the help of the finite element method in cases of nonhomogeneous spacing of points in the layer. The discrete solution leads to velocity profiles which do not appear in the solution of the differential equation of the problem. Thus the introduction of a nonhomogeneous spacing of points introduces a nonhomogenity into the field. One should be careful in evaluating results of the finite element method because the discrete system may have properties different from the continuous field. The numerical model gives good results when one considers typical distributions of velocities in the layer. In cases when higher harmonics are important the influence of discretization of the problem on final results may be significant if the number of points is small.
Acknowledgment The author wishes to thank Professor P. Wilde for his valuable remarks during preparation of the work.
References [l] R.E. Bishop, G.M. Gladwell and S. Michaelson, Matrix Analysis of Vibrations (Wyd. Nauk. Techn., Warszawa, 1972) (in Polish). [2] E. Kreyszig, Advanced Engineering Mathematics (Wiley, New York, 1962). [3] P. Wilde and K. Szmidt, F.E. analysis of waves in a semi-infinite layer of fluid, Proc. 4th Conf. Finite Element Methods in Water Resources, Hannover, 1982. [4] P. Wilde and K. Szmidt, Numerical analysis of waves in a semi-infinite layer of fluid, Comput. Meths. Appl. Mech. Engrg. 36(l) (1983) l-23.