8.1 最小的离散元程序

[Zhao2015] 用 45 行MATLAB的离散元代码模拟了伽利略的比萨斜塔试验，其前处理、求解器和后处理代码 Code 7.1 :

代码 7.1 伽利略的比萨斜塔试验（MATLAB)

% function HelloDEM  % Galileo’s leaning Tower of Pisa experiment 
M_Ball_Steel=10.0; % 10 kg streel ball 
M_Ball_Wood=1.0;   % 1 kg wood ball 
Y0_Ball_Steel = 100.0; % initiall position 
Y0_Ball_Wood  = 100.0; 
dT     = 0.01; % time step 
Nloops = 500; % number of cycles 
G      = 9.8; % gravitational constant 
U_Ball_Steel = 0.0; % initiall displacement 
U_Ball_Wood  = 0.0; 
V_Ball_Steel = 0.0; % initiall velocity 
V_Ball_Wood  = 0.0; 
A_Ball_Stell = 0.0; % accelerated velocity 
A_Ball_Wood  = 0.0; 
Y_Ball_Steel = Y0_Ball_Steel; % initiall position 
Y_Ball_Wood  = Y0_Ball_Wood; 
t   = 0.0; 
t_H = [];% history of time  
Y_H_Ball_Steel = []; 
Y_H_Ball_Wood  = [];
for i=1:NLoops 
    F_Ball_Steel=-G*M_Ball_Steel; %ball force,constitutive model 
    F_Ball_Wood = -G*M_Ball_Wood; 
    A_Ball_Steel= F_Ball_Steel/M_Ball_Steel; %Newton’s second law 
    A_Ball_Wood = F_Ball_Wood /M_Ball_Wood; 
    V_Ball_Steel= V_Ball_Steel+A_Ball_Steel*dT; %update velocity 
    V_Ball_Wood = V_Ball_Wood+A_Ball_Wood*dT; 
    U_Ball_Steel=U_Ball_Steel+V_Ball_Steel*dT; %update displacement 
    U_Ball_Wood  = U_Ball_Wood+ V_Ball_Wood*dT; 
    Y_Ball_Steel = Y0_Ball_Steel + U_Ball_Steel; %update position 
    Y_Ball_Wood  = Y0_Ball_Wood  + U_Ball_Wood; 
    t=t+dT; 
    %record history 
    t_H=[t_H,t]; 
    Y_H_Ball_Steel=[Y_H_Ball_Steel,Y_Ball_Steel]; 
    Y_H_Ball_Wood=[Y_H_Ball_Wood,Y_Ball_Wood]; 
end
plot(t_H,Y_H_Ball_Steel,'k','linewidth',1); 
hold on; 
plot(t_H(1:20:end),Y_H_Ball_Wood(1:20:end),'ks','MarkerFace','y'); 
Ground=zeros(1,length(t_H)); 
plot(t_H,Ground,'r'); 
xlabel('Time [s]') 
ylabel('Height [m]') 
le=legend('Steel Ball','Wood Ball','Ground'); 

比萨斜塔离散元模拟试验中，钢球和木球高度随时间的变化结果见 Fig 7.1 ，可见钢球和木球是同时落地的。

图 8.1 球高度随时间的变化，比萨斜塔离散元模拟试验。

图 8.2 球高度随时间的变化，比萨斜塔离散元模拟试验。

离散元的求解实际上是迭代计算颗粒位移和受力， 可以概括为两部分，如 Fig 7.2 所示。第一步，已知的颗粒所受合力和合力矩，由牛顿第二定律更新每个颗粒的位置（如代码 Code 7.1 中，24~31 行） ；第二步，找到相互接触的颗粒，应用接触力学模型（即力-位移法则，或者称本构模型。如代码 Code 7.1 中，F_Ball_Steel=-G*M_Ball_Steel，该实例中，没有颗粒相互作用，颗粒仅受重力作用，实际这里需要计算颗粒间的相互作用力）计算颗粒间的接触力，进而得到所有颗粒受到的合力与合力矩。反复执行这两个步骤，直到计算结束。

[Zhao2015]

Zhao G-F (2015) High performance computing and the discrete element model: opportunity and challenge. Elsevier,