An Adaptive Mixed Least-squares finite Element Method For .

J. Non-Newtonian Fluid Mech. 159 (2009) 72–80Contents lists available at ScienceDirectJournal of Non-Newtonian Fluid Mechanicsjournal homepage: www.elsevier.com/locate/jnnfmAn adaptive mixed least-squares finite element method forviscoelastic fluids of Oldroyd type夽Z. Cai a , C.R. Westphal b, abDepartment of Mathematics, Purdue University, 150 N. University St., West Lafayette, IN 47907-2067, United StatesDepartment of Mathematics and Computer Science, Wabash College, 301 W. Wabash, Crawfordsville, IN 47933, United Statesa r t i c l ei n f oArticle history:Received 25 June 2008Received in revised form 21 January 2009Accepted 5 February 2009Keywords:Viscoelastic fluidOldroydLeast-squaresRaviart–ThomasFinite elementAdaptive mesh refinementCorner singularities4:1 contractiona b s t r a c tWe present a new least-squares finite element method for the steady Oldroyd type viscoelastic fluids.The overall iterative procedure combines a nonlinear nested iteration, where adaptive mesh refinement isbased on a nonlinear least-squares functional. Each linear step is solved by a least-squares finite elementminimization. The homogeneous least-squares functional is shown to be equivalent to a natural norm,and, under sufficient smoothness assumptions, finite element error bounds are shown to be optimalwhen using conforming piecewise polynomial elements for the velocity, pressure and extra stress, andRaviart–Thomas finite elements for the total stress. In the absence of full regularity, a local weighted-normapproach is used to remove effects of corner singularities. Numerical results are given for an Oldroyd-Bfluid in a 4:1 contraction, showing optimal reduction of the least-squares functional. 2009 Elsevier B.V. All rights reserved.1. IntroductionWhile much progress has been made in recent years toward theaccurate and efficient simulation of viscoelastic fluids under differential constitutive laws, many difficulties persist. We considerthe solution to steady Oldroyd systems, where, in contrast to Newtonian models, the stress cannot be eliminated to form a singlesecond-order equation in terms of the velocity and pressure. Ingeneral, the stress must be directly approximated.As the Weissenberg number increases, the constitutive equationexhibits dominant convective behavior and the nonlinear couplingbetween the unknowns increases. One major difficulty in the simulation of viscoelastic fluids is the failure of the nonlinear iterationto converge at some threshold in the Weissenberg number. Thus,it is critically important to design a method that is robust withrespect to the nonlinear solver. Nonsmooth boundaries can introduce singularities which, if left unaddressed, may cause significantdegradation of numerical methods, or worse, may even contributeto the failure of the nonlinear iteration. In fact, near reentrant cor-夽 This work was sponsored by the National Science Foundation under grantDMS–0511430. Corresponding author.E-mail addresses: zcai@math.purdue.edu (Z. Cai), westphac@wabash.edu,loschados@yahoo.com (C.R. Westphal).0377-0257/ – see front matter 2009 Elsevier B.V. All rights reserved.doi:10.1016/j.jnnfm.2009.02.004ners the extra stress may not be in H 1 ( ), a major difficulty forusing finite element subspaces of H 1 ( ). Finally, the algebraic systems that must be solved from typical discretizations are difficult tohandle and require the use of efficient iterative solvers. The properties of the resulting algebraic systems should be considered closelyin the design of the discretization to ensure an efficient numericalmethod.The analytical and numerical study of viscoelastic fluidshas developed a rich history, distinguished by contributions toindustrial applications as well as providing a generalization to classical fluid mechanics. This range is reflected in works [34,35,2,21,37].In [39] Wang and Carey introduced a least-squares approach toupper-convected Maxwell (UCM) fluids and in [5] Bose and Careygave a more sophisticated approach for UCM and Bingham fluids,utilizing mesh redistribution and polynomial refinement to handleproblems due to boundary singularities. Gerritsma, in [23], useda discontinuous least-squares spectral element approach for UCMfluids. This paper is similar in general approach to these papers, butprovides a novel approach to the nonlinear iteration, treatment ofboundary singularities, and adaptive mesh refinement. Also, to ourknowledge, this paper is the first to prove least-squares functionalellipticity results and rigorous error bounds for linearized viscoelastic fluids. The theory we present here includes the Oldroyd modelsin two or three dimensions and also includes the inertial terms inthe balance of momentum equations.

Z. Cai, C.R. Westphal / J. Non-Newtonian Fluid Mech. 159 (2009) 72–80The least-squares finite element method is based on an explicitminimization principle to form a symmetric variational problem. In recent years, least-squares methods have been developedand refined for many applications in continuum mechanics.See for example general references [26,4,7,10], works on fluids[11,12,28,8,3], or works on solid mechanics [13,14,38,9,32]. Theleast-squares approach displays several attractive features for suchproblems. In its most basic form, a least-squares method is formedby minimizing the norm of the residual of the equations over asolution space compatible with the exact solution. An appropriate formulation of the equations and choice of norms to minimizeover ensures that the error is efficiently reduced in a meaningfulnorm. The locally evaluated least-squares functional can then serveas a reliable and inexpensive error estimator, which can be used inadaptive refinement strategies. Further, a strong advantage of themethod is that it always produces symmetric, positive definite linear systems that can generally be solved with optimal complexitywith standard multigrid methods. It is also worth noting that theleast-squares approach here is not the same as the Galerkin LeastSquares (GLS) approach.Nonlinear PDE problems can be approached by linearizing theequations first and then applying a least-squares method to the linear equations, or by forming a nonlinear variational problem andthen linearizing the resulting equations. In this paper, we use theformer approach and note that it is not clear which is more effectivefor difficult nonlinearities. In this paper, we combine the use of aleast-squares discretization with a nonlinear nested iteration. Thatis, the nonlinear iteration is done in tandem with mesh refinement.At each linearzation step, a nonlinear functional serves as an indicator to decide when to interpolate a current solution to a finer mesh.This not only produces a more efficient iteration by performing theearlier linearization steps with coarser resolution, but it also provides good initial guesses for the finer meshes, producing a morerobust nonlinear iteration.As are many finite element methods, least-squares methods aresensitive to a loss of regularity due to nonsmooth boundaries, thus,care must be taken in this respect. In the presence of such singularities, our approach replaces the L2 norms in the least-squaresfunctional with locally weighted norms, which eliminates the pollution effect and improves discretization accuracy. Details of thisapproach can be seen in [29,30,18,15].The organization of this paper is as follows. Section 2 introduces the notation we use and the equations we consider.Section 3 describes the nonlinear iteration, the function spacesfor the unknowns and the adaptive refinement procedure. Theleast-squares minimization is detailed and the ellipticity of theleast-squares functional is proved in Section 4. In Section 5 thefinite element spaces are discussed and optimal error bounds areproved. Section 6 provides details of the weighted-norm treatment of boundary singularities, and numerical results are given inSection 7.2. Notation and equationsdH k ( )We use standard notation for Sobolev spacesand corresponding norm · k, for k 0. We drop subscript andsuperscript d when the domain and dimension are clear by context. For noninteger k, H k ( ) is the interpolation space betweenH k ( ) and H k ( ) as in [31]. The case of k 0 corresponds tothe Lebesgue measurable space, L2 ( ), in which case we generallydenote the norm and inner product by · and ·, · , respectively.dDefine the subspace of L2 ( ) induced by the divergence of u bydH(div) {u L2 ( ) : · u },73with norm satisfying u 2H(div) u 2 · u 2 .The divergence and trace of d d tensor are given by d d ij / xj and tr( ) . The space of ten( · )i j 1i 1 iisor valued functions with each row in H(div) is denoted byH(div)d . We also define the operator u · duk / xk so that k 1dfor each i, j 1, . . . , d we have (u · u)i du ui / xkk 1 kand(u · )ij u ij / xk .k 1 kThe conservation of mass and momentum equations for timedependent incompressible flow are given by · u 0,and u u · u t · f̃,(1)where is the fluid density, f̃ is an internal body force and theunknowns u and are the velocity vector and the stress tensor,respectively.System (1) must be closed by a constitutive equation. In thispaper we consider a general Oldroyd fluid, where an elastic material of viscosity 1 is dissolved in a viscous solvent of viscosity 2 . The total stress has elastic and solvent components, and s ,respectively, which satisfy s pI,where p is the pressure and I is the d d identity tensor. The solventstress satisfies s 2 2 (u),where (u) (1/2)( u ut ) is the standard strain rate tensor.The symmetric elastic stress satisfies u · ga ( u, ) t 2 1 (u),where is a characteristic stress relaxation time for the fluid andga ( u, ) 1 a1 a(( u) ( u)t ) ( ( u) ( u)t )22is a bilinear tensor-valued function, depending on the parametera [ 1, 1].We nondimensionalize the equations by scaling the lengths by L,the velocities by V, time by L/V , and the stresses by V/L, where 1 2 is the total viscosity. Further, we define R VL/ , W V/L, and ω 1 / [0, 1] as the dimensionless Reynolds number,Weissenberg number and retardation parameter, respectively. Thisgeneral viscoelastic model is characterized by the four parametersR , W , a and ω and can be written as · u 0 u u· u · f̂R t pI 2(1 ω) (u) 0 u · ga ( u, ) 2ω (u) 0 W t u 0in ,in ,in ,(2)in ,on ,where for simplicity of presentation, homogeneous boundary conditions are assumed on . Nonhomogeneous velocity boundaryconditions can be treated analogously, but our present analysisis restricted to pure Dirichlet type boundary conditions. System(2) represents unknowns u, p, , and . In general, least-squaresmethods involve an appropriate first-order system of PDEs, oftenformed by introducing new unknowns to a second-orderPDE.