Wave propagation in hyperelastic media

Sorin Mitran (mitran@unc.edu)

Department of Mathematics, University of North Carolina at Chapel Hill

September 18, 2019

Hyperelastic media allow derivation of a constitutive stress-strain relationship from the strain energy density. This offers the opportunity for automated generation of the source code for a solver for such media using the wave propagation approach of LeVeque, a variant of the finite volume method based upon identification of the eigenmodes of a hyperbolic system of partial differential equations. Momentum conservation for hyperelastic media can be formulated as
(Subscript[(ρ v), t])=div σ(F)
(Subscript[F, t])=div (v⊗δ)
with ρ - mass density, v - velocity, σ Cauchy stress, F - deformation gradient, v⊗δ - exterior product of velocity with identity matrix. The above system can be rewritten as
(Subscript[q, t])+div f(q)=0
with the local linearization
(Subscript[q, t])+A·div(q)=0. 
The eigenmodes of the flux Jacobian A represent the displacement waves of the medium. This notebook defines finite strain tensors, computes the eigendecomposition A R=R Λ , and generates Fortran source code to solve the linearized Riemann problem Rδc=δq. The source code is included by the problem module of BEARCLAW to carry out a numerical solution.

Notation

Tensors are generally denoted by double strike Latin letters such as: 𝕩, 𝔽 , with keyboard shortcuts EscdsxEsc, EscdsFEsc. The scalar components are denoted by the indexed, corresponding normal font Latin letter, 𝕩={(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}. Tensors with standard notation by Greek letters are denoted by the corresponding script letter ℯ, with keyboard shortcut EscdseEsc, and scalar components denoted by the indexed Greek letters, ℯ={{(Subscript[ϵ, 11]),(Subscript[ϵ, 12]),(Subscript[ϵ, 13])},{(Subscript[ϵ, 21]),(Subscript[ϵ, 22]),(Subscript[ϵ, 23])},{(Subscript[ϵ, 31]),(Subscript[ϵ, 32]),(Subscript[ϵ, 33])}}. When the scalar components of the tensor are not used, a tensor can also be denoted by a normal face Greek letter.

Define number of space dimensions

d=3

Define identity matrix, Dirac delta

𝕀=IdentityMatrix[d]; δ=𝕀; MatrixForm[δ]

Define some utility functions for Einstein summation notation

idx2[i_,j_]:=ToString[i]<>ToString[j];
idx3[i_,j_,k_]:=ToString[i]<>ToString[j]<>ToString[k];
idx4[i_,j_,k_,l_]:=ToString[i]<>ToString[j]<>ToString[k]<>ToString[l];

Deformation measures

Displacement

Define reference, current coordinates as independent variables

𝕏=Table[(Subscript[X, i]),{i,d}]; 𝕩=Table[(Subscript[x, i]),{i,d}]; TableForm[{𝕏,𝕩}]

Define current coordinates as dependent variables

𝕩𝕩=Table[𝕏[[i]]+(Subscript[u, i])[t,𝕏],{i,d}]; TableForm[𝕩𝕩]

Define displacements as independent variables

𝕦=Table[(Subscript[u, i]),{i,d}]; TableForm[{𝕦}]

Define displacements as dependent variables w.r.t. reference coordinates

𝕦𝕏=Table[(Subscript[u, i])[t,𝕏],{i,d}]; TableForm[{𝕦𝕏}]

Deformation tensors

Define deformation gradient in terms of displacements taken as independent variables

𝔽=δ+Table[(Subscript[u, i,j]),{i,d},{j,d}]; TableForm[𝔽]

Define deformation gradient with displacements taken as variables that depend on reference configuration

𝔽𝕏=δ+Table[D[𝕦𝕏[[i]],𝕏[[j]]],{i,d},{j,d}]; TableForm[𝔽𝕏]

Define left Cauchy-Green deformation tensor in terms of deformation gradient

𝔹 = 𝔽.Transpose[𝔽]; TableForm[𝔹]

Define right Cauchy-Green deformation tensor in terms of displacement functions

𝔹𝕏 = 𝔽𝕏.Transpose[𝔽𝕏]; TableForm[𝔹𝕏]

Define right Cauchy-Green deformation tensor in terms of deformation gradient

ℂ = Transpose[𝔽].𝔽; TableForm[ℂ]

Define right Cauchy-Green deformation tensor in terms of displacement functions

ℂ𝕏 = Transpose[𝔽𝕏].𝔽𝕏; TableForm[ℂ𝕏]

Strain

Linear engineering strain

In terms of deformation gradient

ℯ=Table[(1/2)((Subscript[u, i,j])+(Subscript[u, j,i])),{i,d},{j,d}]; TableForm[ℯ]

In terms of displacement functions

ℯ𝕏=Table[(1/2)(D[𝕦𝕏[[i]],(Subscript[X, j])]+D[𝕦𝕏[[j]],(Subscript[X, i])]),{i,d},{j,d}]; TableForm[ℯ𝕏]

Green-Lagrange (Green-St. Venant) strain

In terms of deformation gradient

𝔼=Expand[(1/2)(ℂ-𝕀)]; TableForm[𝔼]

In terms of displacement functions

𝔼𝕏=Expand[(1/2)(ℂ𝕏-𝕀)]; TableForm[𝔼𝕏]

Constitutive relations

Linear elastic (Hooke)

(Subscript[ℂ, LE])=Table[λ δ[[i,j]] δ[[k,l]]+μ(δ[[i,k]] δ[[j,l]]+δ[[i,l]] δ[[j,k]]),{i,d},{j,d},{k,d},{l,d}]; TableForm[(Subscript[ℂ, LE])]

Strain energy

Define potential for various types of elastic materials

Linear elastic (Hooke)

(Subscript[ℰ, LE])=(1/2)Sum[(Subscript[ℂ, LE])[[i,j,k,l]] ℯ[[i,j]] ℯ[[k,l]],{i,d},{j,d},{k,d},{l,d}]

(Subscript[𝒰, LE])=Expand[(Subscript[ℰ, LE])]

Isotropic soft solid shear modes

Define invariants of Green-Lagrange strain tensor

(Subscript[I, 2])=Sum[𝔼[[i,j]] 𝔼[[j,i]],{i,d},{j,d}]

(Subscript[𝒰, 2])=Expand[(Subscript[I, 2])]

(Subscript[I, 3])=Sum[𝔼[[i,j]] 𝔼[[j,k]] 𝔼[[k,i]],{i,d},{j,d},{k,d}]

(Subscript[𝒰, 3])=Expand[(Subscript[I, 3])]

(Subscript[ℰ, SS])=μ (Subscript[I, 2])+(1/3)𝒜 (Subscript[I, 3])+𝒟 ((Subscript[I, 2])^2);

(Subscript[𝒰, SS])=Expand[(Subscript[ℰ, SS])]

Cauchy stress

The Cauchy stress is the gradient of the strain energy w.r.t. deformation gradient

σ[𝒰_]:=Monitor[Table[Simplify[D[𝒰,(Subscript[u, i,j])]],{i,d},{j,d}],{i,j}]

Linear elastic

(Subscript[σ, LE])=σ[(Subscript[𝒰, LE])]; TableForm[(Subscript[σ, LE])]

In terms of the deformation gradient components

(Subscript[σ𝔽, LE])=(Subscript[σ, LE]) /. Flatten[Table[(Subscript[u, i,j])->(Subscript[F, idx2[i,j]])-δ[[i,j]],{i,d},{j,d}]];TableForm[(Subscript[σ𝔽, LE])]

sigmaFLE=(Subscript[σ𝔽, LE]);
DumpSave["sigmaF_LE.mx",sigmaFLE];

Isotropic soft solid shear modes

(Subscript[σ, SS])=σ[(Subscript[𝒰, SS])];

The expressions for the stress components are very long. Individual components can be accessed as follows:

(Subscript[σ, SS])[[1,1]] /. {𝒜->0,𝒟->0}

(Subscript[σ𝔽, SS])=(Subscript[σ, SS]) /. Flatten[Table[(Subscript[u, i,j])->(Subscript[F, idx2[i,j]])-δ[[i,j]],{i,d},{j,d}]];

sigmaFSS=(Subscript[σ𝔽, SS]);
DumpSave["sigmaF_SS.mx",sigmaFSS];

Linearly polarized shear wave

To check above computations, extract expressions for a linearly polarized shear wave in which only (Subscript[u, 1,3]) is non-zero

setZero = Flatten[Table[(Subscript[u, i,j])->0,{i,d},{j,d}]]

setZeroNot13 = Drop[setZero,{3,3}]

Simplify[(Subscript[σ, SS]) /. setZero]

(Subscript[σ, SSLP])=Simplify[(Subscript[σ, SS]) /. setZeroNot13]; TableForm[(Subscript[σ, SSLP])[[1,3]]]

Momentum conservation is stated as ρ (Subscript[v, 1,t])=(Subscript[σ, 1j,j])=(∂_3)[ μ (Subscript[u, 1,3])+(μ+(𝒜/2)+𝒟)(Subscript[u, 1,3]^3)+(1/4)(𝒜+6 𝒟)(Subscript[u, 1,3]^5)+(1/2) 𝒟 (Subscript[u, 1,3]^7)], which differs from Zabolotskaya et al (2004) in the terms (Subscript[u, 1,3]^5) and (Subscript[u, 1,3]^7). Verify result against different computational steps.

setZeroD23 = Flatten[Table[(Subscript[u, i,j])->0,{i,d},{j,1,2}]]

setZeroU23=Flatten[Table[(Subscript[u, i,j])->0,{i,2,d},{j,1,d}]]

(Subscript[𝒰, SSLP])=(Subscript[𝒰, SS]) /. Flatten[{setZeroD23,setZeroU23}]

(Subscript[σ, SSLP2])=σ[(Subscript[𝒰, SSLP])]; TableForm[(Subscript[σ, SSLP2])]

Simplify[D[(Subscript[𝒰, SSLP]),(Subscript[u, 1,3])]]

Elliptically polarized shear wave

As a further check, extract expressions for an elliptically polarized shear wave in which only (Subscript[u, 1,3]), (Subscript[u, 2,3]) are non-zero

setZero = Flatten[Table[(Subscript[u, i,j])->0,{i,d},{j,d}]]

setZeroNot1323 = Drop[Drop[setZero,{3,3}],{5,5}]

(Subscript[σ, SSEP])=Simplify[(Subscript[σ, SS]) /. setZeroNot1323]; TableForm[(Subscript[σ, SSEP])[[1;;2,3]]]

The above again differs from Zabolotskaya et al (2004). Alternative computational steps

setZeroD23 = Flatten[Table[(Subscript[u, i,j])->0,{i,d},{j,1,2}]]

setZeroU3=Flatten[Table[(Subscript[u, i,j])->0,{i,3,d},{j,1,d}]]

(Subscript[𝒰, SSEP])=(Subscript[𝒰, SS]) /. Flatten[{setZeroD23,setZeroU3}]

(Subscript[σ, SSEP2])=σ[(Subscript[𝒰, SSEP])]; TableForm[(Subscript[σ, SSEP2])]

Simplify[D[(Subscript[𝒰, SSEP]),(Subscript[u, 1,3])]]

Simplify[D[(Subscript[𝒰, SSEP]),(Subscript[u, 2,3])]]

Tensions

Tensions (forces per unit volume) in the directions of the reference coordinates are obtained by taking divergence of Cauchy stress

Utility functions

Define differentiation of deformation tensor

dF[σF_,j_]:= Simplify[D[σF /. Flatten[Table[ (Subscript[F, idx2[k,l]])->(Subscript[F, idx2[k,l]])[x],{k,d},{l,d}]],x] /. Flatten[{Table[ (Subscript[F, idx2[k,l]])'[x]->(Subscript[F, idx2[k,l],j]),{k,d},{l,d}],Table[ (Subscript[F, idx2[k,l]])[x]->(Subscript[F, idx2[k,l]]),{k,d},{l,d}]}]]

Define divergence of stress tensor expressed in terms of 𝔽

𝕗[σ_]:=Monitor[Table[Sum[dF[σ[[i,j]],j],{j,d}],{i,d}],{i,j}]

Replace 𝔽 and spatial derivative of 𝔽 with 𝕦

𝕦𝔽[σ_] := σ /. Flatten[{
Table[ (Subscript[F, idx2[k,l]])->(Subscript[u, k,l]),{k,d},{l,d}],
Table[ (Subscript[F, idx2[k,l],m])->(Subscript[u, k,idx2[Min[{l,m}],Max[{l,m}]]]),{k,d},{l,d},{m,d}]
}]

Tensions for linear elastic solid

(Subscript[𝕗, LE])=𝕗[(Subscript[σ𝔽, LE])]; TableForm[(Subscript[𝕗, LE])]

fFLE=(Subscript[𝕗, LE]);
DumpSave["fF_LE.mx",fFLE];

Verify against vector formula (λ+μ) ∇ ∇·u + μ (∇^2)u

(Subscript[𝕗𝕦, LE])=Map[𝕦𝔽[#]&,(Subscript[𝕗, LE])];  TableForm[(Subscript[𝕗𝕦, LE])]

form =  Table[Simplify[Sum[ (λ+μ) (Subscript[u, j,idx2[Min[{i,j}],Max[{i,j}]]])+ μ (Subscript[u, i,idx2[j,j]]),{j,d}]],{i,d}]

Simplify[(Subscript[𝕗𝕦, LE])-form ]

Tensions for isotropic soft solid shear modes

(Subscript[𝕗, SS])=𝕗[(Subscript[σ𝔽, SS])];

fFSS=(Subscript[𝕗, SS]);
DumpSave["fF_SS.mx",fFSS];

Flux Jacobians

For an incompressible medium in which all model parameters are scaled by the density and deformations are sufficiently small to enable use of the Cauchy stress tensor (as opposed to the first Piola-Kirchhoff stress), the equations of motion (Subscript[v, i,t])=(Subscript[σ, i j,j]) coupled with rate of deformation equations (Subscript[F, i j,j])=(Subscript[v, i,j]) form a conservative system (Subscript[q, t])-∇·f=0 with  q=(((v  F))^T),f=(((σ  v⊗δ))^T), or in component form (Subscript[v, i,t])=(Subscript[σ, i k,k]), (Subscript[F, i j,t])=(Subscript[((Subscript[v, i]) (Subscript[δ, j k])), ,k]).

To use a wave propagation approach to the numerical solution of the equations of motion, the eigenmodes of the system flux Jacobians are required

f[σ_,v_]:=Transpose[Join[σ,Flatten[Outer[Times,v,δ],1]]];

sig=Table[(Subscript[σ, idx2[i,j]]),{i,d},{j,d}];
vel=Table[(Subscript[v, i]),{i,d}];
flux=f[sig,vel]; TableForm[Transpose[flux]]

Define gradient w.r.t. q

q=Flatten[{Table[(Subscript[v, i]),{i,d}],Table[(Subscript[F, idx2[i,j]]),{i,d},{j,d}]}]; nq=Length[q]; TableForm[q,TableDirections->Row]

gradq[f_] := Monitor[Table[D[f,q[[α]]],{α,nq}],α]

Linear elastic

(Subscript[f, LE])=f[(Subscript[σ𝔽, LE]),vel]; TableForm[(Subscript[f, LE]),TableDirections->Row]

Flux Jacobians

(Subscript[A, LE])=Map[Transpose[gradq[#]]&,(Subscript[f, LE])]; TableForm[(Subscript[A, LE])[[1]]]

TableForm[(Subscript[A, LE])[[2]]]

TableForm[(Subscript[A, LE])[[3]]]

fLE=(Subscript[f, LE]); ALE=(Subscript[A, LE]);
DumpSave["qfA_LE.mx",{q,fLE,ALE}];

Eigensystem - eigenvalues

(Subscript[ΛR, LE])=Map[Eigensystem[#]&,(Subscript[A, LE])] /. {(Sqrt[λ+2 μ])->(Subscript[c, P]),(Sqrt[μ])->(Subscript[c, S])} ;
(Subscript[Λ, LE])=Table[(Subscript[ΛR, LE])[[i,1]],{i,d}];
(Subscript[R, LE])=Table[(Subscript[ΛR, LE])[[i,2]],{i,d}];
TableForm[(Subscript[Λ, LE])]

Eigensystem - eigenvectors (displayed as rows in Mathematica)

TableForm[(Subscript[R, LE])[[1]]]

TableForm[(Subscript[R, LE])[[2]]]

TableForm[(Subscript[R, LE])[[3]]]

Verify eigensystem

TableForm[Table[(Subscript[A, LE])[[i]].(Subscript[R, LE])[[i,α]]-(Subscript[Λ, LE])[[i,α]] (Subscript[R, LE])[[i,α]] /. {(Subscript[c, P])->(Sqrt[λ+2 μ]),(Subscript[c, S])->(Sqrt[μ])},{α,nq},{i,d}] , TableDirections->Row]

Decomposition of δq onto eigenbasis, R c = δq

δ𝕢=Table[(Subscript[δq, α]),{α,nq}]; TableForm[δ𝕢,TableDirections->Row]

Define coefficients c

𝕔=Table[LinearSolve[Transpose[R[[i]]],δ𝕢],{i,d}];

Select propagating waves (place arbitrary positive values for wave velocities into eigenvalues)

Λ12 = Λ/.{(Subscript[c, S])->1,(Subscript[c, P])->2};
iW=Table[Select[Table[α,{α,nq}], Λ12[[i,#]]≠0&],{i,d}]; TableForm[iW]

Number of propagating waves. Should be the same for all d spatial dimensions.

nW=Length[iW[[1]]]

λW=Table[Λ[[i,iW[[i]]]],{i,d}]; TableForm[λW]

W=Table[R[[i,iW[[i]]]],{i,d}]; TableForm[W[[1]]]

TableForm[W[[2]]]

TableForm[W[[3]]]

cW=Table[𝕔[[i,iW[[i]]]],{i,d}]; TableForm[cW]

DumpSave["lcW_LE.mx",{λW,W,cW}];

Isotropic soft solid shear modes

(Subscript[f, SS])=f[(Subscript[σ𝔽, SS]),vel]; TableForm[(Subscript[f, SS]),TableDirections->Row]

Dimensions[(Subscript[f, SS])]

(Subscript[f, SS])[[1,1]]

TableForm[(Subscript[f, SS])[[;;,4;;]],TableDirections->Row]

Flux Jacobians

(Subscript[A, SS])=Map[Transpose[gradq[#]]&,(Subscript[f, SS])];

Dimensions[(Subscript[A, SS])]

TableForm[{ (Subscript[A, SS])[[1,4;;,;;]],(Subscript[A, LE])[[1,4;;,;;]]} ,TableDirections->Row]

TableForm[{ (Subscript[A, SS])[[2,4;;,;;]],(Subscript[A, LE])[[2,4;;,;;]]} ,TableDirections->Row]

TableForm[{ (Subscript[A, SS])[[3,4;;,;;]],(Subscript[A, LE])[[3,4;;,;;]]} ,TableDirections->Row]

Eigensystem - eigenvalues

ΛR=Map[Eigensystem[#]&,A];
Λ=Table[ΛR[[i,1]],{i,d}];
R=Table[ΛR[[i,2]],{i,d}];
TableForm[Λ]

Eigensystem - eigenvectors (displayed as rows in Mathematica)

TableForm[R[[1]]]

TableForm[R[[2]]]

TableForm[R[[3]]]

Verify eigensystem

TableForm[Table[A[[i]].R[[i,α]]-Λ[[i,α]] R[[i,α]] /. {(Subscript[c, P])->(Sqrt[λ+2 μ]),(Subscript[c, S])->(Sqrt[μ])},{α,nq},{i,d}] , TableDirections->Row]

Decomposition of δq onto eigenbasis, R c = δq

δ𝕢=Table[(Subscript[δq, α]),{α,nq}]; TableForm[δ𝕢,TableDirections->Row]

Define coefficients c

Table[(Subscript[c, α]),{α,nq}]; TableForm[𝕔,TableDirections->Row]

𝕔=Table[LinearSolve[Transpose[R[[i]]],δ𝕢],{i,d}];

Select propagating waves (place arbitrary positive values for wave velocities into eigenvalues)

Λ12 = Λ/.{(Subscript[c, S])->1,(Subscript[c, P])->2};
iW=Table[Select[Table[α,{α,nq}], Λ12[[i,#]]≠0&],{i,d}]; TableForm[iW]

Number of propagating waves. Should be the same for all d spatial dimensions.

nW=Length[iW[[1]]]

λW=Table[Λ[[i,iW[[i]]]],{i,d}]; TableForm[λW]

W=Table[R[[i,iW[[i]]]],{i,d}]; TableForm[W[[1]]]

TableForm[W[[2]]]

TableForm[W[[3]]]

cW=Table[𝕔[[i,iW[[i]]]],{i,d}]; TableForm[cW]

Directory[]

Save["linelasWaves",{q,nW,W,cW}];