1
$\begingroup$

I am going to solve this advection equation regarding the flow simulation of an energy tower \begin{equation} Y_t+ v(t) Y_x =0 \end{equation} with the following boundary conditions which depend on the direction of the flow \begin{align} Y(x,0) = Y_0 \end{align} And \begin{align} Y(0,t) = Y_{bot} \text{ if }v(t)>0 \\ Y(1,t) = Y_{top} \text{ if }v(t)<0 \end{align} And for the simulations, I have written a Matlab code in the upwind method with different inputs to $Y_0, Y_{bot}$ and $Y_{top}$. But almost all the solution seems to blow up in time. Can someone tell me what has gone wrong here?

t_end = 10;
n_x = 10;

y0 = 1;  % Inputs for the initial condition Y(x,0)
y_bot = 0; %  Inputs for the boundary condition Y(0,t)
y_top = 1; %  Inputs for the boundary condition Y(1,t)

dx = 1/n_x;

v0 =1; % initial velocity
CFL = 1;
dt = CFL*(dx)/((abs(v0)));

n_t = round(t_end/dt)
x = linspace(0,1,n_x);  % space grid
t = linspace(0,t_end,n_t); % time grid

Y = zeros(n_x, n_t);

Y(:,1)= y0;  % initial condition Y_0
Y(1,:)=y_bot; % boundary condition at the bottom of the tower
% Y(n_x,:)= y_top; % boundary condition at the top of the tower

v = sort(randi([-20 50],1,n_t),'descend'); % generating a random array of velocities in the decending order 
v(1) = v0; % make the initial velocity as v0

for n=2:n_x
    for j=2:n_x-1
        
        if v(n) >= 0
            Y(1,n) = y_bot;
        else
            if v(n_x) < 0
                Y(n_x,n) = y_top;
            end
        end
 

        if v(n) >= 0
            Y(j,n) = Y(j,n-1)-v(n)*dt/dx*(Y(j,n-1)-Y(j-1,n-1));
        else
            if v(n)< 0
                Y(j,n) = Y(j,n-1)-v(n)*dt/dx*(Y(j+1,n-1)-Y(j,n-1));
            end
        end
    end
end

solution for Y(x,t)

$\endgroup$
1
  • $\begingroup$ The line if v(n_x) < 0 is meaningless, always true. Even corrected to if v(n) < 0 it is always true if it is reached. $\endgroup$ Commented Nov 21, 2022 at 7:03

1 Answer 1

2
$\begingroup$

You are using n_x at at least two places as boundary for the time index.

Loops are slow, so avoid them where not necessary. The time loop can be condensed to

for n=2:n_t
    if v(n) >= 0
        Y(1,n) = y_bot;
        Y(2:n_x,n) = Y(2:n_x,n-1)-v(n)*dt/dx*(Y(2:n_x,n-1)-Y(1:n_x-1,n-1));
    else
        Y(n_x,n) = y_top;
        Y(1:n_x-1,n) = Y(1:n_x-1,n-1)-v(n)*dt/dx*(Y(2:n_x,n-1)-Y(1:n_x-1,n-1));
    end
end

Your CFL-dt calculation needs to take the highest velocity into account. You compute it for v0=1, but the actual highest values in magnitude are 50 and -20. Being sneaky and setting v(1)=v0 does not change the general problem, the second largest velocity is still too large.

$\endgroup$
3
  • $\begingroup$ Thank you for the feedback. But to take the highest velocity by taking max(v) to calculate time step $dt$, according to the CFL. how can I define $v = sort(randi([-20 50],1,n_t)$, since $n_t$ depends on $dt$? $\endgroup$
    – TMW
    Commented Nov 21, 2022 at 9:14
  • $\begingroup$ You take the theoretical maximum 50. If you want truth in advertizing, linearly stretch the random array after construction to exactly fit the boundaries. $\endgroup$ Commented Nov 21, 2022 at 9:25
  • $\begingroup$ Thank you. This is a question on a bit improved step, I want to know if my $v(t)$ is also developing inside the loop and how I can get the max velocity.? This happens when $v(t)$ also comes out as a solution of an ODE w.r.t. time with some initial velocity $v_0$ $\endgroup$
    – TMW
    Commented Nov 21, 2022 at 10:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.