Message Boards Message Boards

GROUPS:

Solve analytically 1D transport equation?

Posted 5 months ago
948 Views
|
11 Replies
|
3 Total Likes
|

I'm trying to solve for the analytical solution of 1D transport equation to verify the results of the numerical solution.

eqn = D[u[x, t], t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == 50;
bc1 = u[0, t] == 50;
bc2 = D[u[1, t], x] == 0;
DSolve[{eqn, ic, bc1, bc2}, u[x, t], {x, t}]

However, I obtain the following,

DSolve::deqn: Equation or list of equations expected instead of True in the first argument {(u^(0,1))[x,t]==(u^(2,0))[x,t],u[x,0]==50,u[0,t]==50,True}.

Have I missed any step? I'm looking for the symbolic solution of the PDE with Dirichlet boundary condition at the inlet and Neumann boundary condition at the outlet. Could someone help?

11 Replies

You need to write the no flux BC differently, but your specification returns a trivial result of 50 when corrected.

eqn = D[u[x, t], t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == 50;
bc1 = u[0, t] == 50;
bc2 = (D[u[x, t], x] == 0) /. {x -> 1};
DSolve[{eqn, ic, bc1, bc2}, u[x, t], {x, t}] (* {{u[x, t] -> 50}} *)
Posted 5 months ago

Many thanks for the response.I would like to if it is possible to view the complete expression of the analytical solution. I think the solution that is displayed after making the above-suggested change is the steady-state solution.

I think that you will need to rethink your specification. Right now, your initial condition says that the temperature is 50 everywhere. Your first boundary condition says that the temperature at the boundary is 50. Your second boundary conditions says that the system is insulated at the end. How can the system change from the value of 50 everywhere?

Posted 5 months ago

Solving for analytical solution of 1D transport equation

Posted 5 months ago

Solving for analytical solution of 1D transport equation

Posted 5 months ago

Thanks for the advice. I have changed the left boundary condition

bc1 = u[0, t] == 70;

sol = DSolveValue[{eqn, ic, bc1, bc2}, u[x, t], {x, t}]

gives,

70 + 2 Inactive[Sum][(
   40 E^(-225 \[Pi]^2 t (-1 + 2 K[1])^2)
     Sin[1/2 \[Pi] x (-1 + 2 K[1])])/(\[Pi] - 2 \[Pi] K[1]), {K[1], 
    1, \[Infinity]}]

Since the above solution is an infinite series, I am considering first 300 terms.

asol = sol /. {\[Infinity] -> 300}

Next, I am trying to plot the solution at x = 0.5 for time t= 1:20 using the following command

Plot(sol)

The above syntax isn't complete, I would like to ask for some help in specifying the syntax for plotting the above solution.

Look at the DSolveValue>Applications>Classical Partial Differential Equations for an example how to do this.

eqn = D[u[x, t], t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == 50;
bc1 = u[0, t] == 70;
bc2 = (D[u[x, t], x] == 0) /. {x -> 1};
dsol = DSolveValue[{eqn, ic, bc1, bc2}, 
    u[x, t], {x, t}] /. {K[1] -> m};
asol[x_, t_] = dsol /. {\[Infinity] -> 300} // Activate;
Plot[asol[0.5, t], {t, 1, 20}, PlotRange -> All]

Note that a value of 300 causes a warning about the value being too small.

Posted 5 months ago

Thank you. I made the following changes,

eqn = D[u[x, t], t] == 100*D[u[x, t], {x, 2}];

and

asol[x_, t_] = dsol /. {\[Infinity] -> 30000} // Activate;

Plot[asol[0.5, t], {t, 1, 20}, PlotRange -> All]

However, the plot shows only a straight line at u = 70 (y-axis) from 0 to 20 seconds.I had expected the plot to start from 50 at t = 0 and reach 70 after the initial transient phase.

I'm not sure what's causing the problem. Any suggestions on how to get the right solution?

Hi Natasha,

Your diffusion coefficient is huge and that implies that the system responds very quickly. To see this response, you will need to look at times much closer to 0 and not begin at $t=1$. You probably do not need many terms to capture the response.

eqn = D[u[x, t], t] == 100 D[u[x, t], {x, 2}];
ic = u[x, 0] == 50;
bc1 = u[0, t] == 70;
bc2 = (D[u[x, t], x] == 0) /. {x -> 1};
dsol = DSolveValue[{eqn, ic, bc1, bc2}, 
    u[x, t], {x, t}] /. {K[1] -> m};
asol5[x_, t_] = dsol /. {\[Infinity] -> 5} // Activate;
asol10[x_, t_] = dsol /. {\[Infinity] -> 300} // Activate;
Plot[{asol5[0.5, t], asol10[0.5, t]}, {t, 0, 0.001}, PlotRange -> All]

Plot at smaller time

Posted 5 months ago

Thank you very much. I varied the value diffusion constant and check the response.

I want to check how the analytical solution differs when there is convection along with diffusion.

eqn = D[u[x, t], t] == 100*D[u[x, t], {x, 2}] - 50*D[u[x, t], x] ;

For the same set of initial and boundary conditions, I'm checking for the analytical solution

dsol = DSolveValue[{eqn, ic, bc1, bc2}, 
    u[x, t], {x, t}] /. {K[1] -> m};

The output that I obtain is ,

DSolveValue[{
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, t] == -100 
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, t] + 50 
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"2", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, t], u[x, 0] == 50, u[0, t] == 70, 
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[1, t] == 0}, u[x, t], {x, t}]

The summation terms are missing. Could you please help?

Hi Natasha,

You are probably approaching the limit of DSolve's ability to find an analytical solution (at least directly). You should consider a numerical solution for this problem.

If you are okay with a numerical solution, you can use the code below as a starting point.

tmax = 2/100;
ptmax = tmax/5;
eqn = D[u[x, t], t] == 100*D[u[x, t], {x, 2}] - 50*D[u[x, t], x];
ic = u[x, 0] == 50;
bc1 = u[0, t] == 50 + 20 Tanh[20000 t] ;
bc2 = (D[u[x, t], x] == 0) /. {x -> 1};
ndsol = NDSolveValue[{eqn, ic, bc1, bc2}, 
   u[x, t], {x, 0, 1}, {t, 0, tmax}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 50}}, MaxStepFraction -> 1/1000];
ContourPlot[ndsol, {t, 0.`, ptmax}, {x, 0.`, 1}, 
 PlotPoints -> {200, 200}, MaxRecursion -> 4, PlotTheme -> "Web", 
 ColorFunction -> ColorData["DarkBands"]]

Convection Diffusion Image

Note that I used a Tanh function to provide a rapidly rising "step" function so that the boundary and initial conditions were consistent. If you want to explore the effects of velocity on the solution, you could wrap the above (note that flipped the plot 90 degrees so it would flow left to right) into a Module like so.

condiff[v_] := Module[{tmax, ptmax, eqn, ic, bc1, bc2, ndsol, cp},
  tmax = 2/100;
  ptmax =  2 tmax/5;
  eqn = D[u[x, t], t] == 100*D[u[x, t], {x, 2}] - v*D[u[x, t], x];
  ic = u[x, 0] == 50;
  bc1 = u[0, t] == 50 + 20 Tanh[20000 t] ;
  bc2 = (D[u[x, t], x] == 0) /. {x -> 1};
  ndsol = 
   NDSolveValue[{eqn, ic, bc1, bc2}, u, {x, 0, 1}, {t, 0, tmax}, 
    Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MinPoints" -> 2000}}, MaxStepFraction -> 1/1000];
  cp = ContourPlot[ndsol[x, t], {x, 0.`, 1}, {t, 0.`, ptmax}, 
    PlotRange -> {50, 70}, PlotPoints -> {50, 50}, MaxRecursion -> 4, 
    PlotTheme -> "Web", ColorFunction -> ColorData["DarkBands"], 
    FrameLabel -> Automatic];
  Show[{cp, Graphics[Arrow[{{0, 0.004}, {v/500, 0.004}}]]}]
  ]

Now, you can create an animation of velocities from 0 to 500 (it will take a while!).

imgs = Table[condiff[v], {v, 0, 500, 5}];
ListAnimate[imgs]

Animation of velocity

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract