open TELEMAC-MASCARET The mathematically superior suite of solvers

Dear All,

I am trying to solve advection-diffusion-reaction equation (dC/dt + UdC/dx = [Ex]d2C/dC2 -KC) in a 10.0 KM long channel with depth 10 m and width 50 m .

The channel is having a constant inflow discharge of 15 m3/s and downstream elevation is held constant at 10.0 m. At the upstream end I have a non-conservative tracer with constant conc. of 30 mg/l which is supplied for 6 hrs. I am running the simulation for 6 days, the horizontal diffusivity (COEFFICIENT FOR HORIZONTAL DIFFUSION OF TRACERS)[Ex] for tracer has been specified at 30.0 m2/s . The decay rate[K] for tracer is 1.0 per day. I understand that to code in your decay rate you have to use SOURCE-TRAC subroutine and my code is as follows

K1 = 1 ! rate in per days
K1 = K1/(86400.0)
IF(NTRAC .GE. 1) THEN
DO ITRAC=1,NTRAC
S0TA%ADR(ITRAC)%P%TYPR='Q'
S1TA%ADR(ITRAC)%P%TYPR='Q'
ENDDO
DO IPOIN2 = 1, NPOIN2
DO IPLAN = 1, NPLAN
IPOIN3 = IPOIN2 + (IPLAN-1)*NPOIN2
S0TA%ADR(1)%P%R(IPOIN3) = TA%ADR(1)%P%R(IPOIN3)*K1/2.0D0
S1TA%ADR(1)%P%R(IPOIN3) = K1/2.0D0
END DO

With this code in place I don't catch the peak of the analytical solution, could some one please take a look at my code and let me know if I am doing some thing wrong. I am attaching both numerical and analytical solution. All help is greatly appreciated.

Best regards,
Sumit

Hello,

You can look at subroutine CVDF3D, the equation is programmed as :

DC/DT + advection +diffusion + S1TA*C = S0TA

so in your case you need to have a sign - for S0TA

If you confirm I'll add explanations in subroutine source_trac.f, where no example is given for explicit source terms.

With best regards,

Jean-Michel Hervouet

Hello JMH,

Thanks for your reply, I made the change as suggested by you and the peak concentration goes down but still not same as the analytical sol. Please see the attached picture. The only difference in the code now is

S0TA%ADR(1)%P%R(IPOIN3) = -TA%ADR(1)%P%R(IPOIN3)*K1/2.0D0

I wonder do I need to make S1TA also negative.

Thanks,
Sumit

Hello,

No, normally it is now correct. I understand that you program a semi-implicit form of the decay term. You could try a fully implicit or a fully explicit form (this latter not sure to work), to see the effect of numerical approximation. Or try a smaller time step. As a matter of fact the decay alone has a solution in the form of an exponential whereas numerically speaking it is a first order in time approximation.

Regards,

JMH

Hi JMH,

Once again I tried fully implicit scheme with halving the time step and it did not improve the result. My result are almost same as attached in my previous thread SimRes2.png. For your kind information my time-step right now is 0.5 second and "scheme for advection of tracers" is set at 1, do you think using a different advection scheme might diffuse the result to catch the analytical solution. Also right now I am trying fully explicit scheme. By the way I do understand that decay alone has form of decreasing exponential nature but does that mean that it will be tough to produce the analytical solution ?

All help is greatly appreciated.

Thanks,
Sumit

Hello,

Yes, you should use an advection scheme which is mass-conservative, e.g. 3, also with the option that preserves fluxes at the entrance (TREATMENT OF FLUXES AT THE BOUNDARIES : 2;2;...).

What I do not understand is what is your analytical solution ? It seems to be 1D. What if you have differences on the vertical and then where is taken what you plot ? As there is a velocity profile on the vertical I do not really understand how the concentration could remain constant on the vertical, even if there is only decay, because the travel time is different.

With best regards,

Jean-Michel Hervouet

Hi JMH,

Yes indeed I also thought about that, if the velocity is not same every where in the vertical than neither should be conc. I was trying to reproduce the 3D model validation case provided in, "Three-dimensional numerical modeling of cohesive sediment transport and wind wave impact in a shallow oxbow lake" by Chao et al. published in Advances in Water Resources 31 (2008) 1004–1014 (model CCHE3D).

I described there validation problem in the starting of this thread it's a channel flow with constant velocity of 0.03 m/s having depth of 10m and length of 11500 m and tracer at upstream end supplied for 6 hrs with con. of 30 mg/l. I have used 10 layers to set up the problem in T3D. But now it makes me wonder if there 3D model(CCHE3D) can be used in 1D format or how else do you get completely constant velocity in a three-dimensional domain.

Also strangely, although expected I don't see any variation in the tracer conc. in the vertical while using TELEMAC-3D but the velocity is definitely not constant in vertical, please see the attached pic. The 1D plots that I have been uploading is tracer conc. at x = 2000 m from upstream and at 5 m from the bed and the surface. The width of the channel in the attached pic is 50 m.

Best regards,
Sumit

Hello,

This is indeed very strange. I would need to see your steering file. A simple test case with no flow and constant tracer for initial conditions should show only the solution of an exponential decay, independently of diffusion and hydrodynamics.

Regards,

JMH

Aah I see, so I guess the best way would be to initialize the domain with constant velocity and then deactivate advection as well as diffusion for hydrodynamics but I can still have the the advection and diffusion of tracer switched on right ?
I am also attaching my steering file this one I am simulating with smooth bottom . I am yet to see the tracer result for this case, the picture that was attached before was rough bottom.

Thanks,
Sumit

Hello Sumit,

Would you mind sending us the whole set of files, including your computation of the analytical solution (maybe you do this in Fortran ?). I would like to add this test case to the set of validation test cases.

Thank you.
Sébastien