Stokesian Dynamic Simulation of Sedimenting Spheres in a Newtonian Fluid
Authors: Housam Binous and Ahmed Bellagi
Stokesian Dynamic Simulation of Sedimenting Spheres in a Newtonian Fluid
Authors: Housam Binous and Ahmed Bellagi
Authors: Housam Binous and Ahmed Bellagi
Stokesian Dynamics, a method developed by Brady and Bossis in the 80s, simulate the 3D motion of hydrodynamically interacting spheres at low Reynolds numbers . This program was originally written Professor Ron Phillips during his PhD thesis at MIT in 1989. The code was adapted from Fortran 77 to a Mathematica notebook by Housam Binous . Using Euler' s integration scheme was one simplification of the original code that was performed by Housam Binous . The code calls a function NEWTON2, which needs three functions TEMPM, FTMOB and SSI.
Several simulations are presented here . Two sphere' s falling at low Reynolds numbers can have lateral motion during their fall . Two spheres can also rotate, as they sediment, for certain initial configurations . Four vertically aligned spheres, initially at equal distance from each others, will tend to drift apart . Four and six spheres sedimenting will show periodic arrangements that are very peculiar . This Stokesian Dynamics method was extended to non-Newtonian fluids by Binous and Phillips in 1999 by using FENE dumbbell suspensions.
Several simulations are presented here . Two sphere' s falling at low Reynolds numbers can have lateral motion during their fall . Two spheres can also rotate, as they sediment, for certain initial configurations . Four vertically aligned spheres, initially at equal distance from each others, will tend to drift apart . Four and six spheres sedimenting will show periodic arrangements that are very peculiar . This Stokesian Dynamics method was extended to non-Newtonian fluids by Binous and Phillips in 1999 by using FENE dumbbell suspensions.
Acknowledgement :
Permission was given by Professor Phillips to post this work in the WEB.
The code uses several functions called NEWTON2, TEMPM, FTMOB and SSI as well as a data file called data.m.
1/ NEWTON2 computes far-field and adds to it lubrication; it then computes the velocities of the spheres using the mobility matrix and the applied forces.
2/ TEMPM computes far-field mobility.
3/ SSI creates two-body resistance matrix. To do so, it uses interpolation of values found data.m as well as analytical expressions. The program then subtracts the two-body far-field mobility from the two-body resistance matrix. The result is the lubrication matrix.
4/ FTMOB creates the two-body far-field mobility.
Permission was given by Professor Phillips to post this work in the WEB.
The code uses several functions called NEWTON2, TEMPM, FTMOB and SSI as well as a data file called data.m.
1/ NEWTON2 computes far-field and adds to it lubrication; it then computes the velocities of the spheres using the mobility matrix and the applied forces.
2/ TEMPM computes far-field mobility.
3/ SSI creates two-body resistance matrix. To do so, it uses interpolation of values found data.m as well as analytical expressions. The program then subtracts the two-body far-field mobility from the two-body resistance matrix. The result is the lubrication matrix.
4/ FTMOB creates the two-body far-field mobility.
Two spheres sedimenting in a Newtonian fluid. Since the line of centers of the spheres has a y-component, they present lateral motion as they fall.
In[]:=
FU={0,0,-10,0,0,-10,0,0,0,0,0,0};
In[]:=
X[1]=0;Y[1]=1;Z[1]=0;X[2]=0;Y[2]=0;Z[2]=3;
In[]:=
dt=0.1;
In[]:=
NSPHR=2;NDIM=6 NSPHR;Ng=3 NSPHR;
In[]:=
In[]:=
X[1,0]=X[1];Y[1,0]=Y[1];Z[1,0]=Z[1];
X[2,0]=X[2];Y[2,0]=Y[2];Z[2,0]=Z[2];
X[2,0]=X[2];Y[2,0]=Y[2];Z[2,0]=Z[2];
In[]:=
t=1;
In[]:=
For[t=1,t<100,{NEWTON2[{X[1,t-1],X[2,t-1]},{Y[1,t-1],Y[2,t-1]},{Z[1,t-1],Z[2,t-1]},t];
X[1,t]=X[1,t-1]+UNEW[1]*dt;Y[1,t]=Y[1,t-1]+UNEW[2]*dt;Z[1,t]=Z[1,t-1]+UNEW[3]*dt;
X[2,t]=X[2,t-1]+UNEW[4]*dt;Y[2,t]=Y[2,t-1]+UNEW[5]*dt;Z[2,t]=Z[2,t-1]+UNEW[6]*dt;Print[t];t++}]
X[1,t]=X[1,t-1]+UNEW[1]*dt;Y[1,t]=Y[1,t-1]+UNEW[2]*dt;Z[1,t]=Z[1,t-1]+UNEW[3]*dt;
X[2,t]=X[2,t-1]+UNEW[4]*dt;Y[2,t]=Y[2,t-1]+UNEW[5]*dt;Z[2,t]=Z[2,t-1]+UNEW[6]*dt;Print[t];t++}]
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
In[]:=
Table[{X[1,t],Y[1,t],Z[1,t]},{t,0,99}]
Out[]=
{{0,1,0},{0.,1.04937,-1.39986},{0.,1.09873,-2.79973},{0.,1.1481,-4.19959},{0.,1.19747,-5.59946},{0.,1.24684,-6.99932},{0.,1.2962,-8.39919},{0.,1.34557,-9.79905},{0.,1.39494,-11.1989},{0.,1.4443,-12.5988},{0.,1.49367,-13.9986},{0.,1.54304,-15.3985},{0.,1.59241,-16.7984},{0.,1.64177,-18.1982},{0.,1.69114,-19.5981},{0.,1.74051,-20.998},{0.,1.78988,-22.3978},{0.,1.83924,-23.7977},{0.,1.88861,-25.1976},{0.,1.93798,-26.5974},{0.,1.98734,-27.9973},{0.,2.03671,-29.3972},{0.,2.08608,-30.797},{0.,2.13545,-32.1969},{0.,2.18481,-33.5968},{0.,2.23418,-34.9966},{0.,2.28355,-36.3965},{0.,2.33291,-37.7964},{0.,2.38228,-39.1962},{0.,2.43165,-40.5961},{0.,2.48102,-41.9959},{0.,2.53038,-43.3958},{0.,2.57975,-44.7957},{0.,2.62912,-46.1955},{0.,2.67848,-47.5954},{0.,2.72785,-48.9953},{0.,2.77722,-50.3951},{0.,2.82659,-51.795},{0.,2.87595,-53.1949},{0.,2.92532,-54.5947},{0.,2.97469,-55.9946},{0.,3.02406,-57.3945},{0.,3.07342,-58.7943},{0.,3.12279,-60.1942},{0.,3.17216,-61.5941},{0.,3.22152,-62.9939},{0.,3.27089,-64.3938},{0.,3.32026,-65.7937},{0.,3.36963,-67.1935},{0.,3.41899,-68.5934},{0.,3.46836,-69.9932},{0.,3.51773,-71.3931},{0.,3.56709,-72.793},{0.,3.61646,-74.1928},{0.,3.66583,-75.5927},{0.,3.7152,-76.9926},{0.,3.76456,-78.3924},{0.,3.81393,-79.7923},{0.,3.8633,-81.1922},{0.,3.91266,-82.592},{0.,3.96203,-83.9919},{0.,4.0114,-85.3918},{0.,4.06077,-86.7916},{0.,4.11013,-88.1915},{0.,4.1595,-89.5914},{0.,4.20887,-90.9912},{0.,4.25823,-92.3911},{0.,4.3076,-93.791},{0.,4.35697,-95.1908},{0.,4.40634,-96.5907},{0.,4.4557,-97.9905},{0.,4.50507,-99.3904},{0.,4.55444,-100.79},{0.,4.60381,-102.19},{0.,4.65317,-103.59},{0.,4.70254,-104.99},{0.,4.75191,-106.39},{0.,4.80127,-107.79},{0.,4.85064,-109.189},{0.,4.90001,-110.589},{0.,4.94938,-111.989},{0.,4.99874,-113.389},{0.,5.04811,-114.789},{0.,5.09748,-116.189},{0.,5.14684,-117.589},{0.,5.19621,-118.989},{0.,5.24558,-120.388},{0.,5.29495,-121.788},{0.,5.34431,-123.188},{0.,5.39368,-124.588},{0.,5.44305,-125.988},{0.,5.49241,-127.388},{0.,5.54178,-128.788},{0.,5.59115,-130.187},{0.,5.64052,-131.587},{0.,5.68988,-132.987},{0.,5.73925,-134.387},{0.,5.78862,-135.787},{0.,5.83799,-137.187},{0.,5.88735,-138.587}}
In[]:=
b1=ListPointPlot3D[{Table[{X[1,t],Y[1,t],Z[1,t]},{t,0,99}],Table[{X[2,t],Y[2,t],Z[2,t]},{t,0,99}]},BoxRatios->{1,1,1}];
In[]:=
Show[b1,ViewPoint->{1,1,1}]
Four spheres sedimenting in a Newtonian fluid. Initially the spheres are at equal distance from each others. Due to hydrodynamic interaction, the last sphere lags behind after a certain time.
In[]:=
FU={0,0,-3,0,0,-3,0,0,-3,0,0,-3,0,0,0,0,0,0,0,0,0,0,0,0};
In[]:=
X[1]=0;Y[1]=0;Z[1]=0;X[2]=0;Y[2]=0;Z[2]=2.5;X[3]=0;Y[3]=0;Z[3]=5;X[4]=0;Y[4]=0;Z[4]=7.5;
In[]:=
Show[Graphics3D[{{RGBColor[1,0,1],PointSize[0.05],Point[{X[1],Y[1],Z[1]}]},
{RGBColor[1,0,1],PointSize[0.05],Point[{X[2],Y[2],Z[2]}]},{RGBColor[1,0,1],PointSize[0.05],Point[{X[3],Y[3],Z[3]}]},
{RGBColor[1,0,1],PointSize[0.05],Point[{X[4],Y[4],Z[4]}]}}],BoxRatios->{1,1,1},AxesTrue]
{RGBColor[1,0,1],PointSize[0.05],Point[{X[2],Y[2],Z[2]}]},{RGBColor[1,0,1],PointSize[0.05],Point[{X[3],Y[3],Z[3]}]},
{RGBColor[1,0,1],PointSize[0.05],Point[{X[4],Y[4],Z[4]}]}}],BoxRatios->{1,1,1},AxesTrue]
In[]:=
dt=0.01;
In[]:=
NSPHR=4;NDIM=6 NSPHR;Ng=3 NSPHR;
In[]:=
X[1,0]=X[1];Y[1,0]=Y[1];Z[1,0]=Z[1];
X[2,0]=X[2];Y[2,0]=Y[2];Z[2,0]=Z[2];
X[3,0]=X[3];Y[3,0]=Y[3];Z[3,0]=Z[3];
X[4,0]=X[4];Y[4,0]=Y[4];Z[4,0]=Z[4];
X[2,0]=X[2];Y[2,0]=Y[2];Z[2,0]=Z[2];
X[3,0]=X[3];Y[3,0]=Y[3];Z[3,0]=Z[3];
X[4,0]=X[4];Y[4,0]=Y[4];Z[4,0]=Z[4];
In[]:=
For[t=1,t<500,{NEWTON2[{X[1,t-1],X[2,t-1],X[3,t-1],X[4,t-1]},
{Y[1,t-1],Y[2,t-1],Y[3,t-1],Y[4,t-1]},{Z[1,t-1],Z[2,t-1],Z[3,t-1],Z[4,t-1]},t];
X[1,t]=X[1,t-1]+UNEW[1]*dt;Y[1,t]=Y[1,t-1]+UNEW[2]*dt;Z[1,t]=Z[1,t-1]+UNEW[3]*dt;
X[2,t]=X[2,t-1]+UNEW[4]*dt;Y[2,t]=Y[2,t-1]+UNEW[5]*dt;Z[2,t]=Z[2,t-1]+UNEW[6]*dt;
X[3,t]=X[3,t-1]+UNEW[7]*dt;Y[3,t]=Y[3,t-1]+UNEW[8]*dt;Z[3,t]=Z[3,t-1]+UNEW[9]*dt;
X[4,t]=X[4,t-1]+UNEW[10]*dt;Y[4,t]=Y[4,t-1]+UNEW[11]*dt;Z[4,t]=Z[4,t-1]+UNEW[12]*dt;
Print[t];t++}]
{Y[1,t-1],Y[2,t-1],Y[3,t-1],Y[4,t-1]},{Z[1,t-1],Z[2,t-1],Z[3,t-1],Z[4,t-1]},t];
X[1,t]=X[1,t-1]+UNEW[1]*dt;Y[1,t]=Y[1,t-1]+UNEW[2]*dt;Z[1,t]=Z[1,t-1]+UNEW[3]*dt;
X[2,t]=X[2,t-1]+UNEW[4]*dt;Y[2,t]=Y[2,t-1]+UNEW[5]*dt;Z[2,t]=Z[2,t-1]+UNEW[6]*dt;
X[3,t]=X[3,t-1]+UNEW[7]*dt;Y[3,t]=Y[3,t-1]+UNEW[8]*dt;Z[3,t]=Z[3,t-1]+UNEW[9]*dt;
X[4,t]=X[4,t-1]+UNEW[10]*dt;Y[4,t]=Y[4,t-1]+UNEW[11]*dt;Z[4,t]=Z[4,t-1]+UNEW[12]*dt;
Print[t];t++}]
In[]:=
Show[Graphics3D[{{RGBColor[1,0,1],PointSize[0.05],Point[{X[1,499],Y[1,499],Z[1,499]}]},
{RGBColor[1,0,1],PointSize[0.05],Point[{X[2,499],Y[2,499],Z[2,499]}]},
{RGBColor[1,0,1],PointSize[0.05],Point[{X[3,499],Y[3,499],Z[3,499]}]},
{RGBColor[1,0,1],PointSize[0.05],Point[{X[4,499],Y[4,499],Z[4,499]}]}}],BoxRatios->{1,1,1},AxesTrue]
{RGBColor[1,0,1],PointSize[0.05],Point[{X[2,499],Y[2,499],Z[2,499]}]},
{RGBColor[1,0,1],PointSize[0.05],Point[{X[3,499],Y[3,499],Z[3,499]}]},
{RGBColor[1,0,1],PointSize[0.05],Point[{X[4,499],Y[4,499],Z[4,499]}]}}],BoxRatios->{1,1,1},AxesTrue]
Four spheres sedimenting in a Newtonian fluid. Depending on their initial arrangement, the spheres' trajectories can be very peculiar and show a periodic structure.
In[]:=
FU={0,0,-30,0,0,-30,0,0,-30,0,0,-30,0,0,0,0,0,0,0,0,0,0,0,0};
In[]:=
X[1]=1.1;Y[1]=0;Z[1]=1.5;X[2]=-1.1;Y[2]=0;Z[2]=1.5;X[3]=-1.1;Y[3]=0;Z[3]=-0.7;X[4]=1.1;Y[4]=0;Z[4]=-0.7;
In[]:=
Show[Graphics3D[{{RGBColor[1,0,1],PointSize[0.05],Point[{X[1],Y[1],Z[1]}]},{RGBColor[1,0,1],PointSize[0.05],Point[{X[2],Y[2],Z[2]}]},
{RGBColor[1,0,1],PointSize[0.05],Point[{X[3],Y[3],Z[3]}]},{RGBColor[1,0,1],PointSize[0.05],Point[{X[4],Y[4],Z[4]}]}},BoxRatios->{1,1,1},AxesTrue]]
{RGBColor[1,0,1],PointSize[0.05],Point[{X[3],Y[3],Z[3]}]},{RGBColor[1,0,1],PointSize[0.05],Point[{X[4],Y[4],Z[4]}]}},BoxRatios->{1,1,1},AxesTrue]]
In[]:=
dt=0.01;
In[]:=
NSPHR=4;NDIM=6 NSPHR;Ng=3 NSPHR;
In[]:=
X[1,0]=X[1];Y[1,0]=Y[1];Z[1,0]=Z[1];
X[2,0]=X[2];Y[2,0]=Y[2];Z[2,0]=Z[2];
X[3,0]=X[3];Y[3,0]=Y[3];Z[3,0]=Z[3];
X[4,0]=X[4];Y[4,0]=Y[4];Z[4,0]=Z[4];
X[2,0]=X[2];Y[2,0]=Y[2];Z[2,0]=Z[2];
X[3,0]=X[3];Y[3,0]=Y[3];Z[3,0]=Z[3];
X[4,0]=X[4];Y[4,0]=Y[4];Z[4,0]=Z[4];
In[]:=
For[t=1,t<641,{NEWTON2[{X[1,t-1],X[2,t-1],X[3,t-1],X[4,t-1]},
{Y[1,t-1],Y[2,t-1],Y[3,t-1],Y[4,t-1]},{Z[1,t-1],Z[2,t-1],Z[3,t-1],Z[4,t-1]},t];
X[1,t]=X[1,t-1]+UNEW[1]*dt;Y[1,t]=Y[1,t-1]+UNEW[2]*dt;Z[1,t]=Z[1,t-1]+UNEW[3]*dt;X[2,t]=X[2,t-1]+UNEW[4]*dt;Y[2,t]=Y[2,t-1]+UNEW[5]*dt;Z[2,t]=Z[2,t-1]+UNEW[6]*dt;X[3,t]=X[3,t-1]+UNEW[7]*dt;Y[3,t]=Y[3,t-1]+UNEW[8]*dt;Z[3,t]=Z[3,t-1]+UNEW[9]*dt;X[4,t]=X[4,t-1]+UNEW[10]*dt;Y[4,t]=Y[4,t-1]+UNEW[11]*dt;Z[4,t]=Z[4,t-1]+UNEW[12]*dt;Print[t];t++}]
{Y[1,t-1],Y[2,t-1],Y[3,t-1],Y[4,t-1]},{Z[1,t-1],Z[2,t-1],Z[3,t-1],Z[4,t-1]},t];
X[1,t]=X[1,t-1]+UNEW[1]*dt;Y[1,t]=Y[1,t-1]+UNEW[2]*dt;Z[1,t]=Z[1,t-1]+UNEW[3]*dt;X[2,t]=X[2,t-1]+UNEW[4]*dt;Y[2,t]=Y[2,t-1]+UNEW[5]*dt;Z[2,t]=Z[2,t-1]+UNEW[6]*dt;X[3,t]=X[3,t-1]+UNEW[7]*dt;Y[3,t]=Y[3,t-1]+UNEW[8]*dt;Z[3,t]=Z[3,t-1]+UNEW[9]*dt;X[4,t]=X[4,t-1]+UNEW[10]*dt;Y[4,t]=Y[4,t-1]+UNEW[11]*dt;Z[4,t]=Z[4,t-1]+UNEW[12]*dt;Print[t];t++}]
In[]:=
c1=ListPointPlot3D[{Table[{X[1,t],Y[1,t],Z[1,t]},{t,0,640}],Table[{X[2,t],Y[2,t],Z[2,t]},{t,0,640}],
Table[{X[3,t],Y[3,t],Z[3,t]},{t,0,640}],Table[{X[4,t],Y[4,t],Z[4,t]},{t,0,640}]},BoxRatios->{1,1,1}]
Table[{X[3,t],Y[3,t],Z[3,t]},{t,0,640}],Table[{X[4,t],Y[4,t],Z[4,t]},{t,0,640}]},BoxRatios->{1,1,1}]
Six spheres sedimenting in a Newtonian fluid. The same periodic arrangement of the spheres occurs in this case.
In[]:=
FU={0,0,-30,0,0,-30,0,0,-30,0,0,-30,0,0,-30,0,0,-30,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
In[]:=
X[1]=3 Cos[Pi/3.];Y[1]=3 Sin[Pi/3.];Z[1]=1;
X[2]=3 Cos[2 Pi/3.];Y[2]=3 Sin[2 Pi/3.];Z[2]=0;
X[3]=3 Cos[Pi];Y[3]=3 Sin[Pi];Z[3]=1;
X[4]=3 Cos[4 Pi/3.];Y[4]=3 Sin[4 Pi/3.];Z[4]=0;
X[5]=3 Cos[5 Pi/3.];Y[5]=3 Sin[5 Pi/3.];Z[5]=1;
X[6]=3 Cos[2 Pi];Y[6]=3 Sin[2 Pi];Z[6]=0;
X[2]=3 Cos[2 Pi/3.];Y[2]=3 Sin[2 Pi/3.];Z[2]=0;
X[3]=3 Cos[Pi];Y[3]=3 Sin[Pi];Z[3]=1;
X[4]=3 Cos[4 Pi/3.];Y[4]=3 Sin[4 Pi/3.];Z[4]=0;
X[5]=3 Cos[5 Pi/3.];Y[5]=3 Sin[5 Pi/3.];Z[5]=1;
X[6]=3 Cos[2 Pi];Y[6]=3 Sin[2 Pi];Z[6]=0;
In[]:=
Show[Graphics3D[{{RGBColor[0,0,1],PointSize[0.05],Point[{X[1],Y[1],Z[1]}]},{RGBColor[1,0,1],PointSize[0.05],Point[{X[2],Y[2],Z[2]}]},
{RGBColor[0,0,1],PointSize[0.05],Point[{X[3],Y[3],Z[3]}]},{RGBColor[1,0,1],PointSize[0.05],Point[{X[4],Y[4],Z[4]}]},
{RGBColor[0,0,1],PointSize[0.05],Point[{X[5],Y[5],Z[5]}]},{RGBColor[1,0,1],PointSize[0.05],Point[{X[6],Y[6],Z[6]}]}}],BoxRatios->{1,1,1},AxesTrue]
{RGBColor[0,0,1],PointSize[0.05],Point[{X[3],Y[3],Z[3]}]},{RGBColor[1,0,1],PointSize[0.05],Point[{X[4],Y[4],Z[4]}]},
{RGBColor[0,0,1],PointSize[0.05],Point[{X[5],Y[5],Z[5]}]},{RGBColor[1,0,1],PointSize[0.05],Point[{X[6],Y[6],Z[6]}]}}],BoxRatios->{1,1,1},AxesTrue]
In[]:=
dt=0.01;
In[]:=
NSPHR=6;NDIM=6 NSPHR;Ng=3 NSPHR;
In[]:=
X[1,0]=X[1];Y[1,0]=Y[1];Z[1,0]=Z[1];
X[2,0]=X[2];Y[2,0]=Y[2];Z[2,0]=Z[2];
X[3,0]=X[3];Y[3,0]=Y[3];Z[3,0]=Z[3];
X[4,0]=X[4];Y[4,0]=Y[4];Z[4,0]=Z[4];
X[5,0]=X[5];Y[5,0]=Y[5];Z[5,0]=Z[5];
X[6,0]=X[6];Y[6,0]=Y[6];Z[6,0]=Z[6];
X[2,0]=X[2];Y[2,0]=Y[2];Z[2,0]=Z[2];
X[3,0]=X[3];Y[3,0]=Y[3];Z[3,0]=Z[3];
X[4,0]=X[4];Y[4,0]=Y[4];Z[4,0]=Z[4];
X[5,0]=X[5];Y[5,0]=Y[5];Z[5,0]=Z[5];
X[6,0]=X[6];Y[6,0]=Y[6];Z[6,0]=Z[6];
In[]:=
For[t=1,t<310,{NEWTON2[{X[1,t-1],X[2,t-1],X[3,t-1],X[4,t-1],X[5,t-1],X[6,t-1]},
{Y[1,t-1],Y[2,t-1],Y[3,t-1],Y[4,t-1],Y[5,t-1],Y[6,t-1]},{Z[1,t-1],Z[2,t-1],Z[3,t-1],Z[4,t-1],Z[5,t-1],Z[6,t-1]},t];
X[1,t]=X[1,t-1]+UNEW[1]*dt;Y[1,t]=Y[1,t-1]+UNEW[2]*dt;Z[1,t]=Z[1,t-1]+UNEW[3]*dt;X[2,t]=X[2,t-1]+UNEW[4]*dt;Y[2,t]=Y[2,t-1]+UNEW[5]*dt;Z[2,t]=Z[2,t-1]+UNEW[6]*dt;X[3,t]=X[3,t-1]+UNEW[7]*dt;Y[3,t]=Y[3,t-1]+UNEW[8]*dt;Z[3,t]=Z[3,t-1]+UNEW[9]*dt;X[4,t]=X[4,t-1]+UNEW[10]*dt;Y[4,t]=Y[4,t-1]+UNEW[11]*dt;Z[4,t]=Z[4,t-1]+UNEW[12]*dt;X[5,t]=X[5,t-1]+UNEW[13]*dt;Y[5,t]=Y[5,t-1]+UNEW[14]*dt;Z[5,t]=Z[5,t-1]+UNEW[15]*dt;X[6,t]=X[6,t-1]+UNEW[16]*dt;Y[6,t]=Y[6,t-1]+UNEW[17]*dt;Z[6,t]=Z[6,t-1]+UNEW[18]*dt;Print[t];t++}]
{Y[1,t-1],Y[2,t-1],Y[3,t-1],Y[4,t-1],Y[5,t-1],Y[6,t-1]},{Z[1,t-1],Z[2,t-1],Z[3,t-1],Z[4,t-1],Z[5,t-1],Z[6,t-1]},t];
X[1,t]=X[1,t-1]+UNEW[1]*dt;Y[1,t]=Y[1,t-1]+UNEW[2]*dt;Z[1,t]=Z[1,t-1]+UNEW[3]*dt;X[2,t]=X[2,t-1]+UNEW[4]*dt;Y[2,t]=Y[2,t-1]+UNEW[5]*dt;Z[2,t]=Z[2,t-1]+UNEW[6]*dt;X[3,t]=X[3,t-1]+UNEW[7]*dt;Y[3,t]=Y[3,t-1]+UNEW[8]*dt;Z[3,t]=Z[3,t-1]+UNEW[9]*dt;X[4,t]=X[4,t-1]+UNEW[10]*dt;Y[4,t]=Y[4,t-1]+UNEW[11]*dt;Z[4,t]=Z[4,t-1]+UNEW[12]*dt;X[5,t]=X[5,t-1]+UNEW[13]*dt;Y[5,t]=Y[5,t-1]+UNEW[14]*dt;Z[5,t]=Z[5,t-1]+UNEW[15]*dt;X[6,t]=X[6,t-1]+UNEW[16]*dt;Y[6,t]=Y[6,t-1]+UNEW[17]*dt;Z[6,t]=Z[6,t-1]+UNEW[18]*dt;Print[t];t++}]
In[]:=
d1=ListPointPlot3D[{Table[{X[1,t],Y[1,t],Z[1,t]},{t,0,307}],Table[{X[2,t],Y[2,t],Z[2,t]},{t,0,307}],
Table[{X[3,t],Y[3,t],Z[3,t]},{t,0,307}],Table[{X[4,t],Y[4,t],Z[4,t]},{t,0,307}],
Table[{X[5,t],Y[5,t],Z[5,t]},{t,0,307}],Table[{X[6,t],Chop[Y[6,t]],Z[6,t]},{t,0,307}]},BoxRatios->{1,1,1}]
Table[{X[3,t],Y[3,t],Z[3,t]},{t,0,307}],Table[{X[4,t],Y[4,t],Z[4,t]},{t,0,307}],
Table[{X[5,t],Y[5,t],Z[5,t]},{t,0,307}],Table[{X[6,t],Chop[Y[6,t]],Z[6,t]},{t,0,307}]},BoxRatios->{1,1,1}]
Two spheres sedimenting in a Newtonian fluid. They rotate in opposite directions as they fall.
In fact, we find that the rotational velocities are UNEW[7]=0.2082 and UNEW[10]= -0.2082.
In fact, we find that the rotational velocities are UNEW[7]=0.2082 and UNEW[10]= -0.2082.
In[]:=
FU={0,0,-10,0,0,-10,0,0,0,0,0,0};
In[]:=
X[1]=3;Y[1]=0;Z[1]=0;X[2]=-3;Y[2]=0;Z[2]=0;
In[]:=
dt=0.1;
In[]:=
NSPHR=2;NDIM=6 NSPHR;Ng=3 NSPHR;
In[]:=
X[1,0]=X[1];Y[1,0]=Y[1];Z[1,0]=Z[1];
X[2,0]=X[2];Y[2,0]=Y[2];Z[2,0]=Z[2];
X[2,0]=X[2];Y[2,0]=Y[2];Z[2,0]=Z[2];
In[]:=
For[t=1,t<100,{NEWTON2[{X[1,t-1],X[2,t-1]},{Y[1,t-1],Y[2,t-1]},{Z[1,t-1],Z[2,t-1]},t];
X[1,t]=X[1,t-1]+UNEW[1]*dt;Y[1,t]=Y[1,t-1]+UNEW[2]*dt;Z[1,t]=Z[1,t-1]+UNEW[3]*dt;
X[2,t]=X[2,t-1]+UNEW[4]*dt;Y[2,t]=Y[2,t-1]+UNEW[5]*dt;Z[2,t]=Z[2,t-1]+UNEW[6]*dt;Print[t];t++}]
X[1,t]=X[1,t-1]+UNEW[1]*dt;Y[1,t]=Y[1,t-1]+UNEW[2]*dt;Z[1,t]=Z[1,t-1]+UNEW[3]*dt;
X[2,t]=X[2,t-1]+UNEW[4]*dt;Y[2,t]=Y[2,t-1]+UNEW[5]*dt;Z[2,t]=Z[2,t-1]+UNEW[6]*dt;Print[t];t++}]
In[]:=
Table[UNEW[i],{i,1,12}]