[GeoMechanicsApplication] Add a linear elastic constitutive law for line interfaces #12654

avdg81 · 2024-08-30T07:16:48Z

📝 Description
Added a linear elastic constitutive law that uses an incremental formulation for line interface elements. This new constitutive law will eventually replace the existing classes LinearElastic2DInterfaceLaw and LinearElastic3DInterfaceLaw. So far it aims at 2D models only. In addition to the README file, the unit tests of class GeoIncrementalLinearElasticInterfaceLaw document its intended behavior.

🆕 Changelog
Other changes in this PR include adding two new variables to the GeoMechanicsApplication: INTERFACE_NORMAL_STIFFNESS and INTERFACE_SHEAR_STIFFNESS.

@indigocoral: Could you please read the README file carefully and see whether it makes sense to you? Any suggestions from your side are greatly appreciated. Thanks!

…ement and stiffness values

…lative displacement and previous traction from the constructor to `InitializeMaterial`

…ept by the law itself by using `FinalizeMaterialResponseCauchy`

- When checking the constitutive law, also call the `Check` member of the base class. - Replaced a few intermediate variables by using temporaries (return by other functions).

They now have the same order as in class `ConstitutiveLaw`.

- Renamed a test case. - Set a Z coordinate to zero (since we assume a 2D model). - Introduced a constant for the adopted relative tolerance.

mnabideltares

@avdg81 @WPK4FEM
Thank you for adding constutative law for the interface.
I have few minor comments.

mnabideltares · 2024-08-30T07:31:05Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+## Incremental linear elastic interface law
+The constitutive law for an interface element relates tractions $\tau$ to relative displacement $\Delta u$.
+Relative displacement for interface element is the differential motion between the two sides of the interface. As a
+consequence the relative displacement has unit of length [L] and the stiffness has unit of force over cubic length [F/L^3].


Preferly to be "m" and "N" replacing of L and F.
I know that you don't want to relate it to any specific unit system like SI, and wanted to write it as symbolic. But in the standard units symbols, F is a dependent unit and it hasn't specific symbol. The force then becomes [M L T^-2].

Moreover: The unit is preferly to be written as "$\mathrm{[UNIT]}$" , to account for powers, For example, "$\mathrm{[F/L^3]}$" or "$\mathrm{[N/m^3]}$"

We will make use of the formula editor for correct formatting. However, in engineering it is very common to use the derived unit for Force i.s.o. mass * length * time^-2. To be independent of the system used ( S.I or imperical etc. ) we will express in L, F, t, T (i.e. quantities) i.s.o. mm, kN, days and Fahrenheit (i.e. some units).

mnabideltares · 2024-08-30T07:39:13Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+consequence the relative displacement has unit of length [L] and the stiffness has unit of force over cubic length [F/L^3].
+
+### Relative displacement and traction
+In 2D plane strain y is the opening/closing direction of the interface, while differential motion in the tangential direction


"y" is a variable, so it has to be in italic (equation form). So, i think better to use "$y$" instead.

mnabideltares · 2024-08-30T07:46:00Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+
+$$ \Delta u = \begin{bmatrix} \Delta u_y \\ \Delta u_x \end{bmatrix} $$
+
+$$ \tau = \begin{bmatrix} \tau_{yy} \\ \tau_{xy} \end{bmatrix} $$


I am not sure whether I am correct in the case of interface, but in the normal stress tensors, the normal components are defined by $\sigma_{xx}$. $\sigma_{yy}$ etc. and shears by $\tau_{xy}$. Is it also here applicable?

The name for stress in an interface is traction, and it is denoted by $\tau$. The symbol $\sigma$ is not used in this context at all.

mnabideltares · 2024-08-30T07:50:15Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+$$ C = \begin{bmatrix} C_{yy} & 0     \\
+                       0     & C_{xy} \end{bmatrix}$$


What are the $C_{xx}$ and $C_{xy}$? I know they are stiffness parameters which are mentioned in the keywords, but preferly to be explicitely mentioned what they are, and what are their units.

Added the dimension of the stiffness values to be more explicit about what is expected. The normal and shear stiffness terms are clear from the input parameter names.

mnabideltares · 2024-08-30T09:56:18Z

...ons/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.cpp

+    const auto result = BaseType::Check(rMaterialProperties, rElementGeometry, rCurrentProcessInfo);
+
+    KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(INTERFACE_NORMAL_STIFFNESS))
+        << "No interface normal stiffness defined" << std::endl;


"No interface normal stiffness is defined"

mnabideltares · 2024-08-30T10:49:47Z

...ons/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.cpp

+        << rMaterialProperties[INTERFACE_NORMAL_STIFFNESS] << std::endl;
+
+    KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(INTERFACE_SHEAR_STIFFNESS))
+        << "No interface shear stiffness defined" << std::endl;


rfaasse

Thank you for implementing this constitutive law as a well tested + documented package. I have a few questions and suggestions, but nothing major and mostly just to gain a bit of understanding myself

...cation/tests/cpp_tests/custom_constitutive/test_incremental_linear_elastic_interface_law.cpp

applications/GeoMechanicsApplication/geo_mechanics_application_variables.cpp

rfaasse · 2024-09-02T14:09:19Z

...ons/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.cpp

+
+ConstitutiveLaw::SizeType GeoIncrementalLinearElasticInterfaceLaw::WorkingSpaceDimension()
+{
+    return 2;


Could this ever become 3d (e.g. when using it for a interface plane, where you have multiple 'tangential' directions and a normal)?

I can imagine that at some point we will extend this class towards 3D models. In that case, I would suggest to supply the desired dimension on construction of an instance of this class, and store it with the object. Then this query should return whatever value the corresponding data member contains. Do you think this might work, or do you foresee any issues there?

I don't see any issues in this particular function, I think the challenge is more in the formulation, since then there is more than one 'SHEAR' direction. I don't know enough to see if this way of formulating it, makes it harder to go to 3D at some point.

...ons/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.cpp

rfaasse · 2024-09-02T14:23:02Z

...ons/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.cpp

+Matrix GeoIncrementalLinearElasticInterfaceLaw::MakeConstitutiveMatrix(double NormalStiffness,
+                                                                       double ShearStiffness) const
+{
+    auto result  = Matrix{ZeroMatrix{GetStrainSize(), GetStrainSize()}};
+    result(0, 0) = NormalStiffness;
+    result(1, 1) = ShearStiffness;
+    return result;
+}


This is more of a personal preference, but I'd put this function close to where it's called (i.e. directly below CalculateMaterialResponseCauchy)

The definitions of the member functions in the implementation file have the same ordering as the declarations in the header file. So first the public member functions, then the private member functions. The public members are all overrides of the interface defined by class ConstitutiveLaw, and they are in the same order as in the header file of ConstitutiveLaw. From this point of view, I'm a bit reluctant to rearrange their order. Is that acceptable to you, or would you still insist on having the definition of MakeConstitutiveMatrix close to CalculateMaterialResponseCauchy?

I get that reasoning as well. Since I think it improves readability, I still prefer moving it, but this is your PR and I'm also completely fine with choosing for consistency with the base class and header file.

…incremental nature of the law)

rfaasse

I have one comment left about a render issue in the README.md, other than that it's good to go.

I also added one question, but that's out of curiosity and doesn't at all block the PR

rfaasse · 2024-09-03T09:51:40Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+## Incremental linear elastic interface law
+The constitutive law for an interface element relates tractions $\tau$ to relative displacement $\Delta u$.
+Relative displacement for interface element is the differential motion between the two sides of the interface. As a
+consequence the relative displacement has unit of length [$\mathrm{L}$] and the stiffness has unit of force over cubic length [$\mathrm{F/L^3}$].


These mathrms don't render correctly, I think reversing dollar and [ should work $[\mathrm{F/L^3}]$ (would look like $[\mathrm{F/L^3}]$), adding a space between [ and $ also works

Good catch. Done.

rfaasse · 2024-09-03T09:55:28Z

...ons/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.cpp

+    KRATOS_ERROR_IF_NOT(rMaterialProperties[INTERFACE_NORMAL_STIFFNESS] > 0.0)
+        << "Interface normal stiffness must be positive, but got "
+        << rMaterialProperties[INTERFACE_NORMAL_STIFFNESS] << std::endl;
+
+    KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(INTERFACE_SHEAR_STIFFNESS))
+        << "No interface shear stiffness is defined" << std::endl;
+
+    KRATOS_ERROR_IF_NOT(rMaterialProperties[INTERFACE_SHEAR_STIFFNESS] > 0.0)
+        << "Interface shear stiffness must be positive, but got "
+        << rMaterialProperties[INTERFACE_SHEAR_STIFFNESS] << std::endl;


Just out of curiosity, would it every be reasonable that the stiffnesses depend on something else (and thus become time-dependent)?

Stiffness can depend on many things e.g. time, age, temperature, strain, stress etc. That is what you might use when doing computations that involve improving ground stability by grouting.

indigocoral

Thanks a lot for your work on this and also for adding documentation for it.
I've only reviewed the ReadMe file for now and my comments are more on the style of the documentation.

indigocoral · 2024-09-03T10:30:45Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+
+
+## Incremental linear elastic interface law
+The constitutive law for an interface element relates tractions $\tau$ to relative displacement $\Delta u$.


I would suggest to add the word 'linearly' to the sentence below to be more explicit about the nature of the relation:

Suggested change

The constitutive law for an interface element relates tractions $\tau$ to relative displacement $\Delta u$.

The constitutive law for an interface element linearly relates tractions $\tau$ to relative displacement $\Delta u$.

Yes, being more explicit about the nature of the relation may be clearer. Done.

indigocoral · 2024-09-03T10:37:34Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+
+### Relative displacement and traction
+In 2D plane strain $y$ is the opening/closing direction of the interface, while differential motion in the tangential direction
+gives shear. Like for a continuum, normal behaviour is placed first in the relative displacement and traction vector. Shear 


Suggested change

gives shear. Like for a continuum, normal behaviour is placed first in the relative displacement and traction vector. Shear

gives shear. Similar to the continuum, the normal behaviour is placed first in the relative displacement and traction vector. The shear

indigocoral · 2024-09-03T10:38:44Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+### Relative displacement and traction
+In 2D plane strain $y$ is the opening/closing direction of the interface, while differential motion in the tangential direction
+gives shear. Like for a continuum, normal behaviour is placed first in the relative displacement and traction vector. Shear 
+is placed after the normal motion or traction.


Suggested change

is placed after the normal motion or traction.

behavior follows after the normal motion or traction in these vectors.

indigocoral · 2024-09-03T10:49:21Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+$$ \tau = \begin{bmatrix} \tau_{yy} \\ \tau_{xy} \end{bmatrix} $$
+
+### 2D Elastic constitutive tensor
+


I would suggest to add something here just for a clearer documentation:

The elastic behavior of the interface is characterized by the elastic constitutive tensor (or 2D Elastic constitutive tensor) 𝐶, which relates the traction and relative displacement vectors. The constitutive tensor in 2D is expressed as:

indigocoral · 2024-09-03T10:50:47Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+Where $C_{yy}$ is input as `INTERFACE_NORMAL_STIFFNESS` and $C_{xy}$ is input as `INTERFACE_SHEAR_STIFFNESS`. Both stiffness values have dimension $\mathrm{F/L^3}$.
+
+### Incremental relation
+


Also here, I suggest again to add something like:

The relation between the traction and the relative displacement at a given time increment is given by the following incremental equation:

indigocoral · 2024-09-03T11:04:02Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+
+$$ \Delta u = \begin{bmatrix} \Delta u_y \\ \Delta u_x \end{bmatrix} $$
+
+$$ \tau = \begin{bmatrix} \tau_{yy} \\ \tau_{xy} \end{bmatrix} $$


This is more style and I only suggest it because sometimes you just want to quickly check the formula and what the parameters stand for. So I would suggest adding the following: (if correct to my understanding)

Where:

$\Delta u_𝑦$: Relative displacement in the 𝑦-direction (normal to the interface).

$\Delta u_𝑥$: Relative displacement in the 𝑥-direction (tangential to the interface).

$\tau_{𝑦𝑦}$: Traction in the 𝑦-direction (normal to the interface).

$\tau_{𝑥𝑦}$: Shear traction in the 𝑥-direction (tangential to the interface).

indigocoral · 2024-09-03T11:08:13Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+$$ C = \begin{bmatrix} C_{yy} & 0     \\
+                       0     & C_{xy} \end{bmatrix}$$
+
+Where $C_{yy}$ is input as `INTERFACE_NORMAL_STIFFNESS` and $C_{xy}$ is input as `INTERFACE_SHEAR_STIFFNESS`. Both stiffness values have dimension $\mathrm{F/L^3}$.


Also for this one:

Where:

$C$: Represents the 2D elastic constitutive tensor.

$C_{𝑦𝑦}$: Represents the INTERFACE_NORMAL_STIFFNESS, which characterizes the stiffness in the normal direction.

$C_{𝑥𝑦}$: Represents the INTERFACE_SHEAR_STIFFNESS, which characterizes the stiffness in the tangential direction.

Both stiffness values have dimension $\mathrm{F/L^3}$.

indigocoral · 2024-09-03T11:23:26Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+### Incremental relation
+
+$$ \tau_{t + \Delta t} = \tau_t + C \cdot \Delta u $$
+


Also for this one:

Where:

$\tau_{t + \Delta t}$: The traction at the updated time 𝑡 + Δ𝑡.

$\tau_t$: The traction at the current time 𝑡.

$C$: The elastic constitutive tensor.

$\Delta u$: The incremental relative displacement vector.

indigocoral · 2024-09-03T11:27:25Z

applications/GeoMechanicsApplication/custom_constitutive/README.md

+## Incremental linear elastic interface law
+The constitutive law for an interface element relates tractions $\tau$ to relative displacement $\Delta u$.
+Relative displacement for interface element is the differential motion between the two sides of the interface. As a
+consequence the relative displacement has unit of length [$\mathrm{L}$] and the stiffness has unit of force over cubic length [$\mathrm{F/L^3}$].


Suggested change

consequence the relative displacement has unit of length [$\mathrm{L}$] and the stiffness has unit of force over cubic length [$\mathrm{F/L^3}$].

consequence the relative displacement has unit of length $[\mathrm{L}]$ and the stiffness has unit of force over cubic length $[\mathrm{F/L^3}]$.

avdg81 · 2024-09-03T12:18:56Z

Many thanks to @mnabideltares, @rfaasse, and @indigocoral for your valuable suggestions. We feel that those have improved the quality of this PR.

rfaasse

Thank you for the hard work, to me this looks ready to be merged and used!

avdg81 and others added 21 commits August 27, 2024 13:58

Added an almost empty linear elastic law for line interfaces

af84a00

Added an override for GetStrainSize

e866d4e

Added an override for GetStressMeasure

989f286

Added an override for IsIncremental

a019688

Check whether the normal stiffness is defined in the material properties

13e2a4f

Check whether the shear stiffness is defined in the material properties

6702826

Check whether the normal and shear stiffnesses are positive

6fba772

Check whether given geometry is a line interface geometry

5ba2665

Check whether the computed traction follows from the relative displac…

c348699

…ement and stiffness values

Account for the previous traction and relative displacement

73f50de

Fixed two code smells by moving the initialization of the previous re…

81a9482

…lative displacement and previous traction from the constructor to `InitializeMaterial`

The relative displacement and traction at the end of a step are now k…

3fa1047

…ept by the law itself by using `FinalizeMaterialResponseCauchy`

Instances of the new constitutive law for interfaces can now be cloned

920569a

Added support for serialization of the constitutive law for interfaces

ac22a22

Added support for supplying an initial state

bb182d2

Marked two member functions [[nodiscard]]

ec502cf

Processed some review comments

140330b

- When checking the constitutive law, also call the `Check` member of the base class. - Replaced a few intermediate variables by using temporaries (return by other functions).

Extracted a member that constructs the constitutive matrix

b51f063

First setup README for incremental linear elastic interface material.

b2b1e7e

Reordered the member functions

32a002e

They now have the same order as in class `ConstitutiveLaw`.

A few minor changes to the test code

e563262

- Renamed a test case. - Set a Z coordinate to zero (since we assume a 2D model). - Introduced a constant for the adopted relative tolerance.

avdg81 added the GeoMechanics Issues related to the GeoMechanicsApplication label Aug 30, 2024

avdg81 requested review from rfaasse, indigocoral and mnabideltares August 30, 2024 07:16

avdg81 assigned avdg81 and WPK4FEM Aug 30, 2024

avdg81 linked an issue Aug 30, 2024 that may be closed by this pull request

[GeoMechanicsApplication] Add a linear elastic constitutive law for line interface elements using an incremental formulation #12650

Closed

mnabideltares reviewed Aug 30, 2024

View reviewed changes

Added a few debug statements

15f6f1f

avdg81 added 6 commits September 1, 2024 10:50

Attempt to fix unit test failure on GCC and Clang

c402de4

Second attempt to fix unit test failure on GCC and Clang

588a65a

Check effect of early BOOST_UBLAS_CHECK_FALSE

92a425f

Print value of BOOST_UBLAS_CHECK_ENABLE

41db940

Third attempt to fix unit test failure for GCC and Clang builds

6edfa34

Use ternary operator rather than if else construct

44d0f19

rfaasse reviewed Sep 2, 2024

View reviewed changes

avdg81 added 4 commits September 2, 2024 17:21

Defined unit test wide tolerance values

310d120

Extended a unit test case to include an initial state (to stress the …

babf7d4

…incremental nature of the law)

Processed Mohamed's review comments about the README file

f446f80

Minor improvement of error messages

68f4e2e

avdg81 requested review from rfaasse and mnabideltares September 3, 2024 08:24

rfaasse reviewed Sep 3, 2024

View reviewed changes

Fixed formatting of inline math expression

246f141

indigocoral reviewed Sep 3, 2024

View reviewed changes

Processed Marjan's suggestions

3709231

avdg81 requested a review from rfaasse September 3, 2024 12:20

avdg81 added 2 commits September 3, 2024 14:33

Even more improvements by Wijtze Pieter

d1bbd5d

Beter opletten Van de Graaf ;-)

ee43d50

rfaasse approved these changes Sep 3, 2024

View reviewed changes

avdg81 enabled auto-merge (squash) September 3, 2024 12:54

avdg81 merged commit a6aeead into master Sep 3, 2024
9 of 11 checks passed

avdg81 deleted the geo/12650-add-linear-elastic-constitutive-law-for-line-interfaces branch September 3, 2024 14:19


		$$ \Delta u = \begin{bmatrix} \Delta u_y \\ \Delta u_x \end{bmatrix} $$

		$$ \tau = \begin{bmatrix} \tau_{yy} \\ \tau_{xy} \end{bmatrix} $$

		$$ C = \begin{bmatrix} C_{yy} & 0 \\
		0 & C_{xy} \end{bmatrix}$$



		## Incremental linear elastic interface law
		The constitutive law for an interface element relates tractions $\tau$ to relative displacement $\Delta u$.

	gives shear. Like for a continuum, normal behaviour is placed first in the relative displacement and traction vector. Shear
	gives shear. Similar to the continuum, the normal behaviour is placed first in the relative displacement and traction vector. The shear

	is placed after the normal motion or traction.
	behavior follows after the normal motion or traction in these vectors.

		$$ \tau = \begin{bmatrix} \tau_{yy} \\ \tau_{xy} \end{bmatrix} $$

		### 2D Elastic constitutive tensor

		Where $C_{yy}$ is input as `INTERFACE_NORMAL_STIFFNESS` and $C_{xy}$ is input as `INTERFACE_SHEAR_STIFFNESS`. Both stiffness values have dimension $\mathrm{F/L^3}$.

		### Incremental relation

		### Incremental relation

		$$ \tau_{t + \Delta t} = \tau_t + C \cdot \Delta u $$

	consequence the relative displacement has unit of length [$\mathrm{L}$] and the stiffness has unit of force over cubic length [$\mathrm{F/L^3}$].
	consequence the relative displacement has unit of length $[\mathrm{L}]$ and the stiffness has unit of force over cubic length $[\mathrm{F/L^3}]$.

[GeoMechanicsApplication] Add a linear elastic constitutive law for line interfaces #12654

[GeoMechanicsApplication] Add a linear elastic constitutive law for line interfaces #12654

Conversation

avdg81 commented Aug 30, 2024

mnabideltares left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mnabideltares Aug 30, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rfaasse left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rfaasse left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

indigocoral left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

avdg81 commented Sep 3, 2024

rfaasse left a comment

Choose a reason for hiding this comment

mnabideltares Aug 30, 2024 •

edited

Loading