No plots are attached so I can't comment on them.

If you specified TE10 then that is what you will get (I think) - the other modes are orthogonal so won't be excited, except by inaccuracies in the model, which might be why you get some at high frequency.  Lumps and bumps in the waveguide wall excite higher order modes.  At low frequencies these can't propagate (so look like small lumped inductors or capacitors), but at higher frequencies they can propagate.  They won't appear in your TE10 output (at either end) so add to the loss.

Attenuation constant should be in nepers per metre, so multiply by 8.68 and by length to get dB.

Is your mesh (is there a mesh?) fine enough to give good results at 3 THz.  The loss should be something like 32 times higher than at 3 GHz (same current, skin depth is 1/root(30) so surface resistance is 32 time higher).

I have added the files. You may look into that.

Will specifying only TE10 mode to exist tell me that the value of attenuation that I am getting is true and can be assumed safe for use?

I have used layered impedance boundary condition [design file also attached] with default meshing and other frequency sweep settings.

Results obtained and formulae used are also mentioned.

How long is your waveguide?  I would have said that above about 100 GHz the wave is not propagating at all - S11 is exactly 1, which is total reflection - perhaps your model is not correct above 100 GHz.  If the signal was propagating and then being reflected, or there was any loss or propagation, there should be some attenuation in S11.

90% of your S11 and S21 results look to me like what you would get when the software hasn't worked properly.  Check the discretisation of your model.

Thanks for your comments and suggestions @Malcolm. I too feel that there is some problem with the modelling. I am using 100mm length for the waveguide.Moreover, when I try to run simulation with reduced length than 100 mm, I can see an improved S11 and S21 (which should be there, obviously), but does that actually reflect the condition that I wish to study i.e. the loss experienced by a high frequency wave while propagating through a waveguide with dimensions corresponding to low frequency cutoff?

Also, I can see that attenuation constant plot is still unchanged.

By discretisation, do you mean to increase the resolution of simulation points? I have already attached my simulation file, if you can have a look at the parameters and suggest?

I don't have HFSS so I can't read your file.  Large waveguide should not behave in the way that your predictions are showing.  The attenuation constants look sensible, but S11 and S21 seem to be wrong.  The wavelength at 3 THz is 100 microns.  The model will need to be meshed in 10 micron steps, I would expect (1/10 wavelength at the highest frequency).  That means the model has something like 108 variables, if it is a surface mesh, or 1012 if it is a volume mesh.  Have you used too coarse a mesh (1 mm for example)?

The first three peaks of S21 and the first step are at multiples (1, 2, 4, 8) of 0.1875 THz.  This seems like it should give a clue as to what the code thinks it is modelling.

I understood your point of view, so I tried to do skin depth based mesh or length based meshing by reducing default size of 8mm to um range. But, it seems there is huge memory requirement for that. Maybe I will have to make use of cluster (HPC) service of HFSS for that.

Moreover, can you comment on the mode study? If I wish to restrict my study on the loss offered by metallic walls i.e. only attenuation constant in dB/m, then how to include effect of multiple modes that exist in such high frequency regime? In MATLAB, I can compute the individual mode's attenuation constant, then can I assume this case to be a superposition of such multiple modes? If yes then a simple addition would do or the phase in which they interact need to be considered as well? Can you comment on this method?

Waveguides are usually used in single-mode operation, because most modes travel at different speeds to each other, and have different attenuations,  Any pattern, such as a pulse, or a shape needed to couple to an input or output, that creates several modes, won't usually be recreated at the output end, because the amplitudes and phases of the different modes will no longer add up correctly.  A low order mode in a large waveguide has low attenuation.  It may be fed with a long taper from small single-mode waveguide so that other modes aren't generated.  Lumps in the wall, and curves, generate other modes which are lost power, because they never turn into something you can use again.  If you are not careful they get reflected at the output - because they don't couple to the transition - then bounce around creating a little bit of the original mode at each defect that then arrives at the output at the wrong time.

The prediction code creates one mode at a time, I think, and only propagates this.  I don't know if it absorbs other modes at the output (if it doesn't they will bounce about, and anything that created them will create the original mode from them each time they pass), but if you have a very short waveguide there shouldn't be any other modes generated in significant amounts.

@Dr. James, thanks for your suggestion. On changing the boundary condition to symmetry (Perfect E and Perfect H), I can very well simulate the huge frequency ranges with finer resolution. Now I can see that the structure is trying to be useful (by giving troughs in S21 curve) at the given "large" frequency range but still the base value of S21 in the region lies at ~ -80 dB (with reduced length, please see attached file). Does this mean that insertion loss of the waveguide is such that the wave won't be able to propagate?  

I will share more results soon.

From your latest results, and by the positive values for the S-parameters and the rugged structure, it seems that the meshing did not converge to a useful distribution of tetrahedra or that the solution frequency for the meshing was too low for the chosen sweep. Check the solution frequency of the design, as well as the adaptive number of passes and the criterium of variation of S-parameters for convergence.