Assignment 4 - Physically Based Modeling


Main Page Assignment 1 Assignment 2 Assignment 3



Physically Based Modeling

Two different physically based models have been performed for this assignment as follow:



A Particle-based model using Discontinuous Deformation Analysis



Discontinuous Deformation Analysis (DDA) is a type of discrete element method originally proposed by Shi in 1988 at UC Berkeley [1]. DDA is somewhat similar to the finite element method for solving stress-displacement problems, but accounts for the interaction of independent particles (blocks) along discontinuities in fractured and jointed rock masses. DDA is typically formulated as a work-energy method, and can be derived using the principle of Minimum Potential Energy (e.g., Shi) or by using Hamilton's principle. Once the equations of motion are discretized, a step-wise linear time marching scheme in the Newmark family is used for the solution of the equations of motion. The relation between adjacent particles is governed by equations of contact interpenetration and accounts for friction. DDA adopts a stepwise approach to solve for the large displacements that accompany discontinuous movements between particles.


Rigid particle displacements:

By adopting first order displacement approximations, the displacements (u,v,w) at any point (x,y,z) in a particle i, can be related, in two dimensions, to six displacement variables


where (u0,v0,w0) is the rigid body translation at a specific point (x0,y0,z0) within the particle, rx , ry , rz are the rotation angle of the particle with the rotation centre at (x0,y0,z0). As shown by Shi [1], the complete first order approximation of particle displacements takes the following form:




This equation enables the calculation of the displacements at any point within the particle (in particular, at the corners), when the displacements are given at the centre of rotation. In the three-dimensional formulation of the DDA method, the centre of rotation with coordinates (x0,y0,z0) coincides with the particle centroid.  

Simultaneous equilibrium equations:

In the DDA method, individual particles form a system of particles through contacts among particles and displacement constraints on individual particles. Assuming that N particles are defined in the particle system, Shi [1] showed that the simultaneous equilibrium equations can be written in matrix form as follows:


where each coefficient Kij is defined by the contacts between particles i and j. Since each particle i have six degrees of freedom defined by the components of Di in Eq. (16), each Kij in Eq. (4) is itself a 6x6 sub-matrix. Also, each Fi is a 6 x 1 sub-matrix that represents the loading on particle i. Eq. (4) can also be expressed in a more compact form as KD = F where K is a 6N x 6N stiffness matrix and D and F are 6N x 1 displacement and force matrices, respectively. In total, the number of displacement unknowns is the sum of the degrees of freedom of all the particles. It is noteworthy that the system of Eq. (4) is similar in form to that in finite element problems. The solution to the system of Eq. (4) is constrained by a system of inequalities associated with particle kinematics (e.g. no penetration and no tension between particles) and Coulomb friction for sliding along particle interfaces. The system of Eq. (4) is solved for the displacement variables. The final solution to that system is obtained as follows. First, the solution is checked to see how well the constraints are satisfied. If tension or penetration is found along any contact, the constraints are adjusted by selecting new locks and constraining positions and a modified version of K and F are formed  from which a new solution is obtained. This process is repeated until no tension and no penetration is found along all of the particle contacts. Hence, the final displacement variables for a given time step are actually obtained by an iterative process. The simultaneous Eq. (4) was derived by Shi [1] by minimizing the total potential energy Π of the particle system. The ith row of Eq. (4) consists of six linear equations


where the dri is the deformation variables of particle i. The total potential energy Π is the summation over all the potential energy sources, i.e. individual forces. The potential energy of each force and their derivatives are computed separately. The derivatives


are the coefficients of the unknowns dsj of the equilibrium Eq. (4) for variable dri. All terms of Eq. (6) form a 6 x 6 sub-matrix, which is the sub-matrix Kij in the global Eq. (4). Eq. (6) implies that matrix K in Eq. (4) is symmetric. The derivatives


are the free terms of Eq. (5) which are shifted to the righthand side of Eq. (4). All these terms form a 6 x 1 submatrix, which is added to the sub-matrix Fi. Shi’s thesis covers the details for forming sub-matrices Kij and Fi, for point loads, line loads, volume forces, bolting forces, inertia forces and viscosity forces [1]. Both static and dynamic analyses can be conducted with the DDA method. For static analysis, the velocity of each particle in the system at the beginning of each time step is assumed to be zero. On the other hand, in the case of dynamic analysis, the velocity of the system in the current time step is an accumulation of the velocities in the previous time steps.

collision search:

In the current code the computational domain is divided in uniform cubic cells of side. following. Thus, for a particle located inside a cell, only the interactions with the particles of neighboring cells need to be considered. In this way the number of calculations per time step and, therefore, the computational time diminish considerably, from N2 operations to N logN, N being the number of particles. Only 13 out of 26 possible neighboring cells in 3D are considered when centered on a particular ijk cell. The rest were previously considered when centered on adjacent cells using reverse interactions.


couple of years ago I did write subroutines of DDA for interaction of polyhedral blocks using Matlab. For this assignment same skeleton has been used with some changes inside the code.  the major issue that I faced here was memory allocation, unfortunately controlling the memory usage in Matlab got complicated for this assignment and then I would not be able to run models with more than 2500 particles and also running time tight my hands as well. The best solution for memory management would be rewriting the program with C++ that has more flexibility in term of memory allocation which would me my next step of work.




The following example shows a free fall model with 2500 same sized particles against two vertical glasses. The output file from Matlab code has been converted to "pov" format file using a third program and then POV-Ray (a free licensed renderer software) has been used to render the frames.





A Simple Cloth Simulation





In this tutorial I implemented a simple cloth simulator using Newtons second law, Verlet integration, and iterative constraint satisfaction using C++ and OpenGL.



Here cloth is simulated by particles with mass and interconnections, called constraints or springs. Particles move around in space due to forces (e.g. gravity and springs between particles) that affect them.

Force and position

For this assignment,  Force is translated into acceleration through Newtons second law (acceleration = force/mass) and then acceleration is “translated” into position by numerical integration – since position twice differentiated is equal to acceleration.



in order to constrain the particles, at each step all particles that need to be constrained have been returned to the same relative position, the more rigid the cloth will behave. During the verlet integration, particles move around, resulting in particles that are too far away from each other or too near each other. Therefore a constraint satisfaction has been introduced here to modify the position of the two particles so that their distance is once again "rest distance".

Creating cloth and collision

Gravity force has been added to each particles and collision also has been handled with a moving ball. Adding gravity is very easy, I simply added an acceleration vector pointing down to all particles. The ball is defined by a center and a radius. Detecting collision is done by asking whether the particle is within the radius of the ball, if yes the collision is resolved by moving the particle out of the ball along the vector from center to particle so that the distance is equal to the radius of the ball.




This example simply shows the collision between a moving sphere and a hanged fabric which has been programmed using C++ and also has been rendered by POV-Ray.





  • Matlab code for simulating particle based model.

  • A sample pov file created by Matlab code in particle based model.

  • C++ code for simple cloth simulation.

  • Matlab code for converting the output file to pov format in cloth simulation model.

  • A sample pov file created by C++ code in cloth simulation model.




I would like to thank Prof. Ravi and Jiamin for their help throughout the assignments.




For this assignment I did study several documents and available source codes that some major recourses have been listed below:

  1. Shi, G.-H. (1988). "Discontinuous deformation analysis-a new numerical model for the statics and dynamics of block systems," PhD thesis, Dept. of Civ. Engrg., Univ. of California, Berkeley.