Elastic Simulation of Joints with Particle-Based Fluid

: Skinning, which is used in skeletal simulations to express the human body, has been weighted between bones to enable muscle-like motions. Weighting is not a form of calculating the pressure and density of muscle ﬁbers in the human body. Therefore, it is not possible to express physical changes when external forces are applied. To express a similar behavior, an animator arbitrarily customizes the weight values. In this study, we apply the kernel and pressure-dependent density variations used in particle-based ﬂuid simulations to skinning simulations. As a result, surface tension and elasticity between particles are applied to muscles, indicating realistic human motion. We also propose a tension yield condition that reﬂects Tresca’s yield condition, which can be easily approximated using the difference between the maximum and minimum values of the principal stress to simulate the tension limit of the muscle ﬁber. The density received by particles in the kernel is assumed to be the principal stress. The difference is calculated by approximating the moment of greatest force to the maximum principal stress and the moment of least force to the minimum principal stress. When the density of a particle increases beyond the yield condition, the object is no longer subjected to force. As a result, one can express realistic muscles.


Introduction
Skeletal simulation is one of the central components of digital content, such as movies, computer games, and virtual reality (VR). Particularly in VR, the fingers, hands, and arms handling the controller actually move, and at the same time, those motions must be drawn on the VR screen. Therefore, the deformation of the human body by the skeleton should be quickly rendered. Moreover, since the position of the human body expected by moving the body should match the position of the human body viewed with the eyes on the screen, the expression of a realistic human body is also important. However, owing to the complexity of the human body composed of bones, muscles, and soft tissues, it is still challenging to represent the realistic motions of virtual characters. Early human models consisted of several objects connected with joints. It was intended to describe the motion of the human body by predicting the movement of the joints using kinetic or dynamic methods. Although these methods can express human motion somewhat similarly, it is very difficult to obtain a realistic rendering of the human body. The motion of the human body is dependent on the movement of the skeleton. However, more precisely, the motions observed outside the body are not the movements of the skeleton itself, but the movements of the skin associated with the skeleton. Therefore, the skinning technique, which is based on the movement of the skeleton, was developed and utilized. One such skinning technique is typical representative, which defines the position of the skin according to the distance between the skeleton and the skin as a linear function to obtain skin deformation animation for human behavior. However, skin changes are often not linearly expressed in relation to the movement of the skeleton, so there is a limit to obtaining realistic human appearance motions. Since the late 1990s, attempts have been intensively made to simulate skin movements by simulating human muscles based on anatomical facts. One of these attempts is to calculate the movement of a muscle that is non-linearly deformed according to the movement of the skeleton through physics simulation. Then the skin deformation by the skeleton is defined through the computed result. Since accurate animations of the human body require accurate physical properties such as muscle mass, density, and elasticity to be modeled, the body's muscle modeling and physical simulation techniques are different but closely related. In addition, even though the physical properties of the human body can be accurately modeled, the computational cost is significantly high to represent the behavior of real muscles. Fluid simulation, which uses the physical properties of the fluid to describe its motion, is mainly used in movies and games to express the actual motion of a liquid. Reeves proposed particle systems to represent fluids and smoke or fire [1]. However, fluid data composed of so many particles are difficult to render in real-time even with modern graphics hardware. Generally, the smoothed particle hydrodynamics (SPH) technique is used for the efficient rendering of fluid particles. SPH, developed in astrophysics, is a method of calculating the motion of particles with properties interacting with other particles within the controlled range of a weight and smoothing kernel [2]. SPH techniques can be used to intuitively calculate the motion of a solid or fluid, which varies over time owing to interactions between particles, and are therefore suitable for expressing the motion of an easily deformable solid or fluid.
In this paper, we propose a particle-based skinny simulation method using SPH techniques. We defined the muscle as a particle-based model. The movements of particles that replace muscles are calculated using SPH techniques. These particles express the motion of the human body by the movement of the skeleton. When the human body is subjected to external forces, the physical properties of the particles change as the particles are moved by external forces, which are used to calculate the force exerted on each particle to calculate the motion of the particles. This allows the body to accurately express in real-time the physical motions that occur when it moves under external force.
We also experimented with Tresca's yield condition to express the following limitations of human motions. The muscles and joints of the human body have limited movements. Objects applied with yield conditions will not return to their original form with more stress than a certain amount. We experiment with this based on a three-dimensional (3D) elbow model. As a result, further deformation was prevented if the particles received more than a certain amount of stress. This allowed the human body to no longer deform if its deformation limit has already been reached.
By our method, we can implement the motion of the body in the form of particles in less time. In addition, it was possible to create realistic human body motion with the limits of the human body. This paper is organized as follows. Section 2 introduces related studies. Section 3 deals with existing real-time skinning techniques and problems. Subsequently, Section 4 describes particle-based skinning simulations, Section 5 shows the experimental results of the method proposed in this paper, and the conclusions are given in Section 6.

Related Works
Skinning is a necessary method for the natural movement of the human body or animal. The most frequently used method for implementing skinning in real-time applications is the geometric method. Linear blend skinning (LBS), first proposed by Magnenat-Thalmann et al. [3], is the most widely used method for the visual representation of moving characters. Because this method is intuitive, it can be easily implemented with a low computational cost. LBS defines the motion of an object by setting weights on the transformation matrix applied to the skeleton. The result of the motion by the skeleton to one surface of an object using LBS is obtained by the following computational process: sum the weights for each pre-specified skeleton on the surface and multiply by its transformation matrix. Subsequently, the resulting value of the transformation matrix is multiplied by the location of the surface. This formula must specify how much weight the bone has on each surface. This method defines the deformed position of an object as a linear interpolation of the skeleton, resulting in visual errors caused by the linear interpolation. Typically, there is a candy-wrapper phenomenon in which an object contracts unnaturally like a candy bag by linear interpolation and loses its volume.
A number of nonlinear skinning methods have been proposed to address the limitations of LBS. Cordier and Magnenat-Thalmann proposed a skinning method using log-matrix blending [4]. Kavan and Žára proposed the spherical blending method [5]. Subsequently, dual quaternion skinning (DQS) was proposed using a dual quaternion [6]. The concept of a dual quaternion is a dual number representation of one real number and three imaginary coefficients. A total of eight coefficients are used for one dual quaternion. While a regular quaternion is often used to express rotations only, a dual quaternion allows for translation in addition to rotations.
In skinning, a dual quaternion is used to represent the shift and rotation transformation applied to the skeleton and then weighs the dual quaternion to simulate the motion of the object. This method expresses the motion of an object by adding weights to the dual quaternion, summing them all, and normalizing them. DQS reduces visual errors caused by the linear interpolation of LBS while allowing skinning to be represented in frames slightly lower than LBS and also allowing for the real-time representation of object movements.
Example-based methods have been developed to represent deformable objects. Sloan et al. proposed a method of using various motion examples of character models to model the deformation behavior of character models [7]. Park and Hodgins proposed a way to record the surface movement of the skin to create a database of dynamic skin modifications and to apply dynamic equations to pre-recorded data [8]. Lewis and Anjyo used the blendshape technique to express variations of the human body [9]. Example-based methods use pre-generated basic poses, so artists must create multiple basic poses themselves or use methods such as motion capture to pre-create basic poses.
Physically-based simulation is used as a way of expressing the motion of a transforming object. Capell et al. used the linear finite element method (FEM) for controlling a lattice to express bone-induced skin deformation [10]. Kim and Pollard proposed a method to use nonlinear FEM to obtain a high modeling accuracy of biological tissue [11]. An additional method of using appropriate material properties for physics-based simulation methods is also proposed [12]. Physics-based simulation methods provide visual effects that are difficult to represent in geometric methods. They include collision [13], skin sliding and jiggling [14], and wrinkling [15].
Research on muscles has been conducted in the field of computer graphics to express sophisticated movements of the human body. The critter system developed by Chadwick et al. [16] and the bodybuilder system developed by Shen and Thalmann [17,18] separated the body into the skeletal layer, the skin layer, and the middle layer with volume to obtain a more realistic body shape and movement. Scheepers et al. modeled a human body with shapes similar to the anatomical schema defined a connection between the internal and external structures that cause the appearance change and modeled the muscles of the ellipsoid to anatomically approximate the deformation of the human body [19]. Teran et al. studied how to obtain the muscles of realistic human models from cross-sectional images of real human bodies [20].
Various studies have also been conducted to obtain a realistic appearance in accordance with human motion. Barr animated the entire world and regional transform structure [21]. Sederberg and Parry expressed non-rigid animation by applying a free-form variant [22]. Witkin et al. applied energy-constraining techniques [23], and Lazarus et al. applied axial deformation [24]. Wilhelms and Van Gelder introduced a metaball to express nonlinear animation for the skeletal system [25]. Nedel created a human model and behavior based on physics-based simulations [26], and Matias et al. conducted a study to determine the location of muscle-attaching bones based on landmark points in the bones [27]. Levin et al. evaluated myofiber orientation from diffuse tensor images [28]. Yang et al. proposed a method for estimating centers of rotation (CoRs) based on 15 male finger skeleton computerized tomography scanner data for accurate 3D motion of finger bones [29].
Methods for fluid simulation are classified as Eulerian and Lagrangian methods. The Eulerian method is a technique for creating lattices in Euclidean space and indirectly moving fluid particles using velocity vectors calculated on the lattice [30]. This method can represent realistic atmospheric advection and is mainly used in applications where sophisticated fluid representation is valued.
Foster and Metaxas introduced the first 3D fluid simulation by applying the Navier-Stokes equation to a regularized lattice [31]. Stam expressed a more realistic liquid using semi-Lagrangian techniques [32]. Carlson et al. developed an Eulerian solver to simulate melting objects or very sticky liquids [33], and Goktekin et al. proposed a method of expressing viscoelastic objects in the Eulerian method using an object's rigidity [34]. The Eulerian method is difficult to process in reality owing to its high computation volume to express the advection phenomenon of a fluid. In computer games where real-time processing is essential, the Lagrangian method with less computation is mainly used compared with other methods [35]. The Lagrangian method simulates fluid particles directly by focusing on the particles themselves rather than on the environment around the fluid, such as advection. Tonnesen computed particles by a Lagrangian method to represent deformable objects [36]. Among the Lagrangian methods, many SPH techniques are used to efficiently perform inter-relational operations on all particles [37].
Camporredondo et al. utilized an SPH-based fluid to verify the exact pour simulation in a robotic manipulator [38]. The field value of particles in the kernel, such as density, pressure, and viscosity, is calculated as the average value of the field of particles overlapping with the neighboring kernel. The kernel is mainly used by poly6 [37] and as a visibility kernel or a spiky kernel [39]. Navarro et al. used a particle-based viscoelastic fluid to express the skeletal muscle moving due to electrical activity [40]. Ruffini et al. proposed a model that combines DualSPHysics and CHRONO to measure the impact of dams when flooding occurs [41]. Currently, SPH is used in many fluid fields, including fluid-solid strain simulations [42] and melting fluid simulations [43]. Studies have been conducted to express elastic objects using SPH techniques. Clavet et al. proposed a yield condition in which an object deforms under external forces above its limit to be elastic, or otherwise returns to its original form [44].
By approximating the yield condition proposed by von Mises, the original position remains unchanged if the travel distance of the extruded particle is less than the limit, and if the travel distance of the particle is greater than the limit, the elasticity is expressed by varying the travel distance. Unlike conventional von Mises yield condition calculations, where all the forces an object receives in particle-based simulations must be calculated, this method can express elasticity with a small computation because it expresses the limit condition of the object's deformation only at a distance shifted by the force it receives. However, applying the yield condition through the distance traveled by a particle can make real-time simulation difficult if the number of particles increases, as the same force is applied in all directions.

Problems with Traditional Weight-Based Skinning Techniques
Actual muscles surround bones and soft tissue and are restricted by ligaments and skin. When expressing complex motions, connecting a simple straight line between bones alters a lot of information about the actual muscle. However, there is also a problem that a lot of calculation is needed to reflect the physics of all the torque the muscles receive between the bones. LBS approximates the geometry of muscles as soon as joints move, and is used to compute the transformation matrix and bone-based weight values of each bone to increase accuracy and flexibility.
LBS uses a method of linearly interpolating the transformation matrix of all skeletons which deforms the position of the surface. By using the result, the position of each surface of an object is calculated. Equation (1) is a formula for calculating the location of a deformed surface in the world coordinate system by linearly interpolating the transformation matrix from the local coordinate system to the world coordinate system of each skeleton. For the position of the deformed surface in the world coordinate system, the T j of the transformation matrix of skeleton j is multiplied by the weights w i,j of the skeleton j on the surface i. The m in Equation (1) is the total number of skeletons that affect the position of the surface, and adding all of the transformation matrices multiplied by each of their weights results in a transformation matrix for the motion of the deformed surface in the world coordinate system. Multiplying this matrix by the original position v i of the surface in the local coordinate system, we calculate the deformed motion of that surface in the world coordinate system.
LBS quickly calculates the motion of the surface deformed by the skeleton. LBS is used in various animation businesses owing to its fast calculations and smooth motions, but candy-wrapper artifacts that twist mesh data during rotational motion occur, resulting in joint merging when expressing bending motion. To solve this problem, DQS techniques allow for smoother character movements to be represented. However, the weight-based character movements do not calculate the force on the muscles, so in situations where the character is subjected to external force, additional weight customization changes the mesh data.

Muscle Pressure Control Using SPH Fluid Simulation
Human muscles perform body movements according to contraction and relaxation. The usual muscles allow actin and myosin, which are its components, to maintain ideal distance through the tension between each other and prepare for efficient exercise. During exercise, the muscles contract by an external force, but when the external force is removed from the contracted state, the muscles return to their original ideal distance by elasticity. These muscle properties are similar to those of surface tension in water particles.
In this study, we simulate muscle contraction by external forces by applying a particlebased fluid technique, SPH, to the skinning simulation. Each particle has a kernel to maintain the tension between the other particles coming in the range.
As shown in Figure 1, SPH uses a radial symmetric smoothing kernel, where the field values of the particles in the kernel (such as density, pressure, and viscoelasticity) are interpolated by the field values of the particles overlapping the neighboring kernel. Equation (2) is a formula that approximates an arbitrary scalar value A using SPH. Each of m j , p j , and A j mean the mass, density, and arbitrary scalar value A of particle j within the kernel region. M is the total number of particles in the simulation. At r position, the scalar field value A is used to calculate the force applied to a particular particle by multiplying and adding weights to all particles within the kernel region with radius h. The kernel function W, which is widely used in fluid simulations, uses poly6 [37] and viscosity kernel or spiky [39] kernel.

Animation Restrictions Using Tresca's Yield Condition
Muscle tissue is destroyed internally when a strong external force is applied, so it requires conditions under which destruction occurs. We utilize Tresca's yield condition, the concept of maximum stress theory that yield occurs when the maximum stress that occurs on an object equals the maximum stress of the tensile stress test. Figure 2 shows Tresca's yield condition as a graph in 3-dimensional space.

Animation Restrictions Using Tresca's Yield Condition
Muscle tissue is destroyed internally when a strong external force is applied, so it requires conditions under which destruction occurs. We utilize Tresca's yield condition, the concept of maximum stress theory that yield occurs when the maximum stress that occurs on an object equals the maximum stress of the tensile stress test. Figure 2 shows Tresca's yield condition as a graph in 3-dimensional space.

Animation Restrictions Using Tresca's Yield Condition
Muscle tissue is destroyed internally when a strong external force is applied, so it requires conditions under which destruction occurs. We utilize Tresca's yield condition, the concept of maximum stress theory that yield occurs when the maximum stress that occurs on an object equals the maximum stress of the tensile stress test. Figure 2 shows Tresca's yield condition as a graph in 3-dimensional space.  Tresca's yield condition formula is associated with principal stress, which means that stress occurs owing to external forces inside the object. It is expressed as one-half the difference between the maximum and minimum principal stresses and can represent the yield conditions of various objects [45].
Values of an acceptable range under conditions of ductile materials utilize Tresca's yield condition in Equation (3) because von Mises is stable but computationally complex.
σ 1 , σ 2 , and σ 3 are principal stress components that are orthogonal to each other, and the yield condition is when 1/2 of the value at which the difference between the components of both principal stresses is largest is the same as τ f . When pressure surges owing to external forces, muscle tissue accumulates the force of reaction to return to its original state, resulting in loss of muscle tissue when it is greater than the external force applied.
We assume a muscle tissue that maintains a constant density and define the moment when it fails to resist and a defect occurs when it is disturbed by a very strong external force as the muscle yield condition. A form with minimum principal stress is assumed to be a steady state that maintains incompressibility.
We proposed the yield condition expression for body elasticity in Equation (4). ρ min is the density of an object with minimum principal stress, and ρ max is the density of an object with maximum principal stress. Because each muscle particle has a kernel in a certain range and maintains tension between other particles within the kernel, they can be used as a criterion for yield conditions.
Unlike distance-based approximation, if the same force is received in multiple directions of the object, the total change of the object is calculated, rather than the change of a single particle. In this work, we propose an algorithm for the body's shape change for changes in elasticity when shocked from the outside.
The sequence of the algorithm is as follows. First, the distance between particles is calculated and compared with the initial value. The initial value is the density of the object to be tested by the user. If the distance between particles is less than the initial value, the density is calculated by external forces. If the calculated density is greater than the destruction condition, we do not increase the density by more than that. As a result, the size of the principal stress is assumed to be an ideal form, indicating realistic body movements.

Experimental Results
The experiment was performed on a system with 32 GB of main memory on an Intel Core i5 8500 3.00 GHz central processing unit (CPU). The graphics processing unit (GPU) used was the NVIDIA GeForce RTX 3070 with 8 GB of memory and was implemented in a multiprocessor environment utilizing Unity and data-oriented tech stack (DOTS) scripts. Figure 3 shows the expected outcome of the proposed method when the two skeletons are connected and the elbow is flexed by external forces. Each particle's color has an initial gray value and owing to the density and external force between the particles, it becomes red as the pressure increases. When an external force acts and causes movement in the skeleton, each particle determines its motion according to the magnitude and direction of the external force it receives. Initially, the body is bent in the direction in which the force is applied, and then slowly stretches out due to the SPH smoothing kernel, which tries to maintain the same distance as the force in the range. The corresponding experimental results allow us to simulate changes in the force exerted on the actual human body. The elasticity of SPH's smoothing kernel can confirm changes in a movement when external forces are applied to the body, but the body has complex components to express only as changes in particles. The stronger the applied force, the worse the flexion will be, but the actual person is no longer bent in a particular direction unless there is a defect in the substance due to the properties of the joints and muscles.
To apply this characteristic, we assume a temporary pressure of 180 degrees, which is the maximum bendable angle with respect to the elbow, as the maximum principal stress, and the yield condition for the pressure that can be applied is determined by assuming a minimum of principal stress when no force is applied. As shown in Figure 4, Experimental results indicate a process in which the body is no longer bent when external forces are applied in a particular direction. As a result, we were able to solve the problem caused by rotating the body without considering the physics applied to the body in weight-based skinning algorithms, such as candy-wrapper phenomena or bulging phenomena. The elasticity of SPH's smoothing kernel can confirm changes in a movement when external forces are applied to the body, but the body has complex components to express only as changes in particles. The stronger the applied force, the worse the flexion will be, but the actual person is no longer bent in a particular direction unless there is a defect in the substance due to the properties of the joints and muscles.
To apply this characteristic, we assume a temporary pressure of 180 degrees, which is the maximum bendable angle with respect to the elbow, as the maximum principal stress, and the yield condition for the pressure that can be applied is determined by assuming a minimum of principal stress when no force is applied. As shown in Figure 4, Experimental results indicate a process in which the body is no longer bent when external forces are applied in a particular direction. As a result, we were able to solve the problem caused by rotating the body without considering the physics applied to the body in weight-based skinning algorithms, such as candy-wrapper phenomena or bulging phenomena.  Figure 5 shows the result of a speed comparison between the proposed method for the elbow bending process and the existing skinning technique. In the figure, the proposed method shows slower results than DQS but can express the interactions on external forces on the body compared with weight-based skinning techniques that cannot express the interactions on forces. Furthermore, the results for similar elasticity and bending can be obtained with less computation compared with position-based dynamics, which expresses the interaction by calculating the detailed force the muscle receives when the external force of the elbow is applied.  Figure 5 shows the result of a speed comparison between the proposed method for the elbow bending process and the existing skinning technique. In the figure, the proposed method shows slower results than DQS but can express the interactions on external forces on the body compared with weight-based skinning techniques that cannot express the interactions on forces. Furthermore, the results for similar elasticity and bending can be obtained with less computation compared with position-based dynamics, which expresses the interaction by calculating the detailed force the muscle receives when the external force of the elbow is applied.  Table 1 shows the difference between fps and bending angles when the number of particles is adjusted under the same body part and external force conditions. When external force was applied to the particle based on the elbow model in Figure 3, the result of the proposed yield condition was determined at what angle the elbow would stop. We also wanted to derive an efficient experimental environment by comparing frames. The 500 particles showed a very high speed, but the density was not significant when external forces were applied, so the yield conditions determined by elasticity and density were not properly applied. The 1500 particles were found to be subject to strong elasticity and bending restrictions owing to the high density of the particles when external forces were applied, but had the disadvantage of being slow. On the other hand, 1000 particles were measured with elasticity and yield conditions similar to those of 1500 particle operations, but their speed increased significantly.

Conclusions
In this paper, we propose an efficient way to express the elasticity of the body's forces in skinning simulation. Conventional weight-based skinning simulation is not in the form of calculating the pressure and density of muscle fibers generated by muscle contraction in the body, so it cannot express elasticity when subjected to an external force. Moreover, to express the same motion as a person even in the usual situation where there is no power, one has to arbitrarily customize the weight value depending on the situation. Pressure and density values concentrated on SPH-based kernels are utilized to express elasticity to external forces. In addition, Tresca's yield condition-based algorithm is applied to limit elasticity to areas of the body that are under heavy force when the body is moved by external forces, similar to reality. As a result, it is possible to confirm that conventional weight-based algorithms do not detect rotational twists that usually occur. Compared with the weight-based approximation method, we were able to calculate the exact amount of force a muscle receives and verify its ability to bounce back by elasticity. However, the  Table 1 shows the difference between fps and bending angles when the number of particles is adjusted under the same body part and external force conditions. When external force was applied to the particle based on the elbow model in Figure 3, the result of the proposed yield condition was determined at what angle the elbow would stop. We also wanted to derive an efficient experimental environment by comparing frames. The 500 particles showed a very high speed, but the density was not significant when external forces were applied, so the yield conditions determined by elasticity and density were not properly applied. The 1500 particles were found to be subject to strong elasticity and bending restrictions owing to the high density of the particles when external forces were applied, but had the disadvantage of being slow. On the other hand, 1000 particles were measured with elasticity and yield conditions similar to those of 1500 particle operations, but their speed increased significantly.

Conclusions
In this paper, we propose an efficient way to express the elasticity of the body's forces in skinning simulation. Conventional weight-based skinning simulation is not in the form of calculating the pressure and density of muscle fibers generated by muscle contraction in the body, so it cannot express elasticity when subjected to an external force. Moreover, to express the same motion as a person even in the usual situation where there is no power, one has to arbitrarily customize the weight value depending on the situation. Pressure and density values concentrated on SPH-based kernels are utilized to express elasticity to external forces. In addition, Tresca's yield condition-based algorithm is applied to limit elasticity to areas of the body that are under heavy force when the body is moved by external forces, similar to reality. As a result, it is possible to confirm that conventional weight-based algorithms do not detect rotational twists that usually occur. Compared with the weight-based approximation method, we were able to calculate the exact amount of force a muscle receives and verify its ability to bounce back by elasticity. However, the error values were obtained because they are not simulated by implementing accurate myofibers and the yield conditions for each part of the body are different from muscle to bone.