Discussion :

[roman.nedellec]

Bonjour a tous.

J'ai un probleme de transfert thermiques a resoudre (equation matricielle a 20 dim d'equadiff d'ordre 1) et je souhaiterais resoudre ce probleme avec vba excel et RK4 car c'est le seul outil que j epossede la ou je travaille actuellement.
ou puis-je trouver l'algo clairement explique pour tenter l'implementer sous vba excel ?
Merci d'avance

Roman

[salseropom]

Bonjour a tous.

J'ai un probleme de transfert thermiques a resoudre (equation matricielle a 20 dim d'equadiff d'ordre 1) et je souhaiterais resoudre ce probleme avec vba excel et RK4 car c'est le seul outil que j epossede la ou je travaille actuellement.
ou puis-je trouver l'algo clairement explique pour tenter l'implementer sous vba excel ?
Merci d'avance

Roman

As-tu cherché sur google ?

http://fr.wikipedia.org/wiki/M%C3%A9...de_Runge-Kutta

l'algo est clair. A toi de l'implémenter

[FR119492]

Salut !

J'ai retrouvé ça dans mes archives Fortran. Il ne te restera qu'à traduire.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
 
      Subroutine D006(SM,N,T,DT,Y)
C
C***********************************************************************
C
C     Bibliothèque JMBPACK
C
C     Série D: Equations différentielles
C
C***********************************************************************
C
C     Méthode de Runge-Kutta d'ordre 4 (1/6 + 2/6 + 2/6 + 1/6)
C
C***********************************************************************
C
C     Version 3.0: Jean-Marc Blanc, juillet 2004
C
C **********************************************************************
C
C     SM     Nom d'un sous-programme fourni par l'utilisateur pour le
C            calcul des dérivées à partir des valeurs de la variable
C            indépendante et des variables d'état. Ce sous-programme
C            doit être de la forme
C              Subroutine SM(N,T,Y,DY)
C                DY(1)= ...
C                DY(2)= ...
C                 .
C                 .
C                DY(N)= ...
C                Return
C              End
C            et avoir été déclaré dans une instruction External.
C
C     N      Ordre du système différentiel.
C
C     T      Variable indépendante.
C
C     DT     Pas d'intégration.
C
C     Y      Vecteur des variables d'état.
C
C **********************************************************************
C
      Implicit None
C
      Integer N
      Real*8 DT,T,Y(N)
C
      Integer I
      Real*8 DY(N),K1(N),K2(N),K3(N),K4(N),Z(N)
C
      Call SM(N,T,Y,DY)
      Do I=1,N
        K1(I)=DT*DY(I)
        Z(I)=Y(I)+K1(I)/2.d0
      End Do
C
      Call SM(N,T+DT/2.d0,Z,DY)
      Do I=1,N
        K2(I)=DT*DY(I)
        Z(I)=Y(I)+K2(I)/2.d0
      End Do
C
      Call SM(N,T+DT/2.d0,Z,DY)
      Do I=1,N
        K3(I)=DT*DY(I)
        Z(I)=Y(I)+K3(I)
      End Do
C
      Call SM(N,T+DT,Z,DY)
      Do I=1,N
        K4(I)=DT*DY(I)
        Y(I)=Y(I)+(K1(I)+2.d0*K2(I)+2.d0*K3(I)+K4(I))/6.d0
      End Do
C
      T=T+DT
      Return
C
      End

Bonne chance
Jean-Marc Blanc

[millie]

Aussi des implémentations en Java et C++ :

http://www.developpez.net/forums/sho...d.php?t=370936

Et en Fortran (mais c'est celui qu'a déjà posté FR119492)

[roman.nedellec]

Merci pour vos reponses
J'ai donc implementer le runge kutta mais je n'obtiens rien de bon...
Des la 3 iteration, je suis completement out of bounds...
Pouriez vous jeter un coup d'oeil a mon code svp? Je ne comprend vraiment pas d'ou ca peut venir
Merci d'avance...

Code vb :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
Public Sub LancerCalcul_Click()
 
'Definition des Variables'
 
'Boucle K1'
Dim j1 As Single
Dim l1 As Single
Dim K1(20) As Single
 
'Boucle K2'
Dim j2 As Single
Dim l2 As Single
Dim K2(20) As Single
 
'Boucle K3'
Dim j3 As Single
Dim l3 As Single
Dim K3(20) As Single
 
'Boucle K4'
Dim j4 As Single
Dim l4 As Single
Dim K4(20) As Single
 
Dim p As Single 'Boucle pas a pas'
Dim i As Single 'Boucle initialisation'
Dim j As Single 'Boucle [Ti] a t+1'
 
Dim Fij(20, 20) As Single 'Matrice Resistances thermiques'
Dim Ti(20) As Single      'Vecteur Temperatures'
Dim MiCi(20) As Single    'Vecteur 1/MC'
 
 
 
'Initialisation'
 
For i = 1 To 20
    For j = 1 To 20
        Fij(i, j) = Worksheets("ci").Cells(i, j + 8).Value
        Ti(i) = Worksheets("ci").Cells(i, 32).Value
        MiCi(i) = Worksheets("ci").Cells(i, 5).Value
     Next j
Next i
 
p = 1
h = 0.005
 
'Iterations'
 
 
'Do Until Ti(17) = 25'
 
For p = 1 To 10 'Pas de calcul'
 
           'Calcul [K1]'
 
            For j1 = 1 To 20
 
               'Initialisation'
                K1j = 0
                Rk1 = 0
 
                For l1 = 1 To 20
 
                    Rk1 = MiCi(j1) * (Fij(j1, l1) * Ti(l1))
                    K1j = K1j + Rk1
 
                Next l1
                    K1(j1) = K1j
                    Worksheets("intermediaire").Cells(j1, 2).Value = K1(j1)
 
            Next j1
 
            'Calcul [K2]'
 
            For j2 = 1 To 20
 
                'Initialisation'
                K2j = 0
                Rk2 = 0
 
                For l2 = 1 To 20
 
                    Rk2 = MiCi(j2) * (Fij(j2, l2) * (Ti(l2) + ((h / 2) * K1(l2))))
                    K2j = K2j + Rk2
 
                Next l2
                    K2(j2) = K2j
                    Worksheets("intermediaire").Cells(j2, 4) = K2(j2)
            Next j2
 
            'Calcul [K3]'
 
             For j3 = 1 To 20
 
                'Initialisation'
                K3j = 0
                Rk3 = 0
 
                For l3 = 1 To 20
 
                    Rk3 = MiCi(j3) * (Fij(j3, l3) * (Ti(l3) + (h / 2) * K2(l3)))
                    K3j = K3j + Rk3
 
                Next l3
                    K3(j3) = K3j
                    Worksheets("intermediaire").Cells(j3, 6) = K3(j3)
            Next j3
 
            'Calcul [K4]'
 
            For j4 = 1 To 20
 
                'Initialisation'
                K4j = 0
                Rk4 = 0
 
                For l4 = 1 To 20
 
                    Rk4 = MiCi(j4) * (Fij(j4, l4) * (Ti(l4) + h * K3(l4)))
                    K4j = K4j + Rk4
 
                Next l4
                    K4(j4) = K4j
                    Worksheets("intermediaire").Cells(j4, 8) = K4(j4)
            Next j4
 
            'Creation [T] a t+1'
 
            For j = 1 To 20
            Ti(j) = Ti(j) + (h / 6) * (K1(j) + 2 * K2(j) + 2 * K3(j) + K4(j))
            Ti(1) = 25
            Ti(2) = 25
            Ti(3) = 25
            Ti(4) = 25
            Ti(8) = 25
            Ti(12) = 25
            Ti(16) = 25
            Ti(20) = 25
            Worksheets("intermediaire").Cells(j, 12) = Ti(j)
            Next j
    h = h + 0.005
 
 
Next p
 
''Loop
 
End Sub

Merci d'avance et bonne journee

[FR119492]

Salut !

Je n'ai pas lu ton code parce que je ne connais pas le VBA, mais je peux peut-être quand même t'aider. Si tu vois à l'écran ou en imprimant tes résultats des valeurs qui semblent être en progression géométrique croissant très rapidement, c'est que ton intégration diverge. A cela, il peut y avoir plusieurs causes:

1. Tu t'es tout simplement trompé en écrivant tes équations ou en les programmant. Alors contrôle tout.
2. Tu as choisi un pas d'intégration trop grand. Essaie avec un pas 2, 5 ou 10 fois plus petit. Il est probable que ça divergera encore, mais moins vite. Si ton système est linéaire, le pas devrait être nettement (10 à 20 fois) plus court que la plus petite constante de temps, que tu peux déterminer par les valeurs propres de la matrice.
3. Tu a commis une maladresse en modélisant ton problème physique et tu es tombé sur un système "stiff", c'est-à-dire dont la solution est un mélange de phénomènes très lents et très rapides. Alors, tu dois reprendre ta modélisation pour obtenir un système non "stiff".
4. Ton problème conduit obligatoirement à un problème "stiff". Dans ce cas, tu reviens sur le forum et tu expliques le plus clairement possible ce que tu dois faire et on essaiera d'aller plus loin.

Bonne chance
Jean-Marc Blanc

[roman.nedellec]

J'ai verifie ma modelisation et je ne pense pas avoir fais d'erreur...mais bon, les transferts thermiques n'etaient pas ma matiere de predilection...
Quand j'affiche mes resultats des K, ils sont tres proches les uns des autres pour chaque equation. Alors que d'apres la formule mathematique, ils devraient quand meme etre differents...Je ne comprend pas, pourtant la formule implementee dans vba est bien ecrite. Est-ce un hasard...?

Ou alors tout simplement que comme tu l'as dis, mon probleme est stiff...

Je vais donc expliquer la base de ma reflection.

Je vais fabriquer des capuchons en aluminium+isolant utilise pour refroidir un capteur. Le principe est donc de placer le capuchon au congelateur (-18deg) puis de le mettre sur le capteur (env 80deg) pour le refroidir. Mon but est de savoir en combien de temps l'ensemble aura atteint la temperature d'equilibre.



J'ai donc decompose mon systeme en elements finis et fais l'hypothese que dans chaque element, la temperature est homogene. Ensuite, j'ai applique le 1er principe de la thermodynamique a chaque element en fonction des elements en contact avec celui considere.

En piece jointe la modelisation ( j'ai considere le demi capuchon, car il est symetrique axialialement)

Je me retrouve donc avec un systeme matriciel de 20 equations, avec 20 temperatures qui evoluent et interagissent les unes avec les autres.

Mon runge kutta est donc implemente de maniere a ce que K1 (sous forme de vecteur donc, avec une valeur par element) K2, K3 et K4 sont calcules a chaque pas de calcul. j'ai pris un pas de 0.005 sec. Et je tire donc ensuite a partir de ces vecteur Ki, mon vecteur Temperature, avec la temperature de chacun de mes elements, a t+1...
Au bout de 3 iterations, j'ai des temperatures qui devraient augmenter et qui diminuent, et d'autres qui sont aux alentours de 34291 deg C...

Y vois tu un probleme stiff ? ou une autre maniere de resoudre mon probleme...?

Merci d'avance
romam

[FR119492]

Salut !

Effectivement, ton système est stiff, comme un simple raisonnement physique permet de le voir: lorsque tu mets en contact deux corps portés à des températures différentes, l'interface prend instantanément une température intermédiaire, les écarts par rapport aux températures initiales étant inversement proportionnelles aux conductivités thermiques des matières correspondantes. Verse du café à 50 degrés dans deux gobelets, l'un en mousse de polyuréthane et l'autre en aluminium, et essaie de boire; tu sentiras tout de suite la différence.

Donc, la température varie très vite à certains endroits et beaucoup plus lentement à d'autres, ce qui veut dire que ton système est stiff. Joseph Fourier l'avait déjà découvert en 1822, dans le cas d'un "pavé" trempé brusquement dans un liquide, et c'est pourquoi il a inventé les séries qui portent son nom: les fonctions propres du laplacien sont des combinaisons de sinus et de cosinus. On atténue le caractère stiff en tronquant les séries, de manière à faire disparaître les composantes dont les constantes de temps sont les plus courtes.

Revenons à ton problème. On peux imaginer une foule de méthodes très savantes et très compliquées, comme de passer par un développement en série de Fourier-Bessel, mais je suis persuadé que dans ton cas ce serait une perte de temps. Je te propose donc celle-ci:

Tu as déjà ton système sous forme canonique:
dT/dt = A * T

avec les conditions initiales
T(t=0) = T0

Tu calcules les valeurs et vecteurs propres de A; les valeurs propres sont les inverses des constantes de temps et les vecteurs propres définissent le changement de variables qui diagonalise la matrice. Tu ne prends en considération que les plus longues constantes de temps et les vecteurs propres correspondants. Si tu appliques le même changement de variables à T0, tu obtiens un ensemble d'équations différentielles qui ne sont plus couplées et que tu intègres analytiquement. Pour finir, tu refais le changement de variables inverse.

Bonne chance.
Jean-Marc Blanc

[roman.nedellec]

Merci beaucoup !!!

Cela faisait une semaine que je cherchais a comprendre mon erreur, soit de programmation soit de modelisation du probleme...Je n'aurais pas imaginer un probleme Stiff...
Merci, je vais donc me lancer sur votre proposition

Bonne journee

Roman

[salseropom]

Pour résoudre des problèmes stiff, il y a en matlab ode23tb, ode15s etc...
Pour ma part, j'utilise un solveur écrit en C, CVode, qui est vraiment très bien (j'ai envie de dire mieux que tous les solveurs de matlab car matlab me calculait très mal la solution numérique). Ce code est gratuit, il y a même le code source, et il y a une interace pour Matlab si tu veux.

regarde ici (il y a des exemples pour l'utliser) :

http://www.llnl.gov/CASC/sundials/

la méthode de Jean-Marc Blanc est très bien (je ne la connaissais pas) mais je ne pense pas qu'elle puisse s'appliquer à tous les systèmes raides. Dans ma thèse, j'étudie l'intéraction entre 2 phénomènes physiques : la turbulence et le cinétique chimique (appliquée sur la combustion des moteurs). J'ai donc besoin de ces 2 échelles de temps (dans mon modèle) car lorqu'il y a l'explosion, la température de mon mélange augmente très très rapidement, d'où mon pb stiff.

[FR119492]

Salut!

Tout à fait d'accord avec salseropom! Juste un petit complément d'information: pour les systèmes stiff, les méthodes usuelles s'appellent "Euler implicite" et "Adams-Gear". Je n'ai malheureusement pas les sources sous la main, mais tu peux chercher sous ces désignations.

Bonne chance.
Jean-Marc Blanc

[j.p.mignot]

voir une discussion à ce sujet :
http://www.developpez.net/forums/sho...d.php?t=104554