Project topics and materials
Get Complete Project Material File(s) Now! »
The Extended Finite Element Method
Literature Review :
Several methods are found in the literature to numerically model physical phenomena. The particular physics studied in this project (heat transfer and fluid dynamics) and the continuous nature of the material encourages the use of the finite element method (FEM [66] as a modelling tool [10, 27, 35, 45]. This chapter explores the main approaches used in the finite element method to model phase change problems and handle discontinuous interfaces in general. Section 2.2 gives a general description of these techniques, their advantages and their drawbacks. In sections 2.3 and 2.4, a review of their application to the Stefan and Navier-Stokes problems, respectively, are given as well as advantages and drawbacks associated to the specific problem.
Interface Handling Techniques :
The finite element method has been proven to be robust for a variety of problems, is geometrically flexible and is well understood in the literature. However, particular situations can still cause problems. The finite element method uses a polynomial interpolation within individual elements. Consequently, it can only reproduce a function with discontinuities in the variable (strong discontinuity) or its derivative (weak discontinuity) by separating the domain into submeshes, treated separately. This means that the discontinuity’s location must be known a priori in order to properly adapt the mesh to this additional geometrical constraint.
In this case, the interface moves within the mesh and its management becomes very complex and often very difficult from a practical point of view. To handle such situations, three adaptations of the finite element algorithm have been developed in the literature. In the first, the interface is modelled in a diffused manner, which converts the discontinuity into a steep (but finite) gradient [9, 10, 79]. Secondly, the interface can be defined by a set of nodes on the mesh. Mesh control techniques allow the set of nodes to move within the mesh, based on the desired interface velocity, and reorganize the surrounding nodes to maintain an acceptable mesh distortion. This approach is called the arbitrary Lagrangian Euleurian method (ALE) [37, 39, 51]. Lastly, a technique called the extended finite element method (XFEM) provides the ability to have discontinuous solutions within elements [13, 21, 22, 56]. To do so, the Lagrange polynomial interpolation is enriched with a carefully selected function to introduce a discontinuous behavior in the solution at an arbitrary location, independent of the mesh.
Arbitrary Lagrangian-Eulerian Method :
The arbitrary Lagrangian-Eulerian (or ALE) method [15, 37, 38, 44] offers an interesting alternative to diffused interfaces because it stores the interface on mesh nodes. This results in an exact representation of the discontinuity because the finite element method can include discontinuities in the solution that occur at element boundaries. The main idea behind this method is to combine the advantages of a Lagrangian reference frame, which allows the mesh nodes to move, with those of an Eulerian reference frame with a fixed mesh (better suited for problems involving convection, such as the Navier-Stokes equations). The result is a system of equations in an Eulerian reference which takes into account the movement of the nodes.
The Extended Finite Element Method :
The objective of this chapter is to present the extended finite element method. Section 3.1 gives an overview of the method. Section 3.2 describes the different enrichment schemes used in the work to reproduce the various physical discontinuities of the Stefan problem and Navier-Stokes equations. Section 3.3 explains the penalty and Lagrange multiplier methods which are used to apply the Dirichlet boundary conditions on the interface. Section 3.4 details the integration of the finite element formulation in cut elements, particularly the time derivatives. Section 3.5 gives an overview of the level set method which is used to store and move the interface position.
Interface velocity calculation :
The proper evaluation of the interface velocity is crucial in obtaining a precise and robust model. For this particular problem, the interface velocity is determined by theµ jump in heat flux at the interface, described in (5.2). The use of a Lagrange multiplier to impose the melting temperature allows the evaluation of the jump in heat flux directly from the Lagrange field q.
Flow over a cylinder
To validate the Stokes formulation, the flow over a cylinder was modelled. The domain is 1 m long and 0.5 m wide and the cylinder is centered at (0.4,0.25) m with a radius of 0.05 m. The material properties used are given in table 5.1. Two different boundary conditions were tested. In the first, referred to as case I (see figure 5.1), the top, bottom and right boundaries all have velocity boundary conditions, v1 = (−0.1, 0) m/s (Re =10). The left boundary is an open boundary (p0 = 0 Pa). In the second set of boundary conditions, referred to as case II (see figure 5.2), the top and bottom boundaries are walls (v0 = (0, 0) m/s), the left boundary has a pressure condition p1 = −5 Pa and the right wall is an open boundary (p0 = 0 Pa). The no-slip condition was applied using β = 108 and the convergence criteria for the Newton-Raphson algorithm is 10−6 .
Melting cylinder in a channel
The second problem includes the Stefan formulation to allow the cylinder to melt. The problem setup is as follows. A channel, l = 0.167 m in length and h = 0.025 m in height, contains a solid cylinder of radius 0.005 m. Both phases are initially at the melting temperature Tm = 273 K. The cylinder’s centre is initially at (l4,h2). At t = 0, a pressure difference ∆p = 4 Pa is applied between the channel’s inlet and outlet. The inlet temperature is 274 K. Both top and bottom edges are thermally insulated with a no-slip boundary condition. The pressure difference drives the fluid flow and the buoyancy force was removed from (5.3a). The material properties used are given in table 5.2 and a schematic representation of the problem in figure 5.7. The mesh includes 2904 quadrilateral elements.
Melting of pure tin :
A final problem, the melting of pure tin, was simulated based on the experimental and numerical data found in [81] and the phase change example model found in Comsol [23]. The same simulation was then run in Comsol using a moving mesh algorithm (ALE) and the solution was compared with the solution obtained using the purely XFEM approach. The problem setup is as follows. A square cavity, 0.10 m wide and 0.10 m high, is filled with liquid tin on the left and solid tin on the right. Both phases are initially at the melting temperature Tm = 505K. The initial interface is vertical at x = 0.02 m. At t = 0, the temperature of the left wall is increased to 508 K and the right wall decreased to 503 K, causing the metal to melt. Both top and bottom edges are insulated. The four boundaries are considered walls and a no-slip boundary is applied for the Stokes equations. The no slip condition on the interface is applied using a penalty term with β = 108 . The material properties used are given in table 5.3 and a schematic representation of the problem in figure 5.13. The presence of natural convection changes the heat flux within the melt by increasing the influx of heat near the top of the enclosure and reducing it near the bottom, resulting in an angled interface.
A coupled Stefan and Stokes formulation using the extended finite element method was developed for the resolution of phase change problems involving convection. The Lagrange multiplier technique developed for the diffusive case was successfully applied to the convective-diffusive problem. The temperature and velocity fields obtained using XFEM were compared to the moving mesh algorithm found in Comsol with good results. The XFEM formulation required less degrees of freedom and didn’t cause problems with distorted elements. The simple removal of degrees of freedom with a small contribution to the system for the Q2-Q1 Stokes formulation was shown to produce errors in the velocity field for problematic interface configurations. The same observation for a Q1-Q1 formulation was made in [49]. Future work will be done to include the complete Navier Stokes equations and the application of non-zero velocity boundary conditions on the interface to include density changes between solid and liquid phases.
Modelling of Phase Change with Non-Constant Density using XFEM and a Lagrange Multiplier :
Numerous extended finite element models for the solutions of the classical (diffusive) Stefan problem are found in the literature [14, 22, 46, 54]. More complex models involving convection with constant density have also been developed using different numerical techniques [16, 80, 83]. Particularly, a fully XFEM Stefan/NavierStokes model was used by [52]. Models including the density variation are more uncommon [82], mainly because assuming the density is constant generally has little impact on the interface position. In certain applications however, the conservation A straight-forward strategy to include the non-constant material densities is to use a moving-mesh algorithm such as the one found in the commercial code Comsol. This algorithm defines the phase change interface on a set of nodes, allowing the mass conservation boundary condition to be easily applied. However, the moving mesh adds considerable computational costs caused by the increase in degrees of freedom of the overall problem and the remeshing procedure required when the mesh becomes too distorted. These costs may hinder the use of a moving mesh algorithms in large scale multi-physical simulations which often have a large amount of degrees of freedom.
One Dimensional Phase Change Problem :
The first benchmark problem is inspired by the one dimensional two phase analytical solution of the Stefan problem in a semi-infinite domain (x > 0), taken from [54]. The thermal properties are constant except for the density and are given in table 6.1. These properties were chosen so that the convection term would affect the position of the interface with respect to the same problem without convection. The domain is 1m long and 0.25 m wide and the initial interface is at x = 0.515 m with the liquid phase on the right and solid phase on the left, as shown in figure 6.1. The initial temperature is Tm (see table 6.1). The top and bottom edges are insulated. At t = 0, the temperature on the left edge is lowered to 272 K and the right edge is increased to 275 K. For the Navier-Stokes equations, the right boundary is open (no stress) while for the top and bottom edges the following boundary condition is applied: v · n = 0. The time step is 0.05 sec, β =1 × 108 , Ttol =1 × 108 and the convergence criteria for the Newton-Raphson algorithm is 10−6 for both problem. The mesh contains 180 quadrilateral elements in XFEM and 196 in Comsol.
Conclusion :
As described in chapter 1, the main goal of this thesis was to develop a numerical model to predict the behaviour of the cryolite bath inside the electrolytic cell which could include the jump in density between the liquid and solid phases of the cryolite. The presence of multiple phases leads to various discontinuities in the solution which the model had to include. The model also had to scale easily and efficiently with an increase in degrees of freedom either through mesh size (or an increase in element density), the inclusion of additional physics (electromagnetism, chemistry) or the extension to three dimensions.
Table des matières
1 Introduction
1.1 Motivation
1.2 Thesis Contribution and Novelty
1.3 Thesis Organization
2 Literature Review
2.1 Introduction
2.2 Interface Handling Techniques
2.3 Stefan Problem
2.4 Navier-Stokes Equations
3 The Extended Finite Element Method
3.1 Introduction
3.2 Enrichment Strategies
3.3 Dirichlet Boundary Conditions
3.4 Numerical Integration
3.5 The Level Set Method
4 A XFEM Lagrange Multiplier Technique for Stefan Problems:
4.1 Resumé
4.2 Abstract
4.3 Introduction
4.4 Governing Equations
4.5 Numerical Implementation
4.6 Results
4.7 Conclusion
5 A XFEM Phase Change Model with Convection:
5.1 Resumé
5.2 Abstract
5.3 Introduction
5.4 Governing Equations
5.5 Numerical Implementation
5.6 Results
5.7 Conclusion
6 Modelling of Phase Change with Non-Constant Density using XFEM and a Lagrange Multiplier
6.1 Resumé
6.2 Abstract
6.3 Introduction
6.4 Governing Equations
6.5 Numerical Implementation
6.6 Results
6.7 Conclusion
7 Conclusion
