Dear all,
In my project I am exploring ways to predict risk of blood clotting in patient specific geometries, using FEM.
For this application, I am currently using the Oasis solver, implemented in FEniCS, to model blood flow.
Now, I found an interesting paper on continuum modelling of platelet stress and deformation, which could provide a way to quantify platelet activation throughout the domain. The local platelet activation, in turn, could be associated with risk of thrombosis.
The Oasis solver already has scalar transport implemented in its Incremental Pressure Correction solver.
In the above paper, scalar transport solved for the concentration of activated platelets. As a source term, they add a rate of strain/deformation tensor to the scalar transport equation.
Now, I wish I could say adding a term to the scalar transport equation in this solver would be trivial for me, but unfortunate my experience in this regard is lacking.
The transport equation they use in the paper is as follows, however, they neglect the diffusion term because kappa is small:
With phi the (scalar) platelet activation level. The last term in the equation is the local platelet activation rate (scalar stress), for which they use a linear stress-exposure time model.
The scalar stress is defined as the Frobenius norm of the rate of strain tensor:
In the Oasis solver, to implement the above we would need to add this last term (scalar stress) to the scalar transport equation. (and remove the already implemented diffusion term because we neglect this)
In IPCS_ABCN.py, the optimized incremental pressure correction scheme in the Oasis solver, the scalar transport equation is implemented with:
#M is the mass matrix
M = assemble_matrix(inner(u, v) * dx)
Ta=matrix(M)
a_conv = inner(v, dot(u_ab, nabla_grad(u))) * dx
a_scalar=a_conv
#Everything below is done for every timestep
A = assemble(a_conv, tensor=A)
A *= -0.5 # Negative convection on the rhs
A.axpy(1. / dt, M, True) # Add mass
Ta.zero()
Ta.axpy(1.,A,True) #because same convection used as in tentative velocity
Compute rhs
b[ci].zero()
b[ci].axpy(1., Ta * x_1[ci])
b[ci].axpy(1., b0[ci]) #with b0[ci] the scalar source term “assemble(v * fs[ci] * dx)”
Reset matrix for lhs - Note scalar matrix does not contain diffusion
Ta.axpy(-2., Ta, True)
Ta.axpy(2. / dt, M, True)
bc.apply(Ta, b[ci])
c_sol.solve(Ta, x_[ci], b[ci])
So what I would like to do is add the frobenius norm of the strain field in the scalar equation by replacing the source term:
strain=0.5 * (grad(u_) + grad(u_).T)
stress = sqrt(inner(strain,strain))
a_str = v * stress * dx
Tstr = assemble_matrix(a_str * dx)
Compute rhs for all scalars
Tstr=assemble(a_str)
b[ci].zero()
b[ci].axpy(1., Ta * x_1[ci])
b[ci].axpy(sqrt(2)*4e-3, Tstr) #This is the term above that is replaced.
For some reason I am not getting correct results with this implementation.
Are there any obvious mistakes any of you can point out or any direction you can give me?
I am looking forward to your suggestions.
With kind regards,