Wave propagation in hyperelastic media

Sorin Mitran (mitran@unc.edu)

Department of Mathematics, University of North Carolina at Chapel Hill

October 28, 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 P(F)
(Subscript[F, t])=div (v⊗δ)
with ρ - mass density, v - velocity, P first Piola-Kirchhoff 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
where q=(v,F), f=(P,v⊗δ), with the local linearization
(Subscript[q, t])+A·div(q)=0. 
The eigenmodes of the flux Jacobian A=∂f/∂q, represent the displacement waves of the medium. Hyperelastic media exhibit six propagating eigenmodes corresponding to speeds ±(Subscript[c, 1]),±(Subscript[c, 2])±(Subscript[c, 3]) (e.g., for linear elasticity, these are the P-mode and two S-modes). 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.

SetDirectory[NotebookDirectory[]]

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

m=(d^2)+d

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];

Flux Jacobians

For a hyperelastic medium the flux Jacobian A=∂f/∂q can be expressed in terms of the gradient of the (first Piola-Kirchhoff) P(F) stress with respect to the deformation gradient, P'=∂P/∂F.

Define q=(v,F)

𝕧=Table[(Subscript[v, i]),{i,d}]

𝔽=Table[(Subscript[F, idx2[i,j]]),{i,d},{j,d}]

LinearizeFTensor = Flatten[Table[
(Subscript[F, idx2[i,j]])->(Subscript[F, Position[Flatten[𝔽],(Subscript[F, idx2[i,j]])][[1,1]]]),{i,d},{j,d}
]]

𝕢=Flatten[{𝕧,𝔽}]

Define f(q)=(P(F),v⊗δ)

ℙ=Table[(Subscript[P, idx2[i,j]]),{i,d},{j,d}]; TableForm[ℙ]

LinearizePTensor = Flatten[Table[
(Subscript[P, idx2[i,j]])->(Subscript[P, Position[Flatten[ℙ],(Subscript[P, idx2[i,j]])][[1,1]]]),{i,d},{j,d}
]]

𝕧δ=Flatten[Outer[Times,𝕧,δ],1]; TableForm[𝕧δ]

𝕗=Transpose[Join[ℙ /. LinearizePTensor,𝕧δ,1]]; TableForm[Transpose[𝕗]]

fq[i_]:=Join[Table[0,{j,d},{k,d}],Flatten[
Transpose[Outer[Times,δ,δ[[i]]],{3,1,2}] ,1]  ];
dPdF[i_]:=Join[Table[(Subscript[dP, d(i-1)+j,k]),{j,d},{k,(d^2)}],Table[0,{j,(d^2)},{k,(d^2)}]];
A[i_]:=Join[fq[i],dPdF[i],2]

A[1] // TableForm

A[2] // TableForm

A[3] // TableForm

Eigensystem

Analytical derivation

Define a function to return the speeds and propagating eigenmodes of the conservative system (Subscript[q, t])+(Subscript[A, 1])(Subscript[q, ,1])+(Subscript[A, 2])(Subscript[q, ,2])+(Subscript[A, 3])(Subscript[q, ,3])=0.

PropagatingSpeedCoefs[A_]:=Module[{p,x,y},
p=CharacteristicPolynomial[A,x]/(x^6) /. {(x^2)->y,(x^4)->(y^2),(x^6)->(y^3)};
Return[CoefficientList[p,y]]
]

PropagatingSpeedCoefs[A[1]]

PropagatingSpeedCoefs[A[2]]

PropagatingSpeedCoefs[A[3]]

Eigenvalues[A[1]]

PropagatingModes[A_]:=Module[{eigs,speeds,roots,denom,k,λR,vec,modes={}},
eigs=Eigenvalues[A];
speeds = {Part[eigs[[8]],1]->(c1^2),Part[eigs[[10]],1]->(c2^2),Part[eigs[[12]],1]->(c3^2)};
roots={(Sqrt[(c1^2)])->c1,(Sqrt[(c2^2)])->c2,(Sqrt[(c3^2)])->c3};
λR=Eigensystem[A];
Monitor[
For[k=1,k≤Length[A],k++,
vec= λR[[2,k]] /. speeds /. roots;
AppendTo[modes,Simplify[Denominator[Together[Total[vec]]] vec]]],k];
Return[modes];
]

R=Monitor[Table[PropagatingModes[A[i]],{i,d}],i];

TableForm[R[[1]]]

Verify that the first 6 eigenvectors indeed correspond to null eigenvalues

Table[Simplify[A[i].R[[i,j]]],{i,d},{j,2d}]

TeXRules=Flatten[{c1->(Subscript[c, 1]),c2->(Subscript[c, 2]),c3->(Subscript[c, 3]),Table[(Subscript[dP, i,j])->(Subscript[P, i,j]),{i,(d^2)},{j,(d^2)}]}];

R[[1]] /. TeXRules // TableForm

Export["Eigenmodes.tex",TeXForm[R[[1,1;;6]] /. TeXRules],"TeX"]

Generate Fortran code for the general form of the eigensystem

Define Fortran notation for gradient of first Piola-Kirchhoff stress with respect to deformation gradient

dPdFrules = Flatten[Table[(Subscript[dP, i,j])->dPdF[i,j],{i,(d^2)},{j,(d^2)}]];

Eigenvalues

Write the coefficients for the cubic polynomial giving the elastic wave speeds

CubicCoefs[]:=Module[{a,i,k,n,coef,instr,subr},
subr = "SUBROUTINE CubicCoefs(k,dPdF,a)\n";
subr = subr <> "  INTEGER k\n";
subr = subr <> "  REAL (KIND=8) dPdF(9,9),a(4)\n";
subr = subr <> "  SELECT CASE(k)\n";
For[k=1,k≤d,k++,
subr = subr <> "    CASE(" <> ToString[k] <> ")\n";
a= PropagatingSpeedCoefs[A[k]];
For[i=1,i≤Length[a],i++,
coef=ToString[FortranForm[a[[i]]] /. dPdFrules];
instr = "      a(" <> ToString[i] <> ") = " <>  coef <> "\n";
subr = subr <> instr;
];
];
subr = subr <> "  END SELECT\n";
subr = subr <> "END SUBROUTINE CubicCoefs\n";
Export["SpeedsCubicCoefs.f90",subr,"Text"];
Return[subr]
]

CubicCoefs[]

Eigenvectors

Write the coefficients for the cubic polynomial giving the elastic wave speeds

Eigenmodes[]:=Module[{i,j,W,k,n,coef,instr,subr},
subr = "SUBROUTINE Eigenmodes(k,c1,c2,c3,dPdF,R)\n";
subr = subr <> "  INTEGER k\n";
subr = subr <> "  REAL (KIND=8) c1,c2,c3,dPdF(9,9),R(12,12)\n";
subr = subr <> "  SELECT CASE(k)\n";
For[k=1,k≤d,k++,
subr = subr <> "    CASE(" <> ToString[k] <> ")\n";
W= R[[k]];
For[i=1,i≤Length[W],i++,
For[j=1,j≤Length[W[[i]]],j++,
coef=ToString[FortranForm[W[[j,i]]] /. dPdFrules];
instr = "      R(" <> ToString[j] <> "," <> ToString[i] <> ") = " <>  coef <> "\n";
subr = subr <> instr;
];
];
];
subr = subr <> "  END SELECT\n";
subr = subr <> "END SUBROUTINE Eigenmodes\n";
Export["Eigenmodes.f90",subr,"Text"];
Return[subr]
]

Eigenmodes[];

Linear elastic

Verify above general formulation for linear elastic material

(Subscript[P, EL])=Simplify[λ (Sum[(Subscript[F, idx2[l,l]]),{l,d}]-3)δ+μ Table[(Subscript[F, idx2[i,k]])+(Subscript[F, idx2[k,i]])-2 δ[[i,k]],{i,d},{k,d}]] /. LinearizeFTensor;
TableForm[(Subscript[P, EL])]

ElasticStressGradient=Flatten[Table[(Subscript[dP, d(i-1)+j,k])->D[(Subscript[P, EL])[[i,j]],(Subscript[F, k])],{i,d},{j,d},{k,(d^2)}]]

A[1] /. ElasticStressGradient // TableForm

Eigenvalues[A[1] /. ElasticStressGradient]

Eigenvectors[A[1] /. ElasticStressGradient] /. {(Sqrt[μ])->(Subscript[c, S]),(Sqrt[λ+2 μ])->(Subscript[c, P])}// TableForm

The above are identical to results from StainEnergyPotential.nb

Soft solid

Define potential

Define Green-Lagrange strain

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

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

(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])]

Compute first Piola-Kirchhoff stress

(Subscript[P, SS])=Monitor[Table[Simplify[D[(Subscript[𝒰, SS]),(Subscript[F, idx2[i,j]])]] /. LinearizeFTensor,{i,d},{j,d}],{i,j}];

SoftSolidStressGradient=Monitor[Flatten[Table[(Subscript[dP, d(i-1)+j,k])->D[(Subscript[P, SS])[[i,j]],(Subscript[F, k])],{i,d},{j,d},{k,(d^2)}]],k]

MaxPower=1;
Monitor[For[i=1,i≤Length[SoftSolidStressGradient],i++,
MaxPower = Max[Flatten[{MaxPower,Table[Exponent[SoftSolidStressGradient[[i,2]],(Subscript[F, j])],{j,(d^2)}]}]]],i]

MaxPower

FRules = Flatten[{{𝒜->A,𝒟->D,μ->mu},
Table[(Subscript[F, k]^p)->F[k,p],{k,(d^2)},{p,MaxPower}]}]

dPdFSS[]:=Module[{i,j,W,k,n,coef,instr,subr},
subr = "SUBROUTINE StressGradient(MatParams,F,dPdF)\n";
subr = subr <> "  REAL (KIND=8) MatParams(100),F(" <> ToString[(d^2)]<> ","<> ToString[MaxPower]<> "),dPdF(" <> ToString[(d^2)]<> ","  <> ToString[(d^2)]<>")\n";
subr = subr <> "  INTEGER i,j,MaxPower\n";
subr = subr <> "  REAL (KIND=8) A,D,mu\n";
subr = subr <> "  A=MatParams(1); D=MatParams(2); mu=MatParams(3)\n";
subr = subr <> "  MaxPower=" <> ToString[MaxPower] <>"\n";
subr = subr <> "  DO i=2,MaxPower\n";
subr = subr <> "    DO j=1," <> ToString[(d^2)] <>"\n";
subr = subr <> "       F(j,i) = F(j,i-1)*F(j,1)\n";
subr = subr <> "    ENDDO\n";
subr = subr <> "  ENDDO\n";
For[i=1,i≤(d^2),i++,
For[j=1,j≤(d^2),j++,
coef=ToString[FortranForm[(Subscript[dP, i,j]) /. SoftSolidStressGradient] /. FRules];
instr = "  dPdF(" <> ToString[i] <> "," <> ToString[j] <>  ") = " <>  coef <> "\n";
subr = subr <> instr;
];
];
subr = subr <> "END SUBROUTINE StressGradient\n";
Export["StressGradient.f90",subr,"Text"];
Return[subr]
]

dPdFSS[];