First verification, heat equation

Author: Jukka Aho

Abstract: ?

Problem description

Theory / solution

  • Reference, uncertainty

Implementation in JuliaFEM

  • Mesh, boundary conditions, etc.

In [3]:
using JuliaFEM: Seg2, Quad4, FieldSet, Field, PlaneHeatProblem, DirichletProblem, SimpleSolver, interpolate

# Define Problem 1:
# - Field function: Laplace equation Δu=0 in Ω={u∈R²|(x,y)∈[0,1]×[0,1]}
# - Neumann boundary on Γ₁={0<=x<=1, y=0}, ∂u/∂n=600 on Γ₁

el1 = Quad4([1, 2, 3, 4])
push!(el1, FieldSet("geometry", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])]))
push!(el1, FieldSet("temperature thermal conductivity", [Field(0.0, 6.0)]))

el2 = Seg2([1, 2])
push!(el2, FieldSet("geometry", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0]])]))

# Boundary load, linear ramp 0 -> 600 at time 0 -> 1
load = FieldSet("temperature flux")
push!(load, Field(0.0, 0.0))
push!(load, Field(1.0, 600.0))
push!(el2, load)

problem1 = PlaneHeatProblem()
push!(problem1, el1)
push!(problem1, el2)

# Define Problem 2:
# - Dirichlet boundary Γ₂={0<=x<=1, y=1}, u=0 on Γ₂

el3 = Seg2([3, 4])
push!(el3, FieldSet("geometry", [Field(0.0, Vector[[1.0, 1.0], [0.0, 1.0]])]))

problem2 = DirichletProblem()
push!(problem2, el3)

# Create a solver for a set of problems
solver = SimpleSolver()
push!(solver, problem1)
push!(solver, problem2)

# Solve problem at time t=1.0 and update fields
call(solver, 1.0)

# Postprocess.
# Interpolate temperature field along boundary of Γ₁ at time t=1.0
xi = linspace([-1.0], [1.0], 5)
X = interpolate(el2, "geometry", xi, 1.0)
T = interpolate(el2, "temperature", xi, 1.0)
println(X)
println(T)


20-loka 15:50:52:DEBUG:root:total dofs: 4
Residual norm: 9.845568954283847e-14
Array{T,1}[[0.0,0.0],[0.25,0.0],[0.5,0.0],[0.75,0.0],[1.0,0.0]]
[100.00000000000003,100.00000000000003,100.00000000000003,100.00000000000003,100.00000000000003]

Comparison of results

We should have T=100 in edge.


In [4]:
using FactCheck
@fact mean(T) --> roughly(100.0)


Out[4]:
Success :: (line:-1) :: fact was true
  Expression: mean(T) --> roughly(100.0)
    Expected: 100.0
    Occurred: 100.00000000000003