Zero flux in plasma implantion and TDS

User Support
1

Hi,I’m trying to fit a irradiation of tungsten sample by deuterium plasma and subsequent TDS.But I found that the flux is 0 everywhere.I would like to know why this situation occurred.

Here’s my code:

import festim as F
my_model = F.Simulation()

import numpy as np

vertices = np.concatenate([
    np.linspace(0, 3e-3, num=300),
])

my_model.mesh = F.MeshFromVertices(vertices)

tungsten = F.Material(
    id=1,
    D_0=2.9e-7,  # m2/s
    E_D=0.29,  # eV
)

my_model.materials = tungsten

import sympy as sp

implantation_time = 9000  # s

ion_flux = sp.Piecewise((2.3e22, F.t <= 1800),
                        (1.9e22, (F.t>1800) & (F.t <= 3600)),
                        (1.03e22, (F.t>3600) & (F.t <= 5400)),
                        (6.23e21, (F.t>5400) & (F.t <= 7200)),
                        (2.5e21, (F.t>7200) & (F.t <= 9000)),
                        (0, True))

source_term = F.ImplantationFlux(
    flux=ion_flux,  # H/m2/s
    imp_depth=4e-9,  # m
    width=2e-9,  # m
    volume=1
)

my_model.sources = [source_term]

w_atom_density = 6.3e28  # atom/m3


n1=5e-3
n2=8e-3
n3=5e-2

trap_1 = F.Trap(
        k_0=2.9e12/w_atom_density,
        E_k=0.29,
        p_0=1e13,
        E_p=0.95,
        density=sp.Piecewise((n1*w_atom_density, F.x <= 10e-6),
                             (0,True)),
        materials=tungsten
    )
trap_2 = F.Trap(
        k_0=2.9e12/w_atom_density,
        E_k=0.29,
        p_0=1e13,
        E_p=1.3,
        density=n2*w_atom_density,
        materials=tungsten
    )

trap_3 = F.Trap(
        k_0=2.9e12/w_atom_density,
        E_k=0.29,
        p_0=1e13,
        E_p=1.6,
        density=sp.Piecewise((n3*w_atom_density, F.x <= 1.4e-6),
                             (0,True)),
        materials=tungsten
    )

my_model.traps = [trap_1, trap_2, trap_3]

my_model.boundary_conditions = [
    F.DirichletBC(surfaces=[1, 2], value=0, field=0)
]


implantation_temp = sp.Piecewise(
        (F.t,F.t<=560),
        (560, (F.t>560)&(F.t <= 2260)),
        (560 - (F.t-2260),(F.t >2260) & (F.t<=2283)),
        (527,(F.t>2283) &(F.t<=3500) ),
        (527-(F.t-3500),(F.t>3500)&(F.t<=3532)),
        (495,(F.t>3532) &(F.t<=5300)),
        (495-(F.t-5300), (F.t>5300)&(F.t<=5325)),
        (470,(F.t>5325)&(F.t<=7000)),
        (470-(F.t-7000),(F.t>7000)&(F.t<=7095)),
        (375,(F.t>7095)&(F.t<=8925)),
        (375-(F.t-8925),(F.t>8925)&(F.t<=9000)),
	    (300,True)) # K

temperature_ramp = 1  # K/s

start_tds = implantation_time + 50  # s

my_model.T = sp.Piecewise(
        (implantation_temp, F.t < 9000),
        (300,(F.t >9000) & (F.t<start_tds)),
        (300 + temperature_ramp*(F.t-start_tds), True))

my_model.dt = F.Stepsize(
    initial_value=0.5,
    stepsize_change_ratio=1.2,
    max_stepsize=lambda t: 10 if t > start_tds else None,
    dt_min=1e-06,
    milestones=[implantation_time, start_tds],
)

my_model.settings = F.Settings(
    absolute_tolerance=1e10,
    relative_tolerance=1e-9,
    final_time=15000
)

list_of_derived_quantities = [
        F.TotalVolume("solute", volume=1),
        F.TotalVolume("retention", volume=1),
        F.TotalVolume("1", volume=1),
        F.TotalVolume("2", volume=1),
        F.TotalVolume("3", volume=1),
        F.HydrogenFlux(surface=1),
        F.HydrogenFlux(surface=2)
    ]

derived_quantities = F.DerivedQuantities(
    list_of_derived_quantities,
    filename="tds/derived_quantities.csv",  # optional set a filename to export the data to csv
    show_units=True,
)


my_model.exports = [derived_quantities]

my_model.initialise()
my_model.run()

t = derived_quantities.t
flux_left = derived_quantities.filter(fields="solute", surfaces=1).data
flux_right = derived_quantities.filter(fields="solute", surfaces=2).data

flux_total = -np.array(flux_left) - np.array(flux_right)

import matplotlib.pyplot as plt
plt.plot(t, flux_total, linewidth=3)

plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
plt.xlabel(r"Time (s)")
plt.savefig('pre-W-low.png')

Thank you!

2

Hi @Jellycarrot

Have you tried plotting your concentration fields by exporting to XDMF? this may give you some information about what is going on

1 Like
3

I noticed the mesh you’re using is very coarse. This may be why everything is zero.

The volumetric source your are imposing is a gaussian centred on 4 nm with a stdev of 2nm, so very narrow. With such a coarse mesh, this source is so badly discretised that it’s zero everywhere.

You want to have local refinement as shown in the TDS tutorial of the workshop.

Something like

vertices = np.concatenate([
    np.linspace(0, 30e-9, num=200),
    np.linspace(30e-9, 3e-6, num=300),
    np.linspace(3e-6, 20e-6,num=200),
    np.linspace(20e-6, 3e-3, num=200),
])

my_model.mesh = F.MeshFromVertices(vertices)

Also, consider setting the log level to 20 in order to have more information about the newton iterations:

my_model.log_level = 20
1 Like
4

Appreciate for your quick reply! Problem is solved.

Thanks for your time again！

1 Like