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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
|
.. _sphere_in_vessel:
A rigid sphere floating in an hydrostatic tank
----------------------------------------------
This example demonstrates the API of running a rigid fluid coupling problem in PySPH.
To run it one may do::
$ cd ~/pysph/pysph/examples/rigid_body/
$ python sphere_in_vessel_akinci.py
There are many command line options that this example provides, check them out with::
$ python sphere_in_vessel.py -h
The example source can be seen at `sphere_in_vessel.py
<https://github.com/pypr/pysph/tree/main/pysph/examples/rigid_body/sphere_in_vessel_akinci.py>`_.
This example demonstrates:
* Setting up a simulation involving rigid bodies and fluid
* Discuss mainly about rigid fluid coupling
It is divided in to three parts:
* Create particles
* Create equations
* Run the application
Create particles
~~~~~~~~~~~~~~~~~~~~~~~~~~~
In this example, we have a tank with a resting fluid and a sphere falling into
the tank. Create three particle arrays, ``tank``, ``fluid`` and ``cube``.
``tank`` and ``fluid`` has to obey ``wcsph`` scheme, where as ``cube`` has to obey
rigid body equations.
.. code:: python
def create_particles(self):
# elided
fluid = get_particle_array_wcsph(x=xf, y=yf, h=h, m=m, rho=rho,
name="fluid")
# elided
tank = get_particle_array_wcsph(x=xt, y=yt, h=h, m=m, rho=rho,
rad_s=rad_s, V=V, name="tank")
for name in ['fx', 'fy', 'fz']:
tank.add_property(name)
cube = get_particle_array_rigid_body(x=xc, y=yc, h=h, m=m, rho=rho,
rad_s=rad_s, V=V, cs=cs,
name="cube")
return [fluid, tank, cube]
We will discuss the reason for adding the properties :math:`fx`, :math:`fy`, :math:`fz` to the
``tank`` particle array. The next step is to setup the equations.
Create equations
~~~~~~~~~~~~~~~~~~~~~~~~~~~
.. code:: python
def create_equations(self):
equations = [
Group(equations=[
BodyForce(dest='cube', sources=None, gy=-9.81),
], real=False),
Group(equations=[
SummationDensity(
dest='fluid',
sources=['fluid'], ),
SummationDensityBoundary(
dest='fluid', sources=['tank', 'cube'], fluid_rho=1000.0)
]),
# Tait equation of state
Group(equations=[
TaitEOSHGCorrection(dest='fluid', sources=None, rho0=self.ro,
c0=self.co, gamma=7.0),
], real=False),
Group(equations=[
MomentumEquation(dest='fluid', sources=['fluid'],
alpha=self.alpha, beta=0.0, c0=self.co,
gy=-9.81),
AkinciRigidFluidCoupling(dest='fluid',
sources=['cube', 'tank']),
XSPHCorrection(dest='fluid', sources=['fluid', 'tank']),
]),
Group(equations=[
RigidBodyCollision(dest='cube', sources=['tank'], kn=1e5)
]),
Group(equations=[RigidBodyMoments(dest='cube', sources=None)]),
Group(equations=[RigidBodyMotion(dest='cube', sources=None)]),
]
return equations
A few points to note while dealing with *Akinci* formulation,
1. As a first point, while computing the density of the ``fluid`` due to solid,
make sure to use ``SummationDensityBoundary``, because usual
``SummationDensity`` computes density by considering the mass of the
particle, where as ``SummationDensityBoundary`` will compute it by
considering the volume of the particle. This makes a lot of difference
while dealing with heavy density variation flows.
2. Apply ``TaitEOSHGCorrection`` so that there is no negative pressure.
3. The force from the boundary (here it is tank) on fluid is computed using
``AkinciRigidFluidCoupling`` equation, but in a usual case we do it using the
momentum equation. There are a few advantages by doing this. If we are
computing the boundary force using the momentum equation, then one should
compute the density of the boundary, then compute the pressure. Using such
pressure we will compute the force. But using ``AkinciRigidFluidCoupling`` we
don't need to compute the pressure of the boundary because the force is
dependent only on the fluid particle's pressure.
.. code:: python
def loop(self, d_idx, d_m, d_rho, d_au, d_av, d_aw, d_p,
s_idx, s_V, s_fx, s_fy, s_fz, DWIJ, s_m, s_p, s_rho):
# elide
d_au[d_idx] += -psi * _t1 * DWIJ[0]
d_av[d_idx] += -psi * _t1 * DWIJ[1]
d_aw[d_idx] += -psi * _t1 * DWIJ[2]
s_fx[s_idx] += d_m[d_idx] * psi * _t1 * DWIJ[0]
s_fy[s_idx] += d_m[d_idx] * psi * _t1 * DWIJ[1]
s_fz[s_idx] += d_m[d_idx] * psi * _t1 * DWIJ[2]
Since in ``AkinciRigidFluidCoupling`` (more in next point) we compute both
force on fluid by solid particle and force on solid by fluid particle,
which makes our sources to hold the properties ``fx``, ``fy`` and ``fz``.
4. Here first few equations deal with the simulation of fluid in hydrostatic
tank. The equation dealing with rigid fluid coupling is
``AkinciRigidFluidCoupling`` . *Coupling* equation will deal with forces
exerted by fluid on solid body, and forces exerted by solid on fluid. We
find the force on fluid by solid and force on the solid by fluid in a singe
equation.
Usually in an SPH equation, we tend to change properties only of a destination
particle array, but in this case, both destination and sources properties are
manipulated.
5. The final equations deal with the dynamics of rigid bodies, which are
discussed in other example files.
Run the application
~~~~~~~~~~~~~~~~~~~~~~~~~~~
Finally run the application by
.. code:: python
if __name__ == '__main__':
app = RigidFluidCoupling()
app.run()
|
