udcm-su / NTMpy

Simulation of multi temperature and multi layer heat diffusion
http://udcm.fysik.su.se/code
MIT License
18 stars 11 forks source link

Thermal boundary conductance at two-layer interface #27

Open Radioteddy opened 1 year ago

Radioteddy commented 1 year ago

Hi,

In addition to the requests in previous issues I would also suggest to enable boundary conditions at interface between two layers accounting for temperature jump, such as in this paper: https://iopscience.iop.org/article/10.1088/0953-8984/27/1/015007/pdf Such a request is strongly motivated, for example, by the urgent need to calculate heat dissipation from a laser-excited sample into substrate(s) for the data analysis from various experiments with optical and X-ray free-electron lasers and design of future experiments.

As for previous requests for temperature-dependent thermal parameters, recently we published a paper analyzing the effect of lattice temperature dependence of e-ph coupling. Maybe you will find it useful for testing new features in your code. https://www.mdpi.com/1996-1944/15/15/5193

VaSca92 commented 1 year ago

Ciao Fedor, sorry for the late response, but github did not notify the issue on my email (or maybe the email server put it in spam, I do not know).

I thank you for your suggestion, these requests are helpful to improve the package and, hopefully, I will start working on it soon.

Regarding the previous requests, I am not sure what issue you are talking about, but I implemented temperature dependence (C, k and G depending on both Te and Tl) on a new version I am currently working on. Right now I have not released the new version, but there is an almost-up-to-date version on my github page (https://github.com/VaSca92/NTMpy) which has a slightly different user interface. Let me know if you want to try it: I can give you some help.

I hope I can give you some better news soon.

Radioteddy commented 1 year ago

Hi Valentino,

Thank you for the link to the new version, I already played with it a bit 🙂 I have a couple of notices/suggestions regarding to it. Will it be better to put them here, or open a new issue in your repository?

VaSca92 commented 1 year ago

Ciao Fedor, you can write your suggestions here. The project is the same, and I will merge it as soon as the update is complete.

Radioteddy commented 1 year ago

Hi Valentino,

  1. Function lambdize2 contains dummy = fun(0,0) (used for testing I assume), which leads to incorrect lambdization of functions looking like $f(te, ta) = f_0 \cdot te/ta$. Such a shape is typical, e.g., for electron thermal conductivity in metals.
  2. Code is not suited for calculations with zero thermal conductivities. If you set K[i][k] = 0, Neumann boundary conditions break down.
  3. The commonly accepted in community relation between standard deviation (which is source.time in code) and FWHM is:

$$ \sigma^2 = FWHM^2 / \sqrt{8 \log{2}} $$

Then the temporal shape of pulse will look like:

$$ f(t) = \exp{(-(t-t_0)^2/2 \sigma^2)} = \exp{(-4\log{2} (t-t_0)^2/FWHM^2)} $$

You have 2 instead of 8 under squared root.

Hope you find these suggestions useful 🙂

VaSca92 commented 1 year ago

Ciao Fedor.

  1. You are right, I already noticed this problem as I used exactly the conductivity you wrote in the message. I corrected it as f(300,300) but it I forgot to re-upload the corrected version. Actually, that line is just to check if the input is already a 2-variables function, so any value (excluding the ones where the function is singular) should work.

  2. This point is somewhat more tricky. In the code the neumann condition is imposed as $$k\ \partial_x \varphi = \text{assigned value}$$ rather than the usual $\partial_x \varphi =$ given condition. This is because $k\ \partial_x\varphi$ is the heat flux, so the neumann conditions impose the heat conservation on the boundary. Whenever the conductivity is zero, the Neumann condition is lost, as no heat crosses the surface, no matter how strong is the temperature gradient. So, yes, the neumann condition is neglected intentionally as the material does not perceive the temperature gradient and the neumann condition looked somewhat inconsistent to me (the code uses a normal diffusion equation also on the boundary points when k = 0). Consider I am not an experimentalism myself, so I could have missed some important aspects of the physical system.

    1. The commonly accepted relation you wrote is actually the only one corrected, I had a different one because I computed it manually and I messed up... you saved me a lot of time by reporting this as I was confidently wrong on that formula.
Radioteddy commented 1 year ago

Hi Valentino, sorry for the extremely late answer.

  1. $k=0$ and Neumann boundary conditions puzzled me because pdepe solver in MATLAB preserves temperature gradient in that case. In your code it's enough to set $k$ very small, e.g. 1e-10, in order to get similar results.

  2. I also found out that you've completely rewritten transfer-matrix algorithm. It will be nice if you can provide reference formulas. I just found that Lambert-Beer solution and transfer-matrix solution for one thick layer do not coincide.

  3. Perhaps, you have fresher version of code? I will appreciate if you share it :))

Cheers, Fedor

VaSca92 commented 1 year ago

Hi Fedor. In the last few months I could not work on the code because of other activities.

I admit the new version of the code is not properly tested, so if you have some case with inconsistent results, please send it to me. Even if you are not sure, but it looks odd, please send it anyway.

I do not remember what reference I used for transfer matrices, some formulas were computed by hand maybe. Anyway a reference is needed, so I will look for it and post it here and on the wiki as soon as I have it. The code on TM is completely rewritten because the original was written by a colleague who is currently work on other things. The two methods are slightly different as I neglect the phase propagation (which is usually a small error), I will update this.

I will probably start to work on the code again by the end of february. I saw on your repository that you also worked on a "two temperature model" code using a compiled language (fortran). This could be a useful counterpart of ntmpy, does it work?

Radioteddy commented 1 year ago

Hi Valentino,

I took refraction index $n+ik = 5.04 + 3.94i$ and absorption length $\delta=16.16$ nm (corresponds to 800 nm wavelength and the used refraction index). You can see that actual spatial profiles in two cases are different. Perhaps some typos in your derivation; hard to tell without reference unfortunately :( image

Another big issue I found is crash of solution if I stack two identical layers one to another. Top picture is one 100 nm layer, bottom - two 50 nm. image image

Regarding my code: yeah, it works quite well for one-layer systems. However, its extension to multilayer structures seems to be rather hard task; I feel myself much more confident with python than with fortran and prefer use it wherever it possible.

I have a few ideas/small possible improvements for your code. Please let me know when you start to work on it again, I will be happy to pull-request them.

VaSca92 commented 1 year ago

The difference between Transfer-Matrices and Lambert-Beer could be due to the fact that TM always consider air before the first layer of material. Hence there is some reflection which might give a different amplitude. However, if the amplitude is still wrong, the fresnel equations are likely bugged. I will check their implementation as soon as possible (I am not sure they work properly).

The divergence with two identical layers could be due to a bad choice of the time step (which is chosen by the code itself). I did not know about this problem and it looks like it does not always happen (or at least I could not reproduce it). Could you give me some more info about the test case you used? (you can either post it here or send it to me via email valentino.scalera@uniparthenope.it) I must admit this is a really ugly bug, does it occour also on the old version of the code?

I will surely let you know when I will start to work on the code again and I really appreciate your help.