Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

MATLAB Shallow water dam break problem for a 1D matlab model

  1. Feb 20, 2012 #1
    i am trying to spproximate a PDE in the form below using the lax wendroff 2 step method in matlab coding:

    [h ; hu ] = [ hu ; hu^2 + 1/2gh^2] = [0; -ghbz] (where bz will equal zero)

    i believe this is then the case

    d(h)/dt + d(hu)/dx = 0


    d(hu)/dt + d(hu^2 + 1/2gh^2)/dx = 0

    as the st venant equation was formed by depth deriving the 1D continuity equation

    however i have attempted to program this and i keep produced a column matrix full of "NaN" could some one try to point out where i have been going wrong? thanks kyle

    code for matlab:

    %%%%%%%%%%%%%%%% function compare close all clc clear all

    %intial values ntime = 1000000; dt=0.00050; nx=100; time = 0; a=2;output=0.24; %step size calculation dx= (1/nx); %create size of u_int vector u_int = zeros(nx,2); %create little u_int vector for the initial values of height begining and %end depending which value of A is used u_int(nx,2) =1; u_int(1,1) = 1;

    %loop for the two directions needed for vec = 1:2 % %lax_wendroff [H] = lax_wendroff_2(nx,a,output, dt, ntime, time,u_int, vec,dx); figure(vec) H plot(H(:,vec),'r*');grid on;legend('centered'); hold on end end

    function [H] = lax_wendroff_2(nx,a,output, dt, ntime, time,u_int, vec,dx); %compute original size matracies g = 9.81; uh=zeros(nx,2); update = uh ; H = uh; Hy = uh; UH = uh; if vec == 1; %for i = 1:nx u_int(nx-2:nx,2) = 1; u_int(1:3,1) = 1; %end for i = 1:nx UH(i,vec)=u_int(i,vec); Hy(i,vec)=u_int(i,vec); uh(i,vec)=u_int(i,vec); H(i,vec)=u_int(i,vec); A=a; end else uh(nx,vec)=u_int(nx,vec); h(nx,vec)=u_int(nx,vec); A=-a; end for p = 1:ntime; if(time + dt>output); dt=output-time; end % calculate coefficent coeff=dt/dx; %calculate interim steps for i = 1:(nx-1) % first derivative Hy(i+1,vec)= ((H(i+1,vec)+H(i,vec))/2) - ((coeff*A/2)*((uh(i+1,vec) - uh(i,vec))));

    UH(i+1,vec)=((uh(i,vec)+uh(i+1,vec))/2) ;
    -coeff*A*(((uh(i+1,vec)^2)/(H(i+1,vec)+g/2*H(i+1,vec)^2)) - ((uh(i,vec)^2)/(H(i,vec)+g/2*H(i,vec)^2)));
    %calculate full time steps
    for i = 2:(nx-1)
    update_h(i+1,vec) = H(i,vec) - coeff*A*(UH(i+1,vec) - UH(i,vec));%%%% check i notation and convention may be wrong way
    update_uh(i+1,vec) = (uh(i,vec)) - (coeff*A*(((UH(i+1,vec)^2)/(Hy(i+1,vec)+g/2*Hy(i+1,vec)^2))- ((UH(i,vec)^2)/(Hy(i,vec)+g/2*Hy(i,vec)^2))));
    for i = 3:(nx-1)
    H(i,vec) = update_h(i,vec);
    uh(i,vec) = update_uh(i,vec);
    time = time + dt;
    if time==output
    end end %%%%%%%%%%%%%%%%%
  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted