[shpalman]

Hello all,

I'm trying to write a solver, loosely based on solidDisplacementFoam, which will work on orthotropic solids. The stiffness constants and the density are specified as fields (and there is also a misfit field f, since what I'm actually trying to do involves coherent interfaces between materials with different lattice constants).

Since the system is anisotropic I can't just calculate Lamé parameters, but since OpenFOAM can't handle the full fourth-rank 3x3 stiffness tensor I can't just write sigma = c : epsilon (where sigma is the stress tensor and epsilon is the strain tensor). However, since it's orthotropic I can work with vectors for the pressure and shear terms for both stress and strain, and two 3x3 tensors with the elastic constants in (for c_ij where i or j are 1, 2 or 3; and for c_kk where k is 4, 5 or 6).

rho * fvm::d2dt2(D) == fvc::div(sigma)

epsilon = 0.5*(fvc::grad(D) + fvc::grad(D)().T());

epsilonPressure.component(0)=epsilon.component(0)-f;
epsilonPressure.component(1)=epsilon.component(4)-f;
epsilonPressure.component(2)=epsilon.component(8)-f;

epsilonShear.component(0)=epsilon.component(3);
epsilonShear.component(1)=epsilon.component(2);
epsilonShear.component(2)=epsilon.component(1);

sigmaPressure = cPressure & epsilonPressure;
sigmaShear = cShear & epsilonShear;

sigma.component(0)=sigmaPressure.component(0);
sigma.component(4)=sigmaPressure.component(1);
sigma.component(8)=sigmaPressure.component(2);

sigma.component(5)=sigmaShear.component(0);
sigma.component(7)=sigmaShear.component(0);
sigma.component(2)=sigmaShear.component(1);
sigma.component(6)=sigmaShear.component(1);
sigma.component(1)=sigmaShear.component(2);
sigma.component(3)=sigmaShear.component(2);

I'm calculating the strain tensor from the displacement, setting up two strain vectors, getting two stress vectors from those via the stiffness tensors, and rebuilding the stress tensor which I'd like to put back on the right-hand side of the equation. But is it even possible to solve a system of equations like this, rather than having the spatial derivatives of D explicitly on the right-hand side?

I appreciate any feedback at all...

Best,

Danny.