Table Of ContentNASA/TM–2016–219003
Comparison of NACA 0012 Laminar Flow
Solutions: Structured and Unstructured
Grid Methods
R.C.Swanson
NASALangleyResearchCenter,Hampton,Virginia
S.Langer
DLR,DeutschesZentrumfu¨rLuft-undRaumfahrt
Institutfu¨rAerodynamikundStro¨mungstechnik,
Braunschweig,Germany
January 2016
TheNASASTIProgramOffice... inProfile
Sinceitsfounding,NASAhasbeendedicatedto • CONFERENCEPUBLICATION.
theadvancementofaeronauticsandspace Collectedpapersfromscientificandtechnical
science. TheNASAScientificandTechnical conferences,symposia,seminars,orother
Information(STI)ProgramOfficeplaysakey meetingssponsoredorco-sponsoredby
partinhelpingNASAmaintainthisimportant NASA.
role.
• SPECIALPUBLICATION.Scientific,
TheNASASTIProgramOfficeisoperatedby technical,orhistoricalinformationfrom
LangleyResearchCenter,theleadcenterfor NASAprograms,projects,andmissions,often
NASA’sscientificandtechnicalinformation. The concernedwithsubjectshavingsubstantial
NASASTIProgramOfficeprovidesaccesstothe publicinterest.
NASASTIDatabase,thelargestcollectionof
• TECHNICALTRANSLATION.English-
aeronauticalandspacescienceSTIintheworld.
languagetranslationsofforeignscientificand
TheProgramOfficeisalsoNASA’sinstitutional
technicalmaterialpertinenttoNASA’s
mechanismfordisseminatingtheresultsofits
mission.
researchanddevelopmentactivities. These
resultsarepublishedbyNASAintheNASASTI
SpecializedservicesthatcomplementtheSTI
ReportSeries,whichincludesthefollowing
ProgramOffice’sdiverseofferingsinclude
reporttypes:
creatingcustomthesauri,buildingcustomized
databases,organizingandpublishingresearch
• TECHNICALPUBLICATION.Reportsof
results... evenprovidingvideos.
completedresearchoramajorsignificant
phaseofresearchthatpresenttheresultsof
FormoreinformationabouttheNASASTI
NASAprogramsandincludeextensivedataor
ProgramOffice,seethefollowing:
theoreticalanalysis. Includescompilationsof
significantscientificandtechnicaldataand • AccesstheNASASTIProgramHomePageat
informationdeemedtobeofcontinuing http://www.sti.nasa.gov
referencevalue. NASAcounterpartof
peer-reviewedformalprofessionalpapers,but • E-mailyourquestionviatheInternetto
havinglessstringentlimitationsonmanuscript [email protected]
lengthandextentofgraphicpresentations.
• FaxyourquestiontotheNASASTIHelpDeskat
• TECHNICALMEMORANDUM. (301)621–0134
Scientificandtechnicalfindingsthatare
• PhonetheNASASTIHelpDeskat(301)
preliminaryorofspecializedinterest,e.g.,
621–0390
quickreleasereports,workingpapers,and
bibliographiesthatcontainminimal • Writeto:
annotation. Doesnotcontainextensive NASASTIHelpDesk
analysis. NASACenterforAeroSpaceInformation
7115StandardDrive
• CONTRACTORREPORT.Scientificand
Hanover,MD21076–1320
technicalfindingsbyNASA-sponsored
contractorsandgrantees.
NASA/TM–2016–219003
Comparison of NACA 0012 Laminar Flow
Solutions: Structured and Unstructured
Grid Methods
R.C.Swanson
NASALangleyResearchCenter,Hampton,Virginia
S.Langer
DLR,DeutschesZentrumfu¨rLuft-undRaumfahrt
Institutfu¨rAerodynamikundStro¨mungstechnik,
Braunschweig,Germany
NationalAeronauticsand
SpaceAdministration
January 2016
Theuseoftrademarksornamesofmanufacturersinthisreportisforaccuratereportinganddoesnotconstitutean
officalendorsement,eitherexpressedorimplied,ofsuchproductsormanufacturersbytheNationalAeronauticsand
SpaceAdministration.
Availablefrom:
NASACenterforAeroSpaceInformation(CASI) NationalTechnicalInformationService(NTIS)
7115StandardDrive 5285PortRoyalRoad
Hanover,MD21076–1320 Springfield,VA22161–2171
(301)621–0390 (703)605–6000
Abstract
In this paper we consider the solution of the compressible Navier-Stokes equations for a class of
laminar airfoil flows. The principal objective of this paper is to demonstrate that members of this
classoflaminarflowshavesteady-statesolutions. Theselaminarairfoilflowcasesareoftenusedto
evaluateaccuracy,stabilityandconvergenceofnumericalsolutionalgorithmsfortheNavier-Stokes
equations. In recent years, such flows have also been used as test cases for high-order numerical
schemes. Whilegenerallyconsistentsteady-statesolutionshavebeenobtainedfortheseflowsusing
higher order schemes, a number of results have been published with various solutions, including
unsteadyones. Wedemonstratewithtwodifferentnumericalmethodsandarangeofmesheswith
amaximumdensitythatexceeds8×106 gridpointsthatsteady-statesolutionsareobtained. Fur-
thermore, numerical evidence is presented that even when solving the equations with an unsteady
algorithm,oneobtainssteady-statesolutions.
Contents
1 Introduction 2
2 GoverningEquations 3
2.1 PhysicalBoundaryConditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
3 StructuredRK/ImplicitScheme 4
3.1 ElementsofSolutionAlgorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.2 DiscreteBoundaryConditions: StructuredScheme . . . . . . . . . . . . . . . . . 7
4 UnstructuredRK/ImplicitScheme 8
4.1 DescriptionofSteady-StateSolverforUnstructuredScheme . . . . . . . . . . . . 8
4.2 DescriptionofUnsteadySolverforUnstructuredScheme . . . . . . . . . . . . . . 9
4.3 DiscreteBoundaryConditions: UnstructuredScheme . . . . . . . . . . . . . . . . 11
5 LaminarAirfoilFlowResults 11
5.1 ComputationalGrids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
5.2 Subsonicflows(Cases1–4) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
5.2.1 Numericalsolutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
5.2.2 Possiblesourcesofdifferencesinsolutions . . . . . . . . . . . . . . . . . 14
5.2.3 Effectofouterboundarylocation . . . . . . . . . . . . . . . . . . . . . . 15
5.2.4 Effectofskin-frictionmethods . . . . . . . . . . . . . . . . . . . . . . . . 16
5.2.5 Gridconvergenceofaerodynamiccoefficients . . . . . . . . . . . . . . . . 16
5.3 Transonicflow(Case5) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
6 ConcludingRemarks 18
1
1 Introduction
Inthedevelopmentandevaluationofnumericalschemesitisessentialtohaveacomprehensivesetof
testproblems. Suchproblemscanallowspecificpropertiesofthesolutionalgorithmtobevalidated.
Furthermore,theycanpossiblydelineatedeficienciesinthealgorithmthatmustbeovercome.These
test cases should provide either an analytic solution or a fully converged and resolved numerical
solution. Whileexperimentaldatacanoftenprovideconfirmationforthevalidityofthecomputed
solution,oneshouldnotnecessarilyconcludethateverythingisworkingcorrectlyforasteady-state
problem just because there is reasonable agreement with the data according to some engineering
criterion(suchasasmallvariationintheaerodynamiccoefficients)whentheresidualhasonlybeen
reduced a few orders of magnitude (rather than machine zero). An appropriate evaluation should
involveahierarchyoftestcases. Forexample,ifweconsidersolvingtheNavier-Stokesequations,
thetestcasesshouldincludelaminarflows. Laminarflowtestcaseshavethedistinctadvantagethat
theyeliminatetheneedatsufficientlyhighReynoldsnumberstoeitherresolveturbulenceinaflow
fieldortomodeltheturbulenceappropriately.
SomepossiblelaminarflowtestcasesinvolvingtheNACA0012airfoilwereconsideredbythe
firstauthorofthispaperin1984. Oneofthesetestcases(Machnumber(M)of0.5,angleofattack
(α)of0◦,andReynoldsnumber(Re)equalto5000)wasintroducedinthe1985paperbySwanson
andTurkel[1]toevaluateanalgorithmforsolvingthecompressibleNavier-Stokesequations. This
particularcasewaschosenbecauseithasasmallamountoftrailingedgeseparation(beginningat
approximately the 81% chord location), and thus, represented a good way to check the levels of
dissipationbeingproducedbythenumericalscheme. Forexample,iftheschemeistoodissipative,
theeffectiveReynoldsnumberisreduced, andtheseparationpointmovesdownstream. Sincethat
time,thiscasehasbeenconsideredinnumerousevaluationsofschemes,suchasRefs.[1–7]. Ithas
alsobeenusedinastudyofvortex-airfoilinteractionbySv¨ardetal.[8]. Inaddition,thiscaseand
anotheroneoftheoriginalcases(α = 3◦)havealsobeenusedbyVenditti[9]inevaluatingagrid
adaptiveschemeforfunctionaloutputsofflowsimulations.
In recent years, some additional laminar flow cases for the NACA 0012 airfoil have also been
considered. TheflowconditionsarethesameasfortheoriginalcasesofSwansonexceptα = 1◦
orα=2◦. Thesecases,aswellastheoriginalcases,havebeenconsideredintheEuropeanproject
ADIGMAtodevelopadaptivehigh-ordervariationalmethodsforaerospaceapplications[10]. Due
to the strong interest in these laminar flow cases, there is a need to have documentation of their
solutions.
Somerecenteffortstocomputemembersofthisclassofflowproblemshaveproducedsolutions
thatareinconsistentwiththeratherlargenumberofprevioussolutions(computedwithawiderange
of numerical methods), including those already cited in this paper. For example, Abgrall and De
Santis[11]showunexpectedsolutionsfortheα=0◦case,whichshouldhaveasymmetricsolution
basedonpreviousresults, usingbothsecond-orderandthird-orderdiscretizations. Onepossibility
for such results is the presence of an asymmetry in their solution algorithm. Another example
concerns the unsteady solution for the α = 2◦ case presented in the ADIGMA Project [10] paper
byTaubeetal. startingonpage427. Eventhoughthecomputationalmethodincludedadaptation,
thecalculationwasperformedonarathercoarsemesh(815cells)andtheadaptationseemstohave
beenfocusedinthewake. Thus,therecouldhavebeeninsufficientresolutionintheneighborhood
oftheseparationpoint,whichcanproducearesultthatappearstobeanunsteadysolution. Thereis
alsothepossibilitythatthesolutionalgorithmisnotsufficientlystrong,asmeasuredbythelevelof
implicitness. Inthispaper,wewillshowhowaweaksolutionalgorithmcanproducewhatappears
to be an unsteady solution. Dolejs˘´ı [12] has also obtained an unsteady solution for this laminar
flow case. A nonlinear system of algebraic equations (from discontinuous Galerkin finite element
method) is solved with an inexact Newton-time method. An inner iterative method is required to
solve a linear system at each time level. Since the residual (algebraic error) of the iteration is not
being reduced to approximately zero, there can be sufficient error to produce an unsteady result.
In any of these examples, one cannot dismiss the possible influence of boundary conditions. It is
2
difficulttoassessthisbecausedescriptionsofthediscreteimplementationoftheboundaryconditions
fortheproblemsbeingconsideredarenotgiven.
As a consequence of such examples, a primary objective of this paper is to demonstrate that
there are steady-state solutions, even on high density meshes, to this class of problems. Two ad-
ditional objectives of this paper are as follows. The first of these is to examine and document the
behaviorofasetoflaminarflowproblems,whichincludesthosejustdiscussed,aswellasacom-
monlyconsideredlowerRecasewithahigherangleofattack. Thesecondadditionalobjectiveisto
comparetwomethodsforsolvingthecompressibleNavier-Stokesequations. Onemethod[13,14]
isbasedonstructuredgrids,andtheother[15]isbasedonunstructuredgrids. Bothmethodsapply
afinite-volumeapproachforspatialdiscretization. Thestructuredgridmethodusesacell-centered
formulation,andtheunstructuredgridmethodusesanode-centeredformulation. Eachmethodap-
pliesamatrixformfornumericaldissipation. Inaddition,usingthestructuredgridapproachaRoe
dissipation (which is frequently used in numerical schemes) is also considered. When comparing
whatwecallstructuredandunstructuredschemesitisimportanttoclarifythemeaningofthister-
minologywithrespecttocontrastingtheschemes. Inthepaper,thesamesetofstructuredgrids(i.e.,
quadrilateralelements)isusedforbothmethods. Thus,thecomparisonoftheschemesisbasedon
thefactthattherearedifferencesintheirelementsandimplementation.
For all computations in the paper we solve the full Navier-Stokes equations (i.e., all physical
diffusion terms are retained). The solvers for the two methods are comprised of a Runge-Kutta
(RK) scheme with an implicit preconditioner (herein also designated as RK/Implicit scheme) to
extendstabilityandallowalargeCourant-Friedrichs-Lewy(CFL)number. Amultigridschemeis
employed for convergence acceleration. In the comparison of the methods, we consider accuracy,
stability,andconvergencewhencomputingsolutionsofthelaminarflowcases.
In the first section of the paper, we present the governing flow equations and the boundary
conditions for the continuous problem of flow past an airfoil. The next two sections describe and
discuss the structured and unstructured methods, including the discrete boundary conditions, that
areconsidered. Afterdelineatingthespecificdifferencesinthetwomethods,wethenpresentaset
oflaminarflowresultsthatcomparethemethods. Convergencehistoriesareshownforbothsteady-
stateandunsteadycomputations.Comparisonsaremadebetweentheresultswiththestructuredand
unstructuredmethods. Theeffectsofmeshdensityonthesolutions, whichcontainbetween8,192
and about 8.39×106 cells, as well as the surface pressure and skin-friction distributions for each
laminar flow case are presented and discussed. In all cases, there is flow separation, and both the
separationlocationanddownstreamextentofseparationarecompared. Overallflowdetailsarealso
illustratedwithMachandstreamlineplots. Suchdetailsprovideareferenceforotherinvestigators
tocompareagainst.
2 Governing Equations
Weconsiderthetwo-dimensional(2-D)Navier-Stokesequationsforcompressibleflow. Assuming
avolumefixedinspaceandtime,theintegralformoftheseequationscanbewrittenas
ZZ ∂W Z
dV + F ·ndS =0, (1)
∂t
V S
whereWisthestatevectorofconservativevariables,F isthefluxdensitytensor,andV,S,andn
denotethevolume,surface,andoutwardfacingnormalofthecontrolvolume. Onecansplittheflux
densitytensorintoaconvectivecontributionF andaviscouscontributionF ,whicharegivenby
c v
ρq 0
Fc = ρρuvqq++ppeeyx , Fv = ττ¯¯··eexy (2)
ρHq τ¯·q−Q
3
where q is the velocity vector with Cartesian components (u,v), and the unit vectors (e ,e ) are
x y
associatedwiththeCartesiancoordinates(x,y). Thevariablesρ,p,H representdensity,pressure,
andtotalspecificenthalpy,respectively. Thestresstensorτ¯andtheheatfluxvectorQaregivenby
(cid:20) (cid:21) (cid:20) (cid:21)
τ τ ∂T/∂x
τ¯= xx xy , Q=k (3)
τ τ ∂T/∂y
yx yy
with
(cid:18) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19)
∂u ∂v ∂u ∂u ∂v ∂v ∂u ∂v
τ =λ + +2µ , τ =λ + +2µ , τ =τ =λ + .
xx ∂x ∂y ∂x yy ∂x ∂y ∂y xy yx ∂y ∂x
(4)
Here, the symbol ∂ indicates partial differentiation, µ and λ are the first and second coefficients
ofmolecularviscosity,k denotesthecoefficientofthermalconductivityandT representsthetem-
perature. By the Stokes hypothesis, λ = −2/3µ, and µ is determined with Sutherland’s viscosity
law[16].ThethermalconductivitykisevaluatedwiththeconstantPrandtlnumberassumption(i.e.,
thePrandtlnumberPr =0.72).
InordertoclosethesystemgivenbyEq.(1),weusetheequationofstate
p=ρRT (5)
whereRisthespecificgasconstant.
2.1 PhysicalBoundaryConditions
Inthecontinuumcase, weconsideranexternalflowproblem(i.e., flowpastagivengeometry)on
an infinite domain Ω ⊂ R2. Thus, we need to define appropriate conditions at a wall boundary,
whichweassumetobesolid. Later,inthediscretecase,wedefineafinitedomain. Itwillthenbe
necessarytodefinesuitableinflowandoutflowboundaryconditions.
Forviscousflows,thenon-penetrationconditionandtheno-slipconditionareimposed,
q·n=0, q·t=0, (6a)
where t is the unit vector tangent to the surface of a control volume. In addition, a boundary
conditionisrequiredtodeterminethetemperatureandoneotherthermodynamicvariable. Forthe
temperature,weimposetheadiabaticcondition
Q·n=0, (6b)
where Q is the heat flux vector given in (3). Usually, for the other thermodynamic variable re-
quired, the pressure is selected, which is determined from the normal momentum equation. For
example,considerthedifferentialformofthenormalmomentumequationatasolidandstationary
wallboundary. Thisequationcanbewrittenas
(∇p·n) =[(∇·τ)·n] , (6c)
w w
where∇isthegradientoperator.
3 Structured RK/Implicit Scheme
Weapply afinite-volume approachto discretizethe fluiddynamicequations anduse athree point
second-orderdiscretizationoftheadvectiontermswitheitherthematrix-typeorRoe-typenumerical
dissipation. For details of these forms of numerical dissipation see Swanson and Turkel [17] and
4
Roe[18]. Theviscoustermsarediscretizedwithasecond-ordercentraldifferenceapproximation.
Usingthisdiscretization,wedescribeinthissectionthesolutionalgorithmforthestructuredscheme.
ThenumericalsolvercalledtheRK/Implicitschemeisamemberofageneralclassofschemes
with a nonlinear outer iteration and a linear inner iteration (e.g., see [15,19]). In the outer iter-
ation, a Full Approximation Storage (FAS) multigrid method with a RK smoother such as those
of Ref. [14] is applied. The inner iteration determines an approximate implicit preconditioner for
the RK scheme, resulting in a type of implicit RK smoother. The particular implementations for
the structured grid and unstructured grid methods are described below in Section 3.1 and later in
Section4.1,respectively.
3.1 ElementsofSolutionAlgorithm
A semi-discrete form of Eq. (1) can be written as (omitting subscripts of discrete quantities for
convenience),
dW
+R(W)=0, (7)
dt
definingasystemofordinarydifferentialequations. TheR(W)denotestheresidualfunctiondeter-
minedbythespatialdiscretizationbrieflydescribedpreviously. Applyingathree-stageRKscheme
toupdatethediscretesolution,weobtainontheq-thstage
W(q) =W(0)+δW(q), (8)
wherethechangeinthesolutionvectorWis
∆t
δW(q) =W(q)−W(0) =−α LW(q−1) =−α ∆tR(W), (9)
q V q
and L is the complete difference operator for the system of equations, and V is the volume of the
meshcellbeingconsidered. Here,α istheRKcoefficientoftheq-thstage,and∆tisthelocaltime
q
step. To extend the support of the difference scheme, we consider an implicit residual smoothing
approach[20]. Adaptingthesmoothingtechniquetoasystemofequations,wehavethefollowing:
L δW(q) =δW(q), (10)
i
whereL isanimplicitoperator. ByapproximatelyinvertingtheoperatorL ,weobtain
i i
δW(q) =−α ∆tPLW(q−1) =−α ∆tP X F(q−1)S, (11)
q V q V n
allfaces
whereP isapreconditionerdefinedbytheapproximateinverseLe−1,andthesummationdenotesa
i
fluxbalanceonameshcell. WecanalsoviewEq.(11)asEq.(9)withR(W)replacedbyPR(W),
(q)
whichwewilldenoteasR(W).ThechangeδW replacestheexplicitupdateappearinginEq.(8).
Thus,eachstageintheRKschemeispreconditionedbyanimplicitoperator.
AfirstorderupwindapproximationbasedontheRoeschemeisusedfortheconvectivederiva-
tives in the implicit operator. To derive this operator, one treats the spatial discretization terms in
the flow equations implicitly and applies linearization. For a detailed derivation see Rossow [13].
SubstitutingfortheimplicitoperatorinEq.(10),weobtainfortheq-thstageoftheRKscheme
" #
I+ε∆Vt X AnS δW(q) =−αq∆Vt X F(nq−1)S =−αq∆tRb(q−1), (12)
allfaces allfaces
wherethematrixA isthefluxJacobianassociatedwiththenormalfluxdensityvectorF atacell
n n
face,Rb(q−1)representstheresidualfunctionforthe(q−1)-thstage,andεisarelaxationparameter.
Theparameterεistakentobe0.5.
5
ThematrixA canbedecomposedintoA+andA−,whicharedefinedby
n n n
1 1
A+ = (A +|A |), A− = (A −|A |). (13)
n 2 n n n 2 n n
IfwesubstituteforA inEq.(12)usingthedefinitionsofEq.(13),thentheimplicitschemecanbe
n
writtenas
" #
I+ε∆Vt X A+n S δW(i,qj) =−αq∆tRb(i,qj−1)−ε∆Vt X A−nδW(NqB) S, (14)
allfaces allfaces
where the indices (i,j) indicate the cell of interest, and NB refers to all the direct neighbors of
the cell being considered. Consider traversing the boundary of a mesh cell in a counterclockwise
manner, with the positive surface normal vector always pointing outward from the cell. Then, as
discussedbyRossow[21],thequantityA−δWrepresentsthefluxdensitychangeassociatedwith
n
waves having a negative wave speed (i.e., waves that enter the cell (i,j) from outside). Only the
neighborcellsNB cancontributetothesechangesinfluxdensity. Similarly, thequantityA+δW
n
representsfluxdensitychangesassociatedwithpositivewavespeeds(i.e.,wavesthatleavethecell
(i,j)). Thesefluxdensitychangesaredeterminedonlybyinformationfromwithinthecell(i,j).
In order to efficiently evaluate the Jacobian matrices A+ and A−, we rely upon their forms
n n
when expressed in terms of the cell interface Mach number M . For simplification, we transform
0
Eq.(12)totheprimitivevariablesU=[ρ p u v]T. Thus,
" #
I+ε∆t X P S δU(q) =−∂Uα ∆t X F(q−1)S, (15)
V n ∂W q V n
allfaces allfaces
wherethematrixP , whichistheanalogofthenormalfluxJacobianexpressedinprimitivevari-
n
ables,isgivenby
P = ∂UA ∂W = ∂U (cid:0)A++A−(cid:1)∂W =P++P−. (16)
n ∂W n ∂U ∂W n n ∂U n n
TheJacobian∂U/∂Wontheright-handsideofEq.(15)mustmultiplytheconservativefluxbalance
inordertoensureconservation. UsingthedefinitionsofEq.(13)andthedissipationmatrix,which
isdefinedintheappendixofRef.[14],onecandeterminethematricesP+ andP−,whicharealso
n n
giveninthatappendix. Theresultingmatricescaneasilyberecomputed,onlyrequiringstoragefor
the normal velocity magnitude and the Mach number M . The contributions of the viscous flux
0
Jacobians can be included in a straightforward manner using primitive variables (see Ref. [22]).
Note that δU is simply proportional to a new residual R(U) that requires transforming back to
R(W).
(q)
TosolveEq.(14)forthechangesinconservativevariablesδW ,the4×4matrixontheleft-
i,j
hand side of Eq. (14) must be inverted. It is sufficient to approximate the inverse of the implicit
operator. An adequate approximate inverse is obtained with two symmetric Gauss-Seidel sweeps.
Therelaxationprocedureuseslexicographicorderingofthesolutionpoints.Toinitializetheiterative
process,theunknownsaresettozero.
For the RK/Implicit scheme, the relaxation is a combination of point and line SGS relaxation,
which is the method frequently used with unstructured schemes. With line SGS, we apply the
implicitoperatoracrossthewakeline(thelineemanatingfromthetrailingedgeofanairfoil),which
wecallafullyimplicittreatment. Itisalsopossibletolagtheboundaryconditionatthewakeline
ratherthanextendingthelinesolveacrossthewakeline.Thisapproachcanalsoworkwellprovided
localrelaxationorsubiterationsareperformedtocompensateforthelaggingofthewakeboundary
condition. Althoughthismaybeamoreconvenientoptionfromtheimplementationpointofview
formultiblockstructuredgridmethods,weusethethefullyimplicitmethodinthispapertoobtain
6
