TO THE PREDICTION OF ROLL DAMPING OF A TWO-DIMENSIONAL FLOATING BODY
by
ARNE BRAATHEN
DIVISION OF MARINE HYDRODYNAMICS THE NORWEGIAN INSTITUTE OF TECHNOLOGY THE UNIVERSITY OF TRONDHEIM
CONTENTS
PAGE
Abstract 4
Nomenclature 5
INTRODUCTION 7
PREDICTION OF ROLL AND THE EQUATION OF MOTION 9
TWODIMENSIONAL HYDRODYNAMIC THEORY -RIGID FREE SURFACE CASE
3.1 Theoretical formulation 11
3.2 Method of solutIon 12
3.3 FIgures 15
NUMERICAL ANALYSIS AND TEST RESULTS -RIGID FREE SURFACE CASE
4.1 General about the numerical method 16
4.2 Body modeling 19
4.3 Starting up the time simulation 21
4.4 Shedding of vorticity 22
4.5 Calculation of forces and moments 23
4.6 Test of the computer program when no wake is present 24
4.7 Figures 27
5.1 Characteristics of the vortex shedding process 34
5.2 Dumping 36
5.3 Cutting of the shear layer 41
5.4 Additional numerical problems 44
5.5 FIgures 46
RESULTS - RIGID FREE SURFACE CASE
6.1 Calculation of hydrodynamlc coefficients 63
6.2 Discussion of results 68
6.3 Tables and figures 74
TWODIMENSIONAL HYDRODYNAMIC THEORY -FREE SURFACE CASE
7.1 Theoretical formulation 84
7.2 Method of solutIon 84
7.3 Figures 94
NUMERICAL ANALYSIS AND TEST RESULTS -FREE SURFACE CASE
8.1 Calculation of the free surface profile 95 8.2 Test of the numerical method In the heave case 95
8.3 Test of the computer program in the roll case 98 when no wake is present
8.4 Tables and figures 101
RESULTS
-3-Comparison between results from the rigid free surface model, the free surface model, model tests and analytical methods
9.2 General case 119
Comparison between results from the free surface model and model tests
9.3 Tables and figures 122
10. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER WORK 135
Acknowledgement 137 References 138 Appendix A 141 Appendix B 144 Appendix C 147 AppendixD 154 AppendixE 157
ABSTRACT
A vortex tracking method is used to study the development of thin free shear layers emanating from the sharp corners of a two-dimensional floating body perform-ing forced harmonic roll motion.
The free shear layers carry all vorticity present in the fluid.
Under the assumption that the generation of vorticity, and the development of the free shear layers, is not affected by the wave making due to the body motion, the viscous part of the damping may be calculated separately.
The free surface is treated as a rigid wall in this case.
We found this simple numerical model to be unsatisfactory, except for very low and high roll periods, and in the general case it is necessary to include the free surface in the numerical model.
It is shown that the free surface waves have a noticeable effect on the vortex generation and a strong effect on the roll moment, but the generated free surface waves are on the other hand not significantly influenced by the vortex generation.
The time history of the roll moment is found by integrating the pressure obtained from the Bernouilll equation, and the added moment and damping coefficients are cal-culated.
The results are presented graphically whenever possible, in order to give the reader compact information about the physics of the simulations.
Comparison with other methods and experimental results seems to indicate that our method over estimates the roll damping, and possible deficiencies in our method due to the exclusion of certain physical effects are discussed.
Sources of errors due to simplificatlons and modifications in our numerical method are also pointed out.
A added moment coefficient
Av sum of mass moment of inertia and added moment. A see equation (9.1.14)
B roll damping coefficient
equivalent linear roll damping coefficient B see equation (9.1.14)
C force coefficient CR roll restoring moment SE work done pr. roll period
F velocity potential jump at node no. i on the free shearlayer
F force
KC Keulegan Carpenter number, KC = z2 T
L length scale associated with infinite edge vortex moment about bottom of cross section M vortex moment about a given point P N exciting roll moment
Q source strength
R radius of Curvature for the corner of the cross section Remainder function
S wetted body surface
S control surface at infinity
S small cylindrical control surface around the field point control surface far down in the fluid
control surface consisting of the free surface
control surface including both sides of the Riemanri Cuts Control surface including both sides of the free shearlayers
T roll period
convection velocity for velocity potential jump
V, normal velocity
velocity potential at body element i
Z complex vortex force a radius of circular cylinder
see equation (9.1.9) see equation (9.1.2) see equation (9. 1.11) see equation (9.1.13) see equation (9.1.15)
b (t ) parameter specifying the part of the free surface which is discretized brAN vortex roll damping coefficient calculated by Tanaka's procedure by vortex roll damping coefficient
ds element of surface area
f
empirical function used for calculating the force coefficient C g acceleration due to gravityh draft of cross section unit vector along the x-axis
j
unit vector along the y-axis k unit vector normal to the x-y planek wave number
beam of cross section
unit normal vector to a surface
p pressure
s1 s-coordinate of node no. i on the shearlayer,
measured in a curvelinear coordinate system along the shearlayer element length of shear layer element no. j
time
maximum allowable simulation time in a free surface simulation & time step used in simulation
velocity in transformed barge plane
XE x-coordinate of a fluid particle z complex variable, z = x + iy Zp position of roll center in complex plane
see equation (9.1.14)
12 non dimensional frequency
argument of vortex force 8 internal edge angle
free surface elevation
0 roll angle
roll amplitude
9B see equation (9.1.14)
edge parameter, X = 2
-p density of water the velocity potential
-
) velocity pot. jump across free shearlayer or Riemann cut angle between a given shear layer element and the x-axis circular frequency of oscillation1. INTRODUCTION
In the offshore industry an Important problem is associated with the roll mol ion of large barges used for transporting jackets and other cargoes.
The determination of the roll motion damping of such barges is very important with respect to the problem of cargo fastening and risk assessment of the towing operations. At present the barge motion In waves is computed by combining potential flow theory with empirical coefficients for the viscous part of the damping.
The barge damping is caused by two main energy dissipating effects: Damping due to outgoing waves.
Damping due to vortex shedding.
In addition, water-structure interaction in the splash zone may be an important non linear effect, and have strong influence on the roll motion.
Skin friction should also be taken into account in model scale.
The problem which we have studied is the problem of determining the damping due to vortex shedding, which plays an important part in determining the barge response, especially at resonance. Such barges usually have a rectangular cross section with relative sharp corners or bilges witch determine the separation points.
Traditionally, due to the difficulties in calculating viscous effects numerically. damping coefficients have been obtained from model tests.
The difficulties are twofold : First, it is difficult to formulate a mathematical model which represents the physics in a proper way. Second, the numerical formulation of the mathematical model Is usually very complicated. Computer programs based on such numerical formulations therefore tend to be very time consuming.
In this thesis we have compared our results with model test results obtained by Tanaka (1960) and Vugts (1968).
For low frequency roll we have also compared our results with the results obtained by Graham (1980), who calculated the damping due to vortex shedding ntmerica1ly by using the "Matched isolated edge technique".
In the present thesis the roll damping has been studied by using the vortex track-ing method by Faltinsen and Pettersen (1982).
A fundamental assumption of the method is that vorticity is concentrated in thin free shear layers.
Since we are studying forced roll of a rectangular cross section wIth sharp corners the separation points will be fixed and no boundary layer calculation is needed.
-7-The theory is based on a time step integration method where we in each time step have to solve a potential flow boundary value problem outside the thin free shear layers. Two different cases are treated:
Rigid free surface.
Free surface waves generated by the body are allowed to propagate.
When the free surface is rigid, an image-body moving with opposite phase and with the same amplitude as the real body has to be introduced. At each time step we have to solve a Fredholms integral equation of the second kind.
By dividing both the body surface and the free shear layers into line-elements with dipole and source distribution over the body and dipole distribution over the free shear layers, the integral equation may be represented by a linear equation system where the unknowns are the fluid velocity potentials at the body element midpoints.
In the case of a moving free surface we also have to represent the part of the free surface close to the body by linear elements. In the far field the influence of the body is represented by a horizontal dipole satisfying the linearized time dependent free sur-face condition. The problem Is treated as an Initial value problem with the velocity potential on the free surface given at each time step, and the fluid particle velocity on the free surface as the additional unknown variable In our problem.
The equivalence to the Kutta condition is that we require the potential jump at the separation points to be continuous, and that vorticity is shed parallel to the body sur-face on one of the sides of the free shear layers at the separation points.
In order to start up the time simulation a discrete vortex with a given position and cir-culation Is introduced into the fluid in the first time step. The position and strength of this vortex Is calculated on the basis of a pure potential theory calculation as described
by Rott (1956) , and PullIn (1978).
By using the condition of continuous pressure across the free shear layers in combina-tion with the BernoulIll equacombina-tion we get an expression for the rate of shedding of vorti-city into the fluid. The convection velovorti-city of the potential jump over the free shear layer elements Is also calculated. We thus have sufficient information to step the solu-tion forward in time.
A computer program has been developed based on the method described above. The basic input to the program is: Roll amplitude, roll center, roll period and geometry of the cross section, together with parameters describing the free surface model. The basic output from the program may be divided into three categories:
Information about several important parameters at each time step. Time histories.
-9-2. PREDKON OF ROLL AND THE EQUATION OF MOTION
Traditionally the problem of ship motion has been formulated in terms of linear potential theory.
For practical engineering purposes some form of the strip theory, originally proposed by Korvin-Kroukovsky and Jacobs (1957), has often been used. This method reduces the three-dimensional problem to one of two dimensions Involving the sway, heave and roll modes of motion only.
Experience has shown that the more lightly damped motions such as roll, especially near resonance, will be overpredicted by this method. This effect is exaggerated if the ship has appendages or sharp-edged keels.
One reason for this discrepancy between theoretical results and measurements is due to the fact that unmodified potential theory does not include viscous effects
For barges and ships the natural period in roll is usually in the range from 6 to 15 sec. This means that roll resonance problems are likely to occur in sea states where marine operations take place.
Experimental evidence suggests the existence of a viscous damping moment propor-tional to the square of both the frequency and amplitude of roll. (Himeno. 1981)
A slmpl.iñed model of the roll-motion, corresponding to a single degree of freedom system may be written.
Av O+B+CR
e=N(r)
(2.1)Av is the sum of the moment of Inertia and the added moment.
& is the roll angle, CR Is the coefficient of restoring moment,N(<or) is the moment due to periodic external exciting forces with circular frequency <o, and r is the time. B (e) Is some function of the angular velocity representing the damping of the system. One possible function is.
B(9)=B19+B20 I
I +B3e3 (2.2)The coefficients B1 , B2 and B3 may be found by fitting experimental data to the equa-tion.
However, very often the damping is expressed as some form of linear damping.
B(e)=Bee
(2.3)The linear damping coefficient B will depend on both the frequency and the roll-amplitude, but is otherwise considered to be a constant. The coefficient may be found from a forced roll test and expressed in terms of the amplitude independent coefficients
If the method of equivalent linearization is used we get.
B (2.4)
It is often seen that the last coefficientB3 is ignored. This Is the case in the method used by Salvesen et. al. (1970). Here a linear damping coefficientB Is Introduced into the equations of motion based on an estimated value of the response amplitude . If
the difference between calculated and estimated response is too big, the damping
coefficient is re-calculated and the process repeated until acceptable results are
obtained.
In the equations of motion is is customary to neglect coupling terms for the viscous damping. This means that the viscous damping due to the roll excitation
moment Is associated with the roll motion alone.
Aarsnes (1986) has studied roll motion, using the matched isolated edge technique developed by Graham (1980). This method is described in chapter 9.1 in the present thesis. In addition to the incoming wave and the body roll motion, Aarsnes uses the sway and heave motions in order to find the velocity around the shedding edge and hence the total viscous moment.
His conclusion was that for small beam to draft ratios, the results from simulations using the total relative velocity to predict the viscous moment differed significantly from the results obtained from simulations where only the velocity contribution from the body motion in roll was used. This means that viscous coupling terms have to be Included In the equations of motion.
We therefore see that the problem of finding the proper response of a structure, even in the case where the structure Is excited by regular waves, is complicated. Even if the hydrodynamic coefficients acsoclated with the pure heave, roll and sway motions are found, we still may not have sufficient information for solving the problem.
The final solution for the roll problem would be to perform a full time simulation with all relevant physical effects included, meaning that the development of the shear layers have to be simulated as well during the simulations.
This has however, for the time being, only theoretical interest due to the computer costs involved In such a calculation.
-11-3. TWODIMENSIONAL HYDRODYNAMIC THEORY -RIGID FREE SU(FACE CASE
3.1 Theoretical formulation.
We consider forced roll motion of the rectangular cross-section shown in fig. 3.1 In order to Study the effect of viscous damping separately we adopt a model where the free surface is treated as a rigid wall. The only form of energy dissipation, and conse-quently damping, is in this case due to the generation of vorticity at the two sharp corners of the cross section.
The roll motion Is harmonic, defined by
0= ësln(2f..t)
(3.1.1)Here 0 is the roll angle, is the roll amplitude, T is the roll period and t is the time. The fluid domain Is infinite in extent and of infinite depth.
At the two sharp corners of the body vortex shedding will take place and two shear layers, supposed to be infinitely thin, will develop.
In the fluid domain, outside the free shear layers, the fluid is irrotational and there exists a velocity potential satisfying the Laplace equation.
12th 12th
+=O
(3.1.2)1x2 1y2
When the free surface is treated as a rigid wall, the boundary condition on the free surface is expressed in terms of the normal derivative of the velocity potential as
In (3.1.3)
The boundary condition on the body surface is
In (3.1.4)
where V, is the normal velocity at a typical point p on the body surface.
V
=(xF)
(3.1.5)is the unit normal vector at the point p, defined to be positive into the fluid, k is a unit vector normal to the x-y plane. and is the vector from the center of rotation to the point p.
3.2 Method of solution
The roll motion problem is solved in the rigid free surface case" by introducing an lmage-body" and an 'image-shear layer as shown In fig. 3.1
Let the contour S be defined by
S'=S+S+Sv+SR+Si
(3.2.1)Here S is the wetted body surface, S is a control surface at infinity, S is both sides of the free shear layers, SR is both sides of the Riemann cuts, and S1 is a small cylindrical surface with center at the point (x1,y )in the fluid domain.
By applying Green's theorem to the velocity potential and the function !i defined by
we get
f(.)ds(x,y)=O
. (3.2.3)Both the 'Image-body' , 'Image shear layer' and 'image Riemann cuts" are included in the intergratlons. The contribution to the integral from S,, is zero, and we have
2ir (x1.y1)
= (3.2.4)I
fcry
5+ ,,+S
a ln/(x x1 )2+(y Yi )2
an (x ,y)
(x x
1)2+(y y )2 ôx .y) jds (x ,y) 8n (x ,y)Hereds isa surface element for any one of the surfaces, and i' is a unit vector normal
tods. The positive direction of ' is shown in fig. 3.1 The total circulation around 5+3,, +SR has to be zero.
The concept of a Riemann cut is introduced in order to connect different shear layers and/or discrete vortices emanating from the same separation point.
The Rlemann cut is mathematical and not physical. Across the cut the velocity poten-tial is discontinuous , but all physical quantities are continuous. The potenpoten-tial jump across the cut is constant along the cut.
Across the shear layers the tangential velocity is discontinuous and the normal velo-city is continuous.
Equation (3.2.4) may therefore be rewritten
13
-2r
(x1y1)= (3.2.5)f
!4i(x )a ln-v'(xan(x,y)
1n(xx i)2+(yyi)2
V Ids(x ,y)+
f
I(x,y)
(x,y)
Ia(x_xi)l+(y_yi)2
ds(x,y)
s+v+s+zt J
n(x,y)
Here + and - refer to the two sides of the free shear layer or a Riemann cut.
V and
are the source and dipole strength distributions over the body and (*__) is the dipole strength distribution over the shear layers and Riemann cuts.It is possible to show that an alternative expression, where a vortex distribution is used instead of a dipole distribution, may be derived.
Faltinsen and Pettersen points out that the potential Is expected to be less singular than the velocity at possible singular points. The formulation of the problem in terms of a vortex distribution implies that we have to solve for the velocity instead of the poten-tial. Therefore the numerical problems are expected to be greater with this method than with ours.
Equation (3.2.5) is transformed into a linear algebraic equation and solved for each time step. The body surface and the free shear layers are approximated by straight line elements with Constant source and dipole strength over the body, and linear dipole strength over the free shear layers. The body boundary condition is satisfied at the ele-ment midpoints.
A continuity condition on (' -
) is satisfied at the separation points by requiring that the jump in the velocity potential at the separation points along S is equal to the the jump across the free shear layers at the same separation point. This condition has to be satisfied at each time stepBy using the Bernouilli equation and requiring that the pressure is continuous across the free shear layer we obtain
2 [ax1 ax1 J lx
[+--)
(3.2.6)1 1a
ay Jay1
Physically this equation express the rate of shedding of vorticity into the fluid at the separation point, if the separation point is fixed.
If we integrate the equation we obtain the velocity potential jump at the scparauon points and at the free shear layers for each time step.
We observe from the equation that the velocity potential jump
-
J across the free shear layers Is convected with the velocity1+
ox1 Ox1
T+L
2layi
ay1J I. (3.2.7)Here i and j are unit vectors in the x and y direction respectively.
In order to be able to use equation (3.2.6) we must calculate the derivatives of the potentials involved. This is done by differentiating equation 3.2.5 giving
air
1Ox1 2ir
2(xi
x)ds(x,y)
fd
0
(Xl
_X)ds(xy)
2ir On r2 Andx
_L
f
(+_)8(X1
2)ds(x,y)
r 2)ds(x,y)
-..=1V
.i(y1_y)ds(x.y)...Lf 0(YIV
r 21r On r __.L.(f)---(31'31)ds(xy)
On4 r2 (3.2.8.a) (3.2.8.b)3.3 Figures IM.C.E- WAKE -15
-LUD DOMAN
IMAGE WAKER(C.ID FE
SUPACE D(SCETE VOTE)cFig. 3.1 Flow situation around a two dimensional cross section performing forced roll motion.
4. NUMERICAL ANALYSIS AND TEST RESULTh -RIGID FREE SURFACE CASE.
4.1 General about the numerical method
The general approach described by Faltinsen and Pettersen is followed in the numerical implementation of the method of solution.
The procedure is valid for any number of separation points, but in our case we will always have two, one at each corner of the rectangular cross section.
Let us simplify the description by assuming no detached free shear layers.
The body surface and the free shear layers are approximated by straight line elements. On each body element we write
cb=X1 £ = 1 ,N (4.1.1)
X1 , I = l,N are unknown constants.
On each free shear layer element the velocity potential jump over the respective ele-ment is written as a linear function,
(4.1.2)
S1 S1
where s Is a curvilinear coordinate system along the free shear layer,s1 ands1 are
the s - coordinates at the ends of segment no. I . andF1 is the velocity potential jump
ats =s1.
For each time stepF1are known quantities. The body-wake model is shown in fig. 4.1
By letting (x1
'I
approach the midpoints of each segment of the body we obtain a discretized version of equation (3.2.5).N
+ X1 A1 = B1 + (4.l.3.a)
A.1 =
f±li.
ds (4.l.3.b) B1= 2r1=1
j lnr ds
(4.l.3.c)The integrations are over body (and image body") surface element i, and the field. point is situated at the midpoint of element j
1 M,,
F1F,1
a D--
+£ s - s
mr ds+F1 r_.Linr dsl
+The integrations are over shear layer element no. i. M5 is the no. of wake elcmcnts in free shear layer no. k.
The integrals in (4.1.3) are calculated by Faltinsen and Pettersen (1982).
Due to the continuity condition on the potential jump at the two separation points. we are not able to satisfy equation (4.1.3-a) at more than N-2 midpoints.
Following Faltinsen and Pettersen, we choose to replace the equations associated with two of the elements next to the separation points ( see fig. 4.1) wIth equations express-ing the continuity conditions.
We then have. XNp(j)+1 {
1
- 17-SN.-p(1)+I 1 SVSE?U)+1 + SNp(1)+2SNp(1)+1 I SNSEp(l)+2 -- SN.S.Ep(I) + SNp(i)_1 - 5NSEP(i) = (-
I SVSEP(2)+l XNp(2)+j (1
-
} + { SNp(2)+2 - 5NSEP(2)+1 3NSEP(2)+2 -- 5NSEP(2) SVSEP(2) +XNp(z)_j } + { I + S.VL'J)_ - SN.p(2) Syp(2)_1 - SVSEP(2) = (-
)spJomwo.2 (4. l.4.b)Similar expressions are valid for the "Image body".
The two equations above are based on linear extrapolation of the potential from the body to the separation points. See fig. 4.2 for explanation of symbols.
When the unknows X have been found from the linear equation system, the the solu-tion will be stepped forward in time.
From equations (3.2.7) and (3.2.8) a discretized version of the convection velocity of the free shear layer at the midpoint of wake element j is obtained.
Tc M
IFi F11
((
_sj+i).(X1 _X)()
_"J I
ss1
an r2 (4. l.3.d) SVSEp(l) (4.1 .4.a)+F+lLJ_C12X)ds(x,y)
sim1c(F +F1+1) 17 I N N (Yl)ds(x,y)
+ I ----
2ir.1
xL:(YJ
r2 __.L. {F,_F+1L( s+i)
__Yi y
217 - s +1 + r2 )ds (x y)+F1+lxL__I_( r2
)ds(x.y)
i cosiix(F1 +F+1)M is the total number of free shear layer elements, jt is the angle between element no. J and the x-axis and &. Is the length of element no. j.
The two terms above, containing the sin and cos expressions, emanates from the integrated contribution to the velocity over the element where the field point (x1 .y1) Itself Is situated.
When Vç has been calculated, the Convection velocities at the endpoints of the free shear layer elements are found by linear interpolation. The convection velocities of the element endpoints coinciding with the separation points are extrapolated from veloci-ties on the body.
If the shear layer has a free end, the convection velocity of the free end is found by extrapolating from velocities on the wake.
If the shear layer has a discrete vortex at the end, this vortex will not induce any velo-city on Itself, and It's velovelo-city Is found In a separate calculation.
The endpoints of the free shear layers are now stepped forward in time by using the Euler method.
Faltinsen and Pettersen has reported that thc simple method proved to be stable, giv-ing the same result as the more advanced fourth order Runge Kutta method.
The velocity potential jump at the element endpoints is not changed during the convec-tion, as seen from equation (3.2.7).
One problem Is that the convection velocities are different from element to element. Consequently the lengths of the elements will be quite different after a while, leading to an inadequate numerical model of the wake.
This phenomenon Is described by Fink and Soh (1974).
In order to avoid this problem, the wake is rediscretized at each time step ut such a way that all elements except the ones closest to the separation points are of equal length. The potential jumps at the new element ends are found by linear extrapolation.
-
19-The potential jumps at the separation points themselves are found by using equation (3.2.6).
4.2 Body modeling.
In ordinary potential theory (no generation of vorticity at the sharp corners) and for small roll angles it Is customary to solve the roll-problem by linearizing the problem. This means that we can calculate all variables at the average position of the body.
On the other hand, when we have separation at the corners this approach will lead to wrong results.
From fig. 4.3 we see that the motion of the corner is significant compared to the scale of the generated shear layer. The effect of the corner actually moving away from the generated shear layer elements has to be represented In a realistic way.
When we tried to use the results from a calculation based on the average position of the body the solution got unstable after a few time steps.
The conclusion is that we have to model a body which is actually moving. Therefore both the coefficient matrix and the right hand side of the linear equation system has to be updated at each time step.
In a linearized problem it is sufficient to calculate the coefficient matrix in the first time step, and solve the linear equation system at all the other time steps by multiplying the right hand side by the inverse of the coefficient matrix. This is less time consuming than actually setting up the coefficient matrix, and solving the equation system for all time steps.
In our computer program however, we found that most of the computing time was spent on calculations concerning the wake, and consequently (except when we had very few wake elements) It made no big difference in which way the equation system was solved.
We tried two different ways of modeling the (moving) body. The first example is shown in fig. 4.4
The beam of the cross section is lOm and the draft is 5m. The roll period is 4s and the roll amplitude is 5 degrees.
At the start of the simulation, the first and last node is situated at the free surface. The position of the nodes relative to the body is fixed, and after some time steps we see
that the body model has changed position.
If we used this model, a sudden "jump In the roll moment time history was observed at a certain time step (Fig. 4.5)
A11 =
All physical quantities are given in the S.I. system, and the parameter RO is the density of water.
At the following time steps the time history went back to normal. This unexpected and unphysical behaviour was found to have the following explanation:
The gap between the physical body and the "image body" at the left will only introduce a small error in the moment calculation, but at the critical time step we observed that we had a situation where the control point at one of the midpoints of the
body-elements was situated very near the respective image element.
The terms on the main diagonal of the coefficient matrix in our equation system has the form 1 +
r
ln-.J(x xj )2+(y )2 ds (x .y)n(x,y)
+ j IMAGEalnv'(x x1 )2+(y -y1 )
ds (x ,y)
n(x.y)
The integration over the basic element gives -0.5 whereas the integration over the image element "usually" gives a small contribution, say of order e. Therefore the "usual" situation Is that the main diagonal is of order one. All the off diagonal elements are of order .
In the special situation described above however, the integration over the image ele-ment gives a contribution approximately equal to -0.5 . which means that the
corresponding diagonal element is of order just like the off diagonal elements. All the other diagonal elements are still of order one.
The routine for solving the linear equation system was not able to handle this situa-tion, and a change to another routine did not solve the problem.
The problem could have been avoided by changing the model in such a way that the critical situation would not appear, but in order to make the program general we
decided to change the method of modeling. The new model is shown in fig. 4.6
The first and last nodes are situated at a distance I above the free surface. The distance I has to be so great that the body model will remain "wet" for the maximum rotation angle.
The program will, at all time steps, compute the intersection between the body and the free surface, and place the first and last nodes at the intersection points. We could have chosen the first and last body elements in such a way that we only had to adjust the "dry" node of these elements, thus keeping the number of body elements constant, only changing the first and last body element at each time step. This would however, make it necessary to have quite long elements at the free surface, especially for large (4.2.1)
21
-rotation angles.
We therefore decided to model the body in such a way that several elements would alternate between dry" and "Wet".
The number of body elements, and consequently the number of equations in our linear equation system, will therefore change as the body is rolling.
Tests showed that this did not cause any numerical problems in the simulations, but considerable book-keeping had to be incorporated In the computer program in order to keep track of the separation points etc.
We see (Fig. 4.7) that the image body will now never overlap the real body, and the numerical problem that we had with the first model will not occur.
4.3 Starting up the time simulation
In the first time step, before any vorticity is generated at the sharp corners, the fluid particle velocities at the corners will be unrealisticly high because they arc calcu-lated on the basis of pure potential theory, which predicts infinite velocities at the corners.
In order to obtain a more realistic initial condition, and consequently a more realistic solution as soon as possible, we started up the time simulation by positioning a discrete
vortex with known position and strength just outside the sharp corner, as proposed by Faltinsen and Pettersen. (Fig. 4.8) This starting vortex" simulates the situation at a time t which is small compared with the time scale for which the development of the vortex sheet at the tip is strongly iafluenced by other parts of the sheet. The basis for the calculation Is the Kutta condition, saying that the velocity at the sharp corner
should be finite.
The discrete vortex solution was obtained by Rott (1956)
In our case we first have to do a separate computer run for just one time step in order to determine the potential theory solution. Then this solution is matched with the analytical solution (which conisins an unknown constant) for the flow around an infinite wedge. The formulas given by Pullin (1978) for the strength and position of the vortex can then be used directly.
The time simulation is now started up with one discrete vortex situated near each of the corners.
The discrete vortices are represented in the following way: Two shear layer elements with a dipole distribution of constant strength are created with one end attached to the respective corners and the other ends at the position of the discrete vortexes.
Several computer runs were made with different strength of the vortices in order to test how the generation of different starting vortices, and hence different starting con-ditions, did affect the solution.
When we have a weak starting vortex the wake will not have as strong tendency to roll up as when the starting vortex Is stronger. This Is illustrated in fig. 4.9 and 4.10. In the case of a weak starting vortex, the roll up is due to the circulation of the wake elements alone ("self-roll up"). The starting vortex has In fact drifted away from the center of roll up.
The tests showed that the force and moment time histories were identical after approx-imately the first quarter of a period even if they where quite different at the start. These tests also included a test with zero strength of the vortices.
In fig. 4.11 we see how the roll moment time histories approach one another in this case. The same is observed for the amount of circulation shed pr. time step. (fig. 4.12) Note that when we have a starting vortex of zero strength this unphysical situation is compensated early In the simulation by the shedding of wake- elements carrying strong circulation. The self roll up around these strong wake-elements is not so different from the roll-up around the "correct" starting vortex.
Since, as will be discussed later, the first half period will not give a realistic representa-Uon of the steady state solution for the roll problem, and therefore will be discarded anyway, our conclusion Is that It Is not necessary to do a careful calculation of the starting vortex.
Even with no starting vortex at all the system will "correct" itself and converge to the correct "common" solution.
4.4 Shedding of vorticity
When we calculate the amount of vorticity shed into the fluid at each time step from the separation point we must take Into consideration that the body, and conse-quently the separation point, Is moving.
Looking at fig. 4.13 we see that the uid particle velocity at one side of the shear layer is equal to the the tangential velocity u calculated on the basIs of the extrapolated velocity potential. On the other side of the shear layer the velocity Is equal to the nor-mal component v of the velocity of the body.
During the time step & the "newborn" shear layer element is generated by convection a distance where
The shear layer element has also Increased In length by AL 2 because of the motion of the body.
Al2v
The total length AL of the element is
JAL1+AL2=42-(uv)At
The amount of circulation Ay represented by this element is
A')' = - (u - v) AL
23
-(4.4.2)
(4.4.3)
(4.4.4) The minus sign is due to the fact that the positive direction of the circulation is counter
clockwise. The amount of circulation generated in each time step therefore is,
=_..L(
- v
)2which Is minus half the relative velocity squared.
The tangential shedding direction of vorticity from the body corner is determined by observing where the relative velocity between body and fluid is the greater. When the shedding direction changes, the sign of the generated vorticity in equation (4.4.5) will also change.
4.5 CalculatIon of forces and moments
The force F and the moment N Is calculated by integrating the pressure p obtained from the Bernoullli equation. (Hydrostatic term excluded)
p =-j.-f v'2
=p ñdS
(4.5.1)
(4.5.2)
(4.5.3) p is the density of water.
Because the body is moving, the potential is calculated at a different point in space for each time step.
In order to calculate the term jP. at the position ' at the intermediate time step
at & (4.5.4) We will point out that this calculation of the force in very sensitive to numerical errors in the potential i, since the numerator s a small quantity obtained by subtract-ing two large values.
This is clearly demonstrated in calculations where the simulation is stopped at a given time step In order to perform some approximation, and then restarted. (Further expla-nation Is given in chapters 5.2 and .5.3)
A benefit from this method of calculating the force is that instabilities in the simula-tion are easily discovered.
The terms 1'(
'
,) and t(
' . t ) are not explicit output variables of thecom-puter program. We therefore perform the Taylor expansions.
(4.5.5.a) (4.5.5.b) The terms ( r2 , r2) and '( t1) are known quantities.
It is consistent to replace the terms Vit (,2.t2) and 7c with the term V calculated by taking the mean value of the fluid velocities at t and t2, and we finally
get.
±.55:
ti(2,r2)-(1,t1)
V -atV is the body velocity at the respective point.
-
r2-r1
V&
An alternative expression for the force, calculated by Faltinsen and Pettersen based on Newman's derivation (1977) for the force on a body in non- separated flow is.
FA-P- J
dt s+ ntdS (4.5.8)+s
In the present case, due to the presence of the free surface, the formula will only make sense for the force in the x-direction. A comparison between the two formulas for the force will give us information about the accuracy of the numerical model.
4.6 Test of the computer program when no wake is present
(4.5.6)
25
-In order to construct a reliable computer program the program was tested in several "special cases" where the answers where known from other sources (Theory or computer results from other programs)
The first test consisted in setting the roll motion to zero and allow a uniform stream to come In from the left.
As a help in debugging the program, and in order to get a good physical understanding of the processes Involved in the simulation, a considerable effort was spent in writing graphical routines for several of the output parameters.
A detailed description of the graphical post processing routines is given in appendix C.2
and EltoE.3
The test described above showed that the velocity potential over the body was anti-symmetric, arid that the tangential velocity was symmetric. (Fig. 4.14) The maximum of the tangential velocity was observed at the corners, where the analytical solution based on potential theory will give infinite velocity.
In the x-dlrectlon the force was zero, as expected from D'Alambert's paradox.
The force In the the y direction was different from zero and negative. This is also as expected because there Is a negative contribution to the pressure on the bottom of the body coming from the velocity squared term in the Bernoullil equation.
In fig. 4.15 the velocity field around the body, based on one of the graphical routines, is shown.
When the velocity flow field Is studied it has to be pointed out that velocities calcu-lated close to the body (less than an body- element length) will not give correct results except close to the body element midpoints and nodes. Fig. 4.16
This is due to the fact that the body boundary condition is satisfied at the element mid-points.
We see from the same figure that the net flux over the element is zero.
The second test that we performed was to simulate one period of the roll motion with no generation of vorticity at the corners.
In this case the external flow was set to zero.
In the case of a rigid free surface no wave generation will take place. Hence no energy dissipating mechsnism will be present, and the damping will be zero.
The roll moment time history is shown in fig. 4.17
We observe that the curve Is in phase with the angular acceleration of the body. The damping coefficient is zero, and the added moment coefficient is 0.061 which is close to the added moment coefficient given by Faltincen (1969) for zero angular velo-city. (The added moment coefficient is non dimensionalized by dividing with p /i., where I is the beam of the cross section and h is the draft).
Observation of the velocity potential and the tangential velocity at the first time step showed that they were antisymmetric and symmetric respectively.
Fluid particle velocities in the neighbourhood of one of the corners are shown in fig. 4.18
The velocities at some distance from the body were approximately zero.
We concluded that the computer program so far gave reasonable results when limited to pure potential theory calculations in the rigid free surface case.
-S rsp(i-' -SN5.P(' 5r.SP(H-2 SM5EP(')-1
Fig 4.2 Coordinates used in specifying the continuity condition on the potential jump at the separation point.
Linear extrapolation is used.
Fig. 4.1 Scheme for body/wake model - enlarged view of body close to separation pomt no. 1 included.
layer and corner at the start and after a quarter of a period.
Fig. 4.4 Position of physical body and image body at the start of a simulation and at a typical timestep. Beam - lOm draft - 5m.
m
Fig 4.6 Physical body - new model.
ROLL MOMENT: PM
29
-TIME (SEC.)
Fig 4.5 Unphysical behaviour of the roll momentat four 'critical timestcps.
Fig 4.7 Relative position of physical body and image body for the new model. r'Jo overlap occurs. Beam - lOin draft - 5m.
The starting vortex is shown at the right corner of the cross section together with one generated shear layer element.
Fig 4.9 Situation after a quarter of a period.
A starting vortex with a strength and position as prescribed by Pullin has been used.
Fig 4.10 Situation after a quarter of a period. A very weak starting vortex has been used.
m5 S2 4 4-4-. 431
-Fig 4.12 Circulation shed pr timestep, with and without Rott starting vortex. -I ME
Fig 4.11 Roll moment timehistiories with Rott starting vortex, and with starting vortex of zero strength. TIME
Ki)
DiQ.ECflON OF POi1'IUEC C
UI__A-rioNm2 S m S z 0 a-C)
H
0 z Co z 40 60 60 CO BOD aE4 40 60 80 CCC BOFig. 4.14 Velocity pot. and tangential velocity distribution around a cross section sub-ject to a uniform stream.
a m5 STREAM DIRECTION 33 -BODY ELEMENT
Fig. 4.16 FluId particle velocities very closeto a typical body element. a
a-TIME SEC.) Fig. 4.17 Roll moment timehistory- pure potential theory.
0 F F
5. TIME DEVELOPMENT OF THE FREE SHEAR LAYERS.
5.1 Characteristics of the vortex shedding process
The use of graphical representation of the results are essential both for the under-standing of the physics involved In the simulation and for the technical accomplish-ment of the project, especially with respect to the debugging of the computer program. Important information, such as development of the wake, fluid particle velocities, and (later) generation of surface waves is given an 'analog presentation which helps the user to understand what parameters are the important ones in the problem.
A number of curves displaying time histories, circulation-intensities, velocity- distri-bution along the body etc., are also included in the graphical post-processing program.
In order to explain the physical processes involved, and also some of the numerical dlculties that occurred, a series of plots showing the body and wake at different time steps is shown in fig. 5.1 a - p.
The results are given for a cross section with beam 10 m and draft 5 m, rolling with a period of 4 seconds and having a roll amplitude of 5 degrees. The roll center Is In the origin.
Only details at one of the corners of the cross section are shown. The time step used in the simulation was 0.01 5. The motion is started with the angular velocity at its max-imum.
t-0.Ol5:
The discrete starting vortex at the right corner of the cross section is shown. t-0.5 s:
The free shear layer is rolling up around the starting vortex. The shedding-direction" is tangent to the vertical side of the cross section.
tl.0 s
The cross section has reached its maximum deflection angle. Some time before this position was reached, the vorticity being shed from the corners had changed sign, and it is seen that the free shear layer Is now developing in a direction tangent to the bottom of the cross section.
A small bump Is shown at the start of the shear layer. This is due to the fact that we have had a change in the sign of vorticity along the free shear layer close to this bump, causing the shear layer to roll up in an opposite direction as before.
tl.5 s
Here a fully developed Mouble-vortex" is shown. The change in sign of the vorticity distributed along the free shear layer causes this characteristic shape. If we integrate the velocity along a closed curve, starting and ending at the point on the free shear
35
-layer Indicated by an arrow, and Including the double vortex, the contribution to the circulation will be zero.
tl.62 s
In order to keep the computing time down it is necessary to have as few free shear layer elements as possible. Therefore the double vortex is approximated by twc
discrete vortices having opposite signs.
We have Indicated the strength of the vortices by small circles. The vortex with the greatest strength Is represented by a full circle. The opening in the upper circle shows that this vortex Is somewhat weaker than the lower.
t-3.O s
The cross section has reached its maximum deflection angle at the opposite side as in fig. 5.l.c The vortex pair has moved away from the body and a new vortex has developed at the corner.
t-3.Ol s
Part of the new vortex Is represented by a discrete vortex. It Is seen that this vortex is stronger than the vortices In the vortex pair.
t-4.O
Now a full period has been reached. A new double vortex has been formed. t-4.3 s
When the two vortex pairs have moved away from the body we are allowed to cut the wake at the point where the contribution to the circulation from the rest of the wake is zero. The viscous force on the cross section is dominated by the new vortex developing at the corner.
t-5.Ol
The vortex at the corner is represented by a discrete vortex.
t-6.Os
A new double vortex has developed. t-.6.Ol
The wake has been cut due to the same argument as in fig. 5.l.i
t-7.o
A new vortex is developing at the corner. The end of the cut wake is also rolling up. t-7.Ol s
The new vortex, and the vortex at the end of the wake, are represented by discrete vor-tices.
t-7.5 s
t-8.0 s
The wake has been cut due to the same argument as in fig. 5.1.i and 5.1.1.
The results from the simulations show that after approximately half a period we reach a vortex shedding process that repeats Itself.
The periodic shedding of vortex pairs from sharp edged bodies has been demonstrated in experiments (Singh, 1979).
From the numerical simulation we also observe that these vortex pairs move away from the body. This effect Is also confirmed by the same experiments.
The assumption of periodicity will enable us to compute hydrodynamic coefficients based on a very short portion of the moment time history curve.
It is well known, on the other hand, that for two dimensional circular cylinders, and other cross sections where the separation point Is not fixed, many cycles are needed in order to obtain reasonable coefficient estimates.
Before we end the description of the time development of the free shear layers we will point out the symmetry of the two shear layers emanating from the different corners. (Fig. 5.2).
We will later see that this symmetry property Is due o the present rigid free surface model, and is partly lost when we introduce a moving free surface.
The symmetry property could, in the rigid free surface case, have been used explicit in the equations in order to reduce the computing time.
5.2 Dumping
As the vortex sheet rolls up It becomes increasingly difficult to keep the rolled up region of the vortex sheet stable, even if we use the method of rediscretization proposed by Fink and Soh.
Fig. 5.3 a - d , demonstrates how the double vortices begin to break down during the last part of a one period simulation. (The roll period was 4.s. and a time step of 0.01 s. was used)
Instabilities may be due to real physical processes, but does also take place because the numerical method does not represent the problem properly.
Interaction between and merging of cores with vorticity of the same or opposite strength may occur followed by viscous diffusion and decay. In three dimensions com-plicated longitudinal Instabilities may be triggered even before the merging or decay process has started, resulting in a break- up of the vortex system.
37
-Different numerical methods have been introduced in order to handle the problem of numerical instability, such as introducing a finite region of vorticity in which the velo-city remains finite.
An alternative approach Is to use the so-called vortex-in-cell method where the velo-city field is computed by solving Poisson's equation for the stream function due to a
grid dependent region of distributed vorticity.
The main disadvantage of this method is that the definition of the sheet is lost and that a large number of grid dependent small-scale structures appear to which the relevance
to the large scale roll-up process is difficult to asses. The discrete vortex method handles this problem even worse.
it is possible to show (Appendix A) that the velocity field outside a circle with constant vorticity distribution Is equal to the velocity field induced by a point vortex situated at the center of the circle, and with strength equal to the sum of vorticity dis-tributed on the circle.
As an engineering approximation we studied the effect of representing one or more com-plete turns of the sheet by a discrete vortex placed at the center of vorticity of that part of the sheet. This can be done when the spiral Is nearly circular, and if the dipole strength does not vary too much over one turn of the spiral.
This also holds for double branched spiral vortices.
We Introduce the concept of dumping which means that we during the simulation approximate a part of the wake by one or more discrete vortices.
It would also have been allowable to let the sheet roll up by continuously "feeding" the last element of the sheet into the discrete vortex at the center of roll-up.
Several tests were performed in order to get experience with respect to when so called "dumping" was allowed.
This has to be done very carefully in order not to be in conflict with the physical pro-cess being simulated. Only at certain specific time steps dumping will be allowed. There are two main causes for the necessity of dumping.
First we observe that the computing time will be prohibitive if we are to carry outa simulation for more than one roll period without dumping.
Second we have the stability problem discussed above.
An effective way of getting experience with the effect of dumping is to compare two separate computer runs, one which Includes a dump, and another were no dumping Is performed.
The simulation results from two such runs should show no significant difference with respect to moment time history, generation of vorticity, convection of the part of the wake that has not been dumped and so on.
Fig. 5.4 shows the situation after one period from a simulation where several dumps have been performed. (The parameters of the simulation are the same as those given in chapter 5.1)
If we compare with fig. 5.34 (note that different scales have been used in the plots) we see that the two different free shear layers are almost equal to one another close to the body.
The part of the shear layer connecting the two double vortices contains very little vor-ticity, and has therefore negligible lnuence on the physics of the simulation. This weak part of the shear layer may however Introduce numerical problems if it enters the corner of the cross section.
The main difference between the two simulations is seen if we look at the position of the double vortex at the right. This difference does not affect the forces acting on the body or the shedding of vorticity from the body because the double vortex is far away and represents a circulation close to zero.
It Is very difficult to tell a priori at which time steps during the simulation we will be able to dump.
In principle, dumping and other approximations performed during the time simulation could have been built into the computer program as automatic procedures. In this way an inexperienced user would have been able to use the program.
In practice however, It Is very difficult to build such intelligent procedures into the program, and we have to run the simulation program Interactively.
The graphical post processing program Is used, and the development of the wake is examined closely In order to take action at the proper time steps. The distribution of vorticity along the wake Is very important in this respect.
We will demonstrate this by giving some examples:
In fig. 5.5.a we see a double vortex rolling up at the right corner of the cross section at time step 150 (t - 1.5 s.)
In fig. 5.5.b the 4 inner turns of this double vortex is dumped into a discrete vortex. (time step 151)
Fig. 5.5.c shows the situation some time steps later. (time step 160)
In fig. 5.5.d the upper part of the double vortex has also been dumped into a discrete vortex. (time step 161)
Some elements at the end of the wake have been included in the lower discrete vortex. The vorticity in this part of the wake Is very weak, so the merging of these elements into the discrete vortex will not change this vortex appreciably.
We also did the same run (from time step 150 to 161)with no dumping, and observed that the two wakes were practically identical if we looked at the part of the wake that
39
-had not been dumped, and that the discrete vortices were situated approximately at the centers of the (un-dumped) double vortex.
The moment time history was also practically Identical for the two cases and so were the curves showing the generated vorticity for each case.
Let us consider the potential jump across the wake, and the distribution of vot-ti-city over the wake in the situations described above In order to demonstrate the phy-sirs involved.
Before any dumping has taken place we observe three areas of strong vorticity in the wake. (Fig. 5.6).
The first peak is due to the fact that a new positive vortex Is beginning to develop at the corner of the body. The second peak corresponds to the upper (positive) part of the double vortex. From the information at the end of the free shear layer we see that most of the circulation is contained in the discrete vortex represented by the sudden drop at the end.
In order to get fast and efficient information about the intensity of circulation along the shear layer we divided the circulation intensity into six levels and displayed each level in a different colour.
We also let the computer program show which shear layer elements were the strongest and the weakest. The element where the potential jump changes sign was also found. Unfortunately, we did not have the necessary hardware to make plots with six different colours, but we were able to do these important wake intensity observations from our graphical computer terminal.
From fig. 5.5.a, we were then able to see directly that the potential jump changed sign at the end of element 8.
The region of strongest vorticity was red, and element 35, which was the strongest, was marked separately.
A green colour indicated that the contribution to the circulation was negative, and we observed that the first negative element was no. 84.
At element 124 , we had the element with the strongest negative circulation contribu-tion.
Fig. 5.6 and 5.7 should also be consulted in connection with the description above. The dipole strength along the free shear layer is shown In fig. 5.7 At wake element no. 8 the dipole strength is close to zero. This implies that the total circulation represented by the rest of the wake Is approximately zero.
If we look at a close-up of the wake from element 140 to the end (fig. 5.8), we see that the potential jump over the last element Is constant and drops to zero at the end of the wake.
This is in agreement with the fact that the wake end consists of a discrete vortex. When the double vortex is dumped into two discrete vortices the situation is as shown in figs. 5.9 and 5.10.
The distribution of vorticity over the wake shows that most of the vorticity is concen-trated in the two discrete vortexes.
Looking at the dipole strengths we see that we now have two sudden jumps
corresponding to the two vortices.
At time step 300 a fully developed vortex is situated at the corner of the cross section. Fig. 5.11
We see that the wake has rolled up around the upper discrete vortex. This roll up is very sensitive because the two discrete vortices in fig. 5.5.d are of approximately equal strength, the lower being a little stronger than the upper. Therefore we might think that the wake end will roll up around this lower vortex in stead of the upper.
The tendency to roll up Is however dependent on two factors: The strength and posi-tion of near by discrete vortices, and the sign of the vorticity in the wake at the wake end.
If the sign is positive, the wake will have a tendency to roll up around the positive discrete vortex and vice versa, but if for instance a weak but positive wake end is close to a negative vortex, the roll up will take place around the negative vortex.
This was observed during some of the simulations, and the situation was difficult to handle because of Instabilities occurring after a short time.
Looking at the distribution of vorticity along the wake at time step 300 (fig. 5.12), we see that we have a situation where a strong single vortex close to the corner of the body is connected to a double vortex. The part of the wake connecting the single and double vortex has a very weak vorticity distribution.
We may therefore perform a dumping where one part of this weak wake is included in a discrete single vortex at the corner of the body, and the other part Is included in the upper discrete vortex of the double vortex.
The result is shown in fig. 5.13.
Fig. 5.14 shows the situation at time step 370, where a double vortex has developed at the corner of the cross section.
In fig. 5.15 the double vortex has been dumped into two discrete vortices (time step 371)
We now have generated two vortex pairs from the corner. The dumping procedure can continue along these lines.
41
-An Instability problem will occur when a Riemann cut moves into a corner of the body. This will therefore become a problem if we continue to generate vortex pairs during the simulation. We found however that double vortices far from the body did not influence the flow field at the body, and in later simulations only the "youngest" double vortex was kept.
Notice the paths of the double vortices. They move mainly under the influence of each other, which means that the path should be a straight line if they were of equal but opposite strength. If one of them is stronger than the other the weakest will rotate slowly about the stronger.
This effect is demonstrated in fig. 5.17 , where the shear layer pattern is shown at different time steps during the simulation.
5.3 Cutting of the shear layer
The performance of "dumping" during the time simulation gave us valuable experience in understanding what situations and occurrences was physically important and had to be modelled carefully ), and when we were allowed to introduce approxi-mations.
In addition to dumping, we also have other means of reducing the computing time, and stabilizing the simulation.
In our model, the circulation once being generated has to be conserved. Therefore we wlU not violate the conservation of circulation If we are cutting the wake at a point where the potential jump Is close to zero, such as indicated in fig. 5.5.a.
There are also other integral Invarlants of the vorticity distribution which should be conserved during the simulation (Batchelor - l%7) , but we experienced that the con-servation of circulation alone was sufficient In order to obtain acceptable results. The drastic elimination of a significant part of the wake will of course change the flow field around the cross section, but our conclusion was that under special circumstances cutting was allowed.
We know from experience with the vortex shedding process that we start with generation of negative vorticity during almost the whole first quarter of the roll
period.
Consequently a negative circulation is building up if we integrate over the wake. When the shedding direction of the fluid at the corner changes, a positive contribution to the drculatlon Is Introduced, and at some time step the amount of positive vorticity is equal to the amount of negative vorticity, giving zero circulation as a result.
We observe from the post processing routines that this happens for the first time when we have developed the first double vortex, as shown in fig. 5.5.a.
As a test we cut the wake at time step 151, i.e. just a few time steps after the point of zero circulation was observed. (Fig. 5.18) - The parameters of the simulation are still the same as those given In chapter 5.1.
Since the pointof zero circulation will not in general coincide with a wake node, we Cut the wake at the node were the dipole strength was least, and represented the rest of the wake with a single vortex with a strength equal to this circulation.
This representation of the wake, with a very weak vortex at the end, created a numeri-cal problem after a while because the vortex drifted towards the body. (Fig 5. 19.a) This illustrates one of the major differences in the behaviour between double vortices and single vortices.
Double vortices will tend to move away from the body, and thus create no numerical problems, while single vortices will not induce any velocity on themselves, and may therefore, depending on the circumstances, come too close to the body.
Since the problem we experienced with the discrete vortex is numerical, and does not have anything to do with the physical process, we tried to cut the wake at the same point as before without introducing a discrete vortex at the end. The situation at the same time step as in fig. 5.19.a is shown in fig. 5.l9.b.
This would violate the condition of constant circulation, but since the strength of the discrete vortex was small we hoped that the effect on the simulation would be negligi-ble.
The computer test showed that this seemed to be true, but in order to be sure we finally decided to perform interpolation on the last wake element, and adjust the end node to correspond to the exact point of zero circulation. Fig. 5.19.c shows the situa-tion at the same time step as in fig. 5.l9.a and b.
We found that the cutting In all three cases resulted in a jump in the moment time his-tory curve compared to the same time hishis-tory calculated on the basis of no wake cut-ting. (Fig. 5.20)
The curve showing the amount of vorticity generated at the corner of the body in each time step reflects thS same effect.(Fig. 5.21)
This Is due to the fact that the double vortex is close to the corner of the body, and therefore its effect is significant In spite of the fact that it gives no net contribution to the circulation.
However, we see that both the roll-moment curve and the vorticity-generation curve approach the curves where no cutting was done after some time.
The influence of the cutting on the hydrodynamic coefficients which may be calculated from the moment time history will be small if we consider the example above.
43
-It is however not certain that this will be the case in general, and in order to be on the safe side we decided to Cut the wake when the double vortex had moved further away from the body.
We found that If we Cut the wake at a time step where the point of zero circulation was situated between a new (single) vortex developing at the corner and the double vortex, the discontinuity in the time histories were negligible. In this case the flow near the corner of the body is dominated by this new single vortex.
The dipole strength for this situation is shown In fig. 5.22.
Further tests revealed that it is not necessary to have a fully developed new vortex at the corner in order to get reasonable results. It is sufficient that the vortex is just start-ing to develop under the condition that the double vortex is not to close to the body. There is also another numerical benefit that we get if we are cutting the wake as soon as possible after the new vortex has started its development. This has to do with the stretching and rediscretization of the wake.
A double vortex moving away from the body with great velocity will contribute to a stretching" of the part of the shear layer connecting the double vortex to the rest of the wake.
Because we In each time step are generating wake elements with equal ler'gths, this stretching will contribute to the generation of long wake-elements after the rediscreti-zation. Therefore the new vortex at the corner of the body will not be represented in a proper manner.
One can generally say that It would be ideal if the wake elements and body elements were of approximately the same length.
If we are cutting the wake as soon as possible after the double vortex has started to develop, the wake elements will be much smaller than in the case were stretching due to the presence of the double vortex occurs.
We observed from the moment curve that the moment is much more irregular in the case were we have long wake elements close to the body than in the more favorable case were the elements are shorter,
A more drastic cutting than the one performed at time step 150 is demonstrated in fig. 5.23, showing time steps 342 and 343. Here the double vortex is very close to the body.
We see from the moment and vorticity generation curves (fig. 5.20 and 5.21) that the difference between the curves with and without cutting is greater now than at time step 151. but also in this case we see that the simulation is stable, and the solution becomes identical to the solution without cutting after some time.
the the limitations described above. Since the cutting and dumping will alter the state of the system which is being simulated, instabilities may be occur. The probability of generating such Instabilities will increase when the dumping or cutting significantly changes the physical state of the system.
Computer tests with and without cutting and dumping of the wake were performed systematically in order to assure that the difference in results was within acceptable limits.
5.4 Additional numerical problems
During the time simulation the situation shown in fig. 5.24 a-c may occur. Fig. 5.24 is a close up of the area close to one of the corners of the body. We see that a wake element of considerable length enters the corner of the body. One disadvantage with this situation is that the wake element is long, which is due to the stretching of the shear layer caused by the double vortices, but the main problem is that a wake element entering the corner will give wrong shedding velocity, giving rise to an instability in the solution.
One reason why this situation occurs is due to the rather primitive Euler method used in the convection of the wake.
Tests showed that the Instability appeared later when we used a smaller time step. In order to get rid of this (unphysical) instability we introduced a small "barrier" at the corner of the body (Fig. 5.25) and demanded that no wake element were allowed to cross this barrier.
If this condition is satisfied the wake element will not be able to enter the corner. Another situation that may give rise to numerical difficulties was shown in fig.
5.1 .f.
We see that a part of the shear layer is very close to the body, and In some cases the element may be convected into the body. One of the reasons for this may be that we in our method satisfies the boundary condition on the body midpoints.
Since the wake element crossing the body contour is a non-physical event we are
allowed to push the element outside again.
A separate subroutine was written in order to detect if an element had crossed the body contour.
A third problem, already mentioned in chapter 5.2, occurs if a Riemann cut enters one of the corners of the cross section. (Fig. 5.26)
45
-Our tests showed however that the "oldest double vortices did not affect the simula-tion results.
Therefore only the "youngest double vortex was kept, end the problem with the
5.5Figures