闪电动画模拟（Dielectric Breakdown Model）附源码

　　当两个物体之间存在较大的电势差时会出现放电现象，比如生活中常见的闪电现象，闪电形成的条件就是云层积累了大量负电荷之后与地面之间形成了强大的电势差。目前关于闪电建模的方法比较少，下面介绍一种利用电介击穿模型来模拟闪电的方法，电介击穿模型可以模拟自然界许多现象，该方法通过迭代求解Laplace方程得到放电过程的中间状态。

　　初始电位结构如下图所示，首先在2维栅格正中心的单元放置一个负电荷Ф = 0（灰色），然后在其周围放置一圈正电荷Ф = 1（黑色），而其他栅格单元可以通过求解Laplace方程得到：

　　Laplace方程求解完之后，我们构建一个列表记录负电荷（Ф = 0）周围的栅格单元，并随机选取其中一个单元作为下一个被击穿的单元，被选中的栅格单元设置Ф = 0，同时把它作为下一步迭代计算时的边界条件。

　　相邻栅格单元被击穿的概率与其计算得到的电势相关，其概率如下：

式中i代表第i个与负电荷相邻的栅格单元，n代表与负电荷相邻的栅格单元个数，而η是一个用户设定的参数，通过实验表明η值可以控制电弧的分叉密度，当η = 1时，分叉密集，当η逐渐增大时，分叉密度慢慢减小。

图：不同结构的初始电位，其中灰色点代表负电荷Ф = 0，黑色点代表正电荷Ф = 1

图：初始电位为（b）结构时的模拟结果

左：η = 1，中：η = 2，右：η = 3

% Laplacian Growth Model
clc
clear all
close all

w = 256; % width
h = 256; % height
eta = 3;

% 指定正负电荷
ni = floor(w/2)*h + (1:2)'; % negative charge
pi = ((1:w)' - 1)*h + h; % positive charge

% 构建laplace矩阵
S1 = bsxfun(@plus, ((1:(w-1)) - 1)*h, (1:h)');
D1 = bsxfun(@plus, ((2:w) - 1)*h, (1:h)');
S2 = bsxfun(@plus, ((1:w) - 1)*h, (1:(h-1))');
D2 = bsxfun(@plus, ((1:w) - 1)*h, (2:h)');

S = [S1(:); S2(:)];
D = [D1(:); D2(:)];
A = sparse([S D], [D S], 1);
n = size(A,1);
% Lp = speye(h*w) - bsxfun(@rdivide, A, sum(A,2));
Lp = spdiags(sum(A,2), 0, n, n) - A;
rhs = zeros(n, 1);

% 循环计算下一步状态
[X0, Y0] = ind2sub([h,w], (1:h*w)');
im = zeros(h,w);
iter = 1;
i = 1;
while ~any(ismember(ni, pi))
    phi = solve_equation(Lp, rhs, [ni;pi], [zeros(size(ni));ones(size(pi))]);
    
    [~, adj] = find(A(ni,:));
    adj = unique(adj);
    adj = adj(~ismember(adj, ni));

    k = randsample(adj, 1, true, phi(adj).^eta/sum(phi(adj).^eta));
    ni = [ni; k];

    if mod(iter,10) == 1
        [X, Y] = ind2sub([h,w], ni);
        [~, dis] = knnsearch([X,Y], [X0,Y0]);
        im = reshape(1-dis/3, [h,w]);

        imshow(im, 'Border','tight');
        drawnow;

        Frame(i) = getframe(gcf);
        i = i + 1;
    end
    fprintf('iteration: %d ...
', iter);
    iter = iter + 1;
end

本文为原创，转载请注明出处：http://www.cnblogs.com/shushen。

 

 

参考文献：

[1] Theodore Kim and Ming C. Lin. 2007. Fast Animation of Lightning Using an Adaptive Mesh. IEEE Transactions on Visualization and Computer Graphics 13, 2 (March 2007), 390-402.