Heat equation (FE)

In this tutorial, the heat equation (first steady and then unsteady) is solved using finite-elements.

Theory

This example shows how to solve the heat equation in steady and unsteady formulations. The unsteady heat equation is given by:

\[ \rho C_p \frac{\partial u}{\partial t} - \nabla \cdot ( \lambda \nabla u) = f\]

We shall assume that $f, \, \rho, \, C_p, \, \lambda \, \in L^2(\Omega)$. The weak form of the problem is given by: find $ u \in \tilde{H}^1_0(\Omega)$ (there will be at least one Dirichlet boundary condition) such that:

\[ \forall v \in \tilde{H}^1_0(\Omega), \, \, \, \underbrace{\int_\Omega \rho C_p \frac{\partial u}{\partial t} v dx}_{m(\partial_t u,v)} + \underbrace{\int_\Omega \lambda \nabla u \cdot \nabla v dx}_{a(u,v)} = \underbrace{\int_\Omega f v dx}_{l(v)}\]

To numerically solve this problem we seek an approximate solution using Lagrange $P^2$ elements.

As usual, start by importing the necessary packages.

using Bcube
using BcubeGmsh
using BcubeVTK
using LinearAlgebra

First we define some physical and numerical constants

const q = 1500.0 # heat source
const λ = 100.0 # thermal conductivity
const ρCp = 100.0 * 200.0 # density times specific heat capacity
const degree = 2 # Degree of the Lagrange polynomials (for $$P^2$$ elements)
const outputpath = joinpath(@__DIR__, "..", "..", "myout", "heat_equation/")

Steady case

We will first start by solving the steady case: find $u \in \tilde{H}^1_0(\Omega)$ such that $\forall v \in \tilde{H}^1_0(\Omega), \, \, \, a(u,v) = l(v)$.

println(" ----- Solving steady problem -----")

Read 2D mesh

mesh_path = joinpath(@__DIR__, "..", "..", "input", "mesh", "domainSquare_tri.msh")
mesh = read_mesh(mesh_path)

Build function space and associated Trial and Test FE spaces. We impose a Dirichlet condition with a temperature of 260K on boundary "West"

fs = FunctionSpace(:Lagrange, degree)
U = TrialFESpace(fs, mesh, Dict("West" => 260.0))
V = TestFESpace(U)

Define measures for cell integration

dΩ = Measure(CellDomain(mesh), 2 * degree + 1)

Define bilinear and linear forms. The steady problem is defined by the bilinear form a and the linear form l.

a(u, v) = ∫(λ * ∇(u) ⋅ ∇(v))dΩ
l(v) = ∫(q * v)dΩ

Create an affine FE system and solve it using the AffineFESystem structure. By default, an LU decomposition is used to solve the system. In the present case we know that the matrix A is symmetric. Hence a Cholesky decomposition can be used. The result is a FEFunction (Tn). We can extract its dof values: the result is named Tn_dofs.

Cholesky_linsolve!(y, A, x) = y .= cholesky(Symmetric(A)) \ x
sys = AffineFESystem(a, l, U, V, Cholesky_linsolve!)
Tn = Bcube.solve(sys)
Tn_dofs = get_dof_values(Tn)

Compute analytical solution for comparison. An FEFunction is built using the analytical solution and the dof values are extracted.

T_analytical = PhysicalFunction(x -> 260.0 + (q / λ) * x[1] * (1.0 - 0.5 * x[1]))
Ta = FEFunction(U, mesh, T_analytical)
Ta_dofs = get_dof_values(Ta)

Write both the obtained FE solution and the analytical solution to a vtk file.

mkpath(outputpath)
dict_vars = Dict("Temperature (numerical)" => Tn, "Temperature (analytical)" => Ta)
write_file(outputpath * "result_steady_heat_equation.pvd", mesh, dict_vars)

Compute and display the error, which is computed using the previously extracted dof values.

@show norm(Tn_dofs .- Ta_dofs, Inf) / norm(Ta_dofs, Inf)

Unsteady case

We now solve the unsteady problem: find $ u $ such that for all $v , \, \, \, m(\frac{\partial u}{\partial t}, v) + a(u,v) = l(v)$. The code for the unsteady case if of course very similar to the steady case, at least for the beginning.

println(" ----- Solving unsteady problem -----")

Start by defining two additional parameters:

totalTime = 100.0
Δt = 0.1

Read a slightly different mesh

mesh_path = joinpath(@__DIR__, "..", "..", "input", "mesh", "domainSquare_tri_2.msh")
mesh = read_mesh(mesh_path)

The rest is similar to the steady case

fs = FunctionSpace(:Lagrange, degree)
U = TrialFESpace(fs, mesh, Dict("West" => 260.0))
V = TestFESpace(U)
dΩ = Measure(CellDomain(mesh), 2 * degree + 1)

Compute matrices associated to bilinear and linear forms, and assemble. The bilinear form a and linear form l have already been defined in the steady problem. We now define the bilinear form m.

m(u, v) = ∫(ρCp * u ⋅ v)dΩ

A = assemble_bilinear(a, U, V)
M = assemble_bilinear(m, U, V)
L = assemble_linear(l, V)

Compute a vector of dofs whose values are zeros everywhere except on dofs lying on a Dirichlet boundary, where they take the Dirichlet value

Ud = assemble_dirichlet_vector(U, V, mesh)

Apply lift

L = L - A * Ud

Apply homogeneous dirichlet condition

apply_homogeneous_dirichlet_to_vector!(L, U, V, mesh)
apply_dirichlet_to_matrix!((A, M), U, V, mesh)

Form time iteration matrix (note that this is bad for performance since up to now, M and A are sparse matrices)

Miter = factorize(M + Δt * A)

Init the solution with a constant temperature of 260K

ϕ = FEFunction(U, 260.0)

Write initial solution to a file

mkpath(outputpath)
dict_vars = Dict("Temperature" => ϕ)
write_file(outputpath * "result_unsteady_heat_equation.pvd", mesh, dict_vars, 0, 0.0)

Time loop

itime = 0
t = 0.0
while t <= totalTime
    global t, itime
    t += Δt
    itime = itime + 1
        #! format: off
        @show t, itime
        #! format: on

        # Compute rhs
    rhs = Δt * L + M * (get_dof_values(ϕ) .- Ud)

    # Invert system and apply inverse shift
    set_dof_values!(ϕ, Miter \ rhs .+ Ud)

    # Write solution (every 10 iterations)
    if itime % 10 == 0
        dict_vars = Dict("Temperature" => ϕ)
        write_file(
            outputpath * "result_unsteady_heat_equation.pvd",
            mesh,
            dict_vars,
            itime,
            t;
            collection_append = true,
        )
    end
end

This page was generated using Literate.jl.