Fourier Conduction#

This examples shows how OpenPNM can be used to simulate thermal conduction on a generic grid of nodes. The result obtained from OpenPNM is compared to the analytical result.

As usual, start by importing OpenPNM, and the SciPy library.

[1]:
%matplotlib inline
import numpy as np
import scipy as sp
import openpnm as op
%config InlineBackend.figure_formats = ['svg']
np.random.seed(10)
ws = op.Workspace()
ws.settings["loglevel"] = 40
np.set_printoptions(precision=5)

Generating the Network object#

Next, 2D a Network is generated with dimensions of 10x50 elements. The lattice spacing is given by Lc. Boundaries are added all around the edges of Network object using the add_boundariy_pores method.

[2]:
divs = [10, 50]
Lc = 0.1  # cm
pn = op.network.Cubic(shape=divs, spacing=Lc)
pn.add_boundary_pores(['left', 'right', 'front', 'back'])

Creating a Phase object#

All simulations require a phase object which possess the thermosphysical properties of the system. In this case, we’ll create a generic phase object, call it copper, though it has no properties; we’ll add these by hand later.

[3]:
# Create Phase object and associate with a Physics object
Cu = op.phases.GenericPhase(network=pn)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [3], in <cell line: 2>()
      1 # Create Phase object and associate with a Physics object
----> 2 Cu = op.phases.GenericPhase(network=pn)

AttributeError: module 'openpnm' has no attribute 'phases'

Assigning Thermal Conductance to Copper#

In a proper OpenPNM model we would create a Geometry object to manage all the geometrical properties, and a Physics object to calculate the thermal conductance based on the geometric information and the thermophysical properties of copper. In the present case, however, we’ll just calculate the conductance manually and assign it to Cu.

[4]:
# Add a unit conductance to all connections
Cu['throat.thermal_conductance'] = 1
# Overwrite boundary conductances since those connections are half as long
Ps = pn.pores('*boundary')
Ts = pn.find_neighbor_throats(pores=Ps)
Cu['throat.thermal_conductance'][Ts] = 2
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [4], in <cell line: 2>()
      1 # Add a unit conductance to all connections
----> 2 Cu['throat.thermal_conductance'] = 1
      3 # Overwrite boundary conductances since those connections are half as long
      4 Ps = pn.pores('*boundary')

NameError: name 'Cu' is not defined

Generating the algorithm objects and running the simulation#

The last step in the OpenPNM simulation involves the generation of a Algorithm object and running the simulation.

[5]:
# Setup Algorithm object
alg = op.algorithms.FourierConduction(network=pn, phase=Cu)
inlets = pn.pores('right_boundary')
outlets = pn.pores(['front_boundary', 'back_boundary', 'right_boundary'])
T_in = 30*np.sin(np.pi*pn['pore.coords'][inlets, 1]/5)+50
alg.set_value_BC(values=T_in, pores=inlets)
alg.set_value_BC(values=50, pores=outlets)
alg.run()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [5], in <cell line: 2>()
      1 # Setup Algorithm object
----> 2 alg = op.algorithms.FourierConduction(network=pn, phase=Cu)
      3 inlets = pn.pores('right_boundary')
      4 outlets = pn.pores(['front_boundary', 'back_boundary', 'right_boundary'])

NameError: name 'Cu' is not defined

This is the last step usually required in a OpenPNM simulation. The algorithm was run, and now the simulation data obtained can be analyzed. For illustrative purposes, the results obtained using OpenPNM shall be compared to an analytical solution of the problem in the following.

First let’s rehape the ‘pore.temperature’ array into the shape of the network while also extracting only the internal pores to avoid showing the boundaries.

[6]:
import matplotlib.pyplot as plt
sim = alg['pore.temperature'][pn.pores('internal')]
temp_map = np.reshape(a=sim, newshape=divs)
plt.subplots(1, 1, figsize=(10, 5))
plt.imshow(temp_map, cmap=plt.cm.plasma);
plt.colorbar();
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [6], in <cell line: 2>()
      1 import matplotlib.pyplot as plt
----> 2 sim = alg['pore.temperature'][pn.pores('internal')]
      3 temp_map = np.reshape(a=sim, newshape=divs)
      4 plt.subplots(1, 1, figsize=(10, 5))

NameError: name 'alg' is not defined

Also, let’s take a look at the average temperature:

[7]:
print(f"T_average (numerical): {alg['pore.temperature'][pn.pores('internal')].mean():.5f}")
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [7], in <cell line: 1>()
----> 1 print(f"T_average (numerical): {alg['pore.temperature'][pn.pores('internal')].mean():.5f}")

NameError: name 'alg' is not defined

The analytical solution is computed as well, and the result is the same shape as the network (including the boundary pores).

[8]:
# Calculate analytical solution over the same domain spacing
X = pn['pore.coords'][:, 0]
Y = pn['pore.coords'][:, 1]
soln = 30*np.sinh(np.pi*X/5)/np.sinh(np.pi/5)*np.sin(np.pi*Y/5) + 50
soln = soln[pn.pores('internal')]
soln = np.reshape(soln, (divs[0], divs[1]))
[9]:
plt.subplots(1, 1, figsize=(10, 5))
plt.imshow(soln, cmap=plt.cm.plasma);
plt.colorbar();
../../../_images/examples_simulations_steady_state_continuum_heat_transfer_18_0.svg

Also, let’s take a look at the average temperature:

[10]:
print(f"T_average (analytical): {soln.mean():.5f}")
T_average (analytical): 59.24706

Both the analytical solution and OpenPNM simulation can be subtracted from each other to yield the difference in both values.

[11]:
diff = soln - temp_map
plt.subplots(1, 1, figsize=(10, 5))
plt.imshow(diff, cmap=plt.cm.plasma);
plt.colorbar();
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [11], in <cell line: 1>()
----> 1 diff = soln - temp_map
      2 plt.subplots(1, 1, figsize=(10, 5))
      3 plt.imshow(diff, cmap=plt.cm.plasma);

NameError: name 'temp_map' is not defined
[12]:
print(f"Minimum error: {diff.min():.5f}, maximum error: {diff.max():.5f}")
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [12], in <cell line: 1>()
----> 1 print(f"Minimum error: {diff.min():.5f}, maximum error: {diff.max():.5f}")

NameError: name 'diff' is not defined

The maximum error is 0.01 degrees on a 50 degree profile, which is quite good and thus demonstrates that the OpenPNM finite difference approach is versatile despite being simple.