Hi again,
I created 2d .xdmf files with your help and now am trying to make them work. This is strange, but FESTIM solves temperature field ok, but then fails to correctly solve for solute and retention. Looks like these fields are 0 everywhere except the surface, it gives almost immediate (tenths of second saturation) saturation which probably corresponds to equilibration time between incoming flux and the surface. I am definitely missing something important. Let me try to enter the Python file here:
# Imports:
import numpy as np
import meshio
import gmsh
import festim as F
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.constants import Avogadro
from dolfin import *
from fenics import RectangleMesh, MeshFunction, Point, CompiledSubDomain, plot, IntervalMesh
from mshr import *
# Materials and parameters:
flux_surface_id=32
coolant_surface_id=36
breeder_surface_id=33
tungsten_id = 31
cu_id = 30
cucrzr_id = 29
w_density = 6.3e28
cu_density = 9.2e28
cucrzr_density = 9.2e28
damage_percent = 0.01
#default material definitions, taken from Nucl. Fusion 61 (2021) 036038
tungsten = F.Material(
id=31, D_0=1.5e-7, E_D=0.39, S_0=1.87e24, E_S=1.04,
thermal_cond=173)
cu = F.Material(
id=30, D_0=6.6e-7, E_D=0.39, S_0=3.14e24, E_S=0.57,
thermal_cond=350
)
cucrzr = F.Material(
id=29, D_0=3.9e-7, E_D=0.42, S_0=4.28e23, E_S=0.39,
thermal_cond=320
)
default_materials = [tungsten, cu, cucrzr]
#default traps, taken from Nucl. Fusion 61 (2021) 036038
trap_1 = F.Trap( #a grouped trap containing one in w and one in cu and cucrzr
k_0 = [3.8e-17,6e-17,1.2e-16],
E_k= [0.39,0.39,0.42],
p_0=[8.4e12,8e13,8e13],
E_p=[1.2,0.5,0.50],
density = [w_density*0.0005,cu_density*0.00005,cucrzr_density*0.00005],
materials = [tungsten, cu, cucrzr],
id= 1
)
trap_2 = F.Trap( #a grouped trap containing one in w and one in cucrzr
k_0 = [3.8e-17,1.2e-16],
E_k= [0.39,0.42],
p_0=[8.4e12,8e13],
E_p=[1.4,0.83],
density = [w_density*0.005,cucrzr_density*0.04],
materials = [tungsten, cucrzr],
id = 2
)
damage_trap = F.Trap(k_0=1.5e-7/(1.1e-10**2 * 6 * w_density),
E_k=0.39,
p_0=8.4e12,
E_p=1.51,
density = damage_percent*w_density,
materials = tungsten,
id=3)
default_traps = [trap_1, trap_2, damage_trap]
# 2D simulation:
def run_simulation_2d(abs_tolerance, rel_tolerance, dt_start, folder,
traps = default_traps,
materials = default_materials,
flux = 1e20,
flux_surface = 32,
coolant_surface = 36,
breeder_surface = 33,
other_surfaces = [33,34,35],
mean_depth = 1e-8,
temperatures = [273, 173],
simulation_time = 1e-1,
export_xdmf = False,
mumps_solver = False): #The MUMPS solver is better for finer meshes
"""
Runs 2D simulations
inputs:
abs_tolerance (int) : absolute tolerace for the solver. Values of around 1e8-1e12 tend to be good. Higher is less accurate but compiles better
rel_tolerance (int) : relative tolerance for the solver. Values of 1e-5-1e-12 tend to be good. Larger is less accurate but compiles better
dt_start (float): initial step size for the simulation. Larger/longer simulations can sometimes work better with larger values (up to ~1e5!)
folder (string) : folder where you want the data saved
traps (list) : list of trap objects
materials (list) : list of materials
flux (int) : flux of incoming hydrogen
flux_surface (int) : surface id for the incoming flux
coolant_surface (int) : surface id beside the coolant
other_surfaces (list) : ids for the other surfaces that require BCs (which ones you'd like to hold at 0 flux)
mean_depth (float) : hydrogen implantation depth
temperatures (list) : two values, first is plasma-side temperature, second is coolant-side
simulation_time (float) : total simulation runtime
export_xdmf (bool) : do you want to export an XDMF file alongside the CSV?
mumps_solver (bool) : option to use an alternative solver to the default Newton Solver. Sometimes works better for larger meshes
"""
model = F.Simulation()
model.mesh=F.MeshFromXDMF(volume_file="volume_mesh.xdmf", boundary_file="surface_mesh.xdmf")
model.materials = materials
model.traps = traps
model.T = F.HeatTransferProblem(transient = False)
bcs = [F.ImplantationDirichlet(surfaces=flux_surface,
phi=flux,
R_p=mean_depth,
D_0=1.5e-07,
E_D=0.39),
#F.DirichletBC(
# surfaces = other_surfaces,
# value = 0,
# field = 0),
F.DirichletBC(
surfaces = coolant_surface,
value = 0,
field = 0),
F.DirichletBC(
surfaces = flux_surface,
value = temperatures[0],
field = 'T'),
F.DirichletBC(
surfaces = coolant_surface,
value = temperatures[1],
field = 'T'
)]
model.boundary_conditions=bcs
model.dt = F.Stepsize(initial_value = dt_start, stepsize_change_ratio=1.1, dt_min = 0.0005)
model.settings = F.Settings(
absolute_tolerance=abs_tolerance,
relative_tolerance=rel_tolerance,
final_time = simulation_time,
chemical_pot=True,
traps_element_type="DG",
)
if mumps_solver:
model.settings.linear_solver = "mumps"
retention_1 = F.TotalVolume("retention", volume = 31)
retention_2 = F.TotalVolume("retention", volume = 30)
retention_3 = F.TotalVolume("retention", volume = 29)
solute_1 = F.TotalVolume("solute", volume = 31)
solute_2 = F.TotalVolume("solute", volume = 30)
solute_3 = F.TotalVolume("solute", volume = 29)
surface_flux_coolant = F.SurfaceFlux(field="solute", surface = 36)
surface_flux_breeder = F.SurfaceFlux(field="solute", surface = 33)
derived = [retention_1, retention_2, retention_3, solute_1, solute_2, solute_3]
model.exports = F.Exports([
F.XDMFExport("solute", filename=folder+"/solute.xdmf"),
F.XDMFExport("retention", filename=folder+"/retention.xdmf"),
F.XDMFExport("T", filename=folder+"/temp.xdmf"),
F.DerivedQuantities(derived, filename = folder + "/derived_qs.csv"),
])
model.initialise()
model.run()
return model
# Analysis:
def retention_analysis(folder, dimension, section_area,
tungsten_thickness = 3e-3,
cucrzr_thickness = 11e-3,
coolant_radius = 6e-3,
z_extent = 5e-3,
plot_grams = False,
plot_days = True):
data = np.genfromtxt(fname = folder + "/derived_qs.csv", dtype = 'float', delimiter = ',', skip_header = True)
retention = (np.array(data[:,1]) + np.array(data[:,2]))
if dimension == 1:
area_factor = section_area
elif dimension == 2:
area_factor = section_area/(tungsten_thickness*2+cucrzr_thickness*2+coolant_radius*2)
elif dimension == 3:
area_factor = section_area/((tungsten_thickness+cucrzr_thickness+coolant_radius)*z_extent)
integrated_retention = retention * area_factor / 2 / Avogadro * 3
if plot_days:
time = data[:,0]/3600/24
if plot_grams:
plotted_retention = integrated_retention
plt.plot(time, plotted_retention)
plt.xlabel("Time (days)")
plt.ylabel("Tritium Retention (g)")
plt.show()
else:
plotted_retention = retention
plt.plot(time, retention)
plt.xlabel("Time (days)")
if dimension == 1:
plt.ylabel("Tritium Retention (H/m2)")
elif dimension == 2:
plt.ylabel("Tritium Retention (H/m)")
else:
plt.ylabel("Tritium Retention (H)")
plt.show()
return plotted_retention, time
# Run simulation:
two_d_model = run_simulation_2d(1e16, 1e-8, 1e-6, "2d_monoblock", export_xdmf=True, mumps_solver=True)
# Plot results:
results_2d = retention_analysis("2d_monoblock", 2, 233.93)
