Coronado-Arora-GLS4

13 Pages·2006·1.07 MB·English

J.Non-NewtonianFluidMech.xxx (2006)xxx–xxx
Four-field Galerkin/least-squares formulation for viscoelastic fluids
Oscar M. Coronadoa, Dhruv Aroraa, Marek Behra,b,∗, Matteo Pasqualia,∗
aDepartmentofChemicalandBiomolecularEngineeringandComputerandInformationTechnologyInstitute,RiceUniversity,
MS362,6100MainStreet,Houston,TX77005,USA
bChairforComputationalAnalysisofTechnicalSystems,CCES,RWTHAachenUniversity,52056Aachen,Germany
Received4November2005;receivedinrevisedform14February2006;accepted24March2006
Abstract
AnewGalerkin/Least-Squares(GLS)stabilizedfiniteelementmethodispresentedforcomputingviscoelasticflowsofcomplexfluidsdescribed
bytheconformationtensor;itextendsthewell-establishedGLSmethodforcomputingflowsofincompressibleNewtonianfluids.GLSmethods
areattractiveforlarge-scalecomputationsbecausetheyyieldlinearsystemsthatcanbesolvedeasilywithiterativesolvers(e.g.,theGeneralized
Minimum Residual method) and because they allow simple combinations of interpolation functions that can be conveniently and efficiently
implementedonmoderndistributed-memorycache-basedclusters.
Likeotherstate-of-the-artmethodsforcomputingviscoelasticflows(e.g.,DEVSS-TG/SUPG),thenewGLSmethodintroducesaseparatevariable
torepresentthevelocitygradient;withtheaidofthisvariable,theconservationequationsofmass,momentum,conformation,andthedefinition
ofvelocitygradientareconvertedintoasetoffirst-orderpartialdifferentialequationsinfourunknownfields—pressure,velocity,conformation,
andvelocitygradient.Theunknownfieldsarerepresentedbylow-order(continuouspiecewiselinearorbilinear)finiteelementbasisfunctions.
ThemethodisappliedtotheOldroyd-Bconstitutiveequationandistestedintwobenchmarkproblems—flowinaplanarchannelandflow
pastacylinderinachannel.Resultsshowthat(1)themesh-convergencerateofGLSiscomparabletotheDEVSS-TG/SUPGmethod;(2)theLS
stabilizationpermitsusingequal-orderbasisfunctionsforallfields;(3)GLShandleseffectivelytheadvectivetermsintheevolutionequationof
theconformationtensor;and(4)GLSyieldsaccurateresultsatlowercomputationalcoststhanDEVSS-typemethods.
©2006ElsevierB.V.Allrightsreserved.
Keywords: Stabilizedfiniteelementmethod;Viscoelasticflow;Galerkin/least-squares;Oldroyd-Bfluid;Flowpastacylinderinachannel
1. Introduction of models have been proposed for modeling complex fluids:
fine-grainedmodels[1,2],(e.g.,bead-springorbead-rodmodels
Inthepastdecades,extensiveresearchhasbeendoneonflows ofpolymersolutions),wherethemicrostructureisrepresented
ofliquidswithmicro–macrostructure(alsoknownascomplex bymicromechanicalobjectsgovernedbystochasticdifferential
fluids);thesefluidsarefoundinseveralindustrialandbiological equations, and coarse-grained ones, where the microstructure
applications, e.g., polymer processing, coating of polymer ismodeledbymeansofoneormorecontinuumvariablesrep-
solutions, ink-jet printing, microfluidic devices, and human as resenting the expectation value of microscopic features (e.g.,
well as artificial organs (blood, synovial fluid). Usually these theconformationtensorinmodelsofpolymersolutions)[3–6].
liquids display a viscosity dependent on the rate of straining Fine-grained models incorporate a richer degree of molecular
and the flow kinematics (shear versus extension); they also details, but are still limited to fairly simple flows because of
show elasticity on time scales that overlap with the flow time computationalcost[7–9].
scales. Coarse-grained models represent the liquid microstructure
Realistic models of flowing complex fluids are crucial for intermsofoneormoreconformationtensors;currently,these
understandingandoptimizingflowprocesses.Twomainclasses modelsareconsideredthemostappropriateforlarge-scalesim-
ulationofcomplexflowsofcomplexfluids.Typically,thecon-
formationtensorobeysahyperbolicpartialdifferentialtransport
equation.Inpolymersolutionsandmelts,thistensorrepresents
∗ Correspondingauthors.Tel.:[phone];fax:[phone].
the local expectation value of the polymer stretch and orien-
E-mailaddress:[email protected](M.Pasquali),[email protected]
(M.Behr). tation, e.g., gyration or birefringence tensor. The elastic part
0377-0257/$–seefrontmatter©2006ElsevierB.V.Allrightsreserved.
doi:10.1016/j.jnnfm.2006.03.016
JNNFM-2595; No.ofPages13
2 O.M.Coronadoetal./J.Non-NewtonianFluidMech.xxx (2006)xxx–xxx
of the stress is related to the conformation tensor through an ined.Thismethodhasbeenrefinedandextendedmorerecently
algebraicequation[3–6].Suchmodelsincludemost“classical” toimproveconsistencybyrecoveryofthevelocitygradientas
rate-typestress-baseddifferentialmodels(e.g.,Oldroyd-B,PTT, well as a more appropriate expression of the LS stabilization
Giesekus,etc.)[3–6]. coefficient[29].
Simulations of complex flows of complex fluids require Fanetal.[30]independentlyintroducedanincompleteGLS
solving simultaneously the hyperbolic transport equation of methodforviscoelasticflowandtesteditsperformanceinaflow
conformation (or rate-type equation for the stress) together betweeneccentriccylinders,flowaroundasphereinapipe,and
with the momentum and mass conservation equations; this flowaroundacylinderinachannel.Thismethoddidnotinclude
poses several numerical challenges. In particular, obtaining termsduetotheLSformofthemomentumequation(because
mesh-converged solutions in simple benchmark flows at high it degraded performance) and of the constitutive equation;
Weissenberg number (Wi, the product of characteristic strain therefore,themethodofFanetal.[30]isbettercharacterizedas
rate and fluid relaxation time) is still considered an open apressure-stabilizedSUPGmethod—see[31]foradescription
problem. of pressure-stabilized methods for incompressible Newtonian
The Galerkin method is perhaps the most effective method flows.
forflowswithfreesurfacesanddeformableboundaries.How- ThisarticlepresentsacompleteGLSmethodforcomputing
ever, the Galerkin method is unstable in advection-dominated flows of incompressible viscoelastic liquids modeled by con-
problems, and yields spurious oscillations in the variable formationtensororrate-typeequations.Theflowequationsare
fields. Alternative methods have been developed to handle converted to a set of four first-order partial differential equa-
advection-dominatedaswellaspurelyhyperbolicequations— tionsbyrepresentingexplicitlythevelocitygradienttensor(as
e.g., Streamline upwind/Petrov–Galerkin (SUPG) for high in DEVSS-G). The GLS weighted residual equations include
ReynoldsnumberNewtonianflows[10]andviscoelasticflows naturallytheconsistentstreamlineupwindingfortheadvective
[11], also Discontinuous Galerkin (DG) for viscoelastic flows terms in the conformation evolution equation (and in the mo-
[12]. mentumequation,althoughthepresentationbelowisrestricted
When the Galerkin (or SUPG) method is applied to cou- toinertialessflows).Thechoiceofbasisfunctionsforthefour
pled partial differential equations, the selection of the interpo- unknownfields(velocity,pressure,velocitygradient,andcon-
lating functions for the various unknowns can be restricted by formation) is not restricted by compatibility conditions; here,
compatibilityconditions—e.g.,theBabusˇka-Brezziconditionin theunknownfieldsarerepresentedbythesimplestpossiblefi-
flowsofincompressibleNewtonianfluids[13,14].Somecom- niteelementbasisfunctions—continuouspiecewisebilinearon
patibility conditions between the basis functions of velocity, quadrilateral elements. The method is termed GLS4 to distin-
pressure, velocity gradient, and conformation (or stress) must guishitfromthepreviousGLS3[25,29]method,inwhichthe
still be satisfied [15,16] by current Galerkin-type methods for velocity gradient was not represented explicitly. The accuracy
simulatingviscoelasticflows—e.g.,thestate-of-the-artDEVSS- andstabilityofthemethodisdemonstratedbyusingtwobench-
TG/SUPG,whichevolvedfromsuccessivemodificationsofthe markproblems—theflowinaplanarchannelandtheflowpast
EVSS method [17–22] (see also reviews by Baaijens [23] and acylinderinachannel—foranOldroyd-Bfluid.
OwensandPhillips[24]). Itisworthnotingthatrecentwork[32–34]identifiedanother
Thesetwokeyhurdles(handlingadvection-dominatedprob- sourceofinstabilityinlow-orderfinitedifferenceandfiniteel-
lem and satisfying compatibility conditions) have been over- ementmethodsforcomputingviscoelasticflows—namely,the
come in Newtonian flows by using Galerkin/Least-Squares inability of low-order methods to capture exponentially grow-
(GLS)methods[25–27].WorkonGLSmethodsappliedtoNew- ingprofilesofconformationorelasticstressinregionsofstrong
tonian flows has shown that Streamline-upwind terms appear flow.Suchinstabilitycanbeavoidedbyusingthelogarithmof
naturallyintheGLSform,thatequal-orderbasisfunctionscan the conformation tensor as field variable [32], which has the
beusedforallfields(becausetheLeast-Squares(LS)termsre- additional benefit of ensuring that the conformation tensor is
movethecompatibilitycondition),andthattheresultingnonlin- automaticallypositivedefiniteeverywhereintheflow.Thepro-
earalgebraicequationsyieldaJacobianmatrixthatcanbesolved posedGLS4methoddoesnotaddressthissourceofinstability
moreeasilywithpreconditionedGeneralizedMinimumResid- explicitly. However, as discussed in Ref. [32], the logarithmic
ual method (GMRES) (because the LS terms yield a positive- changeofvariableisgenerallyapplicabletoanyfiniteelement
definiteJacobiancomponent).Moreover,usingequal-orderba- method(see,e.g.[34]);thus,itshouldbepossibletocombinethe
sisfunctionsforallfieldsallows“nodal”(ratherthan“elemen- currentGLS4formulationwiththelog-conformationmethodto
tal”)accounting,whichspeedsupgreatlymatrixoperationson improvethemethodfurther.
distributedmemoryparallelmachines[28].
WeaklyconsistentformsofGLSmethodhavebeenappliedto 2. Governingequations
viscoelasticflows.Behr[25]introducedathree-field(velocity–
pressure–elastic stress) GLS method and studied the flow of The steady flow of an inertialess incompressible viscoelas-
anOldroyd-Bliquidina4-to-1contraction.However,adetailed tic fluid, occupying a spatial domain Ω, with boundary Γ is
comparisonbetweenthismethodandotherpublishedresultswas governedbythemomentumandcontinuityequations:
notperformed,andtheeffectoftheexpressionoftheLSstabi-
lizationcoefficientfortheconstitutiveequationwasnotexam- ∇·T=0, onΩ, (1)
O.M.Coronadoetal./J.Non-NewtonianFluidMech.xxx (2006)xxx–xxx 3
∇·v=0, onΩ, (2) Boundaryconditionsonthemomentumequationareneeded
ontheentireboundaryΓ =Γ ∪Γ .Theessentialandnatural
g h
wherev istheliquidvelocityandTisthestresstensor,which boundaryconditionsarerepresentedas
canbedecomposedintoaconstitutivelyundeterminedisotropic v=g onΓ , (12)
g
contributionrelatedtoincompressibility,andviscousandelastic
contributions: n·T=h onΓh, (13)
T=−pI+τ+σ, (3) where g and h are given functions, and n is the outward unit
vector normal to the boundary. Because the equation of trans-
respectively,wherepisthepressure,Itheidentitytensor,τ = portofconformationishyperbolic,boundaryconditionsonthe
2η Dtheviscousstress(usuallyduetosolventcontribution),η conformationtensor,representedbythetensorG,areimposed
s s
thesolventviscosity,andDistherate-of-straintensor,i.e.,the atinflowboundariesΓ wherev·n<0,
G
symmetricpartofthevelocitygradient.Inordertotransformthe M=G onΓ . (14)
G
equations of motion into a set of first-order partial differential
equations(necessaryfordevelopingaconsistentLSformulation 3. Four-fieldGalerkin/least-squaresformulation(GLS4)
forlow-orderelements),anadditionalvariableLisintroduced
torepresentthevelocitygradient: Inthissection,theGLSformulationofthegoverningequa-
tions(7)–(10)ispresented.ThemethodistermedGLS4because
1
L=∇v− (∇·v)I, (4) theequationsethasfourbasicunknownfields—v,p,LandM.
trI
Thebasis(interpolation)andweightingfunctionspacesare:
wheretrdenotestrace.
ThelastterminEq.(4)ensuresthatLremainstracelesseven Shv ={vh|vh ∈[H1h(Ω)]nsd, vh ≡ghonΓg}, (15)
i(nLt+heLfiTn)/it2e.-pInretchiesiOonldrsooyludt-iBonm[o2d2e]l;,wthietheltahsitsicdsetrfiensistiiosnr,elDate≡d Vhv ={vh|vh ∈[H1h(Ω)]nsd, vh ≡0onΓg}, (16)
to the dimensionless conformation tensor M through a simple Shp =Vhp ={ph|ph ∈H1h(Ω)}, (17)
linearrelationshipσ =G(M−I),whereG=η /λistheelastic
modulus,ηp thepolymercontributiontothevispcosity,andλis ShL =VhL ={Lh|Lh ∈[H1h(Ω)]n2sd}, (18)
therelaxationtime.Theconformationtensorobeysahyperbolic Sh ={Mh|Mh ∈[H1h(Ω)]ntc, Mh ≡G onΓ }, (19)
M G
evolutionequation:
Vh ={Mh|Mh ∈[H1h(Ω)]ntc, Mh ≡0 onΓ }, (20)
∇ M G
λM+(M−I)=0, (5) where H1h represents functions with square integrable first-
∇ order derivatives, nsd the number of spatial dimensions and
whereMdenotesanupper-convectedderivative: n =n (n +1)/2 is the number of independent conforma-
tc sd sd
∇ tion tensor components. Bilinear piecewise continuous func-
M=v·∇M−LT·M−M·L. (6) tionsareusedhereafter.TheGLS4formulationis:Findvh ∈Sh,
v
p(cid:1)h ∈Sh,Lh ∈Sh an(cid:1)dMh ∈Sh suchthat:
Theequationsgoverningtheflowcanberecastindimensionless p L M
formas: ∇wh :ThdΩ+ wh·hhdΓ
∇∗·T∗ =0, (7) Ω ⎡ Γh ⎤
(cid:1)
∇∗·v∗ =0, (8) + τmom⎢⎢⎣∇qh−β∇·(Eh+(Eh)T)− 1−β∇·Sh⎥⎥⎦·
1 Ω (cid:5)Wi(cid:6)(cid:7) (cid:8)
L∗−∇∗v∗+ (∇∗·v∗)I=0, (9) A
trI (cid:1)
∇ [−∇·Th]dΩ+ qh(∇·vh)dΩ
WiM+(M−I)=0, (10) Ω
(cid:1)
wherev∗ =v/vc,p∗ =p/(ηvc/lc)andL∗ =L/(vc/lc)aredi- + τcont(∇·wh)(∇·vh)dΩ
mensionlessvelocity,pressureandinterpolatedtracelessveloc- Ω
(cid:1) (cid:12) (cid:13)
ity gradient tensor, respectively. ∇∗ =∇lc is the dimension- + Eh : Lh−∇vh+ 1 (∇·vh)I dΩ
less gradient operator, vc is a characteristic velocity, and lc is Ω trI
acharacteristiclength.ThedimensionlessWeissenbergnumber (cid:1) (cid:12) (cid:13) (cid:14)
isWi=λ(vc/lc).ThedimensionlessstresstensorT∗is + τgradv Eh−∇wh+ 1 (∇·wh)I : Lh−∇vh
trI
T∗ =−p∗I+β(L∗+L∗T)+ 1−β(M−I), (11) Ω (cid:13) (cid:1) (cid:14)
1
Wi + (∇·vh)I dΩ+ Sh : Wi(vh·∇Mh−(Lh)T·Mh
trI
where β= ηsη+sηp is the viscosity ratio. Hereafter, all variables (cid:15)Ω (cid:1) (cid:14)
aredimensionlessandthe(∗)isomittedforclarity. −Mh·Lh)+(Mh−I) dΩ+ τ Wi(vh·∇Sh
cons
Ω
4 O.M.Coronadoetal./J.Non-NewtonianFluidMech.xxx (2006)xxx–xxx
(cid:15) (cid:14)
−(Lh)T·Sh−Sh·Lh)+Sh : Wi(vh·∇Mh−(Lh)T·Mh
(cid:15)
−Mh·Lh)+(Mh−I) dΩ=0,
∀qh ∈Vh,∀wh ∈Vh,∀Eh ∈Vh,∀Sh ∈Vh , (21)
p v L M
where τ ,τ ,τ and τ are the LS stabilization pa-
mom cont gradv cons
rameters for the momentum, continuity, interpolated traceless
velocity gradient and constitutive equations, respectively. The
underbracedtermAisneglectedatlowWibecausethe(1/Wi)
termgrowslargeasWi→0,causingnumericalproblems.
3.1. Designofthestabilizationcoefficients
Theappropriatedesignofthefourstabilizationparameters—
τ ,τ ,τ andτ —inEq.(21)playsacrucialrolein
mom cont gradv cons
theperformanceofthemethod. Fig.2. Mesh-convergencerateforaplanarchannelflowatdifferentWi.The
The τmom-term stabilizes the Galerkin form in advection- slopeofthecurvesgivestherateofconvergencewithmeshrefinement.
dominatedflows,andalsoremovesthecompatibilitycondition
betweenvelocityandpressurespaces.Theparameterdesigned
specificallyforusewithbilinearinterpolations[31]isadapted
hereforthedimensionlesssystem:
h2
τ = . (22)
mom
4
wherehisthedimensionlesselementlength.
The τ -term improves the convergence of non-linear
cont
solvers in advection-dominated problems. Hereafter, τ =0
cont
Fig.3. Geometryofaflowpastacylinderinahalfchannel.
becauseinertiaisneglected.
The τ -term stabilizes Eq. (4); although the associated
gradv
stabilizationtermisnotstrictlynecessary,τ istakenhereas 1
gradv τ = , (24)
τgradv =1. cons2 Wi(cid:8)Lh(cid:8)
The τ -term is introduced to stabilize the Galerkin form
cons
at high Wi, and to bypass the compatibility conditions be- τ = h . (25)
tweenvelocityandconformationspaces.Nosystematicderiva- cons3 2Wi(cid:8)vh(cid:8)
tion for τ is available in the literature. However, the trans-
cons
τ andτ areimportantinregionsoftheflowwheregen-
portequationofconformationcanbeviewedasanadvection– cons1 cons2
eration is dominant, whereas τ is important in advection-
generation equation, and considerable research has been done cons3
dominatedregions.Thesethreecontributionscanbecombined
on stabilization parameters for a simple advection–diffusion–
generationequation[35–39].Applyingthedefinitionproposed
byFrancaetal.[38],basedontheconvergenceandstabilityanal-
ysis of advection-diffusion-generation equation, and extended
byHauke[39],yields
τ =1, (23)
cons1
Fig.1. Schematicofaflowinaplanarchannelwithw/L=1/4.Thetopwall
iskeptfixed,thebottomwallismovingfromrighttoleftatv0andadifferential Fig.4. Flowpastacylinderinachannel,w/Rc=2:FiniteelementmeshM0
pressureisappliedbetweentheleftandrightwalls. (a)completedomain(b)detailofthemeshfromx=−2tox=2.
O.M.Coronadoetal./J.Non-NewtonianFluidMech.xxx (2006)xxx–xxx 5
Fig.5. (Coloronline)Flowpastacylinderinachannel,w/Rc=2:dragforce Fthieg.s7y.mFmloewtryplaisnteaicnytlhinedwerakineaatcWhain=nel0,.w7./(cid:9)Rcfr=om2:Hσuxlxseonnethtealc.y[l3i4n]d.erandon
onthecylinderversusWi.TheGLS4resultsforthefourmeshes(M1,M2,M3
andM4)arecomparedwiththeresultspresentedbySunetal.[21]andHulsen
etal.[34].Inset:detailofthedragforceathighWi.(cid:1)representsthedragforce
onM4atWi=0.6.AtWi=0.6,theextrapolatedvalueofthedragforceis
117.979,whichiswithin0.2%ofthevaluesreportedinRefs.[30,34,40].
as:
(cid:16) (cid:17)
1 1 1 −1/r
τ = + + , (26)
cons τr τr τr
cons1 cons2 cons3
Hereafter, the switching parameter is set to r =2 (see also
Ref.[29]):
⎡ (cid:18) (cid:19)2⎤−1/2
τ =⎣1+(Wi(cid:8)Lh(cid:8))2+ 2Wi(cid:8)vh(cid:8) ⎦ . (27)
cons
h
4. Numericalresults
The proposed GLS4 formulation is tested in flow in a pla-
Fig.8. Flowpastacylinderinachannel,w/Rc=2:DragforceatWi=0.6for
narchannelandflowpastacylinderinachannel.Ananalytical GLS4andDEVSS-TG/SUPGforallmeshes;dashedlinerepresentsthedrag
solution can be obtained in the former case; in the latter, the forcereportedbyHulsenetal.[34]ontheirfinestmesh.
Fthieg.s6y.mFmloewtryplaisnteaicnytlhinedwerakineaatcWhain=nel0,.w6./(cid:9)Rcfr=om2:Hσuxlxseonnethtealc.y[l3i4n]d.erandon tFhieg.d9ra.gFfloowrcepaasttWaicy=lin0d.6e.rinachannel,w/Rc=2:Mesh-convergencerateof
6 O.M.Coronadoetal./J.Non-NewtonianFluidMech.xxx (2006)xxx–xxx
computation of the relative errors e=|(numericalvalue−
analyticalvalue)/(analyticalvalue)|×100%inregionA.Fig.2
showsthemaximumeinM (whichhasthehighesteamong
yy
allunknownfields)versustheelementsizeforWi=3,5and7.
Fromthethreecurvestherateofmeshconvergenceisestimated
tobe1.73,1.63and1.59,respectively.BecauseincreaseinWi
results in increased generation, subsequently forming steeper
boundary layer close to the channel walls, the rate of conver-
gence is found to decrease. At Wi=3, Xie and Pasquali [41]
reportedarateofconvergenceof1.89usingDEVSS-TG/SUPG
methodwithbi-quadraticinterpolationforvelocity.
4.2. Flowpastacylinderinachannel
TheflowofanOldroyd-Bfluidpastacylinderinarectangu-
larchannelhasbeenusedasastandardbenchmarkproblemto
Fig.10. Flowpastacylinderinachannel,w/Rc=2:Mesh-convergencerate testseveralcomputationalmethods[20–22,29,34,40].Forcom-
ofMxxatapointinthewakeflow(x=2;y=0)atWi=0.6. putationalease,thesymmetryoftheproblemisusedandonly
halfofthechannelissimulated.Fig.3showstheschematicof
numerical results from other state-of-the-art methods are used the problem, where L ,L ,R and w represent the upstream
u d c
for validation [20–22,29,34,40]. The flow past a cylinder in a length,thedownstreamlength,thecylinderradius,andthehalf
channelisastandardbenchmarkproblemwithdesirablechar- channelwidth,respectively.
acteristics of smooth boundaries, and poses several numerical Ano-slipboundaryconditionisimposedonthecylindersur-
challenges at high Wi due to the formation of sharp boundary faceandchannelwalls,andfullydevelopedflowconditionsare
layersonthecylinderandinthewake. assumedattheinflowandoutflowboundaries.Consequently:
(cid:16) (cid:17)
Q y2
4.1. Flowinaplanarchannel v =1.5 1− , v =0, (30)
x w w2 y
Fig.1showsacombinationofPoiseuilleflow(pushingliq- (cid:16) (cid:17)
uidfromlefttoright)andCouetteflow(inducedbythebottom Q y 3Qλy 2
M =−3 λ , M =M =1−2 − ,
walldraggingliquidfromrighttoleftwithvelocityv0)inapla- xx w w2 xy yx w3
nar channel of width w=1 and length L=4w. The flow of
M =1, (31)
anOldroyd-Bfluid(β=0.59)issimulated,andtheresultsare yy
compared with the known analytical solution. The figure also where Q is the flow rate. Whereas the velocity is imposed at
shows velocity profiles at the two open flow boundaries; both both inflow and outflow, the conformation tensor components
right and left ends of the channel have respective inflow and arespecifiedattheinflowonly.Atthesymmetryline,n·T=0
outflowsections.Aregion‘A’(dottedareainFig.1),whichis andv =0,wherenistheunitvectornormaltothesymmetry
y
2winlengthandcentrallyplacedinthechannel,ismonitored line.Thecomputeddragonthecylinderf hasbeentraditionally
d
for comparing numerical results with analytical solution; this usedtocomparenumericalmethods:
(cid:1)
sufficientlyeliminatestheinfluencesduetotheboundarycondi-
tions.Theproblemsetupcloselyfollowsthenumericalexample f =−2 e n:TdS, (32)
d 1
employedbyXieandPasquali[41];theanalyticalsolutionfor S
velocityandconformationfieldsare whereSrepresentsthesurfaceofthecylinder,ntheunitnormal
(cid:12) (cid:12)(cid:20) (cid:21) (cid:13) (cid:13) vector,ande1istheunitvectorinthex-direction.
(cid:8)pw y 2 y y
v = − − + −1 v , v =0, (28)
x 0 y
2 L w w w 4.2.1. Flowpastacylinderinachannel:w/R =2
c
(cid:16) dv (cid:17)2 dv Inthiscase,w=2,Rc =1,Lu =−20,Ld =20,Q=2and
Mxx =1+2 λ x , Mxy =λ x, Myy =1, (29) β=0.59.Fig.4showsthemeshM0fromwhichfoursystemati-
dy dy
callyrefinedmeshesareobtained;inthesemeshes,theelements
where (cid:8)p=50 is the differential pressure between the left areconcentratedonthecylindersurfaceandinthewakealong
and right boundaries. Consequently, Wi=λ[(cid:8)pw/(2L)+ the symmetry line. The M1, M2, M3 and M4 meshes are ob-
1](v /w). The Dirichlet conditions are imposed for ve- tained by dividing every element side of M0 by 3, 4, 5 and 6,
0
locity components on all boundaries, and the conforma- respectively.Thenumberofelementsandcorrespondingnum-
tion tensor components are only specified at the cor- berofunknownsarelistedinTable1.Thisflowproblemposes
responding inflows. The numerical results are obtained numericalchallengesathighWi;therefore,themaximumWiup
on four different uniform meshes — 16×16, 24×24, to which the numerical schemes converge has been employed
32×32 and 64×64 — followed by a node-by-node as a measure of robustness (but not necessarily accuracy). For
O.M.Coronadoetal./J.Non-NewtonianFluidMech.xxx (2006)xxx–xxx 7
Fig.11. (Coloronline)Flowpastacylinderinachannel,w/Rc=2:(a)Mxx,(b)Mxyand(c)MyycontoursatWi=0.7onmeshM2.
example,usingDEVSS-G/SUPG,Sunetal.[21]reportedsolu- Here, a sequence of flow states is computed by first-order
tionsuptoWi=1.85,Fanetal.[30]usinganincompleteGLS arc-lengthcontinuationonWiwithautomaticstepcontrol;the
uptoWi=1.05andHulsenetal.[34]usingthelogconforma- continuationterminateswhentheconformationtensorlosesits
tion up to Wi=2.0; however, the accuracy of the solutions at positive-definiteness, which occurs at Wi∼0.7. The positive-
Wi>0.6wasnotconfirmedintheseworks. definiteness of M was not usually considered in past stud-
Fig.12. Flowpastacylinderinachannel,w/Rc=2:Mxxalonglinex=2on Fig.13. Flowpastacylinderinachannel,w/Rc=2:Mxyalonglinex=2on
meshM3.Inset:detailofMxxnearthecenterline(y=0). meshM3.Inset:detailofMxynearthecenterline(y=0).
8 O.M.Coronadoetal./J.Non-NewtonianFluidMech.xxx (2006)xxx–xxx
Fig.14. Flowpastacylinderinachannel,w/Rc=2:Myyalonglinex=2for Fig.16. (Coloronline)DirectcomparisonofGLS4andDEVSS-TG/SUPGwith
M3.Inset:detailofMyynearthecenterline(y=0). respecttothenumberofelements(bottomaxis)andtothenumberofdegreesof
freedomforconformation(topaxis).Theleftandrightaxesrepresentthetime
ies, with exception of the recent work of Hulsen et al. [34]. perNewtoniteration(s)andmemoryusage(MB),respectively.Afrontalsolver
isusedinbothsimulations.
Fig. 5 shows the drag forces on the meshes M1, M2 and
M3; a good agreement is found up to Wi=0.4 with the re-
sults reported by Hulsen et al. [34] and Sun et al. [21]. Be-
yond that, the three methods show slight differences in the
drag predictions, while following the same trend. Figs. 6 and
7 show σ versus s at Wi=0.6 and 0.7, respectively, where
xx
σ =(η /λ)M andsisdefinedas:0<s<πR onthecylin-
xx p xx c
derandπR <s<πR +L −R inthewakealongthesym-
c c d c
metryline.AtWi<0.6,acompleteoverlapisobservedamong
theresultsonthemeshesM1,M2andM3.In Fig.6,σ pro-
xx
files computed with M1, M2 and M3 are overlapping in the
wake; however, the result from M1 shows underprediction on
thecylinder,implyingthatrefinementofM1isnotsufficientto
capture the steep boundary layer on the cylinder. On the other
hand,inFig.7,differencesareobservednotonlyonthecylin-
der but also in the wake flow. The figures also show the re-
sults reported by Hulsen et al. [34] at the corresponding Wi, Fig.17. Flowpastacylinderinachannel,w/Rc=8:FiniteelementmeshM0
(a)completedomain(b)detailofthemeshfromx=−4tox=4.
Fig.15. Flowpastacylinderinachannel,w/Rc=2:σxxonthecylinderand Fig.18. Flowpastacylinderinachannel,w/Rc=8:Dragforceonthethree
onthesymmetrylineatWi=0.6.TheGLS4andDEVSS-TG/SUPGresultsare meshes.ThedottedcurveisobtainedfromSunetal.[21].Inset:detailofthe
obtainedforM2.(cid:9)fromHulsenetal.[34]. dragforceathighWi.
O.M.Coronadoetal./J.Non-NewtonianFluidMech.xxx (2006)xxx–xxx 9
Fig.19. (Coloronline)Flowpastacylinderinachannel,w/Rc=8:(a)Mxx,(b)Mxyand(c)MyycontoursatWi=2.0onmeshM2.
Table1
Flowpastacylinderinachannel,w/Rc=2:Characteristicsofthefiniteelementmeshes
Mesh Elements Unknowns TimeperNewtoniteration(s) Memoryusage(MB)
GLS4 DEVSS GLS4 DEVSS (cid:10)d GLS4 DEVSS (cid:10)d
M0 532 − − − − − − − −
M1 4788 49960 79102 220 322 31.7 608 1140 46.7
M2 8512 87890 139514 648 934 30.6 1410 2647 46.7
M3 13300 136460 216950 1460 2144 31.9 2720 5122 46.9
M4 19152 195670 311410 2914 − − 4650 − −
(cid:10)disthe%relativedifferenceintherespectivevaluesfromDEVSSandGLS4.
Table2
Flowpastacylinderinachannel,w/Rc=8:Characteristicsofthefiniteelementmeshes
Mesh Elements Unknowns TimeperNewtoniteration(s) Memoryusage(MB)
M0 250 − − −
M1 4000 41690 198 396
M2 6250 64610 440 863
M3 12250 125450 1505 1859
10 O.M.Coronadoetal./J.Non-NewtonianFluidMech.xxx (2006)xxx–xxx
and good agreement is found with results on M3. In previous
works,theconvergenceofthestresseshasbeenshownbycom-
paringstressprofilesobtainedonsystematicallyrefinedmeshes
usingp-andh-refinement[30,34].Whileoverlapoftheresults
demonstratesqualitativelymeshconvergence,here,accuracyis
measuredmorepreciselybyRichardsonextrapolation:
f (h)=f (0)+αhn, (33)
d d
wheref (0)isthedragforaninfinitelyrefinedmesh,ntherateof
d
meshconvergenceandαisaconstant.f (0)isusedtocompute
d
therelativeerrorse=|(f (h)−f (0))/f (0)|×100%inf (h).
d d d d
Fig.8showsf predictionsfromGLS4andDEVSS-TG/SUPG
d
alongwiththeresultspresentedbyHulsenetal.[34]atWi=0.6.
The extrapolated values of f for an infinitely refined mesh
d
are 117.979 and 117.778 for GSL4 and DEVSS-TG/SUPG,
respectively. Thus, the GLS4 extrapolated results are within
0.2% of the values computed by high resolution finite volume Fig.20. Flowpastacylinderinachannel,w/Rc=8:Mxxonthecylinderand
(f =117.79[40]),bypressure-stabilizedfiniteelements(f = alongthesymmetrylineinthewakeatWi=1.5.
d d
117.78[30]),byDEVSS-DG-Logconformationfiniteelements
(f =117.77 [34]), and by DEVSS-TG/SUPG calculations.
d 4.2.2. Flowpastacylinderinachannel:w/R =8
Fig.9showseversush,andameshconvergencerateof2.39is c
Inthiscase,w=8,R =1,L =−40,L =40,Q=8and
observed. c u d
β=0.59.Followingthesameprocedureasinthepreviouscase,
Similarly, Richardson extrapolation analysis is also per-
formed for M at a point in the wake flow (x=2;y=0). the M1, M2 and M3 meshes are obtained by dividing every
xx
element side of the mesh M0 by 4, 5 and 7, respectively. Fig.
The extrapolated value of M is 26.05 and the rate of mesh
xx
17 shows the mesh M0, and details of the subsequent meshes
convergence is 1.74. Fig. 10 shows e versus h for M . In
xx
are listed in Table 2. The drag on the cylinder from the three
all cases, results on M4 are also employed to obtain the ex-
meshesarecomparedwiththeresultsreportedbySunetal.[21]
trapolated values. Fig. 11 shows the conformation contours at
Wi=0.7,inwhich,theformationofsharpboundarylayerson inFig.18.Forthiscase,acompleteoverlapofdragpredictions
fromthethreemeshesisobserved,andagoodagreementwith
thecylinderandalongthesymmetrylineinthewakefloware
resultsofSunetal.[21]isfounduptoWi=2.0.Moreover,the
observed.Theseboundarylayersaredifficulttoresolvenumer-
maximum Wi achieved in this simulation is ∼2.7. In Fig. 19,
ically,andtheonsetofoscillationsintheconformationfieldsis
thecontourplotsfortheMcomponentsareshownatWi=2.0.
observedastheboundarylayersgrowathighWi.Theinfluence
Figs. 20 and 21 show M versus s (defined in Section 4.2.1)
of the sharp boundary layer at high Wi is shown in Figs. 12– xx
14, which plot M components along line x=2. At Wi=0.5 atWi=1.5and2.0,respectively.AtWi=2.0,thestreamwise
normalconformationcomponentM hasnotyetconvergedin
and 0.6 a smooth profile for the M components is observed, xx
whereas at Wi=0.7 oscillations appear towards the symme-
try line (y→0). The maximum Wi attained in these simula-
tions is slightly above 0.7; beyond this, M loses its positive-
definiteness.
A direct comparison of computational cost between GLS4
and DEVSS-TG/SUPG is performed, while keeping the same
numberofdegreesoffreedomforconformation;thelatterem-
ploys biquadratic interpolation functions for velocity, whereas
bilinear for pressure, velocity gradient and conformation. The
resultsatWi=0.6onM2areobtainedfrombothmethodsand
comparable accuracy is observed. Fig. 15 shows σ versus s
xx
along with the results of Hulsen et al. [34]. The GLS4 and
DEVSS-TG/SUPGcharacteristics(numberofunknowns,time
perNewtoniterationandmemoryusage)arelistedinTable1.
Fig. 16 shows a direct comparison of GLS4 and DEVSS-
TG/SUPG with respect to the number of degrees of freedom
forconformation,anditcanbeobservedthat GLS4is∼30%
computationallyfasterandusesonly45%ofthememorycom-
paredtoDEVSS-TG/SUPG(forthesamenumberofdegreesof
Fig.21. Flowpastacylinderinachannel,w/Rc=8.Mxxonthecylinderand
freedomforconformation). alongthesymmetrylineinthewakeatWi=2.0.