I am currently solving a 3D TEG model, 1) I wanted to include resistance as shown below in the model:
currently, am only applying voltage and temperature boundary conditions at the required faces.
Can anyone suggest to me how to map the surface nodes to a single master that, so that I can create the terminal node to connect the external load resistor.
My current fenics code is:
import dolfin as d
import matplotlib.pyplot as plt
import numpy as np
import timeit
import PIL, os
from PIL import Image
#time start
start = timeit.default_timer()
################################## mesh part
# load the mesh data from Gmsh
mesh_file = 'Mesh1.h5'
mesh = d.Mesh()
hdf = d.HDF5File(mesh.mpi_comm(), mesh_file, 'r')
hdf.read(mesh, '/mesh', False)
cells = d.MeshFunction('size_t', mesh, dim=3)
hdf.read(cells, '/subdomains')
facets = d.MeshFunction('size_t', mesh, dim=2)
hdf.read(facets, '/facets')
# parameters for thermal simulations
Kcu = d.Constant(401) # Thermal conductivity [W/m/K]
KBi2Te3 = d.Constant(1.46) # Thermal conductivity [W/m/K]
T_top = d.Constant(452) # Temperature [°C]
T_bot = d.Constant(22) # Temperature [°C]
r = d.Constant(0)
################################### FE Thermal Part ###################################
# function space
W = d.FunctionSpace(mesh, 'CG', 2)
print("Number of DOFs: {}".format(W.dim()))
dx = d.Measure('dx', domain=mesh, subdomain_data=cells)
# JZ: use the Measure function
ds = d.Measure('ds', domain=mesh, subdomain_data=facets)
# Trial and Test functions
T = d.TrialFunction(W)
v = d.TestFunction(W)
# boundary conditions for temperature
# dirichlet condition
T0 = d.DirichletBC(W, T_top, facets, 3)
T1 = d.DirichletBC(W, T_bot, facets, 1)
T2 = d.DirichletBC(W, T_bot, facets, 2)
T_bc = [T0, T1, T2]
# JZ: reformulate for steady state thermal analysis
a = KBi2Te3 * d.dot(d.grad(T), d.grad(v)) * dx(7) + KBi2Te3 * d.dot(d.grad(T), d.grad(v)) * dx(6) + Kcu * d.dot(d.grad(T), d.grad(v)) * dx(8) + Kcu * d.dot(d.grad(T), d.grad(v)) * dx(9) + Kcu * d.dot(d.grad(T), d.grad(v)) * dx(10)
L = r * v * dx
T = d.Function(W) # solution for T
Tproblem = d.LinearVariationalProblem(a, L, T, T_bc)
Tsolver = d.LinearVariationalSolver(Tproblem)
Tsolver.parameters["linear_solver"] = "mumps"
# solve for temperature
Tsolver.solve()
# parameters for electrical simulations
sigmacu = d.Constant(59980806.14) # Electric conductivity [S/m]
sigmaBi2Te3 = d.Constant(60975.60976) # Electric conductivity [S/m]
S_P = d.Constant(0.000187) # Seebeck Coefficient [V/K]
S_N = d.Constant(-0.000187) # Seebeck Coefficient [V/K]
S_cu = d.Constant(0.0000065) # Seebeck Coefficient [V/K]
V_h = d.Constant(0.08) # Reference voltage [V]
V_l = d.Constant(0) # Reference voltage [V]
################################### FE Electrical Part ###################################
# function space
W1 = d.FunctionSpace(mesh, 'CG', 2)
# Trial and Test functions
V = d.TrialFunction(W1)
vtest = d.TestFunction(W1)
# boundary conditions for voltage
# dirichlet condition
#Vh = d.DirichletBC(W1, 0, facets, 4)
Vl = d.DirichletBC(W1, V_l, facets, 5)
V_bc = [Vl]
# JZ: reformulate for steady state electric model
a = sigmaBi2Te3 * d.dot(d.grad(V), d.grad(vtest)) * dx(7) + sigmaBi2Te3 * d.dot(d.grad(V), d.grad(vtest)) * dx(6) + sigmacu * d.dot(d.grad(V), d.grad(vtest)) * dx(8) + sigmacu * d.dot(d.grad(V), d.grad(vtest)) * dx(9) + sigmacu * d.dot(d.grad(V), d.grad(vtest)) * dx(10)
L = -(sigmaBi2Te3 * S_P * d.dot(d.grad(T), d.grad(vtest)) * dx(7)) - sigmaBi2Te3 * S_N * d.dot(d.grad(T), d.grad(vtest)) * dx(6) - sigmacu * S_cu * d.dot(d.grad(T), d.grad(vtest)) * dx(8) - sigmacu * S_cu * d.dot(d.grad(T), d.grad(vtest)) * dx(9) - sigmacu * S_cu * d.dot(d.grad(T), d.grad(vtest)) * dx(10)
V = d.Function(W1)
# solution for V
Vproblem = d.LinearVariationalProblem(a, L, V, V_bc)
Vsolver = d.LinearVariationalSolver(Vproblem)
Vsolver.parameters["linear_solver"] = "mumps"
# solve for temperature
Vsolver.solve()
# Electric Field and current density
W2 = d.VectorFunctionSpace(mesh, 'Lagrange', 2)
E = d.project(-d.grad(V), W2, solver_type="cg", preconditioner_type='hypre_amg')
# Current density
J = d.TrialFunction(W2)
j = d.TestFunction(W2)
a = d.dot(J, j) * dx
L = d.dot(-sigmaBi2Te3 * d.grad(V), j) * dx(7) + d.dot(-sigmaBi2Te3 * d.grad(V), j) * dx(6) + d.dot(-sigmacu*d.grad(V), j) * dx(8) + d.dot(-sigmacu * d.grad(V), j)* dx(9) + d.dot(-sigmacu*d.grad(V), j) * dx(10) + d.dot(-sigmaBi2Te3 * S_P * d.grad(T), j) * dx(6) + d.dot(- sigmaBi2Te3 * S_N * d.grad(T), j) * dx(7) + d.dot(-sigmacu *S_cu* d.grad(T), j) * dx(8) + d.dot(-sigmacu *S_cu* d.grad(T), j)* dx(9) + d.dot(-sigmacu *S_cu* d.grad(T), j) * dx(10)
J = d.Function(W2)
# solution for J
Jproblem = d.LinearVariationalProblem(a, L, J)
Jsolver = d.LinearVariationalSolver(Jproblem)
Jsolver.parameters["linear_solver"] = "mumps"
# solve for current density
Jsolver.solve()
#n = d.FacetNormal(mesh)
n = d.Constant((1.0,0.0,0.0))
currentd = -d.dot(J, n) * ds(5)
t_current = d.assemble(currentd)
print("Current:", t_current)
######################### Post - Processing
#time stop
stop = timeit.default_timer()
# Save solution to file in VTK format
vtkfile = d.File('TEG/Temp.pvd')
vtkfile << T
vtkfile = d.File('TEG/Volt.pvd')
vtkfile << V
#vtkfile = d.File('TEG/Electricfiled.pvd')
#vtkfile << E
vtkfile = d.File('TEG/Currentdensity.pvd')
vtkfile << J