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}}}];
Cijkl
if you want folks to be able to try this. $\endgroup$