Hi all,

i have been recently introduced to CFD. I have written a explicit finite difference code in C for numerically solving the wave equation but the results are not very close to the analytical solution. I want to discuss and solve this problem. Please let me know if anyone has an idea on this.

Kunal

In CFD when people talk of wave equation there can be the first order wave equation or the classical one.Also when you talk of some numerical troubles,you need to be specific with what finite difference/volume/element scheme you employed,what boundary conditions were applied. Also every scheme has dissipation and dispersion.Did you try grid convergence studies ? Just be more specific and people would be happy to help you.

U(i)(n+1) = 2*U(i)(n) - U(i)(n-1)+((c*Î"t/Î"x)^2) * (U(i+1)(n) - 2*U(i)(n) + U(i-1)(n))

i = ith step for spatial increment i=1,2,…..I.

n = nth step for time increment n=1,2……N.

I am using the above finite difference (explicit method) scheme for finding the value of the dependant variable in the 1-dimensional ideal wave equation.

How do I take care of the term U(i)(n-1) at n=2 and U(i+1)(n) at i = I i.e. the last grid point in the spatial direction for every n.

the boundary and initial conditions are

U(x,t=0)= e^t U(x=0,t)= e^(-x/c) c=1

the domain of computation is 0<=x<=1 0<=t<=1.

I think there is a problem with your initial condition and boundary condition.In a wave propagation problem you have some sort of wave propagation introduced into the system.The initial condition says that the value of u=1 for all t and the boundary condition says u=1 for all t.I would expect someting like

u(x,t=0)=e^(-x/c) and u(x=0,t)=e^(-t) as conditions.Check to make sure tat you are using the proper conditions.Also what is the value of the courant number you use.Since the scheme is explicit,you might have stability restrictions.

NO the boundary and initial conditions are correct.

the value of u(0,t) will change with every t at every time step. same way the value of u(x,0) will also change with every spatial grid point as x is different for every grid point.the value of in my case c=1, delta_x = 0.01 and delta_t = 0.005. so the courant number becomes C= c*delta_t/delta_x = 0.5.

the stability requirement for wave equation is C<=1.

Is there any mistake with my finite difference formulation.

i can post the code if u want to see it.

Are you sure the boundary condition at x=0 is e^t and not e^-t.Because e^t grows uncontrolled with time and after sometime it becomes almost infinite.

yes i am 100% sure that the boundary condition is e^t. besides i am not worried about the growing because my grid is from 0 to 1 for both time and space.

My problem is that for the numerical solution the at the last point in the grid i am getting the answer for the dependant variable around 2.7 but in the analytical solution it is 1.0

Thus the numerical solution disagrees with the analytical solution by a large amount.

What is the time interval for which you run the code ?. Do you use characteristic boundary condition at domain exit ?.

the time interval is from 0 to 1 and the step size is 0.005.

Did you plot the field for u after t=1.Do you at least have a qualitative match.Did you try to use a more refined grid and a smaller time stepping.Did the results converge for that.

i did not plot anything. but the next thing i am gonna do is that.

i refined the grid but it does not help.

U(i)(n+1) = 2*U(i)(n) -U(i)(n-1)+((c*Î"t/Î"x)^2) * (U(i+1)(n)- 2*U(i)(n) + U(i-1)(n)).

when i reach the last grid point in the spatial direction for a given time step then how do i take care of the value

U(i+1)(n). I don't have any values over there.

Kunal,

Judging from your difference equation, it seems you are solving the 2nd order wasve equation u,xx = c^2 * u,tt . This calls for TWO ICs (initial conditions), e.g., u(x,t=0)AND u,t(x,t=0) and TWO BCs (boundary conditions), e.g., u(x=0,t) AND u(x=1,t). According to you, there is only ONE IC (U(x,t=0)= e^t) and ONE BC ( U(x=0,t)= e^(-x/c)). Therefore, the problem is ill-posed. Add the missing IC and BC, use CFL number <1 and you'll most likely succeed in getting a good approximation to the expected analytic solution.

I hope this helps, Rami