What's wrong with initial conditions?

Question

I want to solve the nonlinear PDE for the anisotropic fluid flow:

$-\rho*\partial{v_i}/\partial{t} + (C_{ijkl}v_{k,l})_{,j}-p_{,j} = 0$

(The nonlinear term $-v_jv_{i,j}$ will be added later)

Here is the part of the code describing PDE (differential operator, boundary and initial conditions):

Ui[t, x, y, z] = {u[t, x, y, z], v[t, x, y, z], w[t, x, y, z]};
coord = {x, y, z};

(*Operator*)
op = {
   -rho*D[Ui[t, x, y, z][[1]], t]+Sum[D[Sum[Sum[Cijkl[[1, j, k, l]]*D[Ui[t, x, y, z][[k]], coord[[l]]], {l,3}], {k, 3}], coord[[j]]], {j, 3}] - D[p[t, x, y, z], x],
   -rho*D[Ui[t, x, y, z][[2]], t]+Sum[D[Sum[Sum[Cijkl[[2, j, k, l]]*D[Ui[t, x, y, z][[k]], coord[[l]]], {l,3}], {k, 3}], coord[[j]]], {j, 3}] - D[p[t, x, y, z], y],
   -rho*D[Ui[t, x, y, z][[3]], t]+Sum[D[Sum[Sum[Cijkl[[3, j, k, l]]*D[Ui[t, x, y, z][[k]], coord[[l]]], {l,3}], {k, 3}], coord[[j]]], {j, 3}] - D[p[t, x, y, z], z],
   Sum[D[Ui[t, x, y, z][[k]], coord[[k]]], {k, 3}]
   };

(*Boundary conditions*)
bcs = {
   DirichletCondition[p[t, x, y, z] == 100, x == 0],
   DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0, w[t, x, y, z] == 0}, 0 < x < 8],
   DirichletCondition[p[t, x, y, z] == 0, x == 8]
   };

(*Initial conditions*)
ics = {
   u[0, x, y, z] == 0,
   v[0, x, y, z] == 0,
   w[0, x, y, z] == 0
   };

(*Solve*)
{xVel, yVel, zVel, pressure} = NDSolveValue[{op == {0, 0, 0, 0}, ics, bcs}, {u, v, w, p}, {t, 0, 1}, {x, y, z} \[Element] mesh, Method -> {"MethodOfLines","SpatialDiscretization" -> {"FiniteElement","InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, p -> 1}}}];

When I try to solve it, I get the following error message:

NDSolveValue::ivone: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable.

I cant understand what's wrong with the conditions? May there be a syntax error somewhere?

I suppose that boundary conditions are ok, because the stationary version of this equation $(C_{ijkl}v_{k,l})_{,j}-p_{,j} = 0$ is well solved with the following code:

Ui[x, y, z] = {u[x, y, z], v[x, y, z], w[x, y, z]};
coord = {x, y, z};

op = {
   Sum[D[Sum[Sum[Cijkl[[1, j, k, l]]*D[Ui[x, y, z][[k]], coord[[l]]], {l,3}], {k, 3}], coord[[j]]], {j, 3}] - D[p[x, y, z], x](*-Sum[Ui[x,y,z][[j]]*D[Ui[x,y,z][[1]],coord[[j]]], {j,3}]*),
   Sum[D[Sum[Sum[Cijkl[[2, j, k, l]]*D[Ui[x, y, z][[k]], coord[[l]]], {l,3}], {k, 3}], coord[[j]]], {j, 3}] - D[p[x, y, z], y](*-Sum[Ui[x,y,z][[j]]*D[Ui[x,y,z][[2]],coord[[j]]], {j,3}]*),
   Sum[D[Sum[Sum[Cijkl[[3, j, k, l]]*D[Ui[x, y, z][[k]], coord[[l]]], {l,3}], {k, 3}], coord[[j]]], {j, 3}] - D[p[x, y, z], z](*-Sum[Ui[x,y,z][[j]]*D[Ui[x,y,z][[3]],coord[[j]]], {j,3}]*),
   Sum[D[Ui[x, y, z][[k]], coord[[k]]], {k, 3}]
   };

bcs = {
   DirichletCondition[p[x, y, z] == 100, x == 0],
   DirichletCondition[{u[x, y, z] == 0, v[x, y, z] == 0, w[x, y, z] == 0}, 0 < x < 8],
   DirichletCondition[p[x, y, z] == 0, x == 8]
   };

{xVel, yVel, zVel, pressure} = NDSolveValue[{op == {0,0,0,0}, bcs}, {u,v,w,p}, {x,y,z} \[Element] mesh, Method->{"FiniteElement","InterpolationOrder"->{u->2,v->2,w-> 2,p->1}}];

So, is there any mistake in the ICs or solve section?

P.S.: there are 2 versions of whole file for your inspection and execution:

1) Steady flow (it works, evaluate whole notebook): https://1drv.ms/u/s!AvmYXV0MDC0Jh4ZFYTSCBZ9vspyjKw

2) Unsteady flow (doesn't work): https://1drv.ms/u/s!AvmYXV0MDC0Jh4ZGjmCorU_DymHq3A

UPD: ok, there is the whole code with the constants and mesh definition:

Needs["NDSolve`FEM`"]

a2 = { {1, 0, 0}, {0, 0, 0}, {0, 0, 0} };

a4 = Table[a2[[i,j]]*a2[[k,l]],{i,3},{j,3},{k,3},{l,3}];
Np = 0;
Ns = 0;
mu = 1;
rho = 1;

Cijklt = Table[2mu*(KroneckerDelta[i,k]*KroneckerDelta[j,l]+Np*a4[[i,j,k,l]]+Ns*(a2[[l,j]]*KroneckerDelta[i,k]+a2[[i,k]]*KroneckerDelta[j,l])), {i,3},{j,3},{k,3},{l,3}];
Delta = Table[0.5*(KroneckerDelta[i,k]*KroneckerDelta[j,l]+KroneckerDelta[i,l]*KroneckerDelta[j,k]),{i,3},{j,3},{k,3},{l,3}];

Cijkl = Table[Sum[Sum[Cijklt[[i,j,m,n]]*Delta[[m,n,k,l]],{n,3}],{m,3}],{i,3},{j,3},{k,3},{l,3}];

omega = RegionUnion[Cylinder[{{0,0,0},{4,0,0}},2],Cylinder[{{4,0,0},{8,0,0}},1]];
mesh = ToElementMesh[omega, MaxCellMeasure-> 0.005];
RegionPlot3D[omega, PlotPoints->50];
mesh["Wireframe"];

Ui[t,x,y,z] = {u[t,x,y,z], v[t,x,y,z],w[t,x,y,z]};
coord = {x,y,z};

op = {
-rho*D[Ui[t,x,y,z][[1]], t]+Sum[D[Sum[Sum[Cijkl[[1,j,k,l]]*D[Ui[t,x,y,z][[k]],coord[[l]]],{l,3}],{k,3}], coord[[j]]],{j,3}]-D[p[t,x,y,z],x](*-Sum[Ui[t,x,y,z][[j]]*D[Ui[t,x,y,z][[1]], coord[[j]]], {j,3}]*),
-rho*D[Ui[t,x,y,z][[2]], t]+Sum[D[Sum[Sum[Cijkl[[2,j,k,l]]*D[Ui[t,x,y,z][[k]],coord[[l]]],{l,3}],{k,3}], coord[[j]]],{j,3}]-D[p[t,x,y,z],y](*-Sum[Ui[t,x,y,z][[j]]*D[Ui[t,x,y,z][[2]], coord[[j]]], {j,3}]*),
-rho*D[Ui[t,x,y,z][[3]], t]+Sum[D[Sum[Sum[Cijkl[[3,j,k,l]]*D[Ui[t,x,y,z][[k]],coord[[l]]],{l,3}],{k,3}], coord[[j]]],{j,3}]-D[p[t,x,y,z],z](*-Sum[Ui[t,x,y,z][[j]]*D[Ui[t,x,y,z][[3]], coord[[j]]], {j,3}]*),
Sum[D[Ui[t,x,y,z][[k]],coord[[k]]],{k,3}]
};

bcs={
DirichletCondition[p[t,x,y,z]==100,x== 0 ],
DirichletCondition[{u[t,x,y,z]==0,v[t,x,y,z]==0,w[t,x,y,z]==0},0<x<8],
DirichletCondition[p[t,x,y,z]==0,x== 8]
};

ics = {
u[0,x,y,z]==0,
v[0,x,y,z]==0,
w[0,x,y,z]==0
};

{xVel,yVel, zVel,pressure}=NDSolveValue[{op=={0,0,0,0}, ics,bcs},{u,v,w,p},{t, 0,1},{x,y,z}\[Element]mesh,Method->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement","InterpolationOrder"->{u->2,v->2,w-> 2,p->1}}}];

Answer 1

Ok, here is the answer.

It is a little bit unexpected, but NDSolveValue requires initial conditions for all unknown functions of PDE. Even those they have no time derivative.

There is the pressure function p[t,x,y,z] in equations. Usualy, the initial conditions for Navier-Stokes-like equations are the velocity field but Mathematica requires the pressure initial condition too.

For example, there is the working version of initial conditions that allows you to avoid error message:

ics = {
   u[0, x, y, z] == 0,
   v[0, x, y, z] == 0,
   w[0, x, y, z] == 0,
   p[0, x, y, z] == 0
   };