Have a look at https://au.mathworks.com/matlabcentral/fileexchange/73050-lamb-wave-dispersion-curve

https://au.mathworks.com/matlabcentral/fileexchange/65138-masw-dispersion-curve

Muhammad Ali

Thank you for your reply. I have read that program; however, it was developed for real wavenumber k, and I think it can not be used for complex wavenumber k.

Dear Hao,

There are many different approaches to solve this problem. I assume you are trying to solve the analytical dispersion relation by numerical root-finding algorithms? This is indeed quite challenging when you need to find all solutions, particularly when the solutions can be complex-valued.

It is much easier to employ one of the formulations that lead to an eigenvalue problem. This way, all solutions can be found regardless whether they are real or complex. Again, there are many similar approaches; the one I have implemented can be found here

Article A numerical approach for the computation of dispersion relat...

If the complex wavenumber is related to the viscoelastic material, please, have a look at: Article Investigation of guided waves propagation in orthotropic vis...

Here, Legendre polynomial and Newton Raphson methods were selected.

Hauke Gravenkamp Cherif Othmani

Thank you so much for your advice and attached paper link！

In the past few days, I have learned a new method named Spectral Collocation Method (SCM). Using the Chebyshev differential matrix and companion displacement vector field, I have successfully plotted the dispersion curves by solving the eigenvalue problem about angular frequency omega or wavenumber k.

However, I have met new problems:

1. Unlike the root-searching algorithm, I found that it is difficult to separate the cross modes using SCM, as shown in the figure .

2. In order to check my codes, I want to plot the real-part results. Since the obtained k are complex numbers, I use k(real(k)/imag(k)>100) to plot the propagating modes whose imaginary part of k should be zero (or at least much less than the real part during the calculation). But this can lead to additional points or missing points.

What can I do to plot exactly the real part solution?

If there is a better way to verify my code?

3. When using a non-uniform Chebyshev grid of N points, the omega based SCM can give a better result as N increases. However, the k based SCM may give a worse result as N increases which should not happen. It can be seen in attached figures and .

What may have caused this problem?

More details about the k based SCM can be found in the attached PDF.

Thank you for spending your time reading my questions.

I would really appreciate it if you can help me with these problems.

Dear Hao,

1.

The issue of tracing modes across these crossing points is actually a lot more complicated than you would expect -- at least in the general case. See:

Article Analyzing modal behavior of guided waves using high order ei...

Article Osculations of spectral lines in a layered medium

Most poeple use either the so-called modal assurance criterion (MAC), which is based on the similarity of eigenvectors at different frequencies, or eigenvalue derivatives, i.e. the slope of the dipersion curves. Both approaches can fail under certain circumstances.

Article The computation of dispersion relations for three-dimensiona...

If you are dealing only with homogeneous plates, the easiest way is to model only half of the plate and successively apply symmetric and antisymmetric boundary conditions. You will obtain two sets of modes which do not cross and, hence, are trivial to trace. This only works if the structure is symmetric.

2.

I usually simply specify a threshold for the real part of the eigenvalue rather than the ratio. From a physical point of view, note that the imaginary part represents the attenuation. Hence, large imaginary parts will cause very small displacements at a given distance from the source. This consideration can be used to choose the threshold based on the physical setup of the measurement.

Regards,

Hauke

Hauke Gravenkamp

Dear Hauke,

Thank you so much for your detailed reply and suggestions and they are really helpful for me. I would try the half-model method using symmetric and antisymmetric boundary conditions.

Besides, I thought I have found my mistake last night. The original question has been translated to a k based eigenvalue question in the form:

M1*U' = k*M2*U' and U' = [Uhat U]T

where T for transpose and Uhat = kU is the companion displacement vector field.

M2 = [Q2 Z; Z I] and M1 = [-Q1 -Q0 ; I Z] where I is the identity matrix and Z is a matrix of zeros.

Since the results of LAMB modes are good with omega based calculation but became bad with k based calculation. And the curves got worse as N increase. I thought the problem must have been caused by the equation transformation.

Then I found this:

M2 is singular. In matrix M1, Q1 and Q0 are in 1e10 order with lame constants in parameters while I and Z are in 1 order. That would lead to an ill-conditioned matrix.

When I print the condition numbers of M1, they reached as high as about 1e16. As N increase, cond(M1) may get even bigger. So that may be the answer to my question.

After I divided Q2/Q1/Q0 with lame constants, cond(M1) decreased to be in 1e5 order and the results of both LAMB modes and SH modes turned right.

Best Regards.

Hao

https://www.ndt.net/article/IranNDT2017/papers/504.pdf

You can find how to calculate complex wavenumber dispersion curves in my work

DOI: 10.1007/s00707-017-2025-y

For generally anisotropic viscoelastic/elastic plates which are more complicated cases compared to Lamb waves, they can be found in

DOI: 10.1016/j.ijsolstr.2020.05.031

DOI: 10.1016/j.ijsolstr.2018.06.001

Feng Zhu Thank you so much for sharing. I have already fixed this problem; however, I am stuck with another problem and the link was attached below:

https://www.researchgate.net/post/How_to_calculate_the_energy_reflection_coefficients_in_LAMB_wave_edge_reflection_problem

I would appreciate it if you could help me with the edge reflection problem.

These three details seem right to me. I suggest you check the calculation of Power flux. A validation for Power flux is that the energy velocity calculated for Power flux equals to group velocity (dw/dk) in small imaginary or purely real wavenumber. You can find some brief information abouth Power flux in DOI: 10.1016/j.ijsolstr.2018.06.001 where the energy velocity was calculated. More detailed infomration can be obtained form these related references there. Hope it is useful to you.