diff --git a/Base/DynamicalModels/ObjectModels/imstkPbdModel.cpp b/Base/DynamicalModels/ObjectModels/imstkPbdModel.cpp
index df28e32d4c63b64c8aa9f52a726e39dc2ffafe1c..344c8f214dc2080f25f54312042b8498e21357e5 100644
--- a/Base/DynamicalModels/ObjectModels/imstkPbdModel.cpp
+++ b/Base/DynamicalModels/ObjectModels/imstkPbdModel.cpp
@@ -439,6 +439,19 @@ PbdModel::updatePbdStateFromPhysicsGeometry()
     m_currentState->setPositions(m_mesh->getVertexPositions());
 }
 
+void
+PbdModel::setViscousDamping(const double damping)
+{
+    if (damping >= 0 && damping <= 1)
+    {
+        m_viscousDampingCoeff = damping;
+    }
+    else
+    {
+        LOG(WARNING) << "WARNING - PbdModel::setViscousDamping:  Viscous damping coefficients is out of bounds [0, 1]";
+    }
+}
+
 void
 PbdModel::setUniformMass(const double& val)
 {
@@ -493,7 +506,7 @@ PbdModel::integratePosition()
         {
             vel[i] += (accn[i] + m_gravity)*m_dt;
             prevPos[i] = pos[i];
-            pos[i] += vel[i] * m_dt;
+            pos[i] += (1.0 - m_viscousDampingCoeff)*vel[i] * m_dt;
         }
     }
 }
diff --git a/Base/DynamicalModels/ObjectModels/imstkPbdModel.h b/Base/DynamicalModels/ObjectModels/imstkPbdModel.h
index 15deb3e616090d693e21fc0abcfd5a6309cdb770..4deaf1e697552d872c353696713df25958245c8b 100644
--- a/Base/DynamicalModels/ObjectModels/imstkPbdModel.h
+++ b/Base/DynamicalModels/ObjectModels/imstkPbdModel.h
@@ -169,6 +169,12 @@ public:
     ///
     void setGravity(const Vec3d& g) { m_gravity = g; };
 
+    ///
+    /// \brief Set/get viscous damping coefficient. Will be applied globally
+    ///
+    void setViscousDamping(const double damping);
+    double getViscousDamping() {return m_viscousDampingCoeff; }
+
     ///
     /// \brief Set uniform mass to all the nodes
     ///
@@ -209,20 +215,22 @@ protected:
     std::vector<std::shared_ptr<PbdConstraint>> m_constraints; ///> List of pbd constraints
 
     // Lame's constants
-    double m_mu; ///> Lame constant
-    double m_lambda; ///> Lame constant
+    double m_mu;                        ///> Lame constant
+    double m_lambda;                    ///> Lame constant
 
     // Mass properties
-    std::vector<double> m_mass; ///> Mass of nodes
-    std::vector<double> m_invMass; ///> Inverse of mass of nodes
+    std::vector<double> m_mass;         ///> Mass of nodes
+    std::vector<double> m_invMass;      ///> Inverse of mass of nodes
+
+    double m_contactStiffness;            ///> Contact stiffness for collisions
+    Vec3d m_gravity;                      ///> Gravity
 
-    double m_contactStiffness; ///> Contact stiffness for collisions
-    Vec3d m_gravity; ///> Gravity
+    double m_viscousDampingCoeff = 0.01;  ///> Viscous damping coefficient [0, 1]
 
-    unsigned int m_maxIter; ///> Max. pbd iterations
-    double m_proximity; ///> Proximity for collisions
+    unsigned int m_maxIter;               ///> Max. pbd iterations
+    double m_proximity;                   ///> Proximity for collisions
 
-    double m_dt; ///> Time step size
+    double m_dt;                          ///> Time step size
 };
 } // imstk
 
diff --git a/Examples/Sandbox/main.cpp b/Examples/Sandbox/main.cpp
index 839b1ba540d722ac23b1250f7a0b8b463cc24f47..28aed444737dfb51ff5eb1699739a80c11ff22df 100644
--- a/Examples/Sandbox/main.cpp
+++ b/Examples/Sandbox/main.cpp
@@ -3421,8 +3421,8 @@ int main()
     /*------------------
     Test physics
     ------------------*/
-    //testPbdVolume();
-    //testPbdCloth();
+    testPbdVolume();
+    testPbdCloth();
     //testPbdCollision();
     //testPbdFluidBenchmarking();
     //testPbdFluid();
@@ -3456,7 +3456,7 @@ int main()
     //testVectorPlotters();
     //testVirtualCoupling();
     //testBoneDrilling();
-    testVirtualCouplingCylinder();
+    //testVirtualCouplingCylinder();
 
 
     return 0;