Numerical approach
3.2 Discrete boundary value problem
In this section a discretisation procedure for the boundary value problem (2.35) is em-ployed.
3.2.1 Discretisation of the integral equations
For the derivation of a numerical solution of the boundary value problem, one defines N collocation points xcolj , j = 1, . . . , N, where the integral equations are applied. A collocation pointxcolj is defined to be the centre of the panelj which is slightly displaced inside the body. For these collocation points the integral equation(2.35)can be used and the boundary value problem reads as:
1
From the above relation one obtains:
XN
Chapter 3. Numerical approach
are the influence functions that describe the dipole or source influence of the paneli on the pointxcolj . The influence coefficients in Equation(3.5)only depend on the geometry and the choice of the collocation points. In case of flat quadrilateral panels, they can be calculated analytically (s. Appendix C).
Applying Equation (3.4) on N collocation points leads to a linear system of equations that consist ofN +NW columns andN rows:
The discrete source strengths are known from the discrete form of Equation(2.27): σni =Vn∞,i·ni, ∀i= 1, . . . , N. (3.7) It follows that there areN+NW unknowns and onlyN equations. The unknowns are the body and wake dipole strengths. The determination of the wake dipole strengths differs in steady and unsteady computations:
• In steady computations the value of the dipoles on the discrete wake sheet is calcu-lated as follows:
1. On the first row of the trailing wake the value of the dipole strengthsµW,(1,k),
∀k = 1, . . . , Nkis obtained by using the Kutta condition.
2. At the following rows of the trailing wake the dipole values are set to the same value as the dipoles at the leading panel: µW,(l,k) =µW,(1,k),∀l= 2, . . . , NW,l,
∀k = 1, . . . , Nk(s. Figure 3.5).
• In unsteady computations the value of the dipole strengths at the wake panels changes for each time step and each radial stripe. In the first time step the values from the steady computations are used. In the following steps the procedure goes as follows:
Chapter 3. Numerical approach
1. On the first row of the trailing wake the value of the dipole strengthsµnW,(1,k),
∀k= 1, . . . , Nkis obtained by using the Kutta condition.
2. The dipole values for the remaining wake panels are set to the value from the previous time step: µnW,(l,k) = µnW,(l−1−1,k), n = 1, . . ., ∀l = 2, . . . , NW,l,
∀k= 1, . . . , Nk.
The number of unknown dipole strengths is the same in the steady and unsteady case, namelyNk. The numerical implementation of the Kutta condition used in the steady and unsteady computations for the determination of the dipole strengths on the trailing wake is introduced below.
3.2.2 Implementation of the Kutta condition
The aim of the Kutta condition is the determination of unknown wake dipole strengths µnW,(1,k) for each radial stripek such that the force-free condition∆p = p+−p− = 0is fulfilled on the wake sheet. The potential jump∆pis a non-linear function of the dipole strengths and cannot be solved directly. In order to solve this equation numerically, ei-ther a linear approximation or the Kutta condition can be used or the non-linear condition must be solved by an iterative procedure. The linear condition is easier to implement, but is less accurate than the non-linear form of the condition. The iterative solution of the non-linear equation, e.g. by means of Newton’s method, complicates the implemen-tation and increases the compuimplemen-tational effort. InpanMARE both solution approaches are implemented.
Linear form of the Kutta condition The linear form of the Kutta condition reads as:
µnW,(1,k) =µn(Nl,k)−µn(1,k), (3.8) where µn(N
l,k) and µn(1,k) are the dipole strengths on the upper and on the lower side of the trailing edge of the solid body, respectively. This condition means that a body panel will be influenced by the wake panels when it is located at the trailing edge. Now, the following substitution is made:
A∗i,j =Ai,jif the panelidoes not lie at the trailing edge, i.e.i=l+ (k−1)Nl and l 6=Nlandl 6= 1,
A∗i,j =Ai,j±
kNW,l
X
h=1+(k−1)NW,l
Ah,jif the panelilies at the trailing edge of the radial stripek, i.e.i=Nl+ (k−1)Nlori= 1 + (k−1)Nl,
(3.9)
Chapter 3. Numerical approach where the sign±stands for the suction and the pressure side, respectively. By means of this substitution, one obtains a linear system of equations of dimensionN ×N for the unknown body dipole strengthsµi,∀i= 1, . . . N:
Figure 3.5: Declaration of the dipole strengths on the wake sheet
Non-linear form of the Kutta condition
The non-linear form of the Kutta condition at each radial stripe of the lifting body reads as:
∆p(1,k) =p+(1,k)−p−(1,k) = 0, ∀k = 1, . . . , Nk, (3.11) wherep±is the pressure value at the upper and lower side of the wake sheet, respectively.
This non-linear system of equations can be solved iteratively by Newton’s method (s.
Kerwin et al., 1987, p. 111f.). The method computes from a starting solution, a new solution that is the root of the non-linear system and solves the linear system of equations:
J·x(m+1) =−f(m), m= 0,1,2. . . . (3.12) whereJis the Jacobian matrix and the vector functionsxandf are defined by:
x(m+1) =
Chapter 3. Numerical approach The solution algorithm reads as:
(1) The solution from the linear Kutta condition is used to have an initial guess for the potential jump at the trailing edge. The initial guess is denoted by x0 = (µ0W,(1,1), . . . , µ0W,(1,N
k)). For this initial guess the pressure jump ∆p0(1,k) is calcu-lated for allk = 1, . . . , Nk.
(2) The previous approximation is corrected byµ(m+1)W =µ(m)W −J−1·f(m). The Jaco-bian matrixJis given by:
Jjk = ∆pβ(1,k)−∆p0(1,k)
µβW,(1,j)−µ0W,(1,j), j, k = 1, . . . , Nk, (3.14) whereµβW,(1,j) = (1−β)µ0W,(1,j) forj = k and, µβW,(1,j) =µ0W,(1,j) forj 6= k. The constantβis usually set to0.01(see Kerwin et al., 1987, p. 111). In order to avoid the inversion of the Jacobian matrix, the linear system of equations (3.12)is first solved to obtainx(m+1). Then, the following relation is used:
µ(m+1)W =µ(m)W +x(m+1). (3.15)
(3) The iteration procedure is repeated as long as for one iteration stepm= 0,1,2, . . . the relationmaxNk
k=1(|∆p(m)(1,k)|)< εfor a given thresholdεis satisfied.
Newton’s method is based on the assumption that the initial guess is not too far off from the exact solution. The starting point is very important since it determines which root is found by Newton’s method, which it is not always the closest root. If the initial solu-tion found by the linear Kutta condisolu-tion is too far off from the actual solusolu-tion, Newton’s method can diverge.
3.2.3 Initial condition for the unsteady problem
In the unsteady case, an initial condition for the first time step must be specified. An initial condition is required because of the dependency of the governing equations on the time-variablet. As the initial condition the solution from the steady computation is used, i.e. the problem is solved for the steady case first. The solution achieved by the steady computation is then used as the initial value for the unsteady computations.
Chapter 3. Numerical approach
3.2.4 Boundary value problem for multiple bodies in the ow domain
It is possible to compute several bodies submerged in the flow domainΩ. Such a case is given for example when calculating propeller-induced pressure pulses on a ship hull. In this case, two bodies, namely the propeller and the ship hull, are contained in the flow domainΩ. In order to adapt the linear system of equations(3.6) for several bodies, the following additional notation is introduced:
M Number of body surfaces in the fluid domainΩ.
Nm Number of panels on the body surfacem= 1, . . . , M. NWm Number of panels on the wake surfacem= 1, . . . , M. N =PM
m=1Nm Overall number of body panels.
NW =PM
m=1NWm Overall number of wake panels.
Based on the Equation(3.4), the modified boundary integral equations forM bodies in the fluid domain read as:
XM m=1
hXNm
i=1
µniAi,j+
NWm
X
i=1
µnW,iAi,j−
Nm
X
i=1
σinBi,ji
= 0, for all xcolj , j = 1, . . . , N, (3.16) where the influence functions are defined in the same way as in Equation(3.5)for∀i = 1, . . . , N +NW.