Apply viscous forces on gridnodes

Hello everyone,

The model consists of grid nodes and grid-connection, and I am interested to implement the viscous component.
In the contact law, viscosity is implemented in Law2 ScGeom6D CohFrictPhys CohesionMoment only for the shear and twisting components and not for the normal component.
For that reason, I am considering implementing the viscosity by fictitious forces acting on each particle in such a way at each iteration and at each grid-node: F= f +f_viscous and f_ viscous= - c* velocity_of _particles.
By employing this approach, I intend to account for viscosity in both the shear as well as the normal component of particle interactions.

Do you have any idea how can we implement these viscous forces while optimizing computation time?

Thank you in advance!

Here is a simple and representative code:

from yade import pack, plot
from yade.gridpfacet import *
import numpy as np

#########################################
### DEFINE VARIABLES ###
#########################################

# Indicate the physical and mechanical properties of the problem:
# 1) Grains:
grainDensity=1e3 # particle density in [kg/m3]
grainYoung=5e7 # contact stiffness in [Pa]
grainPoisson=0.5 # poisson ratio
damp=0.5
# 2) Cohesion
Normal_cohesion=1.25e6
Shear_cohesion=1.25e6
radius=0.1
color = [255. / 255., 102. / 255., 0. / 255.]

#########################################
### DEFINE MATERIALS ###
#########################################
# 1) Gridnodes:
O.materials.append(CohFrictMat(young= grainYoung,poisson=grainPoisson, frictionAngle=radians(0), normalCohesion=Normal_cohesion, shearCohesion=Shear_cohesion,momentRotationLaw=False, fragile=True ,label='gridNodes'))
# 2) GridConnections:
O.materials.append(FrictMat(young=grainYoung,poisson=grainPoisson,frictionAngle=radians(0),density=grainDensity,label='gridConnections'))

#########################################
### DEFINE ENGINES ###
#########################################

# 1)
newton=NewtonIntegrator(gravity=(0,0,0),damping=damp) # Engine integrating newtonian motion equations.
# 2)
O.engines=[
        ForceResetter(), #Reset all forces stored. Typically, this is the first engine to be run at every step, Id=0:
        InsertionSortCollider([ # Id=1
         Bo1_Sphere_Aabb(),
         Bo1_GridConnection_Aabb(),
        ]),
        InteractionLoop([ # Id=2
         Ig2_GridNode_GridNode_GridNodeGeom6D(),
         Ig2_GridConnection_GridConnection_GridCoGridCoGeom()
        ],
        [
         Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True),
         Ip2_FrictMat_FrictMat_FrictPhys()
        ],
        [
         Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
         Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
         Law2_GridCoGridCoGeom_FrictPhys_CundallStrack()
        ]
        ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), # Id=1
 newton, # Id=2
    PyRunner(command='addPlotData()',iterPeriod=20000), # Id=3
    PyRunner(command='saveData()',iterPeriod=200000), # Id=4
]

#########################################
### Import the microstructure ###
#########################################

nodesIds=[]
connectionIds=[]
nodesIds.append( O.bodies.append(gridNode([1.5,0,1],radius,color=color,wire=False,fixed=True,material='gridNodes') ) ) #node1
nodesIds.append( O.bodies.append(gridNode([0.5,1,1],radius,color=color,wire=False,fixed=False,material='gridNodes') ) ) #node2
nodesIds.append(O.bodies.append(gridNode([1.5,2,1.5], radius, color=color, wire=False, fixed=False, material='gridNodes'))) #node3
nodesIds.append(O.bodies.append(gridNode([0.5,3,2], radius, color=color, wire=False, fixed=False, material='gridNodes'))) #node4
connectionIds.append(O.bodies.append( gridConnection(nodesIds[0],nodesIds[1],radius,wire=False,color=color,material='gridConnections') ))
connectionIds.append(O.bodies.append( gridConnection(nodesIds[1],nodesIds[2],radius,wire=False,color=color,material='gridConnections') ))
connectionIds.append(O.bodies.append( gridConnection(nodesIds[2],nodesIds[3],radius,wire=False,color=color,material='gridConnections') ))
O.bodies[nodesIds[-1]].state.blockedDOFs='yY'
O.bodies[nodesIds[-1]].state.vel=[0,0.0001,0]

#########################################
### SAVE DATA ###
#########################################

def addPlotData():

    xnodeMiddle=0
    ynodeMiddle=0
    ynodeTop=0
    if O.interactions[1,2].isReal==True:
        plot.addData(i=O.iter, t= O.time,
                force_cylinderTop=O.interactions[2,3].phys.normalForce.norm(),
                y_nodeTop=O.bodies[3].state.pos[1],
                y_nodeMiddle=O.bodies[2].state.pos[1],
                x_nodeMiddle=O.bodies[2].state.pos[0])
    else:
        O.interactions.eraseNonReal()
        O.pause()

def saveData():
 plot.saveDataTxt('viscosity_trial1.txt')
 plot.saveGnuplot('viscosity_trial1')

# Define what to plot:
plot.plots={'i':('force_cylinderTop'),'i ':('y_nodeTop','y_nodeMiddle','x_nodeMiddle')}

# Show the plot:
plot.plot()