
As NCCAChris mentioned, you can directly damp velocity every time you have a collision. This will be a duplicate damping depending on your coeffRest.
For more flexibility, you must compute friction/drag in the plane of contact. You already have the projection equation:
// From:
newVel = vel  (1+coeffRest)*n*(n dot vel) // Reflect velocity.
// Velocity in plane:
velInDirOfPlane = n*(n dot vel)
velInPlane = vel  velInDirOfPlane

You can apply simple linear drag as:
newVel = vel  velInPlane*coeffDrag // coeffDrag ranges from 0..1

If you want to implement Coulomb friction, you'll have two coefficients, coeffStaticFriction and coeffDynamicFriction. Here's David Baraff's approximation method, which works OK:
http://graphics.stanford.edu/courses/cs44801spring/papers/baraff.pdf http://www.gamasutra.com/features/20000510/lander_03.htm
Another way is to think about everything in terms of momentum and impulses instead of forces (the velocity reflector code is actually a momentum reflector (when mass = 1)).
Mosts references for friction use:
Fs = muS*FN // static friction force Fk = muK*FN // kinetic (dynamic) friction force
Where FN is the normal force and mu* are coefficients of static and kinetic friction.
If the perpendicular force Fp is muS*FN, we use the kinetic (dynamic) friction equation (slowing down the moving particle).
Given:
F = ma
We can also use:
F = d(mv)/dt // were d means delta, so d(mv) is change in momentum.
If only the velocity is changing:
F = m*dv/dt // remember a = dv/dt
The impulse equation:
(mv)1 + F*dt = (mv)2 // (mv) written this way to show that mass and/or velocity can change.
or
F*dt = (mv)2  (mv)1 // the impulse is a change in momentum.
and from F = m*dv/dt
F*dt = m*dv
In the case when m is always 1 (your particle simulation)
F*dt = dv // impulse.
If you think about everything in terms of momentum and impulses (instead of forces), your simulation work becomes easier.
The basic idea is that you have a routine called pointMomentumInDir(), which computes the momentum of a particle in a direction:
// normalDir must be unit length (normalized).
Vec3 pointMomentumInDir(Vec3 vel,float mass,Vec3 normalDir,float coeffRest) {
return mass*(1+coeffRest)*normalDir*(normalDir dot vel)
} // pointMomentumInDir

You should recognize this equation from the previous example. You can use this equation for collisions with static objects, as well as for friction.
For example:
// particleCoeffRest can be substituted with 0 if you don't want restitution
// to affect friction (or your can use different variables for more tuning ability).
velDirProjectedToGround = normalize(vel  groundNormalDir*(groundNormalDir dot vel))
momentumInDirOfGroundNormal = pointMomentumInDir(particleVel,particleMass,groundNormalDir,particleCoeffRest)
momentumPerpToGroundNormal = pointMomentumInDir(particleVel,particleMass,velDirProjectedToGround,0)
if (length(momentumPerpToGroundNormal) <= length(momentumInDirOfGroundNormal)*coeffStaticFriction) {
particleVel = particleVel  momentumPerpToGroundNormal/particleMass // Stop all motion in the ground plane.
} else { // Apply friction proportional to clamping.
particleVel = particleVel  velDirProjectedToGround*length(momentumInDirOfGroundNormal)*coeffDynamicFriction/particleMass
} // if

You can optimize the previous method using the following routine:
// normalDir must be unit length (normalized).
float pointMomentumInDirLength(Vec3 vel,float mass,Vec3 normalDir,float coeffRest) {
return mass*(1+coeffRest)*(normalDir dot vel)
} // pointMomentumInDirLength

You can optimize further when mass is always equal to 1 (if mass varies, store 1/mass as massInv and multiply instead of divide).
