Python Forum
Heat equation (not stationary)
Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Heat equation (not stationary)
#1
Good morning,

I would be very greatful if someone take the time to help me.
I have a scientific research study to deal with, my topic is : The the evolution of the transfers in a mug with temperature.
I have now to simulate this with the Heat equation in a 1D problem (to simplify) as a wall.

I have already modelised it, (with termal sources so as to simulate the exchange of the mug with the water on the one hand (i=0, water temperature not fixed) and the air (i=N+1, air temperature fixed at Text)).
So I have 3 expressions for the temperature at the moment k+1 that requires the température at the moment k.
I would like to calculate these températures time after time ( incrementation 1 by 1 for k ) and to finally only return the temperature at the interface with the water (i=0).

Nevertheless, my code is not working:
File "F:\TIPE\Modélisation Python\Programme.py", line 71
Temp.append(L[k][N]+dt/(dr*rho*capa(L[k][N],0,0))*(h*(Text-L[k][N])+lambd(L[k][N],0,0)/dr*(L[k][N-1]-L[k][N])))
^
SyntaxError: invalid syntax


I have tried to change a lot of things with my basic ability, but I can't go over this "SyntaxError" and I don't know where could be the error. I'll be greatful for you help.

L a list of list, it contains the list of the temperature at the points i=0 to i=N+1, at each instant k from 0 to M.
capa and lambd are functions that provide the evolution of physic parameters with the temperature.
Interface_eau is the list of temperature at the point i=0 at each moment k from 0 to M.


## Importation des modules

import numpy as np
import matplotlib.pyplot as plt

## Definition des fonctions

# Données fixes:

rho=2500
le=1.4e-3
ce=4186
rhoe=1000
h=10
e=5e-3

#Coefficients dépendances thermiques

lambda0=-11.91
lambda1=87.6e-3
lambda2=-145e-6

c0=-243e3
c1=2e3
c2=-2

#Fonctions dépendances thermiques

def lambd(T,coeff1,coeff2):
      '''O si prends les valeurs exp, 1 si prends les valeurs entrées pour coeff1 et coeff2'''
      if coeff1==0 and coeff2==0:
            coeff1=lambda1
            coeff2=lambda2
      return (lambda0+coeff1*T+coeff2*T*T)

def capa(T,coeff1,coeff2):
      '''O si prends les valeurs exp, 1 si prends les valeurs entrées pour coeff1 et coeff2'''
      if coeff1==0 and coeff2==0:
            coeff1=c1
            coeff2=c2
      return (c0+coeff1*T+coeff2*T*T)

#Paramètres de pas

N=10
dt=10
dr=e/N+1

#Paramètres initiaux

Tdeb=90
Text=20 

#Fonction méthode explicite

def résolution(M):
      '''S'arrete à tmax=M*dt, renvoie la liste des T en r=0 aux différents temps'''
      Interface_eau=[Tdeb]
      Temps=[elt for elt in range(M+1)]
      L=[[Tdeb]+(N)*[Text]]
      for k in range(0,M):
            Temp=[]
            #Calcul interface eau i=0
            Temp.append(L[k][0]+lambd(L[k][0],0,0)*dr*dt/(dr*rho*capa(L[k][0],0,0)-le*rhoe*ce)/dr/dr*(L[k][1]-L[k][0]))
            #Calcul intérieur tasse i entre 1 et N
                  for i in range(1,N):
                  Temp.append((L[k][i]+lambd(L[k][i],0,0)*dr*dt/(dr*rho*capa(L[k][i],0,0))/dr/dr*(L[k][i+1]+L[k][i-1]-2*L[k][i]))
            #Calcul interface air i=N+1
            Temp.append(L[k][N]+dt/(dr*rho*capa(L[k][N],0,0))*(h*(Text-L[k][N])+lambd(L[k][N],0,0)/dr*(L[k][N-1]-L[k][N])))
            L.append(Temp)
            Interface_eau.append(Temp[0])
      plt.plot(Temps, Interface_eau)

#Lançage

résolution(10)
Reply
#2
I do not know if you copied your code incorrectly or that's that way it is, but indent on line 66 is wrong. Then, you have 2 opening braces on line 67.

Besides - writhing veeeeery long calculations in one string without spaces makes an unreadable and unsupportable code. In Python, it is customary to surround operation symbols with spaces, for readability.
Couple of more points:
  • One-letter variable name add to code unreadability. Variable names should be meaningful and recognizable, and yes - if you ask help in an international forum, in English
  • In Python, you don't iterate over indices - you can, but you should not. Below - instant readability improvement, much clearer than L[k][i] **hand**
    for l_row in L:
    ....
        for l_cell in in l_row:
        ......
Test everything in a Python shell (iPython, Azure Notebook, etc.)
  • Someone gave you an advice you liked? Test it - maybe the advice was actually bad.
  • Someone gave you an advice you think is bad? Test it before arguing - maybe it was good.
  • You posted a claim that something you did not test works? Be prepared to eat your hat.
Reply
#3
Ok thank you for your reponse, I will edit my code by the night following your advice. (About the line number : I have just remove some French details , reducing the number of line; I'll change the indent on line 66 too, I have failed while reproducing the indent on the forum)
Reply
#4
One more point - default variable values; instead of this monstrosity
def lambd(T,coeff1,coeff2):
     '''O si prends les valeurs exp, 1 si prends les valeurs entrées pour coeff1 et coeff2'''
     if coeff1==0 and coeff2==0:
           coeff1=lambda1
           coeff2=lambda2
that shorthand
def lambd(T, coeff1=lambda1, coeff2=lambda2):
And lambd(L[k][i],0,0) turns into lambd(L[k][i])

PS lambda being a reserved word in Python, using names very similar to it may be not such a good idea
Test everything in a Python shell (iPython, Azure Notebook, etc.)
  • Someone gave you an advice you liked? Test it - maybe the advice was actually bad.
  • Someone gave you an advice you think is bad? Test it before arguing - maybe it was good.
  • You posted a claim that something you did not test works? Be prepared to eat your hat.
Reply
#5
It is done ! There is the new code, hope you find it more clear.

Still with the message error :
File "F:\TIPE\Modélisation Python\Programme.py", line 65
Temperature.append( l_row[N] + dt / (dr*rho*capa(l_row[N])) * ( h*(Text - l_row[N]) + lambd(l_row[N])/dr * (l_row[N-1] - l_row[N]) ) )
^
SyntaxError: invalid syntax





## TIPE

#Importation

import matplotlib.pyplot as plt

#Numeric values

rho=2500
le=1.4e-3
ce=4186
rhoe=1000
h=10
e=5e-3

#Coefficients for thermal dependency

conduc0=-11.91
conduc1=87.6e-3
conduc2=-145e-6

capa0=-243e3
capa1=2e3
capa2=-2

#Functions for thrmal dependency

def conduc(T, coeff1=conduc1, coeff2=conduc2): #Thermal conductivity
       return (conduc0 + conduc1*T + conduc2*T**2)

def capa(T, coeff1=capa1, coeff2=capa2): #Thermal capacity per mass
       return (capa0 + coeff1*T + coeff2*T**2)

#Parameters of discreditation

N=10
dt=10
dr=e/N+1

#Initial parameters

T_start=90
T_ext=20 

#Function

def résolution(M):
       '''Stop at time_max=M*dt, return the list of Temperatures at the interface with the water at every time'''

       Interface_water=[T_start]
       Time = [elt for elt in range(M+1)]
       L = [ [T_start] + N*[T_ext] ] #List where every temperature calculated will remain

       for l_row in L:
              Temperature=[]

              #Calculation at the interface water (i=0)
              Temperature.append( l_row[0] + lambd(l_row[0])*dr*dt / (dr*rho*capa(l_row[0]) - le*rhoe*ce) /dr/dr * (l_row[1] - l_row[0]) )

              #Calculation in the mug (i from 1 to N)
              for i in range(1, N+1):
                     Temperature.append( (l_row[i] + lambd(l_row[i])*dr*dt / (dr*rho*capa(l_row[i])) /dr/dr * (l_row[i+1] + l_row[i-1] - 2*l_row[i]) )

              #Calculation at the interface air (i=N+1)
              Temperature.append( l_row[N] + dt / (dr*rho*capa(l_row[N])) * ( h*(Text - l_row[N]) + lambd(l_row[N])/dr * (l_row[N-1] - l_row[N]) ) )


              L.append(Temperature)
              Interface_water.append(Temperature[0])

       plt.plot(Time, Interface_water)
       return (Interface_water)


#Execution

résolution(10)
Reply
#6
Is there anyone who has the key to my problem, please ?
Reply
#7
Dear Zulian,

I have to repeat myself - you long formulas are not easy to read. Try to split it into smaller, manageable calculations - or just load your code into any decent IDE. 

Your problem is most probably un-balanced braces.

Even Notepad ++ - if you're on Windows - allows you to check paired braces.
Test everything in a Python shell (iPython, Azure Notebook, etc.)
  • Someone gave you an advice you liked? Test it - maybe the advice was actually bad.
  • Someone gave you an advice you think is bad? Test it before arguing - maybe it was good.
  • You posted a claim that something you did not test works? Be prepared to eat your hat.
Reply
#8
The previous line of code (62 in your last post) is missing a closing parenthesis.

There is no shame (and no performance penalty) in splitting complex expressions and using intermediate variables.
Unless noted otherwise, code in my posts should be understood as "coding suggestions", and its use may require more neurones than the two necessary for Ctrl-C/Ctrl-V.
Your one-stop place for all your GIMP needs: gimp-forum.net
Reply
#9
Thank you ! For having helped me; here is the code corrected and working for those who may need it in a future.
(There were several errors : The "for" had no way to be stopped and two errors on incrementation )

I'm gonna change it now to respond to my problem, nevertheless I believe that I will no more need your help for this, so Thank You !



## TIPE

#Importation

import matplotlib.pyplot as plt

#Numeric values

rho=2500
le=1.4e-3
ce=4186
rhoe=1000
h=10
e=5e-3

#Coefficients for thermal dependency

conduc0=-11.91
conduc1=87.6e-3
conduc2=-145e-6

capa0=-243e3
capa1=2e3
capa2=-2

#Functions for thrmal dependency

def conduc(T, coeff1=conduc1, coeff2=conduc2): #Thermal conductivity
        return (conduc0 + conduc1*T + conduc2*T**2)

def capa(T, coeff3=capa1, coeff4=capa2): #Thermal capacity per mass
       return (capa0 + coeff3*T + coeff4*T**2)

#Parameters of discreditation

N=50
dt=1000000
dr=e/N+1

#Initial parameters

T_start=90
T_ext=20 

#Function

def résolution(M,coeff1,coeff2,coeff3,coeff4):
       '''Stop at time_max=M*dt, return the list of Temperatures at the interface with the water at every time'''

       Interface_water=[T_start]
       Time = [elt for elt in range(M)]
       L = [ [T_start] + N*[T_ext] ] #List where every temperature calculated will remain

       for l_row in L:
              Temperature=[]

              #Calculation at the interface water (i=0)
              Temperature.append( l_row[0] + conduc(l_row[0])*dr*dt / (dr*rho*capa(l_row[0]) - le*rhoe*ce) /dr/dr * (l_row[1] - l_row[0]) )

              #Calculation in the mug (i from 1 to N)
              for i in range(1, N):
                     Num = conduc(l_row[i]) * dt #Numerator
                     Den = rho * capa(l_row[i]) * dr**2 #Denominator
                     Temperature.append( l_row[i] + (l_row[i+1] + l_row[i-1] - 2*l_row[i]) * Num / Den )
                     # T_i + ( T_i+1 + T_i-1 - 2 * T_i ) * Num/Den

              #Calculation at the interface air (i=N+1)
              Temperature.append( l_row[N] + dt / (dr*rho*capa(l_row[N])) * ( h*(T_ext - l_row[N]) + conduc(l_row[N])/dr * (l_row[N-1] - l_row[N]) ) )


              L.append(Temperature)
              Interface_water.append(Temperature[0])

              if len(Interface_water)>=M : #Stop "for"
                     Droite=len(Time)*[55]
                     plt.plot(Time, Droite)
                     plt.plot(Time, Interface_water)
                     return (Interface_water, L)


#Execution

print(résolution(100))
Reply


Possibly Related Threads…
Thread Author Replies Views Last Post
  animating 2d heat map schniefen 0 1,314 Nov-20-2022, 10:00 PM
Last Post: schniefen
  Heat equation HugoL 1 28,373 Mar-07-2017, 05:22 PM
Last Post: micseydel

Forum Jump:

User Panel Messages

Announcements
Announcement #1 8/1/2020
Announcement #2 8/2/2020
Announcement #3 8/6/2020