# differentiate between two material

Asked by Mithushan Soundaranathan

Hi,

I have two materials (same material type, but with different size and density), defined as matPH101 and matCCS. I am simulating the swelling (radius) of these materials. Each material have different swelling rate (see def model in the code) and depending on their parameters (P_PH101 and P_CCS). I am updating each particle radius individually through a for loop with:
for b in O.bodies:
if isinstance(b.shape, Sphere):
My question is, how can I differentiate between each material, by adding an if statement. I would like have something like this:
for b in O.bodies:
if isinstance(b.shape, Sphere):
if materials is matCCS
P=P_CCS #using parameter for CCS
elif materials is matPH101
P=P_PH101 #using parameter for PH1101

I have parts of my code, I was sure which part to add. If you need more info please let me know.

Best regard,
Mithu

from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd
import numpy as np
from PIL import Image
from yade import pack, export
from scipy.interpolate import interp1d
from csv import writer
import os
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import csv
from matplotlib.pyplot import figure
from pylab import *
from scipy.optimize import curve_fit
from yade.params.table import *
import scipy.spatial

o = Omega()
save=save
# Physical parameters
fr = 0.41
rho_PH101 = 1561
rho_CCS =1403
D_PH101 = 7.9e-5
r1_PH101 = D_PH101/2
D_CCS = 5.4e-5
r1_CCS = D_CCS/2
#r2 = Diameter/2
k1 = 10000
kp = 140000
kc = k1 * 0.1
ks = k1 * 0.1
Chi1 = 0.34

if save==0:
g=-9.81
elif save==1:
g=0

o.dt = 1.0e-8
particleMass_PH101 = (4.0/3.0)*math.pi*r1_PH101*r1_PH101*r1_PH101*rho_PH101
particleMass_CCS = (4.0/3.0)*math.pi*r1_CCS*r1_CCS*r1_CCS*rho_CCS
m_tab_PH101=m_tab*1e-3*(1-wCCS)
m_tab_CCS=m_tab*1e-3*wCCS
tab_no_p_PH101=m_tab_PH101/particleMass_PH101
tab_no_p_CCS =m_tab_CCS /particleMass_CCS
r_tab=r_tab*1e-3 #real size
h_tab=h_tab*1e-3
v_tab=math.pi*(r_tab**2)*h_tab
no_p_PH101=(v_1mm/v_tab)*tab_no_p_PH101
no_p_CCS=(v_1mm/v_tab)*tab_no_p_CCS

PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

Cyl_height=0.006

Comp_press_up= comp_press
Comp_force_up=Comp_press_up*cross_area

Comp_press_lp= comp_press
Comp_force_lp=Comp_press_lp*cross_area

#*************************************
compression_data_save=[]
sc_por_15=2
#sc_por_2=2
#sc_por_1=1
rho_mix=((wCCS/rho_CCS)+((1-wCCS)/rho_PH101))**-1
data_to_save=[comp_press/1e6,round(no_p_PH101)+round(no_p_CCS),rho_mix]
compression_data_save.append(data_to_save)

matPH101 = O.materials.append(LudingMat(frictionAngle=fr, density=rho_PH101, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))
matCCS = O.materials.append(LudingMat(frictionAngle=fr, density=rho_CCS, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression and walls
sp=pack.SpherePack()
sp.makeCloud((-8*D_PH101,-8*D_PH101,-35*D_PH101),(8*D_PH101,8*D_PH101,35.0*D_PH101), rMean=r1_PH101,rRelFuzz=0.18,num=round(no_p_PH101))
n1 = len(sp)
sp.makeCloud((-8*D_PH101,-8*D_PH101,-35*D_PH101),(8*D_PH101,8*D_PH101,35.0*D_PH101), rMean=r1_CCS,rRelFuzz=0.15,num=round(no_p_CCS))
for i,(c,r) in enumerate(sp):
mat = matPH101 if i < n1 else matCCS
color = (0,1,1) if i < n1 else (1,0,1)
O.bodies.append(sphere(c,r,material=mat,color=color))

##Single particle swelling model
def model(r,t,Q_max,rho_t,rho_w,r_0,Diff):
Q=((rho_w*(r**3))/(rho_t*(r_0**3)))-(rho_w/rho_t)+1;
drdt =((Diff*rho_t)/(r*rho_w))*((Q_max-Q)/Q);
return drdt
P_PH101=[1.45,rho_PH101,1000,396.39e-12]
P_CCS=[3.16, rho_CCS, 1000, 739.75e-12]

def ParticleSwelling():
time_current=(o.iter-initial_save[0])*o.dt
if wCCS==0.02:
Liq_pos=0.1871*(time_current**0.7339) #mm
elif wCCS==0.05:
Liq_pos=0.1359*(time_current**0.8279) #mm
elif wCCS==0.08:
Liq_pos=0.0968*(time_current**0.8739)
Liq_pos=Liq_pos*1e-3 #convert to m
print(Liq_pos)
print(time_current)
fw=[]
z_min=utils.aabbExtrema()[1][2]
for b in O.bodies:
if isinstance(b.shape, Sphere):
par_pos=initial_save[1]-b.state.pos[2]
k=b.id
if Liq_pos>=par_pos:
r_0=r_save[0][k+1]
if swell_t[0][k]==0:
swell_t[0][k]=time_current
continue
time=time_current-swell_t[0][k]
t = np.linspace(0,time)
r = odeint(model,r_0,t,args=(P[0],P[1],P[2],r_0,fw_D[0][k],))
r_new=float(r[-1])
f=float(r[-1])/r_now
mcurrent=b.state.mass
mnew=mcurrent*(f*f*f)
b.state.mass=mnew
Icurrent=b.state.inertia
Inew=Icurrent*(f*f*f*f*f)
b.state.inertia[0]=Inew[0]
b.state.inertia[1]=Inew[1]
b.state.inertia[2]=Icurrent[2]
z_current=b.state.pos[2]
if z_min>z_current:
z_min=z_current
else:
z_min=z_min
elif Liq_pos<par_pos:

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Jan Stránský (honzik) said on 2022-05-10: #1

Hello,

"mat is b.mat" is what you are looking for.
Alternatively, if the materials O.materials.appended, you can use material ID:
### MWE
mat1 = FrictMat()
mat2 = FrictMat()
b1 = sphere((0,0,0),1,material=mat1)
b2 = sphere((0,0,0),1,material=mat2)
#
print(mat1 is mat1) # True
print(mat2 is mat2) # True
print(mat1 is mat2) # False
print(b1.mat is mat1) # True
print(b1.mat is mat2) # False
print(b2.mat is mat1) # False
print(b2.mat is mat2) # True
#
print(mat1.id) # -1
print(mat2.id) # -1
print(b1.mat.id) # -1
print(b2.mat.id) # -1
#
O.materials.append(mat1) # to assign id
O.materials.append(mat2) # to assign id
print(mat1.id) # 0
print(mat2.id) # 1
print(b1.mat.id) # 0
print(b2.mat.id) # 1
###

Cheers
Jan