Heat Equation with Varying Kappa

#! /usr/bin/env python3

from fenics import *

def heat_steady_01 ( ):

#*****************************************************************************80

heat_steady_01, 2D steady heat equation on a rectangle, constant K.

Discussion:

Del K Del U = 0 in Omega

U = 10 on dOmegaTop

U = 100 on dOmega - dOmegaTop

Omega = rectangle with corners (0.0, 0.0 ) and ( 5.0, 1.0 ).

K = 1 (thermal diffusivity)

Licensing:

This code is distributed under the GNU LGPL license.

Modified:

19 October 2018

Author:

John Burkardt

import matplotlib.pyplot as plt

MESH:

The region is a rectangle with corners (0.0, 0.0 ) and ( 5.0, 1.0 ).

The mesh uses 50 horizontal and 10 vertical divisions.

sw = Point ( 0.0, 0.0 )
ne = Point ( 5.0, 1.0 )
mesh = RectangleMesh ( sw, ne, 50, 10 )

FUNCTION SPACE:

piecewise linear Lagrange functions.

V = FunctionSpace ( mesh, “Lagrange”, 1 )

BOUNDARY CONDITIONS:

BC_TOP: Along the top, we impose U=10.

BC_SIDE: Along the bottom and sides, U=100.

Combine these conditions as “BC”.

y_top = 1.0
u_top = 10.0
def on_top ( x, on_boundary ):
return ( y_top - DOLFIN_EPS <= x[1] )
bc_top = DirichletBC ( V, u_top, on_top )

u_side = 100.0
def on_side ( x, on_boundary ):
return ( on_boundary and x[1] < y_top - DOLFIN_EPS )
bc_side = DirichletBC ( V, u_side, on_side )

bc = [ bc_top, bc_side ]

TRIAL and TEST FUNCTIONS:

u = TrialFunction ( V )
v = TestFunction ( V )

BILINEAR and LINEAR FORMS:

k = 1

def k ( x ):
value = conditional ( le ( x[0], 0.75 ), 1.0, 4.0 )
return value

k = Expression ( ‘1 + x[0]’, degree = 0 )

k = Expression ( ‘1’, degree = 0 )

Auv = k * inner ( grad ( u ), grad ( v ) ) * dx

f = Constant ( 0.0 )
Lv = f * v * dx

SOLVE:

Solve the variational problem a(u,v)=l(v) with boundary conditions.

w = Function ( V )
solve ( Auv == Lv, w, bc )

Plot the solution W.

plot ( w, title = ‘heat_steady_01’ )
plu=plot(w)
plt.colorbar(plu)
filename = ‘heat_steady_01.png’
plt.savefig ( filename )
print ( “Saving graphics in file ‘%s’” % ( filename ) )
plt.close ( )

Terminate.

return

def heat_steady_01_test ( ):

#*****************************************************************************80

heat_steady_01_test tests heat_steady_01.

Licensing:

This code is distributed under the GNU LGPL license.

Modified:

21 October 2018

Author:

John Burkardt

import time

print ( time.ctime ( time.time() ) )

Report level = only warnings or higher.

level = 30
set_log_level ( level )

print ( ‘’ )
print ( ‘heat_steady_01_test:’ )
print ( ’ FENICS/Python version’ )
print ( ’ 2D steady heat equation in a rectangle.’ )

heat_steady_01 ( )

Terminate.

print ( ‘’ )
print ( ‘heat_steady_01_test:’ )
print ( ’ Normal end of execution.’ )
print ( ‘’ )
print ( time.ctime ( time.time() ) )
return

if ( name == ‘main’ ):

heat_steady_01_test ( )

Please format your post with 3x` encapsulation, so that your post becomes readable. Please also add a question, and a potential error message if you obtain one. Follow the guidelines in: Read before posting: How do I get my question answered?