Apply Rayleigh damping in wave propagation problem

I assume you want to introduce C \times v_i as an explicit term in the right hand side?

Thanks, Lutz. Yes, I wish to add the damping force term on the right-hand side. Is there a way to do this?

This is a bit tricky. I assume that you are planning to use lumping. in this case you need to make sure that M in the dumping term becomes in fact the lumped version of the mass matrix. Your equation looks like.

M a + (alpha M + \beta K) v + K u = f

where v and u are from the previous time step. You can rearrange this in the the form
M (a+alpha*v)+ K (u+beta * v) = f

which basically means that you calculate \sigma with (u+\beta*v) instead of just u and
then solve

M w + = f-div(sigma)

and then set

a = w-alpha*v

Hope that makes sense (check signs)!

That's brilliant. But it seems the Newmann boundary is now a trouble. The boundary is subjected to a constant traction. How to adjust this traction in this case?

I'm trying to add the stiffness-related damping force term in Y_i as follows. But it seems not working.

damping = tensor_mult(C,symmetric(grad(beta*vel))) # where C is the elastic modulus tensor
proj = Projector(domain)
pde.setValue(Y=-trace(grad(proj(damping)),1),X=-stress) # where stress tensor is calculated with grad(u) only

Is there any problem in this piece of code? Thanks.

One could make it a bit more efficient but I cannot see anything wrong.

I am not en expert in dumping but it is my understanding that a condition of the for KM^{−1}D = DM^{−1}K (D= (alpha M + \beta K)) needs to be fulfilled to make the dumping effective. I don't think your version meets this requirement.
I can see your argument with boundary condition but I would argue that the consistence condition tells you that you also need to apply dumping the boundary conditions (to deal with reflections?). A comparison will be interesting!?