Weak Form & LiveLink For MATLAB Based Modified Uzawa .

Weak Form & LiveLink for MATLAB Based Modified UzawaMethod for Solving Steady Navier-Stokes EquationHuashan Sheng*1，Shuai Zhu11 Department of Mathematics Shanghai Jiao Tong University* Dongchuan Rd. No. 800, Minhang, Shanghai, China, shs3701001@sjtu.edu.cnAbstract: COMSOL Multiphysics can solvenonlinear PDEs easily by the built-in Newtonmethod, and one typical example is the NavierStokes equation. But, sometimes we can’t figureout a good initial guess for Newton method, ormaybe we have some other better algorithms. Inthese cases, The Laminar Flow interface may notso suitable. So in this paper, we test the modifiedUzawa iterative algorithm[1] for Navier-Stokesequation by using weak form and LiveLink forMATLAB in COMSOL. We get the samesolution as other FEM packages (FreeFEM &FEniCs) even at each iterative step.Keywords ： Navier-Stokes equation, ModifiedUzawa Method, Weak Form PDE, LiveLink forMATLABWe introduce some notation for later use. Letdenote the scalar product over. Givena non-negative integer , letbe the usualSobolev space consisting of all functionswhose generalized derivatives withthe total degree no more than are stillintegrable. With the known number , letwith the standard normand seminorm.Furthermore, letanddefine a trilinear form overas follows:(1-3)The weak form of Navier-Stokes Equation (1-1):Findsuch that1. IntroductionThis paper is concerned with using the modifiedUzawa iterative algorithm for solving the ly in COMSOL Multiphysics.1.1 Navier-Stokes Equation TheorySteady incompressible Navier-Stokes equation:(1-1)Where the domain(or 3) is abounded domain with Lipschitz boundary,means the viscosity coefficient, anddenoting the Reynolds number; the right handside is the prescribed body force; and arethe corresponding velocity field and the pressurefield, respectively. The above PDEs are the wellknown mathematical model describing the steadyflow of an incompressible Newtonian fluid, suchas air or water, which are frequently encounteredin engineering applications, from airplane designto the construction of power stations. Sometimes,to simplify, we impose the non-slip condition forNavier-Stokes Equation:(1-2)To ensure the uniqueness of pressure fieldwe always assume that:(1-4)1.2 Mixed Element MethodWith the partition , we associate a stable pairof finite element spaces, or equivalently,(resp. ) is a subspace of (resp. ), whichadmit inf–sup(B.B) condition. We always use theMINI element, the Girault–Raviart element. Inthis paper we use theelement which is alsocalled the Taylor–Hood element. The mixedelement method for (1-4) is described as follows:Findsuch that(1-5)1.3 The Modified Uzawa MethodGiven two initial functionsand , which aregiven by zeros initial guess or the solution of thediscrete Stokes equations, and given a positivenumber, run the following iteration untilwe get the desired accuracy:Excerpt from the Proceedings of the 2014 COMSOL Conference in Shanghai

(1-6)In paper [1], we have proved that the modifiedUzawa methods converge geometrically with acontraction number independent of the finiteelement mesh size .For understanding convenience, detailedsteps of algorithm are given as follows:Step 0: Set initial value,iterative step n 0, giveStep 1: Backup andStep 2: Solve velocity fieldby (1-6)Step 3: Solve pressure fieldby (1-6)Then add the physics field in Componentbranch, which is the most important. As step 1 inAlg. 1, add Coefficient Form PDE with twodependent variables named, andchoose the shape function type to be quadraticLagrange. Let all the coefficients equal toexcept the absorption coefficientand the source term. Add anotherCoefficient Form PDE with one dependentvariable named, and choose linearLagrange shape function type. Similarly, let allthe coefficients equal to except the absorptioncoefficientand the source term.As step 2, add Weak Form PDE with twodependent variables named, and choosequadratic Lagrange shape function type. Writethe weak expressions as follows:Step 4: Update velocity and pressure fieldStep 5: Break condition;;ifepsbreak;Go to Step 6;else;Go to Step 1;endStep 6: Output solution, plot and error analysisAlgorithm 1. The modified Uzawa Method2. COMSOL Multiphysics SettingsThe Navier-Stokes equation can be solvedconveniently by the laminar flow interface withbuilt-in Newton method, but it’s difficult tochange its built-in algorithms. Fortunately, theWeak Form PDE interface in the Mathematicsbranch can solve this problem perfectly. Thesteps for setting are given as follows (2Dproblems):First, start COMSOL Multiphysics, and builda new mph file with 2D domain. Define someparameters and variables in Global Definitionsbranch, and draw the domain for solving inGeometry branch.After input the weak expressions, add someDirichlet boundary conditions. As step 3, addWeak Form PDE with one dependent variablenamed , and choose linear Lagrange shapefunction type. Write the weak expressionssimilarly as follows:Notice that step 4 and step 1 in Alg. 1 arealmost the same, add Coefficient Form PDE withtwo dependent variables named, andchoose the shape function type to be quadraticLagrange. Let all the coefficients equal toexcept the absorption coefficientand the source term. Add anotherCoefficient Form PDE with one dependentvariable named, and choose linearLagrange shape function type. Similarly, let allthe coefficients equal to except the absorptioncoefficientand the source term.Next, add three stationary study steps, and letfirst step only solve for the first two CoefficientForm PDEs, and second step only solve for thetwo Weak Form PDEs, and last step only solvefor the last two Coefficient Form PDEs as shownin Fig. 1. “Values of variables not solved for”Excerpt from the Proceedings of the 2014 COMSOL Conference in Shanghai

should be checked in all study steps as shown inFig. 2.3. LiveLink for MATLABLiveLink for MATLAB is a powerful function inCOMSOL for designing iterative algorithm. Wewill let MATLAB do the “while” loop in Alg. 1.The way starting COMSOL with MATLAB onWindows / Mac OSX / Linux can be found inhelp document. Once they have been connected,the mph model can be handled in MATLAB soconveniently as using inline function.To realize the “while” loop, step 5 and step 6in Alg. 1, first change the current MATLAB foldto the location of the COMSOL file. Then runthese codes as follows:Figure 1. Details for setting in Study Step 1Figure 2. Details for setting in all Study StepsLastly, partition the solving domain with thesize you need, and save the model withoutcomputing. Notice that it’s no need to computethe model in COMSOL, and once click“compute”, the initial value of and will bechanged. The “tree” structure is shown in Fig. 3model mphload('Uza2.mph');ERROR 1;COUNT 1;ERRUVEC [];ERRUPRE [];while ERROR 1e-6 && COUNT 1000model.sol('sol1').run;Q1 mphmax(model,.'abs(u-uk bk)','surface');Q2 mphmax(model,.'abs(v-vk bk)','surface');Q3 mphmax(model,.'abs(p-pk bk)','surface');B1 mphmax(model,.'abs(uk bk)','surface');B2 mphmax(model,.'abs(vk bk)','surface');B3 mphmax(model,.'abs(pk bk)','surface');ErrU max(Q1/B1,Q2/B2);ErrP );ERROR max(ErrP,ErrU);COUNT COUNT 1;ERRUVEC [ERRUVEC;ErrU];ERRUPRE [ERRUPRE;ErrP];endfprintf('---------OK!-------- \n ');%Plot the solution and Err via ERRUVEC);grid on;subplot(2,2,4);semilogy(ERRUPRE);grid on;Figure 3. “Tree” structure of modelExcerpt from the Proceedings of the 2014 COMSOL Conference in Shanghai

4. Numerical ExperimentThe only other change in different models is theComponent branch. For example, the domain ischanged or the inlet velocity is larger, nothingneeds to do except changing the Geometrybranch and the Dirichlet boundary condition. Inthis part, four typical examples will be tested.4.1 Flows in a Confined CavityWe take the domain as, setthe viscosity coefficientand choose theright side function in (1-1) such that the PDEhas the following solution:Figure 5. Compare with FreeFEM (Right part)HereWe choosethe parameter, and let themesh size. Solution and FEMconvergence order are shown in Fig. 4-6.Figure 6. FEM Convergence Order of P2-P1 Element4.2 Back-Step Problem(a) Streamlines and Pressure ValueThe second classic example is the back-step flow.Let domain,and set the viscosity coefficientandchoose the right side function. The leftside of the “back-step” is the inlet boundary, andthe flow velocity isperpendicular to the inlet boundary. The rightside of the “back-step” is the outlet boundarywith do-nothing condition, and this boundary isNeumann boundary. The up side and down sideof the “back-step” are the Dirichlet boundarywith.The solution and iterative error with meshsizeand parameterare shown inFig. 7.(b) Iterative Error between two stepFigure 4. Solution of Cavity Flow at(a) Velocity magnitudeExcerpt from the Proceedings of the 2014 COMSOL Conference in Shanghai

(b) Streamlines(b) Streamlines(c) Pressure value(c) Velocity magnitude(d) Iterative Error between two stepFigure 7. Solution of back step flow4.3 Channel Flow past a Circular CylinderThe third example is the channel flow past acircular cylinder. Let the solving domain, and set theviscosity coefficientand choose theright side function. The left side of the“channel” is the inlet boundary, and the flowvelocity isperpendicular tothe inlet boundary. The right side is the flow-outboundary with do-nothing condition, and thisboundary is Neumann boundary. The other sidesare the entire Dirichlet boundary with.The solution and iterative error with extremelyfine physics-controlled mesh, and parameterare shown in Fig. 8.(a) Velocity magnitude(d) Iterative Error between two stepFigure 8. Solution of channel flow past a circularcylinder4.4 Aneurysm SimulationThe last example is the channel flow in the bloodvessel with aneurysm. The mesh is contributedby Kent-Andre Mardal, and we change it fromVMTK type into the mphtxt type for COMSOL.The inlet and outlet boundary can be easilydistinguished from “Walls” in Fig 9.Figure 9. Aneurysm SimulationExcerpt from the Proceedings of the 2014 COMSOL Conference in Shanghai

5. ConclusionThe weak form PDE and LiveLink forMATLAB in COMSOL give a perfectperformance during testing the three examplesby modified Uzawa method. This attempt givesalgorithm researchers a hint or a method to testtheir algorithm, and gives engineers another wayto compute their models especially those cannotfind appropriate built-in interface in COMSOL.6. References1.P. Chen, J. Huang, H. Sheng, Some Uzawamethods for steady incompressible Navier–Stokes equations discretized by mixedelement methods, Journal of Computationaland Applied Mathematics, 273 ， 313-325(2015)Excerpt from the Proceedings of the 2014 COMSOL Conference in Shanghai

3. LiveLink for MATLAB . LiveLink for MATLAB is a powerful function in COMSOL for designing iterative algorithm. We will let MATLAB do the “while” loop in Alg. 1. The way starting COMSOL with MATLAB on Windows / Mac OSX / Linux can be found in help document. Once they have been connected, the mph model can be handled in MATLAB so