# Discrete method for thermodynamic model

Hello Juan,
Sorry to bother you again. In the original program, the drift-diffusion model is used for current density. But after I consider non-isothermal simulation. I think the influence of temperature gradient to current should not be ignored. So I want to use the thermodynamic model to calculate the current density:

However, to calculate the current density, the code use Scharfetter-Gummel scheme to discrete the drift-diffusion model:

``````    Jn = "ElectronCharge*{0}*EdgeInverseLength*edgeV_t_T*kahan3(Electrons@n1*Bern01,  Electrons@n1*vdiff,  -Electrons@n0*Bern01)".format(mu_n)
Jp ="-ElectronCharge*{0}*EdgeInverseLength*edgeV_t_T*kahan3(Holes@n1*Bern01, -Holes@n0*Bern01, -Holes@n0*vdiff)".format(mu_p)
``````

But I am not familiar to this method. I read some related paper, but their expression is very different from the above equation. So I don’t know how to discrete the thermodynamic model to make sure that it is consistent with the original drift-diffusion model. The thermodynamic model just add a items after the DD model, and I am trying to just discrete the temperature gradient items independently and add it to the original expression. I am not sure if it is right. So I hope get some suggestion from you.

Hi @ghost

Please see the description in this document:

Starting with the first equation in section 3.1, you should be able to come up with a suitable discretization. It looks like the clip you show me is for using quasi-Fermi levels. You can either switch to those variables, or you can keep the current one in terms of electron and hole densities.

Thanks, Juan, I will try it.

Hello Juan!
I have read the doc you wrote. It is helpful, and I think I can change the coefficient α to consider the temperature gradient items. Its discrete form is similar to the DD model:

However, there is a question that confuse me a long time ago since I learn Devsim. The original code write the DD model as follow:

``````Jn = "ElectronCharge*{0}*EdgeInverseLength*edgeV_t_T*kahan3(Electrons@n1*Bern01,  Electrons@n1*vdiff,  -Electrons@n0*Bern01)".format(mu_n)
``````

It seems that the expression didn’t calculate the logarithmic mean ( β^bar). The expression use the edge value of β to calculate. So it means it use arithmetic mean to approximate the logarithmic mean. But the Dn is related to the mobility. The transformation of the mobility may not be so smooth，will this have a big impact on the final result?

If the temperature is rapidly varying, it may be best to use the logarithmic mean, but I have been doing isothermal simulation. You may need to experiment to see what the best approach is.