Update Codes/Résolution numérique équation de Laplace.py, Codes/Modélisation...

Update Codes/Résolution numérique équation de Laplace.py, Codes/Modélisation des lignes de courant.py, Readme.txt

Codes/Modélisation des lignes de courant.py

+import matplotlib.pyplot as plt
+import numpy as np
+import math
+import random as rd
+#Particule fluide dans un plan (0,x,y)
+def Vr(x,y,v0,R,w):       #Vitesse de la particule de fluide selon er
+    r=math.sqrt(x**2+y**2)
+    costh=x/r
+    return(v0*costh*(1-(R**2)/(r**2)))
+def Vthet(x,y,v0,R,w):    #Vitesse de la particule de fluide selon ethet
+    r=math.sqrt(x**2+y**2)
+    sinth=y/r
+    return(-v0*sinth*(1+(R**2)/(r**2))+((R**2)*w)/r)
+def Vx(x,y,v0,R,w):       #Vitesse particule fluide selon x
+    if x<-influence or x>influence:
+        return(0)   #Vitesse de la particule ne dépend pas de l'Effet magnus si
+                         #trop loin du rotor donc considérée nulle
+    costh=x/math.sqrt(x**2+y**2)
+    sinth=y/math.sqrt(x**2+y**2)
+    return(Vr(x,y,v0,R,w)*costh-Vthet(x,y,v0,R,w)*sinth)
+def Vy(x,y,v0,R,w):       #Vitesse particule fluide selon y
+    if x<-influence or x>influence:
+        return(0)
+    costh=x/math.sqrt(x**2+y**2)
+    sinth=y/math.sqrt(x**2+y**2)
+    return(sinth*Vr(x,y,v0,R,w)+costh*Vthet(x,y,v0,R,w))
+
+
+x0 =-5         #Valeur négative car à gauche du cylindre
+y0 =-2          #Idem
+v0 =2
+R =0.055
+w =0
+influence =0.8
+perturbation=0.1
+
+def trajectoire(x0,y0,v0,R,w):#Courbe qui débute en (x0,y0)
+    i=x0
+    dt=0
+    Lx=[]
+    Ly=[]
+    while i<-x0:
+        Vx1=Vx(i,y0,v0,R,w)
+        if Vx1==0:
+            y1=y0
+        else:
+            dt=(i+1-i)/Vx1   #Mesure du temps
+            y1=(Vy(i,y0,v0,R,w)*dt +y0)
+            y0=y1
+        if i>-perturbation and i<perturbation: #Si la trajectoire est trop proche du cylindre
+                                                 #elle est supprimée
+            if y1>-perturbation and y1<perturbation:
+                return([],[])
+        Lx.append(i)
+        Ly.append(y1)
+        i=i+0.05   #Pas choisi entre deux abscisses x
+    return(Lx,Ly)
+
+#Programme principal : traçé du cercle de rayon R et des lignes de champ confondues avec
+#les trajectoires des particules car écoulement considéré stationnaire
+
+def cercle(x,y,R):
+    i=-R
+    Lx=[]
+    Ly=[]
+    L_y=[]
+    while i<=x+R:
+        j=math.sqrt(R**2-i**2)
+        Lx.append(i+x)
+        Ly.append(j+y)
+        L_y.append(y-j)
+        i+=R/100
+    j=0
+    Lx.append(R)
+    Ly.append(0)
+    L_y.append(0)
+    plt.plot(np.array(Lx),np.array(Ly),"r")
+    plt.plot(np.array(Lx),np.array(L_y),"r")
+
+    y=y0
+    while y<-y0:
+        (Lx,Ly)=trajectoire(x0,y,v0,R,w)
+        Tx=np.array(Lx)
+        Ty=np.array(Ly)
+        plt.plot(Tx,Ty)
+        y=y+0.1 #Pas entre deux particules loin du rotor
+
+cercle(0,0,R)
+plt.axis('equal')
+plt.plot(np.array([0,0]),np.array([y0,-y0]),"black") #Axe des abscisses
+plt.plot(np.array([x0,-x0]),np.array([0,0]),"black") #Axe des ordonnées
+plt.show()
+
+
+
+
+

Codes/Résolution numérique équation de Laplace.py

+import matplotlib.pyplot as plt
+import math
+import numpy as np
+import scipy as scp
+import numpy.random as rd
+import time
+
+#Généralités pour décrire le domaine choisi
+def graphe_bords(B):
+    """Affiche un graphique des bords du domaine.
+    """
+    Nx,Ny = B.shape
+    plt.figure()
+    plt.imshow(B,origin='lower',cmap='binary', interpolation='Nearest')
+    plt.colorbar()
+    plt.title('matrice B des bords du domaine')
+    plt.show()
+
+def graphe_equipot(B,f,a=0,nb_equipot=25):
+    """Affiche les bords du domaine B en noir, et des surfaces équipotentielles de f."""
+    Nx,Ny = B.shape
+    plt.figure()
+    x = np.linspace(0,Nx-1,Nx)        # de 0 à Nx-1 inclus
+    y = np.linspace(0,Ny-1,Ny)
+    plt.imshow(B,origin='lower',cmap='binary',interpolation='nearest')
+    if (a==1):
+        plt.imshow(f,origin='lower')  # colorplot de f
+        plt.colorbar()
+    cont=plt.contour(y,x,f,nb_equipot,colors='k')
+    plt.title('equipotentielles')
+    plt.clabel(cont,fmt='%1.1f') # le 1.1 controle le nombre de chiffres après la virgule sur l'affichage
+    plt.show()
+
+def w_opt(Nx,Ny):
+    return 2./(1.+np.pi/(Nx*Ny)*((Nx**2+Ny**2)/2.)**0.5)
+#
+def Magnus_():
+
+    def initialisation_effet_magnus(B,phi,a,b,R):
+        """Initialise un cadre avec à gauche phi=a, à droite phi=b, en bas et en haut
+        phi qui va linéairement de a jusqu'a b, et un cercle représentant le cylindre, au potentiel a."""
+        # bord bas :
+        B[0,:] = 1
+        phi[0:] = np.linspace(a,b,Ny)
+        # bord haut :
+        B[0,:] = 1
+        phi[0:] = np.linspace(a,b,Ny)
+        # bord gauche :
+        B[:,Ny-1] = 1
+        phi[:,Ny-1] = a
+        # bord droit :
+        B[:,Ny-1] = 1
+        phi[:,Ny-1] = b
+        # Cercle centré au centre de l'image :
+        i0 = Nx//2
+        j0 = Ny//2
+        for i in range(Nx):
+            for j in range(Ny):
+                if (i-i0)**2+(j-j0)**2 <= (R)**2:  # cercle de centre i0, j0 et de rayon R
+                    B[i,j] = 1
+                    phi[i,j] = a
+
+    def une_iteration_GS(B,f):
+        f2 = f.copy()
+        ecart = 0.
+        for i in range(1,Nx-1):  # on va de 1 à N-2 car pas traiter le cadre (Dirichlet)
+            for j in range(1,Ny-1):
+                if B[i,j]==0:    #pas sur un bord
+                    f[i,j] = (1.-w)*f[i,j] + w*(f[i+1,j]+f[i,j+1]+f[i-1,j]+f[i,j-1])/4.
+                    ecart = ecart + (f[i,j] - f2[i,j])**2
+        return (ecart/(Nx*Ny))**0.5
+
+    def iterations_GS(B,f,eps):
+        ecart = eps+1.
+        nb_iterations = 0
+        while ecart > eps:
+            ecart = une_iteration_GS(B,f)
+            nb_iterations += 1
+        return nb_iterations
+
+    def calcul_V(phi):
+        Vx = np.zeros((Nx,Ny))
+        Vy = np.zeros((Nx,Ny))
+        norme_V = np.zeros((Nx,Ny))
+        for i in range(1,Nx-1):
+            for j in range(1,Ny-1):
+                Vx[i,j] = -(phi[i+1,j]-phi[i-1,j])/(2.*delta)
+                Vy[i,j] = -(phi[i,j+1]-phi[i,j-1])/(2.*delta)
+        norme_V = (Vx**2+Vy**2)**0.5
+
+        return Vx, Vy, norme_V
+
+    Nx, Ny = 150,150
+    w = w_opt(Nx,Ny)
+    delta = 3.e-2 # pas de la grille en m
+    B = np.zeros((Nx,Ny),bool)
+    phi = np.zeros((Nx,Ny))
+    initialisation_effet_magnus(B,phi,0,(Nx-1)*delta*1.2,0.055)
+    iterations_GS(B,phi,1.e-3)
+
+    Vx, Vy, norme_V = calcul_V(phi)
+    max_V = np.max(norme_V)
+    print("Valeur maximale de V : " + str(max_V) + "m/s")
+    graphe_equipot(B,phi,1,40)
+

Readme.txt

+Ces deux programmes interviennent dans l'étude numérique de l'effet Magnus sur un rotor de Flettner.
+
+En effet, le premier programme correspond à la modélisation des lignes de courant autour d'un cylindre en rotation dans un fluide. Ainsi, en utilisant les équations obtenues, on peut effectuer une modélisation correcte en prenant en compte les bords ainsi qu'une zone où les lignes de courant disparaissent, à savoir là ou se situe le cylindre.
+
+Le second programme utilise la méthode des différences finies, à savoir une discrétisation du Laplacien via un maillage rectangulaire. On rentre ainsi nos conditions de bords, et on adopte la méthode de Gauss-Seidel. 
+
+Ici, la difficulté informatique et algorithmique est relativement faible puisque la difficulté réside dans la correcte modélisation physique du problème, la prise en compte des bonnes équations et conditions de bords, et surtout l'utilisation du bon modèle.