# Second peak in TDS simulation for a single trap

**URL:** <https://festim.discourse.group/t/second-peak-in-tds-simulation-for-a-single-trap/47>\
**Category:** User Support\
**Created:** [May 25, 2024, 7:37am UTC](https://festim.discourse.group/t/second-peak-in-tds-simulation-for-a-single-trap/47 "2024-05-25T07:37:34Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jellycarrot](https://avatars.discourse-cdn.com/v4/letter/j/7ea924/32.png) [@Jellycarrot](https://festim.discourse.group/u/Jellycarrot)\
**Post date:** [May 25, 2024, 7:37am UTC](https://festim.discourse.group/t/second-peak-in-tds-simulation-for-a-single-trap/47/1 "2024-05-25T07:37:34Z")

</div>

Hi，

I’m working on thermo-desorption spectra of hydrogen in tungsten with FESTIM.However,there’s a question that troubles me:  
Why does FESTIM’s TDS simulation for a single trap have a second peak after the first main peak?

Thanks!

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [May 25, 2024, 10:00am UTC](https://festim.discourse.group/t/second-peak-in-tds-simulation-for-a-single-trap/47/2 "2024-05-25T10:00:01Z")

</div>

Hi @Jellycarrot and welcome to FESTIM!

Do you have an example that you could share with us? A minimal working example that we could run.

A few ideas:

- a peak of mobile particles
- a desorption peak of the right hand side

But again it’s hard to tell without having a script

---

<div class="post-metadata">

**Author:** ![Jellycarrot](https://avatars.discourse-cdn.com/v4/letter/j/7ea924/32.png) [@Jellycarrot](https://festim.discourse.group/u/Jellycarrot)\
**Post date:** [May 25, 2024, 11:34am UTC](https://festim.discourse.group/t/second-peak-in-tds-simulation-for-a-single-trap/47/3 "2024-05-25T11:34:16Z")

</div>

Thank you for your prompt attention!

Here is my code as an example

```auto
import festim as F
my_model = F.Simulation()

import numpy as np

vertices = np.concatenate([
    np.linspace(0, 30e-9, num=200),
    np.linspace(30e-9, 49.97e-6, num=500),
    np.linspace(49.97e-6, 50e-6, num=200),
])

my_model.mesh = F.MeshFromVertices(vertices)

tungsten = F.Material(
    id=1,
    D_0=4.1e-07, # m2/s
    E_D=0.39, # eV
)

my_model.materials = tungsten

import sympy as sp

implantation_time = 400 # s

ion_flux = sp.Piecewise((2.5e19, F.t <= implantation_time), (0, True))

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

my_model.sources = [source_term]

w_atom_density = 6.3e28 # atom/m3

trap_1 = F.Trap(
	k_0=4.1e-7/(1.1e-10**2*6*w_atom_density),
	E_k=0.39,
	p_0=2e13,
        E_p=1.25,
        density=2e-4*w_atom_density,
        materials=tungsten
    )

my_model.traps = [trap_1]

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

#Temperature
implantation_temp = 300 # K
temperature_ramp = 1 # K/s

start_tds = implantation_time + 50 # s

my_model.T = F.Temperature(
    value=sp.Piecewise(
        (implantation_temp, F.t < start_tds),
        (implantation_temp + temperature_ramp*(F.t-start_tds), True))
)

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

#Settings
my_model.settings = F.Settings(
    absolute_tolerance=1e10,
    relative_tolerance=1e-09,
    final_time=1000
)

list_of_derived_quantities = [
        F.TotalVolume("solute", volume=1),
        F.TotalVolume("retention", volume=1),
        F.TotalVolume("1", 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
)

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('1.png')
trap_1 = derived_quantities.filter(fields="1").data

contribution_trap_1 = -np.diff(trap_1)/np.diff(t)

plt.plot(t, flux_total, linewidth=3)
plt.plot(t[1:], contribution_trap_1, linestyle="--", color="grey")
plt.fill_between(t[1:], 0, contribution_trap_1, facecolor='grey', alpha=0.1)

plt.xlim(450, 1000)
plt.ylim(bottom=-0.05e18, top=0.6e18)
plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
plt.xlabel(r"Time (s)")

plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
plt.xlabel(r"Time (s)")
plt.savefig("2.png",dpi=600)

```

In the code above (which is slightly modified from your Task 02 in FESTIM-workshop)，I try to run a system of 50μm with single trap.

Then I’ve got two graphs like this:  
 ![jyz](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/2a50c490f9b25c61d64d2980f6772e8f9b8f1e12.png)

 ![yz](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/5275b3196740b7c88f5e330f61da4ffafa948423.png)

So I wonder how did the second peak appear.

Thanks again for answering me!

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [May 25, 2024, 12:13pm UTC](https://festim.discourse.group/t/second-peak-in-tds-simulation-for-a-single-trap/47/4 "2024-05-25T12:13:53Z")

</div>

Instead of plotting contribution\_trap\_1 can you plot flux\_left and flux\_right instead?

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [May 27, 2024, 8:20am UTC](https://festim.discourse.group/t/second-peak-in-tds-simulation-for-a-single-trap/47/5 "2024-05-27T08:20:14Z")

</div>

@Jellycarrot I confirm this is the desorption from the right hand side

 ![image](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/4ad9694380c00b3865c81896e838d4a0961e79cb.png)

```python
plt.plot(t, flux_total, linewidth=3, color="black")

plt.stackplot(
    t,
    -np.array(flux_left),
    -np.array(flux_right),
    colors=["tab:blue", "tab:orange"],
    labels=["Left", "Right"],
)

```

If you increase the sample size to say 1 mm then this peak disappears:

 ![image](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/3faf3150dcdcdf0a3ad9263de2e14aa4838939d3.png)

This is because when you have a thin sample, after detrapping, hydrogen diffuses to both sides and there is a small delay after which it has reached the right side. This also happens with thicker samples except there is a much larger distance to diffuse through and therefore H only desorbs from the left.

---

<div class="post-metadata">

**Author:** ![Jellycarrot](https://avatars.discourse-cdn.com/v4/letter/j/7ea924/32.png) [@Jellycarrot](https://festim.discourse.group/u/Jellycarrot)\
**Post date:** [May 27, 2024, 9:53am UTC](https://festim.discourse.group/t/second-peak-in-tds-simulation-for-a-single-trap/47/6 "2024-05-27T09:53:35Z")

</div>

Thank you，It really helps！

I plot too：

 ![022 (1)](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/1d2b955b2d275bdebb2fc71c6c0f5b771e095fd4.png)

My problem is solved.

Wish you have a good day！