Stress gradient: Soft solid medium, 2D

Define P(F), (Subscript[∇, F])P

SetDirectory[NotebookDirectory[]]

Import 2d flux definitions

NotebookEvaluate["/home/student/projects/ShearWaveAIR/NumericalModels/bearclaw/FluxDefinitions2d.nb"];

Strain energy

Define deformation tensors

Define Green-Lagrange strain and invariants

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

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

FtoGraduRules=Map[Reverse[#]&,LinearizeFTensor] /. Flatten[Table[(Subscript[F, idx2[i,j]])->δ[[i,j]]+(Subscript[u, i,j]),{i,d},{j,d}]]

Expand[𝔼 /. FtoGraduRules] // 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])]

Define zero-deformation state

ZeroDeformRules=Flatten[{Table[{(Subscript[F, 1]^k)->1,(Subscript[F, 2]^k)->0,(Subscript[F, 3]^k)->0,(Subscript[F, 4]^k)->1},{k,8}] }];

(Subscript[I, 2]) /. ZeroDeformRules

(Subscript[I, 3]) /. ZeroDeformRules

Define strain energy in terms of deformation invariants

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

ℰ /.  ZeroDeformRules

𝒰=Expand[ℰ]

𝒰 /. ZeroDeformRules

Gradient of strain energy (First Piola-Kirchhoff stress)

Pℰ=Table[Simplify[D[ℰ,(Subscript[F, i])]],{i,(d^2)}];

Manipulate[Pℰ[[i]],{i,1,(d^2),1}]

Manipulate[Pℰ[[i]] /. {𝒜->0,𝒟->0},{i,1,(d^2),1}]

Manipulate[Expand[Pℰ[[i]]/μ /. {𝒜->0,𝒟->0} /. FtoGraduRules],{i,1,(d^2),1}]

Linearly polarized wave (in 2D)

LinPolRules={(Subscript[F, 1])->1, (Subscript[F, 2])->0,(Subscript[F, 4])->1};
Manipulate[Expand[Pℰ[[i]]] /. LinPolRules,{i,1,(d^2),1}]

Zero deformation stress

Simplify[Pℰ /. ZeroDeformRules]

Hessian of strain energy (gradient of first Piola-Kirchhoff stress)

Hℰ=Table[Simplify[D[D[ℰ,(Subscript[F, i])],(Subscript[F, j])]],{i,(d^2)},{j,(d^2)}];

Manipulate[Hℰ[[i,j]],{i,1,(d^2),1},{j,1,(d^2),1}]

Zero deformation stress gradient

Simplify[Hℰ /. ZeroDeformRules]

Average Hessian at interface between left/right states

Define integration path between left and right states

FpathRules=Table[(Subscript[F, i])->(Subscript[F, L,i])+((Subscript[F, R,i])-(Subscript[F, L,i]))ξ,{i,(d^2)}]

Also define a parametrization for variation of material properties

MaterialPathRules={𝒜->AL+(AR-AL)ξ,𝒟->DL+(DR-DL)ξ,μ->muL+(muR-muL)ξ}

Hℰpath[ξ_] =Table[Expand[ Hℰ[[i,j]] /. FpathRules  /. MaterialPathRules],{i,2d},{j,2d}];

Long expressions result:

Manipulate[Hℰpath[ξ][[i,j]],{i,1,(d^2),1},{j,1,(d^2),1}]

Evaluate integral of Hessian elements under parametrization F(ξ)=(SuperMinus[F])+((SuperPlus[F])-(SuperMinus[F]))ξ, analogous for material properties, and store the result.

If[FileExistsQ["Hav.mx"],Print["Reading Hav.mx ..."]; <<Hav.mx; Print["Done."],Hav=Monitor[ Table[Simplify[Integrate[Hℰpath[ξ][[i,j]],{ξ,0,1}] ],{i,2d},{j,2d}],{i,j}];DumpSave["Hav.mx",Hav]];

Manipulate[Hav[[i,j]],{i,1,(d^2),1},{j,1,(d^2),1}]

Investigate  zero-deformation state

ZeroDeformRulesLR=Flatten[{Table[{(Subscript[F, L,1]^k)->1,(Subscript[F, L,2]^k)->0,(Subscript[F, L,3]^k)->0,(Subscript[F, L,4]^k)->1,(Subscript[F, R,1]^k)->1,(Subscript[F, R,2]^k)->0,(Subscript[F, R,3]^k)->0,(Subscript[F, R,4]^k)->1},{k,8}] }];

Hav0 =Monitor[Table[Simplify[Hav[[i,j]]/. ZeroDeformRulesLR],{i,(d^2)},{j,(d^2)}],{i,j}]

SameMaterialRules = {AL->𝒜,DL->𝒟,muL->μ,AR->𝒜,DR->𝒟,muR->μ};

Hav0 /. SameMaterialRules

M=Table[H[[k]].U[[k]],{k,d}];
Map[TableForm[#]&,M]

Prules = Flatten[Table[(Subscript[P, i,j])->Hav[[i,j]],{i,(d^2)},{j,(d^2)}]];

MP = Monitor[Table[M[[k,i,j]] /. Prules,{k,1,d,1},{i,1,d,1},{j,1,d,1}],{k,i,j}];

Manipulate[MP[[k,i,j]] ,{k,1,d,1},{i,1,d,1},{j,1,d,1}]

MP0  = Simplify[ MP /. ZeroDeformRulesLR]

Manipulate[Simplify[M[[k,i,j]] /. Prules /. ZeroDeformRulesLR],{k,1,d,1},{i,1,d,1},{j,1,d,1}]

Generate Fortran code

Define additional, material-dependent replacement rules

Frules =Flatten[{ Table[(Subscript[F, L,i]^k)->FL[i,k],{i,(d^2)},{k,8}],Table[(Subscript[F, R,i]^k)->FR[i,k],{i,(d^2)},{k,8}]}]

Modify StressGradient Fortran template to take into account specific material properties

dPdF[]:=Module[{i,j,k,n,coef,instr,subr},subr="SUBROUTINE StressGradient(k,MatParams,F,M)\n";
subr=subr<>"  IMPLICIT NONE\n";
subr=subr<>"  INTEGER,PARAMETER::MaxParams=200,L=1,R=2\n";
subr=subr<>"  INTEGER k,n\n";
subr=subr<>"  REAL (KIND=8) MatParams(MaxParams,2),F(4,2),M(2,2),FL(4,8),FR(4,8)\n";
subr=subr<>"  INTEGER i,j\n";
subr=subr<>"  REAL (KIND=8) AL,AR,DL,DR,muL,muR\n";
subr=subr<>"  AL=MatParams(1,1); DL=MatParams(2,1); muL=MatParams(3,1)\n";
subr=subr<>"  AR=MatParams(1,2); DR=MatParams(2,2); muR=MatParams(3,2)\n";
subr=subr<>"  FL(1:4,1) = F(1:4,1)\n";
subr=subr<>"  FR(1:4,1) = F(1:4,2)\n";
subr=subr<>"  DO n=2,8\n";
subr=subr<>"    FL(1:4,n) = FL(1:4,n-1)*FL(1:4,1)\n";
subr=subr<>"    FR(1:4,n) = FR(1:4,n-1)*FR(1:4,1)\n";
subr=subr<>"  END DO\n";
subr=subr<>"  SELECT CASE(k)\n";
For[k=1,k≤d,k++,
subr=subr<>"    CASE(" <> ToString[k] <>")\n";
For[i=1,i≤d,i++,For[j=i,j≤d,j++,coef=ToString[FortranForm[M[[k,i,j]]/.Prules/.Frules]];
instr="      M("<>ToString[i]<>","<>ToString[j]<>") = "<>coef<>"\n";
subr=subr<>instr;];
For[j=1,j<i,j++,
subr=subr<> "      M("<>ToString[i]<>","<>ToString[j]<>") = M(" <>ToString[j]<>","<>ToString[i]<>")\n"];
];
];
subr=subr<>" END SELECT\n";
subr=subr<>"END SUBROUTINE StressGradient\n";
Export["StressGradient.f90",subr,"Text"];
Return[subr]]

dPdF[];