左志華(在讀碩士) zuo.zhihua@qq.com
@深圳市羅湖區(qū) 2021.1.12
ThinkPad S2 2018 Windows10 Home
關(guān)鍵詞：SPH法；一維帝嗡；啟發(fā)幸海；B樣條核函數(shù)勺疼；Octave 6.1
1. 背景
在5月31號(hào)仗扬，我的“傻弟弟”室友不太會(huì)編程粥诫，有一道SPH習(xí)題不會(huì)做幔荒，當(dāng)時(shí)我給他寫(xiě)了一份Matlab代碼朝抖。
如今正式自己學(xué)SPH時(shí)發(fā)現(xiàn)啥箭，網(wǎng)上這方面的基礎(chǔ)資料還是挺匱乏的，所以找到了這份代碼槽棍，希望對(duì)初學(xué)者有啟發(fā)效果捉蚤。
2. 題目
拉格朗日型無(wú)網(wǎng)格粒子法在描述結(jié)構(gòu)毀傷行為時(shí)具有先天的優(yōu)勢(shì)。假定單位長(zhǎng)度細(xì)長(zhǎng)梁結(jié)構(gòu)如下圖所示炼七，梁的截面內(nèi)被離散為單個(gè)粒子，在初始狀態(tài)下梁內(nèi)變量布持，請(qǐng)選用三次樣條核函數(shù)豌拙，光滑長(zhǎng)度
，通過(guò)SPH方法中的粒子近似原理編程求解每個(gè)粒子處的變量值题暖。
（1）粒子分辨率和
情況按傅；
參考答案：（1）對(duì)于函數(shù)的粒子近似存在以下表達(dá)形式：
三次樣條函數(shù)為：
其中，
胧卤。
①粒子分辨率為情況：
支持域半徑為唯绍，粒子從左到右編號(hào)分別為1~6：
粒子1：與其相互作用的粒子編號(hào)：1，2枝誊，3
因此
粒子2：與其相互作用的粒子編號(hào)：1况芒，2，3叶撒，4
因此
粒子3：與其相互作用的粒子編號(hào)：1绝骚，2，3祠够，4压汪，5
因此
3. 代碼
% func B-Spline_w 是B-樣條函數(shù)(光滑核函數(shù)).
function w = B_Spline_w(r,h)
r = abs(r);
alpha = 1/h;
q = r/h;
if q>=0 & q<1
w = alpha*(2/3-q^2+0.5*q^3);
elseif q>=1 & q<2
w = alpha*(1/6*(2-q)^3);
elseif q>=2
w = 0;
else
fprintf("q = %d 是否小于了0; 請(qǐng)檢查! ",q);
end
end
% func partcleApproximation 粒子近似值函數(shù).
function f = particleApproximation_f(i,tol,h,dx)
A1 = 1; A2 = 1; A3 = 1; A4 = 1;
if i == 1
A1 = 0; A2 = 0;
elseif i == 2
A1 = 0;
elseif i == tol
A3 = 0; A4 = 0;
elseif i == tol-1
A4 = 0;
else
% 此處不執(zhí)行任何語(yǔ)句.
end
f = A1*dx*(i-3)*B_Spline_w(2*dx,h)*dx+A2*dx*(i-2)*B_Spline_w(dx,h)*dx+...
dx*(i-1)*B_Spline_w(0,h)*dx+A3*dx*i*B_Spline_w(-dx,h)*dx+...
A4*dx*(i+1)*B_Spline_w(-2*dx,h)*dx;
end
% SPH(Smoothed Particle Hydrodynamics)法 是光滑粒子流體動(dòng)力學(xué)方法的縮寫(xiě).
% 描述：
% 一維無(wú)網(wǎng)格實(shí)例
%
% 方法：
% B樣條核函數(shù)
%
% 代碼所有者：
% 左志華[zoziha] zuo.zhihua@qq.com
%
% 版本 日期 注釋
% ---- -------- -------
% V1.0 2020-6-1 原始版本 左志華
%
% 軟件標(biāo)準(zhǔn)：Octave 6.1
%
close all; % 關(guān)閉所有圖像.
% float dx 是粒子分辨率; 下語(yǔ)句 并且賦值 dx=0.2. [2020-6-1 wrote; 可修改 dx 值]
dx = 0.2;
% float tol 是粒子總數(shù). [2020-6-1 wrote; 可修改 tol 值]
tol = 6;
% float h 是光滑長(zhǎng)度.
h = 1.3*dx;
% float r 是域半徑.
r = 2*h;
% 下語(yǔ)句 在屏幕打印出 f
for i = 1:tol
f(i) = particleApproximation_f(i,tol,h,dx);
end
disp(f);
plot(f,"rp");
grid on;
4. 解析
可以發(fā)現(xiàn)SPH法在邊界處的精度不高，在求解域內(nèi)精度較高古瓤。
如果要開(kāi)發(fā)一個(gè)核心的SPH計(jì)算程序止剖，則需要包含三個(gè)模塊：
- 核函數(shù)模塊；
- 粒子近似模塊落君；
- 實(shí)例求解模塊穿香。
如果需要配合核心的SPH計(jì)算程序，還需要： - 建模的前處理模塊叽奥；
- 數(shù)據(jù)的可視化模塊扔水。
（挖坑一個(gè)，有空整理和修改朝氓。）
5. 思考
網(wǎng)上有挺多的SPH法的代碼魔市，還有基于分子動(dòng)力學(xué)的MPS法代碼和程序主届，如LAMMPS。要我這數(shù)學(xué)基礎(chǔ)弱的工科人來(lái)用別人的代碼待德，看懂都費(fèi)勁君丁；自己編寫(xiě)也不現(xiàn)實(shí)，需要花費(fèi)很多時(shí)間和精力将宪，而且有可能失敗绘闷。難熬~
- 現(xiàn)有的常用的建模工具都是基于有限元的，少有粒子建模较坛；
- 我估計(jì)還得學(xué)一個(gè)粒子可視化軟件印蔗。
- SPHysics是有Fortran代碼的，我可以參考丑勤，我學(xué)Fortran的华嘹；
- SPHysics依靠的建模與可視化工具是Blender軟件；
- LAMMPS偏分子學(xué)科法竞，完成度很高耙厚，基于C++，我不太想多學(xué)這么多語(yǔ)言岔霸；
- LAMMPS文檔挺多薛躬，資料也全，但是我一個(gè)學(xué)船的挺難看進(jìn)去分子層面的編程的呆细；
- 自己編程一般幾年時(shí)間內(nèi)能編寫(xiě)和維護(hù)的代碼量?jī)H為5千~5萬(wàn)的中小型程序型宝，還存在不確定性；
- 我現(xiàn)在正在請(qǐng)教師兄侦鹏，看看能給我點(diǎn)意見(jiàn)不诡曙。之后找找老師。希望能給我指條明路略水。
前路漫漫价卤，上下求索~
