[CIG-SHORT] Kelvin-Voigt viscoelasticity

[okuyama.satoshi]

Dear Developers,

TL;DR: How can I model kelvin-voigt material when GenMaxwellPlaneStrain/
       GenMaxwellIsotropic3D does not allow shear-ratio-1 being huge value?

I recently started analyzing kelvin-voigt viscoelastic deformation using 
PyLith. Since PyLith does not support kelvin-voigt model directly, I tried 
to do this using Generalized Maxwell material.

Since kelvin-voigt model consists of single spring and single dashpot in 
parallel, I have to allow only spring0 and dashpot1 to be effective.


  +--------- spring0 (mu0) ------------+   strain0, stress0
  |                                    |
  +- spring 1(mu1) -- dashpot1 (eta1) -+   strain1, stress1
  |                                    |
--+                                    +-- total strain, total stress
  |                                    |
  +- spring2 (mu2) -- dashpot2 (eta2) -+   strain2, stress2
  |                                    |
  +- spring3 (mu3) -- dashpot3 (eta3) -+   strain3, stress3


Each maxwell model (and a spring) is connected in parallel, so

(total strain) = (strain0) = (strain1) = (strain2) = (strain3)
(total stress) = (stress0) + (stress1) + (stress2) + (stress3)

To disable spring2, 3 and dashpot2, 3, I set mu and eta for corresponding 
spring/dashpot to zero. This should make strain2 and strain3 always zero as 
required.

Next, spring1 and dashpot1 are connected in series, so

(strain1) = (spring1 strain) + (dashpot1 strain)
(stress1) = (spring1 stress) = (dashpot1 stress)

To disable spring1, I want to set mu1 to infinity so that (spring1 strain) 
becomes zero.

(mu0, mu1, mu2, mu3) = (mu, inf, 0, 0)
(eta1, eta2, eta3) = (eta, 0, 0)

Since PyLith defines fraction of each shear modulus rather than shear modulus 
itself, I set shear-ratio-1 to 1e+10.

---------------------------- part of spatialdb ----------------------------
  value-names =  density vs vp shear-ratio-1 shear-ratio-2 shear-ratio-3 viscosity-1 viscosity-2 viscosity-3
  value-units =  kg/m**3  m/s  m/s None None None Pa*s Pa*s Pa*s
                   .
                   .
                   .

0 0 2700 3333 5800 1e+10 0 0 1e+18 0 0
---------------------------------------------------------------------------

However, PyLith does not allow sum of fractions greater than 1!

----------------------- part of the error message -------------------------
RuntimeError: Shear modulus ratios sum to a value greater than 1.0 for
Generalized Maxwell model.
Ratio 1: 0
Ratio 2: 0
Ratio 3: 0
Total:   1e+10
---------------------------------------------------------------------------

I looked into the source and confirmed this. It just adds up the shear ratio
and throws runtime_error when the sum is greater than 1.

--------------------- part of GenMaxwellPlaneStrain.cc --------------------
272  PylithScalar visFrac = 0.0;
273  for (int imodel = 0; imodel < numMaxwellModels; ++imodel) 
274    visFrac += dbValues[db_shearRatio + imodel];
275  if (visFrac > 1.0) {
276    std::ostringstream msg;
277    msg << "Shear modulus ratios sum to a value greater than 1.0 for\n"
278	<< "Generalized Maxwell model.\n"
279	<< "Ratio 1: " << propValues[db_shearRatio  ] << "\n"
280	<< "Ratio 2: " << propValues[db_shearRatio+1] << "\n"
281	<< "Ratio 3: " << propValues[db_shearRatio+2] << "\n"
282	<< "Total:   " << visFrac << "\n";
283    throw std::runtime_error(msg.str());
284  } // if
---------------------------------------------------------------------------

Now, how can I configure kelvin-voigt visco-elasticity using GenMaxwellPlane
Strain or GenMaxwellIsotropic3D?

I would appreciate any advice. Thank you.



# I think "propValues" in line 279-281 should be "dbValues".

---
Satoshi Okuyama