<jupyter_start><jupyter_text>École des Ponts ParisTech  SPH pour l'hydraulique   Projet SPH GMM 2021  Parois mobiles   Auteurs : <NAME> & <NAME>  Encadrants : <NAME> et <NAME>   3. Etude d'un écoulement suite à la chute d'un solide.  3.1 Explication de la situation étudiée. Si dans les deux premières parties du projet on imposait le mouvement des parois, l'objectif de cette dernière partie est de programmer les fonctions adéquates afin de modéliser le mouvement d'un solide. On choisit de simuler un écoulement induit par la chute d'un solide dans un bassin.Précisons que l'on ne résout les équations de la dynamique que pour le déplacement selon la verticale. En effet, nous n'avons pas résolu les équations du moment.   3.2 Etude numérique  Import des fichiers et des modules utiles pour programmer un code SPH.<jupyter_code>%load_ext autoreload
%autoreload 2
import numpy as np
from sys import exit
import os.path
from os import path
import csv
import time
import matplotlib
import matplotlib.pyplot as plt
#matplotlib.rcParams['text.usetex'] = True
from src.spyh import *
from src.checkTD3 import *
from src.sphvar import *
from src.plotParticles import *
from src.state import *
from src.contrib import *
from src.analytical_solutions import *
from src.solidStuffManagement import *<jupyter_output><empty_output><jupyter_text>Définition des grandeurs du problèmes et du fluide<jupyter_code>#FLUID PARAMETERS 
e = 0.5 #half width in meters
U = 1 # Upper bound velocity
#FLUID PARAMETERS
Re = 1 # Reynolds number
rhoF = 1000
mu = 200 #(Pa.s)
print(mu)
grav = np.array([0.0,-9.81])#np.array([0.0,-9.81])#np.array([0.0,0.0]) #m.s-2#OTHER FLUID PARAMETERS
c0 = 10*U
gamma = 7
B = rhoF*c0**2/gamma 
#DENSITY & SHEPARD THRESHOLDS : 
shepardMin = 10**(-6)
rhoMin = 0.5*rhoF
rhoMax = 1.5*rhoF<jupyter_output>200<jupyter_text>Définition des propriétés des particules<jupyter_code>#GEOMETRY
Lx = 2    #Cavity length (m)
Ly = 2    #Cavity height (m)
lx = 0.5  #Solid length  (m)
ly=1      #Solid height  (m)
#PARTICLES & SPACES PARAMETERS : 
dr = ly/50
h = smthfc*dr
m=dr*dr*rhoF
lspace = 2*h
#SOLID PARAMETERS
ls = 0.4        
ws = 0.2
rhoS = 100
ms = dr*dr*rhoS<jupyter_output><empty_output><jupyter_text>Création du domaine de calcul<jupyter_code>#COMPUTATION DOMAIN : 
xOrigin = -nBound*dr
yOrigin = -nBound*dr
xSize = Lx+2*nBound*dr
ySize = Ly+2*nBound*dr
xMax = xOrigin+xSize
yMax = yOrigin+ySize
#INIT SPACES & PART:
posSpace,neibSpace,partSpace,listNeibSpace = \
init_spaces(xOrigin,yOrigin,xSize,ySize,lspace,dr)
part = init_particles()
#INIT PART 
part = addBox(part,[Lx,Ly/2],FLUID,dr,rhoF)   #Fluid domain
part = addBox(part,[-nBound*dr,-nBound*dr,Lx+2*nBound*dr,nBound*dr],BOUND,dr,rhoF) #Lower boundary domain
part = addBox(part,[-nBound*dr,0,nBound*dr,Ly],BOUND,dr,rhoF) #Left boundary domain
part = addBox(part,[Lx,0,nBound*dr,Ly],BOUND,dr,rhoF) #Right boundary domain
numberOfNonSolidParticles = len(part)
part = addBox(part,[(Lx-ws)/2, (Ly-ls), ws, ls], MOBILESOLID, dr, rhoS)  #Mobile solid domain
part,partSpace = sortPart(part,posSpace,partSpace,xOrigin,yOrigin,xSize,ySize,lspace,dr)
listNeibSpace= getListNeib(partSpace,neibSpace,listNeibSpace)
numberOfSolidPart = len(part)-numberOfNonSolidParticles<jupyter_output><empty_output><jupyter_text>Affichage du domaine<jupyter_code>%matplotlib notebook
Umax= 1
tabUx = part[:,VEL[0]]
domain = [xOrigin,xMax,yOrigin,yMax,0,Umax]
plotPropertiesWithBound(part,tabUx,r'$U_x$',domain,dr,1)<jupyter_output><empty_output><jupyter_text>Affichage des espaces<jupyter_code>plotSpaces(posSpace,'k',lspace,1)
neibSpace[2][neibSpace[2]>-1]
spacesOutline(posSpace[neibSpace[2][neibSpace[2]>-1]],'r',lspace,1)<jupyter_output><empty_output><jupyter_text>Initialisation de la pression<jupyter_code>def initWaterPressure(part):
    infoTab = part[:,INFO]
    y_i = part[infoTab == FLUID] [:,POS[1]]
    p_i = np.linalg.norm(grav)*rhoF*(Ly/2-y_i)
    rho_i = density(p_i,B,rhoF,gamma)
    part[infoTab == FLUID, RHO] = rho_i
    return part<jupyter_output><empty_output><jupyter_text>Affichage du domaine après initialisation de la pression<jupyter_code>part = initWaterPressure(part)
%matplotlib notebook
Pmax= rhoF*np.linalg.norm(grav)*Ly/2
tabPres = pressure(part[:,RHO],B,rhoF,gamma)
domain = [xOrigin,xMax,yOrigin,yMax,0,Pmax]
plotPropertiesWithBound(part,tabPres,r'Pressure',domain,dr,1)<jupyter_output><empty_output><jupyter_text>Détermination du centre de gravité, initialisation de la vitesse et de l'accéleration<jupyter_code>infoTab = part[:,INFO]
OG = np.array([1/numberOfSolidPart*np.sum(part[infoTab == MOBILESOLID] [:,POS[0]]),1/numberOfSolidPart*np.sum(part[infoTab == MOBILESOLID] [:,POS[1]])])
V_OG = np.array([0,-10])
A_OG = np.array([0,0])<jupyter_output><empty_output><jupyter_text>Gestion des enregistrements<jupyter_code>current_directory = os.getcwd()
case_directory = os.path.join(current_directory, r'Results/Floating_solid_'+time.strftime("%Y%m%d_%H%M%S"))
os.mkdir(case_directory)
data_directory = os.path.join(case_directory,r'Data')
figures_directory = os.path.join(case_directory,r'Figures')
os.mkdir(data_directory)
os.mkdir(figures_directory)<jupyter_output><empty_output><jupyter_text>Gestion du pas de temps et des fréquences de sorties des figures<jupyter_code># Here we specify the output frequencies
dt_figure = 0.005#0.025*e**2*rhoF/mu
t_print = 0
#final time :
t_end = 3
t=0
it=0
im_count=0<jupyter_output><empty_output><jupyter_text>3.3 Simulation SPH<jupyter_code>%matplotlib notebook
while t<t_end:
    #STEP1 : Calcul de la CFL
    dt = CFLConditions(part[:,VEL],h,c0,grav,rhoF,mu)
    #STEP2 : Interpolation des conditions au bord
    part[:,RHO],part[:,VEL] = interpolateBoundary((part[:,INFO]==BOUND),\
                                     part[:,SPID],\
                                     part[:,POS],\
                                     part[:,VEL],\
                                     part[:,RHO],\
                                     listNeibSpace,\
                                 aW,h,m,B,rhoF,gamma,grav,shepardMin)
    #STEP2BIS : Interpolation des conditions au bord sur les mobiles bound
    part[:,RHO], part[:, VEL] = interpolateMobileSolidBoundary((part[:,INFO]==MOBILESOLID),\
                                     part[:,SPID],\
                                     part[:,POS],\
                                     part[:,VEL],\
                                     part[:,RHO],\
                                     listNeibSpace,\
                                 aW,h,m,B,rhoF,gamma,grav,V_OG,A_OG,shepardMin)
    
    #STEP3 : Calcul des forces et des termes de densité
    part[:,FORCES],part[:,DRHODT] = computeForcesMorris((part[:,INFO]==FLUID),\
                                                 part[:,SPID],\
                                                 part[:,POS],\
                                                 part[:,VEL],\
                                                 part[:,RHO],\
                                                 listNeibSpace,\
                                                 aW,h,m,B,rhoF,gamma,grav,mu)
    #STEP4 : Integration en temps
    part[:,POS],part[:,VEL],part[:,RHO] = integrationStep((part[:,INFO]==FLUID),\
                                                       part[:,POS],\
                                                       part[:,VEL],\
                                                       part[:,RHO],\
                                                       part[:,FORCES],\
                                                       part[:,DRHODT],\
                                                       dt)
    #STEP4BIS : integration en temps des particules du solid
    OG, V_OG,A_OG = IntegrateCenterOfMassMovement((part[:,INFO]==MOBILESOLID),\
                                            part[:,SPID],\
                                            part[:,POS],\
                                            part[:,VEL],\
                                            part[:,RHO],\
                                            listNeibSpace,\
                                            aW,h,m,ms,B,rhoF,rhoS,gamma,grav,mu,OG,V_OG,A_OG,part,numberOfSolidPart,dt)

    part[:,POS],part[:,VEL] = MoveSolidParticles((part[:,INFO]==MOBILESOLID),\
                                                part[:,POS],\
                                                part[:,VEL],\
                                                OG,V_OG)
    print("dOG, V_OG :")
    print(OG, V_OG)
    
    #STEP5 : Corriger densité trop basse
    part[:,RHO] = checkDensity(part[:,RHO],rhoMin,rhoMax)
    #STEP6 : Mise à jour des voisins (pas forcément à tous les pas de temps)
    part,partSpace = sortPart(part,posSpace,partSpace,xOrigin,yOrigin,xSize,ySize,lspace,dr)
    listNeibSpace= getListNeib(partSpace,neibSpace,listNeibSpace)
    t +=dt
    it +=1
    if t>=t_print:
        fig = plt.figure(1)
        plt.clf()
        plt.title(r'$t\nu/e^2 = %2.2f$'%(t*mu/rhoF/e**2))
        velMagn = (part[:,VEL[0]]*part[:,VEL[0]]+part[:,VEL[1]]*part[:,VEL[1]])**.5
        domain = [xOrigin,xMax,yOrigin,yMax,0,1]
        plotPropertiesWithBound(part,velMagn,r'$U_{mag}$',domain,dr,1)
        figname = os.path.join(figures_directory,r'vel_%06d.png'%im_count)
        fig.savefig(figname,bbox_inches='tight')
        fig.canvas.draw()
        plt.pause(0.01)        
        plt.pause(0.01)
        im_count = im_count+1
        t_print +=dt_figure<jupyter_output>dOG, V_OG :
[ 0.         -0.00400157] [  0.       -10.003924]
{'fraction_non_alphanumeric': 0.1266476853316461, 'fraction_numerical': 0.017715912685858905, 'mean_word_length': 2.604713036868111, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 54, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 54, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '431237', 'license': 'MIT', 'n_tokens_mistral': 3338, 'n_tokens_neox': 3273, 'n_words': 607, 'path': '3Floating_solid.ipynb', 'repo_name': 'rubenpersicot/mouving_boundaries', 'size': '9515'}