A Method for Accurate Force Calculations in Potential Flow
Heinrich Söding, Instit ut für Schiffbau Hamburg1 TECHISCIE UNIVERSIT1T
Laboratorium voor
ScheepShYdrOm0i Archief
i
Introduction
MekeiWeg 2, 2628 CD DeiftThree-dimensional potential flows are traditionally calculated using thIpt nrtnoi 01
Hess and Smith (1962) or similar methods. Foi- a few hundred or a few thousand panels on one half of the wetted ship huh surface, such methods give velocities and pressures with a relative accuracy of some 10_2 to some iU. To determine the wave resistance of merchant ships. however, we can tolerate relative inaccuracies of 1 10 at most because the wave resistance is
a very small quantity compared to stagnation pressure times wetted surface. The inaccuracies of the usual panel methods are obvious when the wave resistance coefficient C,, is plotted over speed: Instead of the correct limiting value O for small speed. one often finds a limiting value
in the range of +2 i0
to ±iO. That was one reason for .Jcnsen (1988) to introduce thespliei-e niethod'; hut also with this method noti-zero C values for small speed show that the accuracy of this method is still tiot fully sufficient for resistance prediction.
The following remedies to this problem could be conceived:
- to subtract the resistance coefficient found for small speed from the results for arbitrary speed;
- to use Lagahli's theorem instead of pressure integrations to calculate the resistance; - to cleternune the resistance from tile far-fleld waves.
The fit-st method is at least dubious auch probably inaccurate for higher speed. It eliminates also a very useful accuracy check. Tile secoiid niethod cannot account for the actual inunersion of the ship, at least not exactly. Further, this method avoids oniy one part (maybe one half) of the inaccuracy. The third method requires a larger area of the free surface around the ship to be discretized than would be necessary if tile resistance were derived froni integrating the pressure over the hull, and it causes additional coinplexJty of the program.
Another application of resistance values in potential flow is connected to Navier-Stokes flow computations foi- tile afterbody only. To determine the viscous resistance of a total ship
from such a calculation, one has to add the resistance of the forebody to the Navier-Stokes resistance' of tile afterbody. Because flow separation effects are very small in the forebod. the forebody resistance can be estimated as a sum of two parts: (1) the viscous resistance
of the forebody determined by boundary-layer calculations or simply from the plate friction coefficient together with a foi-m factor; (2) the pressure resistance on the forebody determined, e.g., by a. panel method. (The pressure resistance of a closed body, but not of a half-body, is zero in potential flow.) Foi- this application we need the potential pressure resistance of the forebody with an accuracy of about 1O times stagnation pressure times vetted hull surface.
The method described hiere produces such accurate predictions of forces on half-bodies or full bodies in unbounded ideal fluid. Applications to wave-resistance calculations have not yet been tested. The purpose of the method is not to increase the accuracy of computed velocity and pressure at points on the hull or within tite fluid; iristea.d it aims to accurately determine the velocity and pressure avelage over surface patches. Surface patches are small parts of the body surface which together cover the body surface without overlaps or gaps; in the present miplemen tation they are curved triangles or quadrilaterals.
High accuracy of velocity and pressure at arbitrary points on tile hull would require very fine patches and extreme care in reproducing tile curvature of tite hull in the computational model. By restricting tile accuracy requirements on the practically interesting patch averages.
'Li.mmersie.th 90, D 22305 Hamburg, Germany
a simple, fast, second-order accurate method is found. For moderate patch numbers the error of hull forces is decreased to less than 10% compared to other known methods.
2
Considerations about panel methods for exterior incompressible flow
Panel methods and other potential flow calculation methods approximate the potential () at a point .ïofa flow around a smooth body by that of tite incoming how (here -Ux for uniform flow of velocity U in negative r direction) plus a linear superposition of functions pj():
= -Ux +>
ii().
(1)The functions j(.i) satisfy the differential equation = O (which corresponds to the con-dition of incompressibility) and the boundary concon-dition in infinity V( ) = (U, 0, 0) (vanishing disturbance of the incoming flow) exactly. The accuracy is thus determined by the errors in satisfying the body boundary condition which requires that there is no flow through the body surface.
In so-called collocation methods this condition is satisfied at a number of points (collocation points) distributed on the body surface. The methods of Hess orni Smith (1962) and Jensen
(1988), but not the present method, are such collocation methods. To obtain small errors of the body boundary condition also between the collocation points, collocation methods have to use sufficiently smooth functions (Z). In Hess and Smith's method, the potential of sources distributed with constant source density over plane panels is used as p(.Z); the panels approximate the surface patches. Other methods arrange the panels not at the body surface,
but within the body ( Webster 1975), or they use curved panels with linearly varying source
density. These modifications smooth the functions for points ou the body surface and thus improve the accuracy for a giveil number of panels.
Better accuracy can he obtained by satisfying the body boundary condition in the average
over each surface patch. That means: the total flux through each surface patch is made to
vanish. Such niethods have been used before: but in connection with the functions pj(.) corre-sponding to sources distributed over panels, they require ittiinerical integrations to determine
the total flux due to each (.), and these integrations introduce additional errors and cost computer time. Difficulties arise also from the fact that the gradient of the panel potential tends to infinity at the panel edges; this makes integrations Rear tite panel edges inaccurate or impossible.
3
PrincipIes of the present method
The present method satisfies the condition of no flux through each surface patch, i.e. it satisfies tite body boundary condition in the average over each surface patch. Contrary to the usual panel methods. it uses no source panels (distributed sources) hut the potential of point sources a.s functions p(F). Therefore it is called hiere a. patch method, not a panel muelliod. The point sources may be placed at the approximate centres of tite surface patches. The flux through a polygon induced by a poiut source caii be determined exactly; thus numerical integrations are avoided if the surface patches are bounded by polygons. Oit a generally curved surface tite patch boundaries are curves; they are approximated simply by straight hues. The error introduced if the edges of a surface patch are curves instead of straight lines are conservative: a positive error in one patch is compensated by a negative error of equal magnitude in the adjacent patch. Titus errors introduced by curved patch boundaries may he expected to have little effect. The area enclosed between a curved panel side and tite straight-line connection of the end points of the panel side is proportional to the third lower of the panel length. Thus the error decreases strongly with decreasing panel size. The flow through titis area could he
app roxiina.ted if titis correction were found wort Itwhtile.
a
Fig. 1. Source point S and patch ABC
The functions () used here. i.e. the potential of point sources, may he considered less
smooth than typical panel source potentials. However, smooth p() are required only by col-location methods. not necessarily by methods satisfying tile boundary condition in theaverage over patches. In fact, averaging the boundary condition over a patch has a similar effect as averaging the source location over a panel. In spite of this, smoothing of j(.) is advantageous also here. I do this he shifting the sources from the body surface somewhat into the interior of the body (which is often called desingularisation').
Panel methods have to avoid evaluating V near to the panel boundaries because VØ is not smooth (frequently infinite) there. The present method must avoid to evaluate V at or near the patch centres because there, near the point sources, corresponding difficulties occur.
My program uses triangles and quadrilaterals as surface patches. The curvature of the
patch sides is neglected. For determining velocities and pressurès on the hull, quadrilaterals are divided by a diagona.l into two triangles. For non-plane quadrilaterals, of the two possible subdivisions one produces a concave. the other a convex patch: the convex subdivision is chosen. The average velocity V within each triangle is determined not as an analytical derivative of i), but from the difference of the potential evaluated a.t the three corner points of the triangle. Errors are again conservative. i.e. an error of the poteittial in one corner point produces errors of different signs in the triangles around this point.
Tile l)1'eSSU1e depends quadratically on the absolute value of the velocity. Using the average
velocity in Bernoulli's equation would result in a pressure average being too high in cases
of non-constant absolute velocity. To decrease tlus error, the pressure average is determined from integrating a quadratical approximation of the velocity derived froni the potential and its gradient at the three corners of each triangle.
4 Details
The body boundary condition is
(2)
where 7 is the inward normal vector ou the hull and n1 is its z component. p() is 1/distance between . and the ith source point; q1 is the negative source strength of this source divided by
4ir.
Fig. I shows a triangular patch ABC and a source poilit S. Integrating (2) over tile area. ABC gives
(1(aX
2 q1c1 = O. (3)
Tile first terni is the volume flow through ABC due to the uniform flow; tite index i indicates
the z component of the vector product of two side vectors of the triangle (Fig. I). The ¿th
term of the sum in (3) 15 tile volume flow through ABC induced by the ith source
is tite solid angle in which ABC is seen from 5. This is easily verified: Tile tota.l volume flow
--approaching the ith sink (source with negative source strength) is 4irq. As the flow is directed radially to S. the flow passing the triangle ABC is proportional to the solid angle a in wiuch the triangle is seen froui S. lt would amount to the total flow 4rq for a solid angle of 47r (all directions).
The rules of spherical geometry give as the sum of the angles between each pair of planes 5.-IB. SBC, and SCA minus 7r:
= 3SAB,SBC + I3SBC,SCA + J3SC,,SAB - ir (4) where. e.g..
[(.x
Thx (ilx)].
/3S.4B,SEÌC = arctan . (.5)
(Ax B)(BxC)IBI
Here
.,
ê, (7 are the vectors pointing froni the source point S to the panel corners A, B, C. Further details referring to the non-uniqueness of arctan and of aj (Q'j +47r is equivalent to a)may be seen in the program in the Appendix.
A simpler approximate formula is used if the distance d between patch vertices and source point exceeds a given limit. typically .5 times the maximum panel length. In this case the solid angle is approximated as
= f/d
(6)where f is the patch area projected on a plane normal to the direction from the source to the patch centre:
1A+B+C
d=kTI;
f=x.
(7)For a c1uadrilateral patch, results for two triangles are added.
Source points are arranged on the straight lines normal to the surface through the patch
centres. Patch centres are approximated as the average of the 3 or 4 patch vertices. The procedure for determining a gives results changing froni 2 to 2ir if the source point moves
from inside the body through the patch plane to outside. To avoid difficulties following from this discontinuity and prolleiìs of numerical accuracy, and to smooth tite functions j, the
sources are pla.ced inside the body in a user-specified distance from the body surface. If this distance brings a source closer to time centre plane than 1/4 of the body breadth (in y direction) a.t the panel centre, the distance is decreased for those panels so that the sources are placed midway between body surface and centre plane.
To account for port-starboard synimetry, patches are assumed on port side only. but sources of equal strength are assumed both ou port and starboard side when determining the sum in (3). Other symmetry planes (e.g. the undisturbed water plane) are handled correspondingly.
Equations (3) for each patch form a. linear system froni wInch the unknowns qj are deter-mined. I use an elimination of the most important elements below the main diagonal iii the coefficient matrix followed by a. Cauß-Seiclel iteration. hut other methods may be used as well. With known qj one can determine o and V at all patch corners using (1) and its gradient
= (U.0,0)+
qV(F).
(S) From the values at corners .4, B, C of a triangle, the average velocity within the triangle is found as= Vo =
(& òc,) -.
- AB +(B )
-
'i_mc -b114ßvit h
= and n .,-
-
c - --- b.
With known and the corner velocities Q, ÎB,VC, tile pressure force on the tríaiigle can be
determined:
- I
[P
'
f = n
J
pda = nJ
(U - v) da.
(11)To evaluate this expression. tile velocity within the triangle is approximated as
(12)
r is the 'triangle coordinate' directed to A: r is
1 at A and O at the line BC.
s and t
are the corresponding 'triangle coordinates' directed to R and C respectively. Using this tri-quadratic f formula, the integral of i over the triangle area in (il) is found after some algebraic nianipula.tions as
f2da
=f da.
+VB_V)2+(Vc_U)2
(13)
5
Results
As a. first test, the flow around a deeply submerged sphere in uniform inflow in negative
s (lirectioti was investigated. For sphere radius 1 and inflow velocity 1, the force on 1/8 of
the sphere having positive s, y and negative z coordinates (corresponding to tile port lorebody of a half-spherical ship) was determined. Exact values are
= ir/64 = 0.04909 and
f, = -f. = lLir/128 = 0.26998. Only triangular patches were used (Fig. 2). Tests showed that coiubuung triangles in pairs to quadrilaterals gave practically tile same accuracy with a little more than half tile number of panels. Table i gives numerical results of my method. It indicates that the force errors decrease about proportionally to the square of the patch length, i.e. inversely to the panel number. In this example, errors depend substantially on the distance of tile sources from tite sphere surface which was chosen as 10% of tile sphere radius for the data of Table 1. If a source distance proportional to the panel size is used, the errors decrease approximately proportional to the 1.5th power of the patch length.
Table 1. Numerically computed forces on 1/8 of a sphere in uniform flow ;
Fig. 2. Two discretiza-tions of 1/4 sphere
i "J
The 'HSVA tanker' wa.s used as a second test. Mirror images reflected at tIte waterline and a.t tile centre plane were used. No exact results are known with one noii-trivial exception:
180 Schiffstechnik Bd. 40 - 1993 / ShipTechnology Research VoI. 40 - 1993
110. of patches 011 1/8 sphere f. error of L
f
error of h4 -0.0400 0.0091 0.1143 0.1557
16 -0.0427 0.0064 0.2200 0.0500
64 -0.0481 0.0010 0.2581 0.0119
the longitudinal force on the whole ship is zero. For a numerical method based on pressure integration this is non-trivial because of the extreme differences in forehody and afterbody shape (Fig. 3). Results for the resistance coefficient
U. (.ietz, Institut für Schiffbau, prtvate communication
Schiffstechnik Bd. 40 - 1993 I Ship Technology Research Vol. 40 - 1993 181
Fig. 3. Three patch arrangements of HSVA tanker
Pine: Arrar.gemenc I
?ane Arrangemen; 2
?ane Arrangement
Compared to the coarsest panel arrangement 1, the finer subdivision of the section con-tour used for arrangement 2 was no improvement, whereas increasing the number of sections (arrangement 3) was helpful.
The above results of my method were produced using a source distance from the hull surface of 0.5m. Table 2 shows the effect of varying this quantity for patch arrangement 1. W'eb.tei's
C = _f1/(LT2S)
(14)were determined using the methods of Hess and Smith (1962,) and Jensen (1988) and the
present method:
Jensen's method with patch arrangement i (Fig. 3; 421 patches): C = +1.54V i0
Present method with patch arrangement 1: C7, = -0.16. i0-Jensen's method with patch arrangement 2 (Fig. 4; 818 patches):
C = +1.90 iO
Present method with patch arrangement 2: = -0.16. iO
Jensen's method with patch arrangement :3: (Fig. 5: 780 patches):
C = +0.75
i0-Present method with patch arrangement 3:
C = +0.03 iO
i0-(1975) recommendation of using about O.5.minimum body radius of curvature is found here to give optimum results (the ship has a bilge radius of 1.63m). Smaller or zero curvature radius at
the centre plane (e.g. sharp waterline ends) need not be considered in this connetion because of the smaller source distance used for sources approaching the centre plane. Only for very fine panel arrangements the source distance should he decreased to improve tile condition nuniber of the equation system (3).
Table 2. Longitudinal forces in kN for U = i on the HSVA tanker hull
(L B - T = 253 - 38.33. 14.2 m) depending on the distance of sources from the hull surface in ni
Lagalli's equation f1 = U -47rq gives a longitudinal force coefficient ill my met hod of 0.00? iü for all patch arrangements when using four-byte arithmetic for real numbers:
without machine roundoff it would give exactly zero in this method. As explained above,
however, this is not very relevant in practice.
Naturally the longitudinal force on the fore- or afterbody alone need not be zero. However, for the IISVA tanker this force turned out to he smaller than the accuracy obtained with all panel arrangements tested here; presumably jC for a half-body of the lIS VA tanker is smaller
than 0.05 . i0.
Computing times of the patch and the sphere method are similar: about 10 minutes on
a VAX 6310 (ranging in speed between a PC and a modern workstation) for the 780 patch arrangement.
References
JENSEN, G. and SÖDENG, H. (1988), Ship wave resistance calculations, Notes on Numerical Fluid
Mech. 25, Vieweg, Braunschweig
HESS, .1.L. and SMITH, AMO. (1062), Calculation of non-lifting potential flow about arbitrary three-dimensional bodies, Douglas Aircraft Division Report E.S.40622
WEBSTER, WC. (1975), The flow about arbitrary three-dimensional smooth bodies, J. Ship Res. 19
(4)
Appendix: Fortran-77 program
C Patch method of 3D potential floc using triangles and quadrilaterals C Mirror planes y=O and zZSP
PARAMETER (NEMAX= 1000, NPMAX= 1000)
DIMENSION X(3,NEMAX),NRPE(0:4,NPMAX),Z(3),PZ(3,NPMAX),A(3),B(3),
C(3) ,AA(NPMAX,NPMAX),R(NPMAX) ,xQ(3) ,Q(NPMAX) ,xA(3) ,XB(3) ,xC(3),
& VNAC(3),VNAB(3),PHI(NEMAX),XBA(3),XCA(3),F(3),UVW(3,NENAX),
Sc VQ(3),VA(3),VB(3),VC(3)
COMMON RLIM,IDREI,FL(3,2*NPMAX) C Reads data of NE points and NP patches
READ(7,*)NE,((X(ID,IE),1D1,3),IE=1,NE)
READ(7,*)NP, (NRPE(0,IP) ,(NRPE(IE,IP),IE=1 ,NRPE(0,IP)) ,IP=1 ,NP)
IF(NE.GT.NEMAX.OR.NP.GT.NPMAX)STOP'Too many points or panels'
WRITE(*
182 Schiffstechnik Bd. 40 - 1993 / Ship Technology Research Vol. 40- 1993
distance afterbody force forebodv force total force
0.2 -0.24 2.74 2.50
0.4 0.08 1.75 1.83
0.6 0.08 0.99 1.06
0.8 0.24 0.43 0.67
'Coordinates of point seeing all patches from interior side?'
READ(* ,*)Z
WRITE(*,*)'Mirror plane height? x midship section?' WRITE(*,*)'Minimum distance for simplified formula?' WRITE(*,*)'Sources inside body?'
READ(z ,*)zsp ,XM,RLIM,VMAX
RLIM= RL IM* *2
C For all patches DO 30 IP=l,NP C Patch centres
DO lo ID=l,3
PZ(ID,IP)=X(ID,NRPE(1,IPfl+x(ID,NRPE(2,Ipfl+x(ID,NRpE(3,Ip))
IF(NRPE(O 1?) .E .4)PZ(ID ,IP)=PZ(ID ,IP)+x(ID ,NRPE(4,IP))
10 PZ(ID, IP)=PZ(ID, IP)/NRPE(O, IP)
C Clockwize renumbering of patch vertices seen from outside DO 20 ID=l,3 A(ID)=X(ID,NRPE(3,IPfl-X(ID,NRPE(l,IP)) 20 B(ID)=X(ID,NRPE(NRPE(0,IP),IPfl-X(ID,NRPE(2,IP)) C(1)=A(2)*B(3)-A(3)*B(2) c(2)=A(3)*B(l)-A(1)*B(3) C(3)=A(1)*B(2)-A(2)*B(1) IF(C(l)*(PZ(l,IP)-Z(l))+C(2)*(PZ(2,IP)-Z(2fl+C(3)*(PZ(3,IP)-Z(3)) & .LT.0)THEN SAVENRPE(NRPE(0,IP) ,IP) NRPE(NRPE(0,IP) ,IP)=NRPE(2,IP) NRPE(2,IP)=SAVE C(l)=-C(l) C(2)=-c(2) c(3)=-C(3) ENDIF
C Renumbering for convex triangular subdivision of quadrilateral panels
IF(NRPE(0,IP) .EQ.4)THEN IF(SOLIDA(X(l,NRPE(l,IP)),X(l,NRPE(2,IP)),X(1,NRPE(3,IP)), & X(1,NRPE(4,IPfl).LT.0)THEN ISAVE=NRPE(l ,IP) NRPE(l ,IP)=NRPE(2,IP) NRPE(2,IP)NRPE(3,IP) NRPE(3,IP)=NRPE(4,IP) NRPE(4,IP)=ISAVE ENDIF ENDIF
C Determination of source point coordinates FAK=VMAX/SQRT(C(1)**2+C(2)**2+C(3)**2)
IF(PZ(2,IP)-FAK*C(2).LT.0.5*PZ(2,IP))FAKO.S*PZ(2,IP)/C(2)
PZ(i ,IP)=PZ(l ,IP)-FAX*C(l) PZ(2,IP)=PZ(2, IP)-FAK*C(2)
Pz(3 ,IP)=PZ(3 , IP)-FAK*C(3) 30 CONTINUE
C For i or 2 triangles of every patch. Vertices 11,12,13, area vector FL IDREI=0 DO 45 IP1,NP R(IP)=O DO 45 IDR=3,NRPE(0,IP) IDREI=IDREI+ 1 Ii=NRPE(l IP) 12=NRPE(IDR-1 ,IP) I3NRPE(IDR, IP) DO 38 ID'1,3 A(ID)=X(ID ,13)-X(ID ,12)
38 B(ID)=X(ID,12)-X(ID,I1)
FL(1 ,IDREI)=(A(2)*B(3)-A(3)*B(2fl/2
FL(2,IDREI)=(A(3)*B(1)-A(1)*B(3fl/2
FL(3 ,IDREI)=(A(1)*B(2)-A(2)*B(1))/2
C Coefficient matrix AA and inhomogeneous terms R for speed i
iii
direction -xR(IP)R(IP)-FL(i, IDREI)
DO 40 JP=i,NP
IF(IDR.EQ. 3)AA(JP ,IP)=0
DO 40 IYSYM=-i,1,2 DO 40 IZSYM=O,1
C XQ source point; function SOLIDA gives solid angle alpha_i
XQ(i)=PZ(1 ,JP) XQ(2)=PZ(2,JP)*IYSYM XQ(3)=PZ(3,JP) IF(IZSYM.NE.0)XQ(3)=2*ZSP-Xq(3) AA(JP,IP)=AA(JP,IP)+SOLIDA(X(1,Ii),X(1,12),X(1,13),XQ) 40 CONTINUE 45 CONTINUE
C Solving linear equation system A Q = R for source strengths Q. C Other methods (e.g. Gauss elimination) will work, too.
CALL SINQIT(AA,R,NP,NPMAX,1.E-S,SQRT(1./NP),i.E-6,Q)
C Sum of source strengths
QSUM=0
DO 50 IP=1,NP
50 QSUMQSUM+Q(IP)
WRITE(6,*)'Sum of source strengths',QSUM
C Determination of the potential and
its gradient at all vertices
DO 75 IE=1,NE PHI(IE)=-X(1 ,IE) Uvw(1 ,IE)=-1 U1.rw(2,IE)=0 UVW(3,IE)=0 DO 70 JP=1,NP DO 70 IYSYM=-1,1,2 DO 70 IZSYM=0,1 XQ(i)=PZ(i ,JP) XQ(2)=PZ(2,JP)*IYSYM XQ(3)=PZ(3 3?) IF(IZSYM.NE. 0)XQ(3)2*ZSP-XQ(3) RR=SQRT((X(1,IE)-XQ(1))**2(X(2,IE)-XQ(2fl**2(X(3,IE)_xQ(3fl**2) PHI(IE)=PHI(IE)-Q(JP)/RR Uvw(1,IE)=Uvw(1,IE)+Q(Jp)/RRX*3*(x(1,IE)-xQ(1))Uvw(2 ,IE)=UVW(2,IE)+Q(JP)/RR**3*(X(2, IE)-XQ(2))
70 UVW(3, IE)=UVW(3,IE)Q(JP)/RR**3*(X(3, IE)-XQ(3))
75 CONTINUE
C Force integration. FE (FV) = long, force ort afterbody (forebody)
F(1)=0
F(2)=0
F(3)=0
FH=0 FV=0
C For one or two triangles of all panels
DO 150 IP=1,NP DO 150 IDR=3,NRPE(0,IP) I1=NRPE(i,IP) 12=NRPE(IDR-1 1?) 13=NRPE(IDR, IP) DO 120 10=1,3 VA(ID)=UVW(ID Ii)
VB(ID)=UVW(ID 12) vc(ID)=Uvw(ID,13) XBA(ID)X(ID,12)-X(ID,I1) 120 XCA(ID)=X(ID,13)-X(ID,I1) C Average velocity VQ SBACA=XBA(1)*XCA(1)+XBA(2)*XCA(2)+XBA(3)*XCA(3)
SBABA=XBA( 1) **2-FXBA( 2) **2+XBA (3)
SCACA=XCA(1)**2+XCA(2)**2+XCA(3)**2 DO 130 10=1,3 VNAB(ID)=XCA(ID)-SBACA/SBABA*XBA(ID) 130 VNAC(ID)=XBA(ID)-SBACA/SCACA*XCA(ID) XBANAC=XBA(1)*VNAC(1)+XBA(2)*VNAC(2)+XBA(3)*VNAC(3) XCANAB=XCA(1)*VNAB(1)+XCA(2)*VNAB(2)XCA(3)*VNAB(3) DO 140 ID=1,3 140 VQ(ID)=(PHI(I2)-PHI(Ilfl*VNAC(ID)/XBANAC (PHI(13)-PHI(Ilfl*vNAB(ID)/xCANAB C Integration of force on triangle
vQvq=v(1)**2+vQ(2)**2v(3)**2 VAVQ(VA(1)-VQ(1))**2+(VA(2)-VQ(2))**2+(VA(3)-VQ(3))**2 VBV=(VB(1)-VQ(1))**2+(VB(2)-VQ(2D**2+(VB(3)-VQ(3))**2 vcv=(vC(1)-vQ(1fl**2(vc(2)-vQ(2))**2+(vC(3)-VQ(3))**2 VAVB=(VA(1)-VQ(1))*(VB(1)-VQ(1fl+(VA(2)-VQ(2fl*(VB(2)-VQ(2)) & +(VA(3)-VQ(3))*(vB(3)-vQ(3)) vBvc=(VB(1)-vQ(1))*(vc(1)-v(1fl+(vB(2)-vQ(2fl*(vc(2)-vQ(2)) & (VB(3)-vq(3))*(VC(3)--VQ(3)) VCVA=(VC(1)-VQ(1))*(VA(1)-VQ(1))+(VC(2)-VQ(2fl*(VA(2)-VQ(2)) & (vc(3)-VQ(3))*(VA(3)--vq(3)) RINT=VQV(VAVQVBVQ+VCVQ)/30.-(VAVB+VBVC+VCVA)/9O. p=o.s*(1-RINT) F(1)=F(1)+P*(-XCA(2)*XBA(3)+XCA(3)*XBA(2)) F(2)=F(2)+P*(-XCA(3)*XBA(1)+XCA(1)*XBA(3)) F(3)F(3)+P*(-XCA(1)*XBA(2)+XCA(2)*XBA(1)) IF(PZ(1 1?). GT.XM)THEN FV=FV+P*(-XCA(2)*XBA(3)+XCA(3)*XBA(2)) ELSE
FH=FE+P* (-XCA (2) *XBA (3) +XCA( 3) (2)
ENDIF
160 CONTINUE
WRITE(6,*)Total force (3 components) for U=1',F
RITE(6,*)'Long. force on afterbody and forebody,FH,FV END
FUNCTION SOLIDA(AO,B0,CO,Z)
C Solid angle between vectors A0-Z, BO-Z, C0-Z. >0 if, seen from Z, ABC C turn clockwize. Not applicable for Z on plane ABC inside triangle ABC.
PARAMETER (NEMAX= 1000, NPMAX= 1000)
DIMENSION A(3),B(3),C(3),ABZN(3),BCZN(3),CAZN(3),Z(3), & A0(3),BO(3),CO(3),ABC(3)
COMMON RLIM,IDREI ,FL(3,2*NPMAX)
DO 10 10=1,3 A(ID)=A0(ID)-Z(ID) B(ID)=B0(ID)-Z(ID) 10 C(ID)=CO(ID)-Z(ID) IF(A(1)**2+A(2)**2+A(3)**2.GT.RLIM)THEN ABC( 1) =A( 1)+B( 1) +C( 1) ABC ( 2) =A (2) +B (2) +C (2) ABC (3) A (3) +3 (3) +C (3) SOLIDA-9*(ABC(1)*FL(1,IDREI)+ABC(2)*FL(2,IDREI)
& +ABC(3)*FL(3,IDREI)) /SQRT((ABC(1)**2+ABC(2)**2+ABC(3)**2)**3)
ELS E ABZN(1)=A(2)*B(3)-A(3)*B(2) ABZN(2)=A(3)*B(1)-A(1)*B(3) ABZN(3)=A(1)*B(2)-A(2)*E(1) BCZN(1)=B(2)*C(3)-B(3)*C(2) BCZN(2)=B(3)*C(1)-B(1)*C(3) BCZN(3)=B(1)*C(2)-B(2)*C(1) CAZN(1)=C(2)*A(3)-C(3)*A(2) CAZN(2)=C(3)*A(1)-C(1)*A(3) cAZN(3)=c(1)*A(2)-C(2)*A(1) CALL SVPR(ABZN,BCZNE,CO1,SI1) CALL SVPR(BCzM,CAzN,C,CO2,s12) CALL SVPR(CAZN,ABZN,A,CO3,S13) SOLIDA=O. FAKSI1*SI2*SI3 IF(FAK.EQ .0)RETTJRN SIJLIDA=ATAN2(SI1,CO1)+ATAN2(S12,CO2)+ATAN2(S13,CO3)-3. 1415926 IF(FAK . LT. O)SOLIDA=SOLIDA+2*3. 1415926 ENDIF END SUBROUTINE SVPR(A,B,E,CO,SI) DIMENSION A(3),B(3),E(3) CO=-A(1)*B(1)-A(2)*B(2)-A(3)*B(3)
SI( (A(2)*B(3)-A(3)*B(2) )*E(1)+(A(3)*B(1)-A(1)*B(3))*E(2) (A(1)*B(2)-A(2)*B(1))*E(3)) /SQRT(E(1)**2+E(2)**2+E(3)**2)
END
From the editors' software collection
Minimization of a Linear Function Subject to Linear Constraints
(rest, continued from issues i and 2 Vol. 40)
C Search for line
AM1 .E30 DC 60 1=4,11+3 TOT(I ,Ni) IF(T0.GT. 1.E-20)THEN AM1T(I ,NN+1)/T0 ELSEIF(TQ.GT.-1.E-20.OR.IFIX(ABS(T(I,NN+3fl).GT.NN)TREN AM11 .520 ELSE AM1=(T(I ,NN+1)-1j/T0 ENDIF IF(AM1 LT. AM)THEN 11=1 AM AMi T2=T0 ENDIF 60 CONTINUE
C Change to the other boundary of a variable if necessary J=ABS(T(II5,N1))
IF(J.LE.NN.AND.AM.GE.1.)THEN
C Change to the other boundary in column nl
DO 61 1=3,11+5 T(I ,N1)=-T(I ,N1) T(I,NN1)=T(I,NN+l)+T(I Ml) 61 CONTINUE GOTO 45 ENDIF IF(T2.LT.0.)TEEN
C Change to the other boundary for basis variable corr. to line il
DO 62 N=1,NN+3
IF(N.NE. NN+2)T(I1 ,N)=-T(I1 ,N) 62 CONTINUE T(I1 ,NN±1)T(I1 ,NN+1)+1 T2=-T2 ENDIF C Change of basis SAVE=T(I1 ,NN+3) T(I1 ,NN+3)=T(II+5,N1) T(II+5,N1)SAVE T(I1 ,N1)=-1 DO 70 I=K+2,II+4 70 T(I,N1)=-T(I,N1)/T2 DO 80 Nl,NN+1 IF(N.EQ.N1.OR.T(II+5,N).EQ.0.)GOTO 80 T0=T(I1 ,N) DC 75 1K+2,II+4 75 T(I,N)=T(I,N)+T0*T(I,N1) T(I1,N)=T(I1 ,N)-T0 80 CONTINUE IF(K.E.2)GOTO 45
C Prepare equations and inequalities
J=ABS(SAVE)-NN+3
IF(J.LT.4.OR.T(J,NN+2).LT.0.)GOTO 90
IF(T(J,NN+2) .EQ.0. )THEN
T(II+S ,N1)=O. ELSE DO 85 1=3,11+4 85 T(I,N1)=-T(I,N1) T(J ,NN+2)=-T(J ,NN+2) T(II+4,N1)T(II+4,N1)+1. END IF GOTO 90 C Determination of optimum x 150 DO 160 N=1,NN N1=T(II+5 ,N) N2=IAES(N1)
IF(N2. GT.NN. OR.N2.E .0)GOTO 160
T(3,N2)=T(1,N2)+WENN(N1.LT.O,T(2,N2),0.) 160 CONTINUE DO 170 1=4,11+3 N1=T(I ,NN+3) N2=IASS(N1) IF(N2.GT.NN)GOTO 170 T(3,N2)=T(1,N2)+(WENN(N1.LT.0,1.,0.)+SIGN( & T(I,NN+1),FLOAT(N1)))*T(2,N2) 170 CONTINUE T(3 ,NN+1)-T(II+4,NN+1) RETURN END FUNCTION WENN(L,A,B) WENN=B IF (L) WENN=A END
Jahrbücher
beliebt
-
begehrt- bewährt
Internationales Jahrbuch der
Luft- und Raumfahrt
Redaktion: 1-jans M. Namislo 184 Seiten. Format DEN AS. 150 s/w- und 16
Farb-Abbildun-gen. Broschur. mit Preisrätsel,
DM 19,80 Best.-Nr. 3904.
Der neue Fliegerkalender bietet
wieder eine Reihe von
interessan-ten Beiträgen. hier die Auswahl:
. ATTAS, der fliegende Simulator
. Clement Ader: Die Fledermaus
war sein Modell
Galileo unterwegs zum Jupiter . Ariane starlet mit Gewinn . Luftfahrt auf Postkarten
Hubschrauber im modernen
Rettungswesen
Das deutsche Jahrbuch der
Seefahrt seit 1)I
224 Seiten. Format DEN AS. 120 slw- und 30 Farb-Abbil-dungen. Broschur.
DM 19,80, mit Preisrätsel
Best.-Nr. 5-141.
Im Jahre 1901 ist erstmals der
Flottenkalen-der erschienen. Vom Gründungsjahr bis
heute durchiebte dieses Jahrbuch zahlrei-che Höhen und Tiefen, die auf zahlreizahlrei-chen
Seiten als deutsche Seefahrtgeschichte dokumentiert sind. Hier ist eine ausge-zeichnete Möglichkeit, längst
'verschol-Iene Bücher wieder zu erhalten - und damit
die eigene Bibliothek zu bereichern. Alle folgenden Bände können im Abonnement
bezogen werden.
Lieferbar sind jetzt die Ausgaben von 1901 und 1902
Jedes Jahr - beginnend ab 901 - erscheint eine weitere Reprint-Ausgabe der bei allen
Freunden der Seefahrt so beliebten Buchreihe.
INENATOAL.m JAH6UC14
OLUFT UND UUUMFUHRT
Der Flottenkalender ist ein Spiegel der internationalen Seefahrt und berichtet
seit Jahrzehnten in bewährter Form über ein breites Spektrum maritimer
Themen.
Hier eine kleine Auswahl von interes-santen Beiträgen:
Seenafen Benin, Hocnseeschiffe für die Bin nenstadt
Auf Einsatztahrt mit Seenotkreuzer-Neubau NIS RANDERS
Pillau-Baltisi<
Vorn Dampfschlepper zum
Motorgüter-erscheinende
FlortenkaIend bietet für 19.80 Mark viel
LeSeStOft aus denunterschd
lichsten Bereichen von
Schiffahrt und Marine.Meeres forschu9.
Weefl
und HafenWirt5 Das seit
1901 erscheinende Jahrbuch
ist nunmehr
wer Glück hat, gewinntmit dem
Floffenka1ee
preisratse'
Reprints
-
für Sammler eine »Raritat«
Ca 240 Seiten, Format 14 5 X 21 cm. viele Abbil-dungen gebunden. ca. DM 39,80. Best.-Nr 5361.
Aus dem inhalt:
Kaiser Wilhelm lin Admiraisuniform Die bedeutenden Reedereien Deutschlands Ein Wal ais Schiffszerstörer
Humoristisches
Soeben
erschiet!f1
Das Jahrbuch «Fliegerkalender 1993» ist spannend zu lesen und bestens illustriert, bietet informa-tive Daten. Fakten, Berichte und lebendige Schilderungen über
Flugzeuge. die Fliegerei und die
R a um fahrt.
Ubereinstimmendes Urteil aus
dem Spektrum der Rezensionen:
Dieses empfehlenswerte
Jahr-buch hat deshalb einen so großen
Leserkreis gewonnen, weil es
wertfrei, informativ und unterhal-tend berichtet. Besonders
Samm-1er freuen sich über dieses
«Nach-schlagewerk'. in dem sie immer wieder interessante Ereignisse, Ergebnisse und Daten nachlesen
können.
Der alljährlich erscheinende
Flot-tenkalender, eine
«Pflichttektü-re», die bereits seit 1901 als das beliebte und bekannte Jahrbuch
der deutschen Seefahrt
er-scheint. ist unterhaltsam und in-formativ gestaltet.
Fotos und Abbildungen bieten
Jahr für Jahr zusammen mit der Flottenrevue einen interesssan-ten Uberblick über die Handels-schiffahrt und Marinen der Welt.
32052 Herford
REPRINT Köhlers Flottenkalender
Kalender des Deutschen Flotten-Vereins für 1901
216 Seilen, Frmai 4 5 X 21 cm. roil 4 Farn- und 130 historiscnen Schwarzweiß.Abbildungen bzw. Tabellen, Broschur. DM 39,80, Best.-Nr, 4901
REPRINT Köhlers Flottenkalender
