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 waxycrude oils, simulation of lava flows, optimization of processes in industrial bakery, production and packing of jellies and honeybased 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 (nonlocal parameters). Following a similar perspective, in [11] the authors extend this analysis to the steady flow of the HerschellBulkley model for a specific shearthinning 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 mixturestype 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 fluidheat 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 NavierStokes 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 heatfluid 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 fluidheat 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 (crossgrid ) elements. Next, we discuss the time advancing technique based on the application of a semiimplicit BDF2 method. Section 4 is devoted to the construction of the BDF2SSN 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 spacetime sets and .
We are concerned with the nonisothermal 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
() 
() 
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
() 
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
(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
(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
() 
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 .
Let us recall that we use the notation for the spaces .
2.2 Multiplier Approach
The use of tensor multipliertype 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 nonisothermal 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.
() 
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
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 biviscosity 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.
() 
Here the regularized active and inactive sets are defined by
Note that, for given a.e. in , we have that a.e. in .
Given , we usually rewrite the equation for as follows
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 nonisothermal flow under study.
3 SpaceTime Discretization
In this section, we propose a spacetime 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 firstorder finite element approximation for the space variable with (crossgrid ) elements and the semiimplicit 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 semiimplicit BDF2 is an advancing scheme that leads us to convectionindependent 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 firstorder finite element method for the space variable and a BDF2 method for the time variable. In such a way, we also obtain a convectionindependent 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 (crossgrid ) finite elements, , and , as follows.
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 (crossgrid ) elements satisfy the LadyzhenskayaBabuskaBrezzi (LBB) or infsup condition and, therefore, lead to a stable approximation of the NavierStokeslike 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
where is defined above. The degrees of freedom are shown in Figure 1, left.
By using the classical Galerkin approach, we obtain the following semidiscrete approximation for the coupled system ().
(3) 
(4) 
Here,
denotes a vector whose components are the values of
in the center of gravity of each . , , , , and are the timedependent vectors of coefficients in the finite element representation of the 4tuple , 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
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
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
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
Since we are working with affine function on (see (2)), we can decompose any term in the integral above as follows.
and
Therefore, by following the classical Galerkin method, the discretization of this terms reads as follows
where
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
where
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.
3.2 Time advancing method
Usual time advancing techniques, such as the onestep method (see [12]), applied to Navier–Stokestype 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 Stokestype matrices with no convective term active in every time step (see [2, 5, 4]). In this article, we focus on a semiimplicit method proposed in [5] for the Bingham flow, based on the secondorder backward differentiation formulae (BDF2) and on the introduction of a lagoperator. This approximation enjoys the same kind of property: a system whose associated matrix is a Stokestype 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
(5) 
where represents the approximation of . The right hand side is given by
Here, stands for the lag operator and it is defined by
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 CrankNicholson 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 .
(6) 
where represents the approximation of , and the right hand side is given by
The initialization follows from the discretized initial condition and the calculation of by using a CrankNicholson 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 BDF2SSN Algorithm
In this section, we discuss the combined BDF2Semismooth 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 BDF2SSN algorithm for the nonisothermal Bingham flow with temperature dependent parameters in the time interval .
Algorithm 4.1.
(nonisothermal BDF2SSN)

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

For do

Flow update: Given , obtain by applying the SSN algorithm to the following system
(7) 
Temperature update: Use the calculated to obtain by solving
(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 (
Comments
There are no comments yet.