yade-users team mailing list archive

[Mithushan Soundaranathan <question701747@xxxxxxxxxxxxxxxxxxxxx>]

New question #701747 on Yade:
https://answers.launchpad.net/yade/+question/701747

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):
            b.shape.radius=model (P)
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
          b.shape.radius=model (P)

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
readParamsFromTable(comp_press=0.782e8,h_tab=2.04,m_tab=0.2102, r_tab=5.02,wCCS=0.08,tab_porosity=15,tab_height=1,save=0)
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 
Tab_rad=0.001
r_tab=r_tab*1e-3 #real size
h_tab=h_tab*1e-3
v_tab=math.pi*(r_tab**2)*h_tab
v_1mm=math.pi*(Tab_rad**2)*(1e-3)
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
cross_area=math.pi*(Tab_rad**2)


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)

# Add material
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)
    radius=[]
    fw=[]
    z_min=utils.aabbExtrema()[1][2]
    radius.append(time_current) 
    for b in O.bodies:
        if isinstance(b.shape, Sphere):
            par_pos=initial_save[1]-b.state.pos[2]
            k=b.id
            r_now=b.shape.radius
            if Liq_pos>=par_pos:
                r_0=r_save[0][k+1]
                if swell_t[0][k]==0:
                    swell_t[0][k]=time_current
                    radius.append(b.shape.radius-r_save[0][k+1])
                    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])
                b.shape.radius = float(r[-1])
                radius.append(b.shape.radius-r_save[0][k+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:
                radius.append(b.shape.radius-r_save[0][k+1])

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.