In this section, we describe our full technique on liquid simulation in a virtual lab. First, in sections 3.1 and 3.2, we address the problem of simulating liquids inside various vessels of complex structure by leveraging advanced analytical geometry. Thereafter, in section ..., building on that liquid simulation method, we propose another method, which deals with the simulation of liquid transfer between laboratory vessels of various shapes. The results from both of those simulation methods’ implementation under Unity/C# are presented in section 4.
3.1. Description of Pipes and Bottles with Circular Cross-Sections and Liquid Volume Representation
During the experimental procedure for determining the rate constant of methyl acetate hydrolysis in the presence of an acid as well as in all experimental procedures in Biology and Chemistry laboratories, bottles of various sizes and shapes are extensively used. Instruments with tubes, which facilitate the flow of liquid from one bottle to another, are also employed.
In the experiment examined in this study, there is extensive use of and transfers between conical, cylindrical/volumetric bottles and burettes. Specifically, during the use of an automatic burette, the following processes occur:
Liquid flow through a tube, bottle emptying, and volumetric tube filling with free-flow during the filling of the automatic burette’s volumetric tube.
Liquid flow through a tube and bottle refilling with free-flow from the automatic burette during the zeroing of the volumetric tube’s measurement.
Emptying of the volumetric tube with liquid flow into the spout tube of the automatic burette’s measurement tap.
From the above, the utility and significance of simulating liquid flows, as well as the filling and emptying of bottles with liquids, become evident for enhancing the sense of realism experienced by the student/user. To achieve this, the movement of the liquid as a function of time must be simulated, ensuring this function closely mirrors the actual quantitative and temporal behavior of liquid flow within a tube or bottle.
3.1.1. Types and Classification of Liquid Movement
Liquid movement can be broadly classified, based on the applications of interest, into the following categories:
Steady liquid flow in a tube
Transient liquid flow in a tube
Filling and emptying of a bottle
This work seeks a unified mathematical approach to describing pipes and bottles, facilitating seamless implementation in Unity/C#. Specifically, the goal was to parameterize the description of an entire composite pipe, enabling straightforward modeling and simulation across consecutive pipe segments.
3.1.2. Simplifying Assumptions
The following simplifying assumptions were made for describing pipes and bottles:
Pipes and bottles are assumed to have circular cross-sections only.
A single parameter, s, is used along the entire length of the composite pipe to represent the current length of the line formed by the centers of the circular cross-sections (central line).
The diameter is continuous along each pipe/bottle segment and remains consistent during transitions between segments (i.e., the diameter at the end of one segment equals the diameter at the start of the next).
The slope of the central line relative to the horizontal plane is continuous at every point along the pipe and during transitions between segments.
The diameter of an individual pipe/bottle remains constant or varies linearly as a function of its length.
For pipes, the liquid surface level is considered perpendicular to the central line at every point.
Only vertical bottles are considered.
A common curvature parameter, κ, is used for each pipe segment.
The radius of curvature of the central line must be greater than or at least equal to the pipe’s cross-sectional radius.
The central line of the composite pipe lies within a plane.
Essentially, the aim is to determine the Equations describing the coordinates of any point along the central line in the plane and the pipe’s slope at that point. This enables drawing the corresponding liquid level at that specific point.
3.1.3. Central Line Coordinate Equations for Curved Pipes with Constant Diameter and Straight Pipes with Constant or Linearly Varying Diameter
Recall that a composite pipe consists of various segments, some straight and some curved. To position each segment of the composite pipe in the plane, we initially identify a single point—the starting point of the segment’s central line. Naturally, the endpoint of one segment serves as the starting point for the next, so only the initial points of consecutive segments need to be specified. For each starting point of a pipe segment’s central line, the parameters that must be measured or determined include its horizontal and vertical coordinates and the slope (angle) of the central line relative to the horizontal plane.
Once these three values (horizontal coordinate, vertical coordinate, and slope) are defined for the start of a pipe segment, they can be used to calculate corresponding values for any arbitrary point along that segment. This is the basis for fully defining the positioning of the segment in the plane.
Let’s assume that the horizontal coordinate, vertical coordinate, and slope at point A of the central line of a curved segment of the burette’s composite pipe are known. To calculate the corresponding values at another arbitrary point B along the same segment, refer to Figure 7 for illustration.
Figure 7.
(a) Initial point A and arbitrary point B of the central line of a curved segment of a composite pipe. (b) Diagram for deriving the coordinate Equations of arbitrary point B.
Figure 7.
(a) Initial point A and arbitrary point B of the central line of a curved segment of a composite pipe. (b) Diagram for deriving the coordinate Equations of arbitrary point B.
According to
Figure 7(b), the angle, horizontal coordinate y, and vertical coordinate of point B are given by the following Equations:
The angle ψ is defined such that its sides are perpendicular to those of the angle
θs/2, leading to:
For the angle θ
s, considering that the curvature κ equals 1/R, we have:
As for the length m, it is calculated as follows:
Substituting ψ, θ
s and m from Equations (4), (5) and (6) into Equations (1), (2) and (3), respectively, we derive the following expressions for the coordinates and angle of point B:
For straight pipe segments, the above Equations degenerate into the following:
Equations (10), (11) and (12) for straight pipe segments apply not only to the straight sections of the composite pipe of the burette but also to volumetric cylinders and conical flasks, whether wide-necked or not (Figures Figure 3 and Figure 4).
3.1.4. Geometric and Mathematical Description of the Coordinates of the Perimeter of a Pipe or Vessel with an Inclined Circular Cross-Section
Initially, we will describe the mathematical coordinates of a random point C on the liquid’s perimeter in a pipe with an inclined circular1 cross-section (Figure 8). Knowing the coordinates of the points on the circular perimeter, we can later in Unity draw a polygonal approximation of the latter. By successively drawing the perimeters corresponding to each point on the central line, we have the full liquid lateral surface inside the pipe.
Figure 8.
The circular cross-section of an inclined pipe with the horizontal level.
Figure 8.
The circular cross-section of an inclined pipe with the horizontal level.
In order for the calculation to be done, the diagram in Figure 8 is deployed into the ones of Figure 9.
Figure 9.
(a) Cross-sectional view at a point on the central line with a horizontal orientation. (b) Central-line plane of a pipe segment.
Figure 9.
(a) Cross-sectional view at a point on the central line with a horizontal orientation. (b) Central-line plane of a pipe segment.
From
Figure 9(a), we have:
Now, according to
Figure 9(b), when the cross-section is inclined by angle θ,
can be analyzed into d
y and d
z as follows:
Combining Equations (14) and (15) with the z-coordinate of the point on the central line of the cross-section, and Equations (14) and (16) with the y-coordinate of the same point on the central line of the cross-section, we derive the general Equations for the points of the circular perimeter in 3D space:
We remind that xc is already given in Equation (13).
3.2. Discretization of Pipes and Bottles’ Cross-Sections
In sections 3.1.3 and 3.1.4, we calculated the Equations for the co-ordinates of any point on a central line as well as the Equations for the co-ordinates of the respective points on the central line’s circular cross-section with the horizontal level. Having done so, we can dynamically create in Unity the meshes of the interior liquids of the various vessels of our simulation.
However, the circles upon which the liquid mesh creation is made cannot obviously be drawn as such. Instead, we approximate those with polygons through a process called discretization. The latter depends heavily on the number of vertices we decide our polygons to have (Figure 10).
Figure 10.
An example of 8-point discretization.
Figure 10.
An example of 8-point discretization.
Figure 11 shows the notation for points in each cross-sectional circle.
Figure 11.
Discretization notation in cross-sectional circles.
Figure 11.
Discretization notation in cross-sectional circles.
Specifically, in Figure 11, δ represents the cross-section number, ranging from 0 to the total number of cross-sections minus 1, while k denotes the sequential number of the point within each cross-section, ranging from 1 to the total number n of discretization intervals.
The number n, chosen by the developer of the simulation, applies to all cross-sections and determines the fidelity of the circle representation. In the example above, the value of n equals 4.
The coordinates of the aforementioned points will be calculated based on angular intervals of 2π/n at angles ϕk=k⋅2π/n. In this specific example, the angles are π/2, π, 3π/2 and 2π. For these angles, using the radius of the circle of the current cross-section, the inclination of the cross-section, the coordinates of the circle’s center (the central line point), and Equations (13), (17) and (18), the coordinates of each discretized point of the circle/cross-section are computed.
This method of notation—continuous sequential numbering of the discretized points of cross-sectional circles—was mandated by the way points are defined in Unity’s method for creating triangular finite surface elements. In this method, the sequential number of each point is declared in a one-dimensional array whose elements are the above 3D points (Vector3 variable type).
3.2.1. Definition of Triangular Surface Finite Elements Using Cross-Section Points
The following figures illustrate possible combinations of points forming surface finite triangles. Connections are chosen only between neighboring points within the same cross-section and points at adjacent angles between neighboring cross-sections, as shown in Figure 12.
Figure 12.
Connections Between Points for Optimal Resolution.
Figure 12.
Connections Between Points for Optimal Resolution.
Each point on the last cross-section can participate in two to four triangular elements.
Each line connecting points of adjacent cross-sections at the same angle forms the edge of exactly two triangular elements. Consequently, the total number of triangular elements equals 2n.
Each such line can participate in three different triangular element combinations; that is, one edge of each triangle always starts in one cross-section and ends in the neighboring cross-section and at a neighboring angle:
Triangles with one edge in the previous cross-section.
Triangles with one edge in the last cross-section.
Triangles with edges in different cross-sections.
Figure 13 illustrates the triangle combinations for each edge between cross-sections at the same angle.
Figure 13.
Triangle combinations for each edge between cross-sections at the same angle.
Figure 13.
Triangle combinations for each edge between cross-sections at the same angle.
To ensure the visible side of the triangle (the outer surface in this application) appears in Unity, the order of point declaration for the triangle must follow a clockwise sequence (Figure 14(a)).
Figure 14.
Orders of point declaration for visible surface side: (a) Clockwise. (b) Counter-clockwise.
Figure 14.
Orders of point declaration for visible surface side: (a) Clockwise. (b) Counter-clockwise.
3.2.2. Mathematical Handling and Implementation Specifics of Triangles Between Cross-Sections
The diameters (
d) of pipes/bottles (half of which is the radius
r used in the circle coordinate Equations for inclined cross-sections) as well as the lengths (
l) and the curvatures (
κ) of the central lines are stored as data in a description matrix.
To account for potential tapering in the first section, the initial diameter must also be declared. For subsequent sections, since continuity between sections is assumed, the initial diameter of a section is determined by the final diameter of the previous section. For instance, for the second section, the initial diameter is d1 and the final diameter is d2.
Using the initial diameter of the composite pipe, the initial slope, and the starting coordinates of the central line, discretization points are generated in the initial cross-section to form the first triangular surface finite elements using the discretization points of the initial cross-section (excluding the start point).
The surface between two points on the central line is created through the following steps:
Assuming a sufficiently small step length for the central line and incrementally increasing the current length in a loop.
Calculating the coordinates of the central line point corresponding to the current length at each iteration.
Based on the above coordinates, iteratively calculating the discretization points of the cross-section’s perimeter using angular intervals (as determined by the chosen resolution). These points are recorded sequentially in a one-dimensional array.
Using the discretization points of the previous and current cross-sections, triangular surface finite elements are created by recording the sequential numbers of each point in the previous array. Each set of three points, in order, forms a triangular element. The two arrays are assigned as corresponding vertex and edge attributes in a mesh structure.
To minimize iterations and algorithm execution time (thereby reducing time between iterations, improving dynamic simulation of level variations, and increasing surface accuracy), the grouping of vertices for triangle creation is further analyzed. That is, the existence of a possible pattern for the set of each cross-section is investigated so that the implementation can be included in the existing point/peak generation loop.
For each cross-Section , starting from the first defined point , the following two triangular finite elements can be declared using points from the previous cross-section:
and
For the next point , combinations with points incrementally higher than those in the previous sets are tested:
and
Following this approach, triangular finite elements are defined sequentially for the last point in the current cross-section:
and
Figure 15 illustrates a graphical depiction of the sequence for recording/representing the triangular finite elements for each point of a cross-section with 4 discretization points.
Figure 15.
Creation of triangles from the first point of the new cross-section (point 5).
Figure 15.
Creation of triangles from the first point of the new cross-section (point 5).
In Figure 15, the triangles that can be defined after creating the first point of the new cross-section (point 5) are shown. Thereby the two created triangles for point 5 are:
and
Next, the triangles for the remaining three points of the cross-section are created (at this point, in addition to the points from the previous cross-section, the neighboring points of the new cross-section that have already been created are available).
Thus, the triangles for point 6 (Figure 16) are:
and
Figure 16.
Creation of triangles from the second point of the new cross-section (point 6).
Figure 16.
Creation of triangles from the second point of the new cross-section (point 6).
Next come the triangles for point 7 (Figure 17) which are:
and
Figure 17.
Creation of triangles from the penultimate point of the cross-section (point 7).
Figure 17.
Creation of triangles from the penultimate point of the cross-section (point 7).
Finally, the triangles for point 8, which is the last point of the cross-section (Figure 18) are:
and
Figure 18.
Creation of triangles from the last point of the cross-section (point 8).
Figure 18.
Creation of triangles from the last point of the cross-section (point 8).
In total, we have the following pairs of triangles for each of the points 5, 6, 7 and 8:
and
and
and
and
It is easily observed that for a point with a total serial number and being the number of discretization points in each cross-section, a numerical pattern is followed for the point’s two triangles:
and
For the next point , the numbers in the above triplets forming the finite triangular elements are increased by 1 (that is, the next point from the corresponding triplets of the previous point is recorded):
and
This pattern can be directly utilized in implementation to generate the finite surfaces simultaneously, without requiring a new loop, during the creation of discretization points for each cross-section.
3.3. Simulation of Dynamic Flow and Liquid Level Variation in Bottles and Tubes
As mentioned in section 3.1, the developed methodology for modeling the surface of tubes and liquids aims to define the surface as a function of a single variable, that of the current local length of the bottle or tube segment. By knowing the current local length—representing the liquid level for bottles or the volume in general for the container (bottle or tube)—and the characteristic data of the container (from the corresponding registry), the current diameter and local slope of the container at every point along its centerline are calculated accordingly. This approach enables the generation of the corresponding surface from the base to the current length of the container.
The goal of this modeling effort has been, from the outset, the dynamic simulation of liquid level variation during filling or removal from the container. The developed methodology successfully separates the implementation of surface creation from the mathematical modeling of the one-dimensional dynamic liquid flow (which dynamically provides the current liquid volume in each container).
The remaining task is to mathematically calculate the local liquid level (the current length) in each circular cross-section container, assuming that the simplifying assumptions are outlined in section 3.1.2 are adhered to. The resulting function for calculating the local current length and slope is used as input to the implementation described in section 3.1 to design the corresponding surface as its outcome.
The mathematical modeling of liquid flow is achieved using fluid mechanics Equations. Subsequently, a mathematical approach is presented for modeling the filling of a bottle under a known liquid volume flow rate.
In cases where the filling occurs with a known liquid volume flow rate over time, the bottle is filled by transferring liquid via free flow (without losses and at atmospheric pressure) from a different container.
To calculate the instantaneous liquid volume flow rate into the target bottle, it suffices to know the outflow velocity, the cross-sectional area of the outflow opening of the supply container, and the proportion of volumetric outflow directed into the target bottle. The instantaneous flow rate is given as the product of these three factors, which can be determined by solving a system of fluid mechanics Equations for the flow within the supply container.
3.3.1. Geometric and Mathematical Description of Bottle Levels with Respect to a Known Liquid Flow Rate
Given a liquid flow rate—either constant or expressed as a time-dependent function—and knowing the time step of each iteration of the simulation loop (handled by a C# Unity function), the additional volume added to the tube between the previous and current iterations can be calculated. In this sense, the volume is treated as a known quantity.
Figure 19 shows a schematic representation of a tube with an increasing diameter cross-section.
Figure 19.
Schematic Representation of Liquid Volume in a Diverging Frustum Cone as a Function of Its Level.
Figure 19.
Schematic Representation of Liquid Volume in a Diverging Frustum Cone as a Function of Its Level.
The volume V between two levels, s
1 and s
2, in
Figure 19 is given by:
To relate parameter H with the levels and diameters, the tangent of angle ψ is used:
The linear function of diameter with respect to the level/point on the central axis of the tube is expressed as:
where the parameter
equals the tangent of the cone angle
:
In Equation (19), we substitute H with its value from Equation (20); we then substitute d1 and d2 with their values from Equation (21) for s being s1 and s2, respectively; and finally, we rearrange all terms to one side of the Equation, simplifying into terms of s2 (the level of interest).
The resulting Equation becomes:
Equation (23) is a complete cubic polynomial of the form:
The discriminant Δ is calculated as:
Substituting the coefficients, the discriminant becomes:
Since Δ is always negative, the Equation has one real solution and two conjugate complex roots. This result is desirable because, in other cases, there could be either three distinct real roots or one double real root and one single real root. The real root provides the exact analytical expression for the new liquid level based on the given data and the previous known level.
The general solution for Equation (24) is:
Substituting the coefficients from Equation (23), the real solution simplifies to:
Verification of this solution by substitution into Equation (23) shows the left-hand side becomes zero, confirming its correctness. Simplifying by cubing both sides and rearranging terms verifies equivalence to Equation (23).
We will now consider the case of the truncated conical tube of decreasing diameter, according to the schematic representation of Figure 20.
Figure 20.
Schematic Representation of Liquid Volume in a Converging Frustum Cone as a Func-tion of its Level.
Figure 20.
Schematic Representation of Liquid Volume in a Converging Frustum Cone as a Func-tion of its Level.
The liquid volume
V in
Figure 20 is given by:
Using the tangent of angle
ψ to relate
h to the other parameters, we get:
The same process applied to Equation (28) results in a cubic polynomial identical to Equation (23). Hence, its solution is also identical to Equation (27).
For completeness, the case of a straight tube with constant cross-sectional diameter is included. For this case, therefore, we have:
3.3.2. Key Points and Implementation Examples for Bottle Filling in Unity
The modeling and implementation of surface creation accept the following as input parameters for the initial cross-section characteristics:
coordinates of the central line point of the initial cross-section
slope of the initial cross-section
diameter of the initial cross-section
In the implementation, it is necessary to retrieve the time interval Δt that has elapsed since the previous iteration of the execution loop in the corresponding Unity C# script.
The implementation requires tracking and monitoring of the following variables:
The current volume V of liquid in the bottle,
The cumulative volume Vtot,i that each segment of the bottle can contain.
This information is critical to determine the segment of the bottle where the liquid level will rise after adding the new volume to the existing volume of liquid. Specifically, if the flow rate is sufficiently high (or the segment volume is sufficiently small) so that within the time interval between two loop iterations the liquid level can move through one or more segments, the surface can be rendered for all such segments up to the level of the current segment.
At the beginning of the implementation code (during the initialization function) and before the execution loop starts, the cumulative volumes stored from the bottom of the bottle to each segment Vtot,i are calculated and recorded.
The volume of a frustum of a cone with
da and
dτ being its bases and
l its height (
Figure 19) is given by the formula:
If, in Equation (32), dα and dτ are replaced by d (for a bottle segment with a constant cross-sectional diameter) and the length l is replaced by s2−s1, then the equation becomes identical to Equation (30).
Thus, Equation (32) can be used to calculate the individual volumes of all types of bottle segments examined in this study (constant-diameter and linearly varying diameter shapes).
We will now express the volume of a curved pipe section as a function of its length
We know from calculus [
10] that the volume V of a solid of revolution, where the cross-sectional area is bounded by the functions
R(x)=y1 and
r(x)=y2, is given for the segment of the solid between heights α and β by the following integral:
In our case, we are interested in the volume of liquid within a segment of circular cross-section, as illustrated in Figure 21.
Figure 21.
Liquid volume within a segment of circular cross-section.
Figure 21.
Liquid volume within a segment of circular cross-section.
The upper semicircle in the said figure is expressed as:
while the lower semicircle is expressed as:
For rotation through an angle
θs, the volume
Vs is expressed as:
By substituting
s=θsRs, we get:
If we further replace s with s2−s1 and r with d/2, the equations remain consistent with Equations (30) and (32), which also hold for cylindrical volumes.
Before the execution loop begins, the points of the pipe’s initial cross-section are created, and if an initial liquid level exists, the points of the cross-section at the initial liquid level and the surface between them are also created. Finally, the initial volume is calculated (this functionality has been implemented for an initial level in the first segment). During the execution loop, based on the given liquid flow rate, the new volume of liquid to be added to the pipe for each time interval between iterations is calculated as ΔV=Q⋅Δt. The new total current volume is then compared with the cumulative volumes up to each segment to determine in which segment the new liquid level will be located:
Two cases are distinguished: if the new volume remains in the same segment as the previous one, the new liquid level is calculated using
V=Vloc=ΔV and either Equation (27) or (31), depending on whether the segment’s diameter is variable or constant. Otherwise, if the new volume spans multiple segments, and if the new volume exceeds the total volume of the composite bottle, the entire bottle surface is rendered, and the loop execution is terminated. If not, the surfaces up to the previous segment are rendered, and the local volume for the current segment is calculated as:
This value is then substituted back into either Equation (27) or Equation (31), depending on whether the segment’s diameter is variable or constant (where s1, the initial level, is replaced with zero).