三维海浪曲面动态模拟matlab仿真,包含仿真操作录像,代码注释

clc; clear; close all; warning off; addpath(genpath(pwd)); rng('default') % 水滴模拟的参数设置 n = 64; %网格大小 g = 9.8; % 重力常数 dt = 0.01; % 硬连线时间步长 dx = 1.0; dy = 1.0; nplotstep = 8; %x和y方向上的网格大小 ndrops = 5; %最大滴数 dropstep = 500; % 最大滴数 D = droplet(1.5,21); %下降间隔 %创建水滴初始条件 %初始化模拟的图形 [surfplot,top,start,stop] = initgraphics(n); %模拟将继续进行，直到满足特定的停止标准。 while get(stop,'value') == 0 set(start,'value',0) %初始化高度字段和流量字段的数组。 H = ones(n+2,n+2); U = zeros(n+2,n+2); V = zeros(n+2,n+2); Hx = zeros(n+1,n+1); Ux = zeros(n+1,n+1); Vx = zeros(n+1,n+1); Hy = zeros(n+1,n+1); Uy = zeros(n+1,n+1); Vy = zeros(n+1,n+1); ndrop = ceil(rand*ndrops); nstep = 0; %计算将添加到系统中的液滴数。 %内环，用于时间步进。模拟将继续进行，直到满足特定的停止标准。 while get(start,'value')==0 && get(stop,'value')==0 nstep = nstep + 1; %向系统中随机添加水滴。这是在每一个“dropstep”时间步长完成的，并且受到“ndrops”的限制。 if mod(nstep,dropstep) == 0 && nstep <= ndrop*dropstep w = size(D,1);%获取初始条件数组的大小 i = ceil(rand*(n-w))+(1:w);%随机选择放置滴剂的位置 j = ceil(rand*(n-w))+(1:w);%随机选择放置滴剂的位置 H(i,j) = H(i,j) + rand*D;%在IC确定的某个初始高度处，将落差添加到该位置的高度字段 end %反射边界条件在此适用于所有字段。这意味着边界处的场在乘以-1之后被复制到相反的边界。 H(:,1) = H(:,2); U(:,1) = U(:,2); V(:,1) = -V(:,2); H(:,n+2) = H(:,n+1); U(:,n+2) = U(:,n+1); V(:,n+2) = -V(:,n+1); H(1,:) = H(2,:); U(1,:) = -U(2,:); V(1,:) = V(2,:); H(n+2,:) = H(n+1,:); U(n+2,:) = -U(n+1,:); V(n+2,:) = V(n+1,:); %前半步 %x方向 i = 1:n+1; j = 1:n; %高度 Hx(i,j) = (H(i+1,j+1)+H(i,j+1))/2 - dt/(2*dx)*(U(i+1,j+1)-U(i,j+1)); %x动量 Ux(i,j) = (U(i+1,j+1)+U(i,j+1))/2 - ... dt/(2*dx)*((U(i+1,j+1).^2./H(i+1,j+1) + g/2*H(i+1,j+1).^2) - ... (U(i,j+1).^2./H(i,j+1) + g/2*H(i,j+1).^2)); %y动量 Vx(i,j) = (V(i+1,j+1)+V(i,j+1))/2 - ... dt/(2*dx)*((U(i+1,j+1).*V(i+1,j+1)./H(i+1,j+1)) - ... (U(i,j+1).*V(i,j+1)./H(i,j+1))); %y方向 i = 1:n; j = 1:n+1; %高度 Hy(i,j) = (H(i+1,j+1)+H(i+1,j))/2 - dt/(2*dy)*(V(i+1,j+1)-V(i+1,j)); %x动量 Uy(i,j) = (U(i+1,j+1)+U(i+1,j))/2 - ... dt/(2*dy)*((V(i+1,j+1).*U(i+1,j+1)./H(i+1,j+1)) - ... (V(i+1,j).*U(i+1,j)./H(i+1,j))); %y动量 Vy(i,j) = (V(i+1,j+1)+V(i+1,j))/2 - ... dt/(2*dy)*((V(i+1,j+1).^2./H(i+1,j+1) + g/2*H(i+1,j+1).^2) - ... (V(i+1,j).^2./H(i+1,j) + g/2*H(i+1,j).^2)); %后半步 i = 2:n+1; j = 2:n+1; %高度 H(i,j) = H(i,j) - (dt/dx)*(Ux(i,j-1)-Ux(i-1,j-1)) - ... (dt/dy)*(Vy(i-1,j)-Vy(i-1,j-1)); %x动量 U(i,j) = U(i,j) - (dt/dx)*((Ux(i,j-1).^2./Hx(i,j-1) + g/2*Hx(i,j-1).^2) - ... (Ux(i-1,j-1).^2./Hx(i-1,j-1) + g/2*Hx(i-1,j-1).^2)) ... - (dt/dy)*((Vy(i-1,j).*Uy(i-1,j)./Hy(i-1,j)) - ... (Vy(i-1,j-1).*Uy(i-1,j-1)./Hy(i-1,j-1))); %y动量 V(i,j) = V(i,j) - (dt/dx)*((Ux(i,j-1).*Vx(i,j-1)./Hx(i,j-1)) - ... (Ux(i-1,j-1).*Vx(i-1,j-1)./Hx(i-1,j-1))) ... - (dt/dy)*((Vy(i-1,j).^2./Hy(i-1,j) + g/2*Hy(i-1,j).^2) - ... (Vy(i-1,j-1).^2./Hy(i-1,j-1) + g/2*Hy(i-1,j-1).^2)); %更新绘图 if mod(nstep,nplotstep) == 0 C = abs(U(i,j)) + abs(V(i,j)); % Color shows momemtum t = nstep*dt; tv = norm(C,'fro'); set(surfplot,'zdata',H(i,j),'cdata',C); set(top,'string',sprintf('t = %6.2f, tv = %6.2f',t,tv)) drawnow end if all(all(isnan(H))), break, end % Unstable, restart end end close(gcf)