A BDF2-Semismooth Newton Algorithm for the Numerical Solution of the Bingham Flow with Temperature Dependent Parameters

This paper is devoted to the numerical solution of the non-isothermal instationary Bingham flow with temperature dependent parameters by semismooth Newton methods. We discuss the main theoretical aspects regarding this problem. Mainly, we focus on existence of solutions and a multiplier formulation which leads us to a coupled system of PDEs involving a Navier-Stokes type equation and a parabolic energy PDE. Further, we propose a Huber regularization for this coupled system of partial differential equations, and we briefly discuss the well posedness of these regularized problems. A detailed finite element discretization, based on the so called (cross-grid P_1) - Q_0 elements, is proposed for the space variable, involving weighted stiffness and mass matrices. After discretization in space, a second order BDF method is used as a time advancing technique, leading, in each time iteration, to a nonsmooth system of equations, which is suitable to be solved by a semismooth Newton algorithm. Therefore, we propose and discuss the main properties of a SSN algorithm, including the convergence properties. The paper finishes with two computational experiment that exhibit the main properties of the numerical approach.

Authors

• 3 publications
10/21/2019

Finite Element Error Analysis of Surface Stokes Equations in Stream Function Formulation

We consider a surface Stokes problem in stream function formulation on a...
01/08/2021

Parallel Newton-Krylov-BDDC and FETI-DP deluxe solvers for implicit time discretizations of the cardiac Bidomain equations

Two novel parallel Newton-Krylov Balancing Domain Decomposition by Const...
01/25/2022

A Discontinuous Galerkin and Semismooth Newton Approach for the Numerical Solution of the Nonhomogeneous Bingham Flow

This paper is devoted to the study of non-homogeneous Bingham flows. We ...
04/09/2021

A Dual-Mixed Approximation for a Huber Regularization of the Herschel-Bulkey Flow Problem

In this paper, we extend a dual-mixed formulation for a nonlinear genera...
05/19/2021

Coupled-Cluster Theory Revisited

We propose a comprehensive mathematical framework for Coupled-Cluster-ty...
10/27/2021

Numerical simulation of multiscale fault systems with rate- and state-dependent friction

We consider the deformation of a geological structure with non-intersect...
03/03/2016

Using Newton's method to model a spatial light distribution of a LED with attached secondary optics

In design of optical systems based on LED (Light emitting diode) technol...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The numerical simulation of complex fluids has received an increasing amount of attention because of the wide variety of fields where these materials play central roles. In particular, viscoplastic materials are widely used in food industry (production of sauces and pastes), geophysics and other important fields of application. Recently, the focus on the heat transfer in these materials and its effects on the different material regions are of great interest. The development of theoretical and numerical tools to understand and simulate the effect of temperature fields on the evolution of rigid zones of the material has a particular impact in the understanding of several processes like transportation and exploitation of waxy-crude oils, simulation of lava flows, optimization of processes in industrial bakery, production and packing of jellies and honey-based products, etc. From the modelling perspective, this study is carried out by coupling the viscoplastic Bingham flow model with energy equations for the temperature fields.

In this interesting context, the analysis and simulation of viscoplastic flows with temperature dependent parameters represent a challenging and interesting field of research. From the analytical point of view, the viscoplastic model needs to be coupled with energy equations which are nonlinear and complex to analyze. Further, the temperature field should be inserted in the differential operators of the flow model, as functional weights. These models allow us to directly affect the mechanical properties of the material by changing the viscosity and the yield stress due to the action of temperature. This constitutes a first key step in the design of a control strategy for viscoplastic materials subject to heating from controlled heat sources.

There are several contribution on the theoretical analysis of the problem, starting with the classical work [7], where the problem was introduced and studied from the variational point of view. Recently, the steady flow of the Bingham model with temperature dependent parameters has been deeply analyzed in [3]. Besides a detailed discussion regarding the existence and uniqueness of solutions for the problem, the authors analyze the asymptotic behavior of the velocity field when the thermal conductivity of the material tends to infinity. This regime can be understood as a superfluid model, where the parameters stop depending on the geometry (non-local parameters). Following a similar perspective, in [11] the authors extend this analysis to the steady flow of the Herschell-Bulkley model for a specific shear-thinning setting.

In contrast with the analytical work developed, the numerical solution of these models is, to the best of our knowledge, scarce. In [14], the authors propose a numerical solution of this problem with the background of waxy crude oil flows. The authors consider that this kind of crude oils behaves as a viscoplastic material (Bingham model), and they numerically simulate the flow in a geometry representing a pipe. They propose a discretization based on a finite volume method. However, they dismiss the effect of the dissipation term in the energy equation, following mechanical properties of the flow under analysis. Finally, they proposed an Augmented Lagrangian method with an Uzawa type algorithm to solve the discretized system.

In [8], the authors introduce the so called Houska model, which represents a generalized approach to this kind of coupled problems. In this model, the energy equations is replaced by a transport equation and the flow is coupled with a structural parameter function which can be seen as temperature, concentration of certain components in mixtures-type materials, etc.

In this article, we study the Bingham flow with temperature dependent parameters, by analyzing the instationary Bingham model coupled with a parabolic partial differential equation for the temperature field. We first briefly discuss the variational formulation of this coupled fluid-heat system. Next, we discuss the existence of solutions for the coupled model and propose a multiplier formulation, which yields a system of equations formed by a Navier-Stokes type system (the Bingham model) and a parabolic PDE, coupled through the convection terms in both equations and a dissipation term in the energy equation. The analytical part of the article closes with the introduction of a smoothing procedure for the coupled multiplier system, based on a Huber approach. We propose such approximation to the problem since it has proven to be an efficient way to obtain fast convergent algorithms, maintaining an accurate approximation of the mechanical behavior of the flow (see [5]).

For the numerical solution of the problem, we aim at extending the approach developed in [6] to the case of the coupled heat-fluid flow under analysis. We start by discussing in detail a finite element discretization for the space variable. This discussion includes the construction of weighted stiffness and mass matrices for both the flow and the energy partial differential equations. Further, we discuss the time advancing scheme, considering a second order BDF time method, combined with a lag operator, as proposed in [2, 5]. In such a way, we obtain a fully discretized coupled fluid-heat system to be solved. The main characteristic of this system is that, because of the discretization approach, the system matrices are convective free. This means that the convective matrices can be constructed only with information of the velocity and temperature fields in previous time steps. Further, the coupled system, involves, at each time step, a semismooth system of equations related to the regularized Bingham model. The solution of this nonsmooth system is carried out with a semismooth Newton method, wich we propose and discuss. The main result is that the superlinear convergence of the algorithm holds in this case.

The paper is organized as follows. In Section 2

we introduce the problem given by a coupled system of nonlinear PDEs. Next, we briefly discuss the variational formulation of the PDEs system, which yields a coupled system of a variational inequality of the second kind for the velocity field and a parabolic nonlinear energy PDE for the temperature field. Finally, we characterize the solutions for this coupled system by the introduction of a tensor multiplier and propose a local smoothing procedure on the multiplier system based on the Huber approach. In Section

3, we discuss the fully discretization of the regularized system. First, we propose a stable finite element discretization for the space variable based on the so called (cross-grid )- elements. Next, we discuss the time advancing technique based on the application of a semi-implicit BDF2 method. Section 4 is devoted to the construction of the BDF2-SSN algorithm to numerically solve the coupled problem. Also, we construct the semismooth Newton inner algorithm to solve the flow equations. Finally, in Section 5, several numerical experiments are carried out in order to verify the theoretical properties of the proposed approach.

2 Problem Statement

Let , , be an open and bounded set with Lipschitz boundary . Let us assume that there exist , such that Given a final time , we define the following space-time sets and .

We are concerned with the non-isothermal flow of a Bingham fluid, considering that both the viscosity and the plasticity threshold depend on temperature. We assume the existence of volume forces acting on the fluid. The constitutive equations for this phenomenon are given by the following problem: find a velocity field , a pressure field and a temperature field such that

 ∂tu+(u⋅∇)u−∇⋅τ(θ)+∇p=f,in Ω∇⋅u=0,in Ωτ(θ)=μ(θ)Eu+g(θ)Eu∥Eu∥if Eu≠0,% ∥τ(θ)∥≤g(θ)if Eu=0,u=0,on Σ,u(x,0)=u0. (B)
 ∂tθ+u⋅∇θ−κΔθ=τ(θ):Eu−αθ,in Ω.∂θ∂n=0,on Γ0×(0,T)κ∂θ∂n+βθ=0,on Γ×(0,T)θ(x,0)=θ0. (E)

System () represents the Bingham flow model. However, in this case, both the viscosity and the yield stress depend on a temperature which is provided as the solution of the the energy equation (). The analysis and numerical solution of the coupled system ()-() is the main goal of this work. Here, stands for the thermal conductivity and for the volume forces. and are given parameters and, as usual, stands for the unit outward normal to . Notice that, in addition to the action provided by the dissipation energy factor , we admit the effect of a possible external heat source proportional to the temperature, if . Finally, we assume only nonslip boundary conditions for the flow velocity, while we impose a Robin boundary condition for the temperature in with given coefficient and Neumann homogeneous boundary conditions in and given initial conditions and for the velocity and temperature, respectively.

In the following section we analyze the variational formulation of the coupled system ()-(). Further, we discuss existence of solutions and propose a mixed formulation, based on a regularization approach, for this system.

2.1 Variational Formulation

We consider the Bingham incompressible flow model. Therefore, we introduce the divergence free spaces and . Next, the variational formulation of the flow model is given by the problem: find , a.e. in such that

 (∂tu(t),v−u)+∫Ωμ(θ)Eu(t):E(v−u(t))dx+∫Ω⟨(u(t)⋅∇)u(t),v−u(t)⟩dx+∫Ωg(θ)∥Ev∥dx−∫Ωg(θ)∥Eu(t)∥dx≥∫Ω(f,v−u(t))dx,∀v∈Vu(0)=u0 (VB)

One important issue regarding this formulation is the definition of functions and , which measure the dependence of the viscosity and the plasticity threshold with the temperature, respectively. These functions must be defined in such a way that the integrals in () are well posed. Usually, they are required to satisfy the following hypotheses [3, 7, 11]: , and there exist , and and such that

 μ0≤μ(λ)≤μ1<+∞,∀λ∈Rg0≤g(λ)≤g1<+∞,∀λ∈R. (1)

In our approach, we propose to follow a similar setting as the one given by the Houska model (see [8]). This model is a generalization of the coupled system studied here, where the flow parameters depend on a structure parameter which can be seen as several magnitudes, including the temperature. The main difference is that the Houska model does not include a diffusion term for the energy equation, which is replaced by a transport equation. In that model, the viscosity and the plasticity threshold are supposed to be affine functions of this structure parameter. Further, it is accurate for us to suppose that , since the fluid is suppose to be a Bingham material which will be altered by the action of temperature. Summarizing, we propose to analyse the problem considering the following functions

 g(λ):=⎧⎪⎨⎪⎩g0if λ≤0(g1−g0)λ+g0if 0<λ<1g1if λ≥1.andμ(λ):=⎧⎪⎨⎪⎩μ0if λ≤0(μ1−μ0)λ+μ0if 0<λ<1μ1if λ≥1. (2)

Considering these choices for and , it is clear that all the integrals in () are well posed, so the variational formulation for the flow system holds.

Now, let us focus on the variational formulation of the energy equation (). Let . Then, the variational formulation of the energy equation is given by: find a.e. in such that

 ∫Ω∂tθ(t)ϕdx+κ∫Ω(∇θ(t),∇ϕ)dx+β∫Γθ(t)ϕdx+∫Ω(u(t)⋅∇θ)ϕdx=∫Ω[μ(θ(t))∥Eu(t)∥2+g(θ(t))∥Eu(t)∥]ϕdx−α∫Ωθϕdx,∀ϕ∈W1,q′(Ω),θ(0)=θ0, (VE)

where . Note that the dissipated energy term makes sense only if , a.e. in . Therefore, we use the corresponding form of in the energy equation. Furthermore, the associated integral term in the variational formulation is well posed since for , we have that which implies that is continuously embedded in .

Theorem 2.1.

Let . Let us assume that , and . Then, there exist such that , and solutions of ()-().

Proof.

See [7, Th. 3.1] and [5, Th. 2.1]. ∎

Let us recall that we use the notation for the spaces .

2.2 Multiplier Approach

The use of tensor multiplier-type functions provide a versatile characterization for the solutions of problems involving variational inequalities. This approach has been used in previous contributions focused on Bingham flow (see [5, 4]). With such a characterization, a partial differential equation involving a multiplier, together with additional complementarity relations, is obtained. Further, in the case of non-isothermal flow, [7, Th. 3.2] has established the existence of such a multiplier for (). The resulting system, in variational form, reads as follows: find , and a.e. in such that.

 ∫Ω(∂tu(t),v)dx+∫Ω⟨(u(t)⋅∇)u(t),v⟩dx+∫Ωμ(θ)(Eu:Ev)dx+∫Ωg(θ)(q:Ev)dx−∫Ωp(t)∇⋅vdx=∫Ω(f(t),v)dx,∀v∈H10(Ω)∫Ωr∇⋅u(t)=0dx,∀r∈L20(Ω)∥q(x,t)∥≤g(θ),a.e. in Q(q(x,t):Eu(x,t))=g(θ)∥Eu(x,t)∥,a.e. in Qu(x,0)=u0. (MB)

In this system, it is clear that the pressure is recovered by a direct application of the D’Rahm’s Theorem, so the variational formulation is constructed by using the Sobolev space as the test space.

The active and inactive sets for the flow are defined, respectively, by

 A:={(x,t)∈Q:∥Eu(x,t)∥≠0}% and I:=Q∖A.

Further, since the multiplier is undetermined in the solid regions and not necessarily unique (see [7]), the computational approach depends on projection techniques (see [8, 14]). Moreover, instability in the numerical methods may occur due to the nonuniqueness of the multiplier. In this contribution, we propose to extend the approach based on local regularization techniques (see [4, 5]) which leads us to second order methods of semismooth Newton type.

2.3 Huber Regularization Approach

In this section, we introduce a family of regularized problems to approximate the solution of ()-(). This regularization approach is based on the so called Huber local smoothing of the stress tensor and it has proved to be efficient in the numerical solution of the Bingham flow in the stationary as well as in the instationary cases (see [6, 5]). Further, we have used this approach in the numerical solution of the convective flow of Bingham fluids in the Bousinessq paradigm [4]. The Huber regularization performs in a similar way as the bi-viscosity approach. Here, the active and inactive sets are approximated by sets depending on a regularization parameter. As the parameter grows, the regularized inactive set tends to the actual inactive set, which makes that the viscosity of the flow, in the solid regions, tends to infinity, which is the expected mechanical behavior.

For a parameter , the regularized problem consists in: find , , and a.e. in such that.

 ∫Ω(∂tuγ(t),v)dx+∫Ω⟨(uγ(t)⋅∇)uγ(t),v⟩dx+∫Ωμ(θ)(Euγ:Ev)dx+∫Ωg(θ)(qγ:Ev)dx−∫Ωpγ(t)∇⋅vdx=∫Ω(f(t),v)dx,∀v∈H10(Ω)∫Ωr∇⋅uγ(t)=0dx,∀r∈L20(Ω)qγ:=⎧⎨⎩g(θ(x,t))Euγ(x,t)∥Euγ(x,t)∥,a.% e. in AγγEuγ(x,t),a.e. in Iγuγ(x,0)=u0∫Ω∂tθ(t)vdx+∫Ω(uγ(t)⋅∇θ)vdx+κ∫Ω(∇θ(t),∇v)dx+β∫Γθ(t)vdx=∫Ω[μ(θ(t))∥Euγ(t)∥2+g(θ(t))∥Euγ(t)∥]vdx−α∫Ωθvdx,∀v∈W1,q′(Ω),θ(0)=θ0, (RB)

Here the regularized active and inactive sets are defined by

 Aγ:={(x,t)∈Q:γ∥Eu(x,t)∥≥g(θ(x,t))} and Iγ:=Q∖Aγ.

Note that, for given a.e. in , we have that a.e. in .

Given , we usually rewrite the equation for as follows

 max(g(θ(x,t)),γ∥Eu(x,t)∥)qγ(x,t)=γg(θ(x,t))Euγ(x,t),a.e. in% Q.
Remark 2.2.

Considering that existence of solutions for the energy equation is established, the existence of solutions for () follows from [5, Th. 3.1]. Furthermore, [5, Th. 3.3] guarantee that for a given , strongly in and weakly in . This convergence result shows that our approach is consistent and produce reliable approximations for the non-isothermal flow under study.

3 Space-Time Discretization

In this section, we propose a space-time discretization scheme for both the flow equation and the energy equation in the coupled system (). For the flow equation we use the scheme proposed in [5], which is based on a combination of a first-order finite element approximation for the space variable with (cross-grid )- elements and the semi-implicit BDF2 scheme for the time variable. This class of finite elements allows us to use the same test functions for the velocity gradient and the dual variable. In such a way, a direct relation between these two variables is obtained, and an accurate determination of active and inactive sets is achieved. On the other hand, the semi-implicit BDF2 is an advancing scheme that leads us to convection-independent systems in each time step. This is particularly useful to reduce the computational cost in every step of the Newton iteration.

For the energy equation we propose a first-order finite element method for the space variable and a BDF2 method for the time variable. In such a way, we also obtain a convection-independent system. In this case, the main advantage is that the fully discretized equation becomes linear and, consequently, computationally cheap to be solved. The temperature is discretized in the same nodes as the velocity, and a simple restriction method (weighting) is used to obtain the value of the temperature at each triangle.

3.1 Finite element discretization

We start the discussion of the space discretization by introducing the finite dimensional spaces corresponding to the so called (cross-grid )- finite elements, , and , as follows.

 Vh:=(Vh∩H1(Ω))2,% whereVh:={vh∈C(Ω):vh|T∈Π1, for all % T∈Th},Uh:=Rh∩L2(Ω),whereRh:={rh∈C(Ω):rh|Q∈Π0, for all Q∈Qh},Wh:={(qh1,qh2,qh3,qh4)⊤∈(L2(Ω))4:qhj|T∈Π0, for all T∈Th}.

Here, , and with , respectively. is a regular quadrangulation of , and is the regular triangulation obtained by dividing any square in by using its two main diagonals [12, Sec. 6]. To simplify the analysis, we assume that

has a polygonal boundary. The degrees of freedom associated to these elements are depicted in Figure

1.

It is remarkable that the (cross-grid )- elements satisfy the Ladyzhenskaya-Babuska-Brezzi (LBB) or inf-sup condition and, therefore, lead to a stable approximation of the Navier-Stokes-like systems, such as the Bingham model. (see [12, p. 435])

Further, for the energy equation, we discretize the temperature in the following finite dimensional space

 Xh:=W1,q(Ω)∩Vh and Xh′:=W1,q′(Ω)∩Vh

where is defined above. The degrees of freedom are shown in Figure 1, left.

By using the classical Galerkin approach, we obtain the following semi-discrete approximation for the coupled system ().

 Mh∂t→u(t)+Ch(→u(t))→u(t)+Ahμ(→θ(t))→u(t)+Bh→p(t)+Qhg(→θ(t))→q(t)=→f(t)−(Bh)⊤→u(t)=0,max[GT(→θ(t)),γNh(Eh→u(t))]⋆→q(t)=γdiag(GT(→θ(t)))Eh→u(t),a.e. in [0,Tf] and ∀T∈Th→u(0)=→u0. (3)
 Mh∂t→θ(t)+Ch(→u(t))→θ(t)+κAh→θ(t)+αMh→θ(t)+βMhΓ→θ(t)=[Khμ(→u(t))+Khg(→u(t))]→θ(t)→θ(0)=→θ0. (4)

Here,

denotes a vector whose components are the values of

in the center of gravity of each . , , , , and are the time-dependent vectors of coefficients in the finite element representation of the 4-tuple , and the initial conditions and , respectively. and are the mass matrices for and , respectively, while stands for the stiffness matrix associated to . is the boundary mass matrix constructed for the Robin boundary condition of the energy equation (see [10, Sec. 4.6.2]). Matrix is obtained in the usual way from the bilinear form .

The discretized convective matrices and are given by

 Ch(→u(t))ij:=2d∑k=1uk∫Ω⟨(φk⋅∇)φj,φi⟩dx and Ch(→u(t))ij:=2d∑k=1uk∫Ω(φk⋅∇ϕj)ϕ′idx,

where , are the basis functions of , , are the basis functions of and , are the basis functions of , respectively. The discrete approximation, , for the deformation tensor and the right hand side are constructed by using the basis functions , (see [6, Sec. 5]). Finally, the function is defined by

 Nh(q)i=Nh(q)i+m=⋯=Nh(q)i+4m:=|(qi,qi+m,...,qi+4m)|,

for and . The values of and are given in the gravity centers of each (see Figure 1).

Let us now explain the matrices and . These matrices are defined as follows

 Ahμ(→θ)ij:=∑T∈Th∫TμT(→θ(t))(Eφi:Eφj)dx and Qhg(→θ)ij:=∑T∈Th∫TgT(→θ(t))(ψi:Eφj)dx,

where, as before, , are the basis functions of and , , are the basis functions of . Here, represents the application of a simple weighting operator to obtain the value of function in the gravity center of each triangle from the values of at each vertex of , a.e. in . The same applies for . This is depicted in Figure 2.

Finally, let us discuss the space discretization for the dissipation term

 ∫Ω[μ(θ(t))∥Euγ(t)∥2+g(θ(t))∥Euγ(t)∥]ϕdx.

Since we are working with affine function on (see (2)), we can decompose any term in the integral above as follows.

 ∫Ωμ(θ(t))∥Euγ(t)∥2vdx=(μ1−μ0)∫Ω∥Euγ(t)∥2θ(t)vdx+μ0∫Ω∥Euγ(t)∥2vdx,

and

 ∫Ωg(θ(t))∥Euγ(t)∥vdx=(g1−g0)∫Ω∥Euγ(t)∥θ(t)vdx+g0∫Ω∥Euγ(t)∥vdx.

Therefore, by following the classical Galerkin method, the discretization of this terms reads as follows

 Khμ(→u(t)):=(μ1−μ0)Mh,2(→u(t))→θ(t)+μ0Θh,2(→u(t)),

where

 Mh,2(→u(t))i,j:=∑T∈Th∫T[Nh(Eh→u(t))]2T(ϕiϕ′j)dx, and (Θh,2)i:=16|T|[Nh(Eh→u(t))]2T,

for all . Here stands for approximated value of at each , and is the measure of the given triangle. The components of vector are constructed by following the approximation given in [1, Sec. 6].

By using the same argumentation, we have that

 Khg(→u(t)):=(g1−g0)Mh,1(→u(t))→θ(t)+g0Θh,1(→u(t)),

where

 Mh,1(→u(t))i,j:=∑T∈Th∫T[Nh(Eh→u(t))]T(ϕiϕ′j)dx, and (Θh,1)i:=16|T|[Nh(Eh→u(t))]T,

for all . Here stands for approximated value of at each .

Remark 3.1.

By construction, matrices , , and are weighted stiffness and mass matrices, respectively. In the particular case of the dissipation term, it is possible to obtain the respective weighted mass matrices because of the specific form of the functions and . Other class of functions need a further analysis.

Usual time advancing techniques, such as the one-step -method (see [12]), applied to Navier–Stokes-type equations lead to the numerical solution of nonlinear and convective systems of algebraic equations, which change in every time step. This fact provokes an increase of the computational cost. One approach to avoid this issue is to use operator splitting techniques (see [13]). However, the use of such methods needs the solution of several other systems to construct the solution of the problem. Another approach is the use of semi implicit methods. One important characteristic of some of these methods is that they lead to Stokes-type matrices with no convective term active in every time step (see [2, 5, 4]). In this article, we focus on a semi-implicit method proposed in [5] for the Bingham flow, based on the second-order backward differentiation formulae (BDF2) and on the introduction of a lag-operator. This approximation enjoys the same kind of property: a system whose associated matrix is a Stokes-type one and does not change in every time step.

We discuss separately the time advancing for the flow and the energy equation. Lest us start by the flow equation. By following [5, Sec. 4.B.], we formulate the BDF2 approximation for (3), as follows: given a vector , at each time level , for , solve the system

 (32δtMh+Ahμ(→θ))→uk+1+Bh→pk+2+Qhg(→θ)→qk+2=→Fk+2−(Bh)⊤→uk+2=0max[GT(→θ),γNh(Eh→uk+2)]⋆→qk+2=γdiag(GT(→θ))Eh→uk+2, (5)

where represents the approximation of . The right hand side is given by

 →F:=→fk+2−Ch(Λ(→uk))Λ(→uk)+2δtMh→uk+1−12δtMh→uk.

Here, stands for the lag operator and it is defined by

 Λ(→uk):=2→uk+1−→uk.

The BDF2 scheme combined with the lag operator allows us to approximate the convection matrix with information given by the function in the two previous time steps, which implies that this matrix can be moved to the right hand side of the system.

Next, regarding the initialization of this scheme, having the discretized initial condition , we calculate two intermediate steps and , by using a Crank-Nicholson type scheme. Next, we define . This approach guarantees that the second order approximation in time holds. For further details, we recommend the reader to see [5, Sec. 4.B] and [2, p. 371].

Let us now focus on the energy equation. Given , we use the same approach based on the BDF2 method and the lag operator to obtain the following system at each time level , for .

 (32δtMh+κAh)θk+2+αMhθk+2+βMhΓθk+2−[(μ1−μ0)Mh,2(→u)+(g1−g0)Mh,1(→u)]→θk+2=→Fk+2, (6)

where represents the approximation of , and the right hand side is given by

 →Fk+2:=Ch(Λ(→uk))Λ(→θk)+2δtMhθk+1−12δtMh→θk+μ0Θh,2(→u)+g0Θh,1(→u)

The initialization follows from the discretized initial condition and the calculation of by using a Crank-Nicholson scheme. Due to the use of the lag operator, this discretized equation is a linear equation, which does not need the application of Newton type methods.

4 Combined BDF2-SSN Algorithm

In this section, we discuss the combined BDF2-Semismooth Newton Algorithm to solve the system (5)-(6). We propose a sequential algorithm, which means that we solve the flow equation with a given temperature field, and then we update the temperature with the resulting velocity field. Due to the discretization scheme, the fully discretized energy equation (6) is now a linear equation, which implies that its solution depends only on the nonsigularity of the system matrix. On the other hand, the flow equation requires a nonlinear algorithm to be solved. As stated before, we propose a semismooth Newton algorithm to find the numerical solution of system (5).

Let us show the proposed combined BDF2-SSN algorithm for the non-isothermal Bingham flow with temperature dependent parameters in the time interval .

Algorithm 4.1.

(non-isothermal BDF2-SSN)

1. Initialization: Given and , calculate . Introduce in (6), calculate , and set .

2. For do

1. Flow update: Given , obtain by applying the SSN algorithm to the following system

 (32δtMh+Ahμ(→θk+1))→uk+2+Bh→pk+2+Qhg(→θk+1)→qk+2=→Fk+2−(Bh)⊤→uk+2=0max[GT(→θk+1),γNh(Eh→uk+2)]⋆→qk+2=γdiag(GT(→θk+1))Eh→uk+2, (7)
2. Temperature update: Use the calculated to obtain by solving

 (32δtMh+κAh+αMh+βMhΓ)θk+2−[(μ1−μ0)Mh,2(→uk+2)+(g1−g0)Mh,1(→uk+2)]→θk+2=→Fk+2, (8)
Remark 4.2.

We have rewritten the fully discretized version of (6) to emphasize the fact that the solution of this system depends only on the non singularity of the system matrices. This fact directly follows from the positive definiteness of mass and stiffness matrices [12, p.148]. Further, this property is shared by the weighted matrices and , since the weights are positive numbers which are constant in every .

Note that, the Remark 4.2 implies that the success of Algorithm 4.1 depends only on the convergence of the inner SSN algorithm developed to numerically solve (7).

4.1 Semismooth Newton Algorithm

The main difficulty regarding system (