Since we have the Navier-Stoke equation for simulating fluid, we can break it down to simulate the movement of individual fluid particles, for a particle , we have
Smoothed-particle hydrodynamics (SPH) is a computational method used for simulating the mechanics of continuum media, such as solid mechanics and fluid flows. It can be used to solve the equation (281) for any quantity of one particle. It is a meshfree Lagrangian method(where the co-ordinates move with the fluid), and the resolution of the method can easilybe adjusted with respect to variables such as density.
The key assumption of the SPH method is that calculations are based on a weighted average (by ) of field values. will give more strength to points that are closer and make points that are further away have a weaker influence. At points that are more than a certain distance away, will become zero, which means those points would not have any influence at all.
importnumpyasnpimportmatplotlib.pyplotaspltN=1000# number of points to plot (ignore)h=16.0ri_minus_rj=np.linspace(-h,h,N)defW(ri_minus_rj,h):return315/(64*np.pi*h**9)*(h**2-ri_minus_rj**2)**3plt.figure()plt.plot(ri_minus_rj,W(ri_minus_rj,h))plt.xlabel('$r_i$ - $r_j$')plt.ylabel('W($r_i$-$r_j$,h)')plt.title('Smoothing Kernel')plt.show()
The general expression for any quantity is given below:
where is the location of a neighboring point and is the location of our point; is the interaction radius, the points that are more than away will stop interacting with our point.
So for density, the discrete sum of approximation is
But this equation does not obey the conservation of momentum because it does not have the right symmetry (force derived from this equation felt by particle A from B does not equal to the force felt by B from A). To conserve quantities, we can apply Lagrangian, then the equation becomes:
So instead of the approximation from equation (284), we could have the following expression by plugging the and into equation (285), keeping the momentum a constant:
# Open Anaconda Prompt (anaconda 3), then type#pip install ipyvolume#pip install bqplot#jupyter nbextension enable –py bqplotimportipyvolumeasipvimportbqplot.scalesimportnumpyasnpfromnumpyimportrandomimportipywidgetsaswidgets# https://www.youtube.com/watch?v=hOKa8klJPyoimportmatplotlib.cmascm# for coloring temperature
# some constantsG=np.array([0.0,-9.8,0.0]);# external (gravitational) forcesREST_DENS=1000.0;# rest density of liquidGAS_CONST=2000.0;# const for equation of stateH=16.0;# kernel radiusHSQ=H*H;# radius^2 for optimizationT0=1.0# initial temperature of the fluidT_bottom=1.5# bottom heating temperatureT_ceiling=0.5# upper cooling temperaturealpha_v=0.0002# thermal expansion coefficientk=0.0003# thermal conductivity of the particlek_wall=100.0# thermal conductivity of the wallMASS=65.0;# assume all particles have the same massVISC=250.0;# viscosity constantdt=0.0008;# integration timestep# smoothing kernels defined in Müller and their gradientsPOLY6=315.0/(65.0*np.pi*H**9);SPIKY_GRAD=-45.0/(np.pi*H**6);VISC_LAP=45.0/(np.pi*H**6);# simulation parametersEPS=H;# boundary epsilonBOUND_DAMPING_NORMAL=-0.0;# particle does not bounce back verticallyBOUND_DAMPING_HORIZONTAL=0.0;# no slip-boundary# solver dataFluid_Particles_list=[];X=[[],[],[]]# positions for all particlesTs=[]# temperature of all particles# interactionMAX_PARTICLES=3000;# rendering projection parametersWINDOW_WIDTH=425;WINDOW_HEIGHT=150;WINDOW_DEPTH=300VIEW_WIDTH=1.5*WINDOW_WIDTH;VIEW_HEIGHT=1.5*WINDOW_HEIGHT;VIEW_DEPTH=1.5*WINDOW_DEPTH;
# particle data structure# stores position, velocity, and force for integration# stores density (rho) and pressure values for SPHclassParticle:def__init__(self,X,V,F,rho,p,T):self.X=X# x position vectorself.V=V# velocity vectorself.F=F# force vectorself.rho=0.0# density of the fluid at the position of the particleself.p=p# pressure of the fluid at the position of the particleself.T=T0# temperature of the fluid at the position of the particle
defInitSPH():globalXglobalfigglobalscatterglobalipvglobalTs# for 3d simulation, un-comment the following 2 lines and indent the codes below#for zi in range(int(2*EPS),int(WINDOW_DEPTH-EPS*2.0),int(H)):#zi += Hforyiinrange(int(1*EPS),int(WINDOW_HEIGHT-EPS*1.0),int(H*1/16)):yi+=H#+random.rand()*H/1000forxiinrange(0,int(WINDOW_WIDTH),int(H)):iflen(Fluid_Particles_list)<MAX_PARTICLES:xi+=H#+random.rand()*H/1000Vx=0.0# gives the particle a random tiny initial velocity in the x-directionVy=0.0# gives the particle a random tiny initial velocity in the y-directionVz=0.0V=np.array([Vx,Vy,Vz])Fx=0.0# gives the particle a random tiny initial force in the x-directionFy=0.0# gives the particle a random tiny initial force in the y-directionFz=0.0F=np.array([Fx,Fy,Fz])rho=0.0p=0.0T=T0# for 3d simulation, delete the line belowzi=2*EPS# make position in z-direction the same because its a 2d problemXi=np.array([xi+np.random.random_sample(),yi+np.random.random_sample(),zi])# + np.random.random_sample(); for 3d simulationFluid_Particles_list.append(Particle(Xi,V,F,rho,p,T))# for 3d simulation, indent the codes aboveforpiinFluid_Particles_list:X[0].append(pi.X[0])X[1].append(pi.X[1])X[2].append(pi.X[2])Ts.append(np.array([pi.T,pi.T,pi.T]))x=np.array(X[0])y=np.array(X[1])z=np.array(X[2])# the scale does not seem to work...scales={'x':bqplot.scales.LinearScale(min=0,max=VIEW_WIDTH),'y':bqplot.scales.LinearScale(min=0,max=VIEW_HEIGHT),'z':bqplot.scales.LinearScale(min=0,max=VIEW_DEPTH),}fig=ipv.figure(scales=scales)Ts=np.array(Ts)Ts-=Ts.min()Ts/=Ts.max()colors=Tsscatter=ipv.scatter(x,y,z,marker='sphere',color=colors)scatter.size=3ipv.show()print("initializing mantle convection with ")print(len(x))print("particles")
In order to compute the at every point for equation (223), we need to calculate the following equations respectively:
where is a constant, is the resting density (density of fluid at equilibrium).
defComputeDensityPressureTemperature():forpiinFluid_Particles_list:pi.rho=0.0dT=0.0# the change in temperatureforpjinFluid_Particles_list:h=np.linalg.norm(pi.X-pj.X)if(h<H):pi.T+=(pj.T-pi.T)*k*dt*POLY6*(HSQ-h**2)**3;pi.rho+=MASS*POLY6*(HSQ-h**2)**3;pi.rho-=pi.rho*alpha_v*(pi.T-T0)pi.p=GAS_CONST*(pi.rho-REST_DENS)*pi.T;
Then we can compute the pressure and viscosity term
After equation (290) to (293) are computed, we can put them together into N-S equation to obtain acceleration.
defComputeForces():forpiinFluid_Particles_list:fpress=0.0fvisc=0.0forpjinFluid_Particles_list:d=pj.X-pi.Xh=np.linalg.norm(d)r=d/hif(pi==pj):continue;if(h<H):# compute pressure force contributionfpress+=-r*MASS*(pi.p+pj.p)/(2*pj.rho)*SPIKY_GRAD*(H-h)**2# compute viscosity force contributionfvisc+=VISC*MASS*(pj.V-pi.V)/(pj.rho)*VISC_LAP*(H-h);fgrav=G*pi.rhopi.F=fpress+fvisc+fgrav
And then we can use a simple time stepping numerical integration scheme to advance the velocity and position of particles.
defIntegrate():forpinFluid_Particles_list:# forward Euler integrationp.V+=dt*(p.F)/p.rho;p.V[2]=0.0p.X+=dt*p.V;if(p.X[0]-EPS<0.0):# left boundary conditionp.V[0]*=1.0#BOUND_DAMPING_NORMALp.X[0]=WINDOW_WIDTH-EPS# coming out from the right boundaryp.V[1]*=BOUND_DAMPING_HORIZONTAL#p.Vz *= BOUND_DAMPING_HORIZONTALif(p.X[0]+EPS>WINDOW_WIDTH):# right boundary conditionp.V[0]*=1.0#BOUND_DAMPING_NORMALp.X[0]=EPS# coming out from the left boundaryp.V[1]*=BOUND_DAMPING_HORIZONTAL#p.Vz *= BOUND_DAMPING_HORIZONTALif(p.X[1]-EPS<0.0):# bottom boundary conditionp.V[1]*=BOUND_DAMPING_NORMALp.X[1]=EPSp.V[0]*=BOUND_DAMPING_HORIZONTALp.T+=(T_bottom-p.T)*k_wall*dt# heating from the bottom#p.Vz *= BOUND_DAMPING_HORIZONTALif(p.X[1]+EPS>WINDOW_HEIGHT):# upper boundary conditionp.V[1]*=BOUND_DAMPING_NORMALp.X[1]=WINDOW_HEIGHT-EPSp.V[0]*=BOUND_DAMPING_HORIZONTALp.T+=(T_ceiling-p.T)*k_wall*dt# cooling from the ceiling#p.Vz *= BOUND_DAMPING_HORIZONTAL# only for 3D boundary condition (in the z-direction)if(p.X[2]-EPS<0.0):# front boundary conditionp.V[2]*=BOUND_DAMPING_NORMALp.X[2]=EPS#p.Vx *= BOUND_DAMPING_HORIZONTAL#p.Vy *= BOUND_DAMPING_HORIZONTALif(p.X[2]+EPS>WINDOW_DEPTH):# back boundary conditionp.V[2]*=BOUND_DAMPING_NORMALp.X[2]=WINDOW_DEPTH-EPS#p.Vx *= BOUND_DAMPING_HORIZONTAL#p.Vy *= BOUND_DAMPING_HORIZONTAL
C:\Users\JianouJiang\anaconda3\lib\site-packages\ipykernel_launcher.py:58: RuntimeWarning: invalid value encountered in true_divide
initializing mantle convection with
3000
particles
# Run this code to update the mantle convectionTotalTime=1000iterations=int(TotalTime/dt)print("start iterating")fortimeinrange(iterations):Update()print("iterations")print(time)
start iterating
start updating
C:\Users\JianouJiang\anaconda3\lib\site-packages\ipykernel_launcher.py:12: RuntimeWarning: invalid value encountered in true_divide
if sys.path[0] == '':
C:\Users\JianouJiang\anaconda3\lib\site-packages\ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in true_divide
after removing the cwd from sys.path.
The simulation takes hours because each iteration is time consuming. The figure below shows the result after many iterations. One problem of this simulation is the boundary condition, where particles are sticked to the wall when they are in contact with it. This means that after many iterations, all particles will be sticked to the wall and the fluid dynamics would be poorly simulated.
Figure 1: Two-dimensional cellular convection in a fluid layer heated from below simulation shows similar character as the figure 6.38 in the book called , [Turcotte_D.L.,_Schubert_G.].
One way of resolving the problem is called the ghost particle. These virtual particle (ghost particle)will interact with the fluid particles through both pressure and shear stress. If no slip boundariesare required, then the reflected particles are given a velocity that has the same magnitude, butopposite direction to that of the matching particle, while for full slip boundary conditions thereflected particle has the same velocity.