Régime stationnaire

unsteady.py

@@ -1,5 +1,7 @@
 # -*- coding: utf-8 -*-
 
+import numpy as np
+import matplotlib.pyplot as plt
 
 class Spreader_2D_adim:
     """Spreader 2D adimensionné"""
@@ -8,10 +10,10 @@ class Spreader_2D_adim:
     K = 1
     S = 1
     B = 1
-    Q = 1
+    Q = 0
 
     # Initialisation
-    def __init__(self, K, S, B, Q):
+    def __init__(self, K, S, B, Q=0):
         """Constructeur de Spreader_2D_adim"""
         self.K = K
         self.S = S
@@ -44,11 +46,48 @@ class Spreader_2D_adim:
         return self.Q
 
     # Méthodes privées
-    
+    def _beta_old(self, n, Y):
+        F = self.K * self.S
+        S = self.S
+        B = self.B
+        P = B/F
+        PI = np.pi
+
+        return 2*np.sinc(n/S)*((1/S - P/(n*PI))*np.exp(-(2-Y)*F*n*PI/S)
+                                + (1/S + P/(n*PI))*np.exp(-Y*F*n*PI/S))\
+            / (F*(P + n*PI*np.tanh(F*n*PI/S)/S)
+                * (1 + np.exp(-2*F*n*PI/S)))
+
+    def _beta(self, n, Y):
+        Knpi = self.K*n*np.pi
+        return 2*np.sinc(n/self.S)*((Knpi - self.B)*np.exp(-Knpi*(2-Y))
+                                    + (Knpi + self.B)*np.exp(-Knpi*Y))\
+            / (self.S*Knpi*(self.B + Knpi*np.tanh(Knpi))
+                * (1+np.exp(-2*Knpi)))
 
     # Méthodes publiques
-    def theta_ss(self, X, Y):
+    def theta_ss(self, X, Y, N=100):
         """Température en régime permanent"""
-        return
+        f = (1 - Y + 1/self.B)/self.S
+        g = self.Q*(1/self.B + (1 - 0.5*Y**2))
+        return f + g + sum([self._beta_old(n, Y)*np.cos(n*np.pi*X) for n in range(1, N+1)])
+
+
+# Execution en temps que main
+if __name__ == "__main__":
+
+    K = 0.025
+    S = 4
+    B = 0.005
+    Q = 0
+
+    spreader = Spreader_2D_adim(K, S, B, Q)
+    
+    X = np.linspace(0, 1, 20)
+    Y = np.linspace(0, 1, 20)
     
+    XX, YY = np.meshgrid(X, Y)
     
+    plt.figure()
+    plt.contourf(XX, YY, spreader.theta_ss(XX, YY))
+    plt.colorbar()