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. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp