Convection solve

1

I solve diffusion-convection equation, for problem:

Cube
Cube1058×824 44.9 KB

Solve done, but i have T>1, my field in paraview:
image
image1920×1080 271 KB

Max T = 1.4, but i expect to receive T from 0 to 1 only.
I think that i did’t right chose my FunctionSpace. My code:
I read mesh from gmsh and:

T = FunctionSpace(mesh, 'P', 3)
V = VectorFunctionSpace(mesh, 'P', 2)

bc_inlet = DirichletBC(T, Constant(1), facet_domains, 7)
bct = [bc_inlet]

u = Function(V)
u = project(Constant((1,1)), V)
plot(u)

t_n = Function(T)
t_n = project(Constant(0), T)

t = TrialFunction(T)
ti = TestFunction(T)

Time = 10.0
num_times = 1000
dt = Time/num_times
k = Constant(dt)

D = Constant(0)
F4 = dot((t - t_n) / k, ti)*dx \
    + div( t * u)*ti*dx \
    - D * inner( grad(t), grad(ti))*dx

a4 = lhs(F4)
L4 = rhs(F4)

A4 = assemble(a4)
[bc.apply(A4) for bc in bct]

import numpy as np

time_s = 0
t_ = Function(T)
#num_steps
for n in range(num_times):
    time_s += dt
    
    b4 = assemble(L4)
    [bc.apply(b4) for bc in bct]
    solve(A4, t_.vector(), b4)
    
    t_n.assign(t_)
    max_t = np.abs(np.array(t_n.vector())).max()
    print('t = %.2f: max_f = %.10f' % (time_s, max_t))

What type of elements should i choose?

2

Please have a look at the advection diffusion demo: Bitbucket

2 Likes