# yade-users team mailing list archive

## [Question #708272]: i.geom.penetrationDepth = nan value when saving ViscElCap interactions

 Thread Previous • Date Previous • Date Next • Thread Next

```New question #708272 on Yade:

Hello, I need to save then reload simulation with ViscElCap interaction inside but when I reload it the penetration Depth of ViscElCap interaction is equal to nan. This  nan value contaminates position and velocity values making disappear particles.

I softly fix that by erasing and recreate each ViscElCap interaction, but obviously it didn't conserve friction history which is important to me.

Maybe someone know why and how to fix this ?

I'm using yade 2022.01a

You can test that with :

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# import
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

import sys
import numpy as np
from yade import plot, qt

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#DEF
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

s_name = sys.argv[0] [0:-3]

#sim param --------------------------------------------------------------------

p_speed = 0.01

tc_frac = 0.0002

#Mechanics --------------------------------------------------------------------

#Value in order :

R1 = 19 * 1e-6
R2 = 27.5 * 1e-6 #35 #32.5 #27.5

#radius = 2 * R1 * R2 /(R1 + R2)
radius = 0.0002 #0.002381 #2e-4

#Law {NoCoh,LiCoh,SoCoh}
Law = 'ViscElCap' #CohFrictMat

#particles
p_density = 2650 #2650
f_density = 1000

#Modulus
E_p = 1e8 #1e8

#poisson
P_p = 0.2

#friction angle
fa_p =  0.5

#restit
rcn_p = 0.7
rct_p = 0.7

#cohesion ---------------------------------------------------------------------

#contact time
tc = 0.001

#surface tension
Gamma = 72.8 * 1e-3 #27 #24 # 28

#Wetting angle [°]
Theta = 10  #0

#Liquid bridge volume def

#[%] of the Void volume filled with water
water_content = 0.1 #from [0.0001 ; 0.1]

Nl = 5.00481359595265 #6 #to precise

#VB = 36 * 1e-19 #2 #12 #36
#VB = 2 * water_content * Vviod / Nl
VB = water_content * np.pi * p_density * radius**3 / ( 3 * Nl * f_density )

#seed type=int (no seed = -1)
seed_gen = 1

#Liquid Coh
CapModel = ['Willett_numeric',
'Willett_analytic',
'Weigert',
'Rabinovich',
'Lambert',
'Soulie',]

SolidModel = ['']

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#Pre compute
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#create result folder if no existing

#print
print()

print('Bridge volume : ' + str(VB) )

print('\n\n\n\n\n\n\n\n\n\n')

#r.seed(1)

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#Material def
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#------------------------------------------------------------------------------
#Cohesive (Liquid)
#------------------------------------------------------------------------------

for i in range(len(CapModel)) :
O.materials.append(ViscElCapMat(label = CapModel[i], frictionAngle = fa_p, mR = 0.1, mRtype = 1, en = rcn_p, et = rcn_p, tc = tc, Capillar = True, gamma = Gamma, Vb = VB, theta = 10, CapillarType = CapModel[i] ))

nbCm = len(CapModel)
Models = CapModel

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#Engine def
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#Engine def -------------------------------------------------------------------

O.engines = [ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),

#Collision detection
InteractionLoop([Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],

[Ip2_FrictMat_FrictMat_FrictPhys(),

#Liquid
Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys()],

[Law2_ScGeom_FrictPhys_CundallStrack(),

#liquid
Law2_ScGeom_ViscElCapPhys_Basic()]
#[Law2_ScGeom_CapillaryPhys_Capillarity],

),

#GlobalStiffnessTimeStepper(),
NewtonIntegrator(gravity=(0, 0, 0), damping = 0),
]

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#Cohesion law test
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#Init scene -------------------------------------------------------------------

sep_exp = 3 * radius

x_og_pos = 2 * radius
y_exp_pos = 0

fp_id = []
mp_id = []

for cm in range(nbCm) :

#1
fp_id.append( O.bodies.append([sphere(center=(0, y_exp_pos, 0),
material = Models[cm],
fixed = True)])[0] )

#2

mp_id.append( O.bodies.append([sphere(center=(x_og_pos, y_exp_pos, 0),
material = Models[cm],
fixed = True)])[0] )

O.bodies[mp_id[cm]].state.vel = [p_speed,0,0]

y_exp_pos += sep_exp

#Define Time-step
O.dt = tc_frac * tc

#Energy tracking
O.trackEnergy = True

#Init Var ---------------------------------------------------------------------

ntt_iter = 0
last_time = 0

time = np.full([1],0)

#O.run(1,True)

first = True

#open graphique ctrl
qt.View()

plot.plots = {'Displacement' : ( Models[0], Models[1], Models[2], Models[3], Models[4], Models[5] )}

# show the plot on the screen, and update while the simulation runs
plot.plot()

###############################################################################
################################### Compute ###################################
###############################################################################

O.run(1000, True)

O.save('./test')