3D Navier-Stokes + Heat Transfer

[bragostin]

Dear Fenics Community,

I am trying to simulate the classical laminar flow with heat transfer in a square duct with fixed temperature at the wall, which should yield a Nusselt number of 2.98, but I cannot make it work. It converges, but the values are clearly wrong. The code I am using is below:

from dolfin import *
from msh2xdmf import import_mesh_from_xdmf, msh2xdmf

# MAKE CFD
def sim_flow(u_0, nu, cp, k, p_0, T_0, fileName, Le, He, We):
    
    # LOAD MESH
    
    mesh, boundaries, association_table = import_mesh_from_xdmf(prefix=fileName, dim=3)
    
    # Build function space
    P2 = VectorElement('Lagrange', mesh.ufl_cell() , 2)
    P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
    element = MixedElement([P2, P1, P1])
    W = FunctionSpace(mesh, element)
    #define test functions
    (v, q, s) = TestFunctions(W)
    #split functions
    upT = Function(W)
    (u, p, T) = split(upT)

    n = FacetNormal(mesh)
    
    # Define boundary conditions
    # FLOW
    # Define inflow profile from Shah & London 1976 (velocity profile in a rectangular duct)
    alpha = min(We/He, He/We) # <1
    mv = 1.7 + 0.5 * alpha**-1.4
    if alpha<=1./3.:
        nv = 2
    else:
        nv = 2 + 0.3 * (alpha - 1./3.)
    Umax = u_0 * (mv+1)/mv * (nv+1)/nv
    print("alpha, m, n, Umax = ", alpha, mv, nv, Umax)
    inflow_profile = ('0', '0', 'Umax * (1. - pow(abs(x[0]/H*2), n)) * (1. -  pow(abs(x[1]/W*2), m))')
    inflow_profile = Expression(inflow_profile, Umax=Constant(Umax), H=Constant(He), W=Constant(We), m=Constant(mv), n=Constant(nv), degree=2)
    bcu_inflow = DirichletBC(W.sub(0), inflow_profile, boundaries, association_table["inlet"])
    bcu_noslip = DirichletBC(W.sub(0), Constant((0, 0 ,0 )), boundaries, association_table["noslip"])
    bcu_outflow = DirichletBC(W.sub(1), Constant(p_0), boundaries, association_table["outlet"])
    bcu = [bcu_inflow, bcu_noslip, bcu_outflow]
    
    # ENERGY
    bcT_inflow = DirichletBC(W.sub(2), Constant(T_0), boundaries, association_table["inlet"])
    bcT_noslip = DirichletBC(W.sub(2), Constant(T_0+10.), boundaries, association_table["noslip"])
    bcT = [bcT_inflow, bcT_noslip]
    
    bcs = bcu + bcT
    
    # DEFINE VARIATIONAL FORM
     
    F1 = (inner(grad(u)*u, v)*dx +                 # Momentum advection term
        nu*inner(grad(u), grad(v))*dx -          # Momentum diffusion term
        inner(p, div(v))*dx +                    # Pressure term
        div(u)*q*dx                            # Divergence term
    ) 
    F2 = (cp*inner(dot(grad(T), u), s)*dx - # Energy advection term
        k*inner(grad(T), grad(s))*dx # Energy diffusion term
    )
    F = F1 + F2
    
    # Create VTK files for visualization output
    vtkfile_u = File('results/u.pvd')
    vtkfile_p = File('results/p.pvd')
    vtkfile_T = File('results/T.pvd')
    
    J = derivative(F, upT) # Jacobian

    solve(F == 0, upT, bcs)
    
    # Save solution to file (VTK)
    (u, p, T) = upT.split(deepcopy=True)
    vtkfile_u << u
    vtkfile_p << p
    vtkfile_T << T
    
    # POST-PROCESSING
    ds_bc = ds(subdomain_data=boundaries)
    
    T_in_avg = assemble(T*ds_bc(association_table["inlet"])) / (We * He)
    T_out_avg = assemble(T*ds_bc(association_table["outlet"])) / (We * He)
    DT_avg = T_out_avg - T_in_avg
    Heat_Load = u_0 * We * He * rho * cp * DT_avg
    
    htc = project(k*grad(T), VectorFunctionSpace(mesh, 'Lagrange', 1))
    s_htc = 2 * (We + He) * Le
    Tbulk = (T_0 + DT_avg) / 2.
    htc_avg = assemble(dot(n, htc)/(T-Tbulk)*ds_bc(association_table["noslip"])) / s_htc
    
    Nu = htc_avg * Dh / k
    print("htc_avg = ", htc_avg)
    print("Nu = ", Nu)
    
    return Nu


# Dimensions
Le = 55e-3
He = We = Dh = 2.8e-3

# Fluid properties
rho = 1.127 # air density at 20 degC, 1 atm, [kg/m3]
nu  = 16.92E-6 # air kinematic viscosity at 20 degC, 1 atm, [m2/s]
cp = 1008. # air heat capacity @ 40°C (J/kg K) 
k = 27.35e-3 # air thermal conductivity @40°C (W/m/K) 
p_0 = 0. # outlet air pressure (atmospheric pressure), normalized
T_0 = 0. # Inlet temperature (K) 
    
fileName = "SquareDuct"

u_0 = 5.67 # Inlet velocity (m/s)
Re  = u_0 * Dh / nu
print("Dh = ", Dh)
print("Re = ", Re)

res = sim_flow(u_0, nu, cp, k, p_0, T_0, fileName, Le, He, We)

The velocity profile looks good with 0 at the wall and maximum at center, but the temperature profile below does not look physical. The mesh file is here https://github.com/bragostin/Fenics. Is there anything obviously wrong in this code?

[kamensky]

One thing that stands out when skimming the code is that the advection–diffusion problem for T is advection-dominated for the given values of cp and k. The basic Galerkin method is unstable in that regime. (Typically, you will get bad results unless the Péclet number is small, taking the element diameter in the mesh as a length scale. Thus, for a sufficiently-fine mesh, you can get still good results, but it’s often not practical.) A standard way to stabilize advection–diffusion is to use a formulation called the streamline-upwind Petrov–Galerkin (SUPG) method, which is implemented for advection–diffusion in an undocumented demo distributed with DOLFIN. (Given that nu is also very small, I would also recommend stabilizing advection in the Navier–Stokes subproblem. For Galerkin’s method to reliably produce good results for Navier–Stokes, you want the Reynolds number to be small when taking the element diameter as a length scale, which is also impractical for many common flows.)

The original reference on SUPG from the 1980s includes the extension to incompressible Navier–Stokes. I cover the theory behind SUPG from the more recent variational multiscale (VMS) perspective in Chapter 6 of the lecture notes I’m writing for my graduate course.

[bragostin]

Thank you very much, I am going to try SUPG and post the results here. And thank you for the link to the course, I lack the theoretical background to fully understand this method, so this is perfect!

[bragostin]

If I understand correctly the SUPG method artificially increases the diffusion part to stabilize the solution, and since this is multiplied by the residual it should disappear as it converges to the stable solution. Following the available examples I have started to implement SUPG for the velocity part of my problem by adding to code above:

    a = Constant((0,0,u_0))
    k = 2
    # Formulation:
    res_strong = dot(a,grad(v)) - div(nu*grad(v))
    h = CellDiameter(mesh)
    Cinv = Constant(16*k**2)
    vnorm = sqrt(dot(a, a))
    tau = Min(h**2/(Cinv*nu), h/(2*vnorm))
    res_SUPG = tau*inner(dot(a, grad(u)), res_strong)*dx
    F += res_SUPG

I have also initialized the velocity with

    # Define initial conditions
    e_u0 = Expression(('0.', '0.', '5.67'), degree=1)
    e_p0 = Expression('0.', degree=1)
    e_T0 = Expression('0.', degree=1)
    u0 = interpolate(e_u0, W.sub(0).collapse())
    p0 = interpolate(e_p0, W.sub(1).collapse())
    T0 = interpolate(e_T0, W.sub(2).collapse())
    assign(upT, [u0, p0, T0])

u_0=5.67 is my inlet velocity, which in z direction only.

SUPG actually makes diverging what was converging before, so I am doing it wrong I guess. I tried changing Cinv value but it does not make it converge. Any clue?

[kamensky]

There are a number of issues with your SUPG term, although, looking at the original code again, I think I may have spotted a more significant typo: the sign of the diffusion term in F2 is incorrect, and corresponds to anti-diffusion, which is inherently unstable.

It is still a good idea to stabilize advective problems, though. Some notes on your attempt at the SUPG term:

• The advection velocity a should be the (unknown) fluid velocity, for both the Navier–Stokes and advection–diffusion subproblems.
• The strong residuals for these subproblems should also be written in terms of the unknown solutions, i.e., \mathbf{u}\cdot\nabla\mathbf{u} - \nabla\cdot(2\nu\nabla^s\mathbf{u} - p\mathbf{I}) - \mathbf{f} for the Navier–Stokes subproblem and c\rho\mathbf{u}\cdot\nabla T - \kappa\nabla^2T - Q for the advection–diffusion problem (where \mathbf{f} and Q are source terms I’ve added for completeness, even if they’re zero in your problem).
• Even with an inf-sup-stable pressure–velocity pair, it is a good idea to use PSPG in conjunction with SUPG for Navier–Stokes, so there is a coercive \nabla p–\nabla q term to “hide” non-coercive \nabla p–(other stuff) cross-terms behind, at least if one follows Franca’s original error analysis. (Actually, I was involved in a collaboration recently on a method which lets you keep (q,\nabla\cdot\mathbf{u}) = 0~~\forall q with PSPG-like stabilization, but at the price of introducing a new unknown field.)
• I have SUPG/PSPG for Navier–Stokes implemented using FEniCS for this example from my class, albeit in the more complicated setting of fluid–structure interaction (and without dividing Navier–Stokes through by \rho).

[bragostin]

Thank you very much for spending the time to explain these things and giving expert advice! I have worked on the problem and tried to implement it as close as possible to the latest example of your class.

1. Setting the sign of the diffusion term in F2 to + solved the issue with the temperature distribution which is now physical, I should have spotted this one, thank you! My system is now:

    F1 = (rho*inner(grad(u)*u, v)*dx +                 # Momentum advection term
        mu*inner(grad(u), grad(v))*dx -          # Momentum diffusion term
        inner(p, div(v))*dx +                    # Pressure term
        div(u)*q*dx                            # Divergence term
    ) 
    F2 = (k*inner(grad(T), grad(s))*dx + # Energy advection term
        rho*cp*inner(dot(grad(T), u), s)*dx # Energy diffusion term       
    )
    F = F1 + F2

2. Indeed with higher velocity the solver diverges so that SUPG/PSPG is needed. I came with the following implementation, using the unknown velocity u:

    sigma = 2.*mu*sym(grad(u)) - p*Identity(len(u))
    res_strong = rho*dot(u, grad(u)) - div(sigma)
    k = 2
    Cinv = Constant(16*k**2)
    vnorm = sqrt(dot(u, u)) 
    tau_SUPG = Min(h**2/(Cinv*nu), h/(2*vnorm))
    F_SUPG = inner(tau_SUPG*res_strong, rho*dot(grad(v),u))*dx
    tau_PSPG = tau_SUPG
    F_PSPG = tau_PSPG*inner(grad(q), res_strong)*dx
    F = F + F_SUPG + F_PSPG

Then the solver becomes extremely slow, diverges and I get this message below.
So I guess I am doing something wrong still.

  WARNING: The number of integration points for each cell will be: 343
           Consider using the option 'quadrature_degree' to reduce the number of points
Calling FFC just-in-time (JIT) compiler, this may take some time.
  WARNING: The number of integration points for each cell will be: 3375
           Consider using the option 'quadrature_degree' to reduce the number of points.
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 4.958e+03 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 7.215e-01 (tol = 1.000e-10) r (rel) = 1.455e-04 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 6.061e+05 (tol = 1.000e-10) r (rel) = 1.222e+02 (tol = 1.000e-09)```

[kamensky]

Indeed with higher velocity the solver diverges so that SUPG/PSPG is needed.

SUPG/PSPG won’t necessarily help with nonlinear convergence, but it may help with the accuracy of the resulting converged solution. For difficult nonlinear problems, you might look into using the PETSc SNES solver with a backtracking line search. This can be accessed through FEniCS’s Python API as shown in this discussion.

Then the solver becomes extremely slow, diverges and I get this message below.

By default, FEniCS will try to estimate an appropriate polynomial degree to integrate exactly with Gauss quadrature. However, for complicated variational forms (such as those obtained from stabilized methods), this estimate is usually overkill. For complicated problems, it’s usually a good idea to manually set the degree of numerical quadrature, e.g.,

dx = dx(metadata={"quadrature_degree":2*k})

(and similar for ds and dS, if applicable) which is what the warning message is suggesting to do.

Also: Another type of stabilization that can add robustness for incompressible Navier–Stokes is “least squares on the incompressibility constraint” (LSIC; called “grad-div” in some references), which is also defined in my notes. It essentially penalizes the deviations from \nabla\cdot\mathbf{u} = 0 that can occur at individual points when incompressibility is only enforced weakly with respect to some finite-dimensional test space.

[bragostin]

So I think I have it rounded now thanks to your advice:

1. I added LSIC and to my surprise it made the solver converge at higher velocites (2x5.67 m/s):

    tau_LSIC = rho * 2*nu/3
    F_LSIC = tau_LSIC*div(u)*div(v)*dx2
    F = F + F_LSIC

2. In order to make SUPG/PSPG converge I had to use this formula with Reynolds number for Cinv:

    Cinv = Constant(16*Re)

It seems that higher values of Cinv are needed at high velocites in this problem. It yields Nu=4.88 instead of Nu=4.74, so the accuracy is indeed changed.

3. The flow is hydrodynamically developped (fully developped velocity profile imposed at inlet) and thermally developping. The script yields a Nusselt number of 4.88 compared to 6.44 in Shah & London, Laminar Flow Forced Convection in Ducts, 1978. This difference may come from the meshing of the square duct: a better meshing at the wall where the temperature gradient is calculated is probably needed.

4. If it can be useful to anybody, the whole script and mesh are available here https://github.com/bragostin/Fenics

Thanks again for your help, I am glad it works!

[kamensky]

I added LSIC and to my surprise it made the solver converge at higher velocites (2x5.67 m/s)

The LSIC constant, often denoted \tau_\text{C}, is typically selected as

\tau_{\text{C}} = h^2/\tau_{\text{M}}\text{ ,}

which is motivated by the original error analysis for a linearized model problem.

In order to make SUPG/PSPG converge I had to use this formula with Reynolds number for Cinv

There should be some velocity-independent value for C_\text{inv}. (In principle, you can compute it exactly by solving a local eigenvalue problem in each element, as discussed in this paper, but, in practice, everyone just uses some heuristic formula in terms of the element degree). For purposes of nonlinear convergence, it can help to use “smoothed” versions of the \text{min} formula that comes from the analysis, e.g.,

\text{min}\,\{a,b\}~~~\to~~~\left(a^{-2} + b^{-2}\right)^{-1/2}\text{ ,}

and generalizations of “h” based on the full Jacobian of the mapping from the parent element to physical space (as in the FSI example I linked above), but this is, again, a separate issue from accuracy of the resulting converged algebraic solution.

I will mention, though, that when I’ve previously tried using an SUPG/PSPG(-ish) formulation with Taylor–Hood elements, I needed larger values of C_{\text{inv}} than I was used to using. (This type of stabilization is mostly used with equal-order interpolation of velocity and pressure, because PSPG stabilization means you no longer need to worry about inf-sup stability of the velocity–pressure pair.)

[bragostin]

1. With tau_LSIC = h**2/tau_SUPG it works as well, thank you.
2. In Cinv = Constant(16*Re) is actually not function of the unknown velocity u, but of u_0 the inlet velocity because here Re = u_0 * D / nu
3. I had to use dx = dx(metadata={"quadrature_degree":2*3}) because with 2*2 the solver spitted out nan

[kamensky]

@bragostin

In Cinv = Constant(16*Re) is actually not function of the unknown velocity u , but of u_0 the inlet velocity because here Re = u_0 * D / nu

It should still, in principle, be independent of any problem data. In the analysis, C_\text{inv} comes from a so-called “inverse estimate” (which is what the subscript “\text{inv}” stands for). Inverse estimates are properties associated with the finite element function spaces themselves, and should therefore only depend on things like polynomial degree.

[bragostin]

If I take tau_SUPG = Min(h**2/(Cinv*nu), h/(2.*vnorm)) and substitute with Cinv=16*u*D/nu there I get Min(h/(8D)*h/(2.*vnorm), h/(2.*vnorm)), so it’s essentially like having the second term h/(2.*vnorm) multiplied with h/(8D). I don’t know if that makes any sense though.

[bragostin]

I worked further on implementing the stabilization scheme for the energy problem and came to this:

    res_strong = rho*cp*dot(u, grad(T)) - k*div(grad(T))
    tau_E = h**2 * (rho*cp/k) / 1000.
    F_E = inner(tau_E*res_strong, rho*cp*dot(u,grad(s)))*dx2
    F = F + F_E

Setting tau_E=tau_SUPG made the solver return Nan.
It’s converging with and without the stabilization, but to very different temperature profiles:
Without stabilization, yielding Nu=6.42:

With stabilization, yielding Nu=18.2:

with

    T_in_avg = assemble(T*ds_bc(association_table["inlet"])) / (We * He)
    T_out_avg = assemble(T*ds_bc(association_table["outlet"])) / (We * He)
    DT_avg = T_out_avg - T_in_avg
    Heat_Load = u_0 * We * He * rho * cp * DT_avg
    LMDT = ((Tw - T_in_avg) - (Tw - T_out_avg)) / ln((Tw - T_in_avg) / (Tw - T_out_avg))
    htc_avg = Heat_Load / (Area * LMDT)
    Nu = htc_avg * Dh / k

One of the 2 solutions is not correct. The one without stabilization at Nu=6.42 is correct I would think.

[kamensky]

The choice of \tau_\text{E} looks off to me. (As a sanity check, look at the units of the SUPG term, compared with the rest of the formulation.) I would expect something more like

\tau_\text{E}\sim\operatorname{min}\left\{\frac{h}{2\rho c_p\vert\mathbf{u}\vert},\frac{h^2}{C_{\text{inv}}\kappa}\right\}\text{ ,}

by analogoy to advection–diffusion.

The nan result is likely due to divide-by-zero when \vert\mathbf{u}\vert = 0. This is avoided in the “smoothed” \text{min} I mentioned above, but could also be fixed with a +DOLFIN_EPS in the denominator.

[bragostin]

Yes this makes sense of course, and it reduces a lot the difference between the solution with and without stabilization. Then the final solution is sensitive to the value of Cinv. Is it to be expected?

[kamensky]

If you let C_\text{inv}\to\infty, it essentially switches off stabilization, since you’ll revert to Galerkin’s method. If you follow the analysis of Franca et al. in detail, C_\text{inv} would correspond to 2/m_k in their notation, so you should use C_\text{inv}=6 for P1 elements (because \tilde{C}_k can be arbitrarily large for P1 elements; my current notes are a bit sloppy in letting C_\text{inv} = \tilde{C}_k^{-1}, which can run into problems when taking \tilde{C}_k as large as possible for low-order elements).

[bragostin]

Evidently I need to study more, thank you for the references and your help! I have a working script now. It is a very interesting topic.