3
$\begingroup$

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}}}];
$\endgroup$
  • 3
    $\begingroup$ you should give Cijkl if you want folks to be able to try this. $\endgroup$ – george2079 May 7 '18 at 15:04
  • $\begingroup$ They are just numbers. You can find them in attached files. The whole code is too long to copy it into the post. The problem does not realy depend on which numbers they are. $\endgroup$ – Данил Семёнов May 7 '18 at 15:48
  • 1
    $\begingroup$ The thing is that it's better if you put dummy values by yourself in the original question rather than asking every user who is trying to help you to choose dummy values themselves. In short: help people help you :) $\endgroup$ – anderstood May 7 '18 at 20:05
  • $\begingroup$ Is it better to type the whole huge code in the question body or just attach a fite to press a 1 button to execute? =) $\endgroup$ – Данил Семёнов May 8 '18 at 9:08
  • $\begingroup$ UPD: ok, the whole code has been added, as you asked $\endgroup$ – Данил Семёнов May 8 '18 at 10:31
4
$\begingroup$

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
   };
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.