Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Model involving SOD of two variables

  1. Sep 1, 2011 #1
    1. The problem statement, all variables and given/known data
    I'm going to model the velocity of a river stream v(y,z) as a function of the location (y,z) in the cross-section of the river. Here y represents the distance from one of the shores and z the depth of the river.
    I have done measurements of the depth z on different locations y, so I can find an approximation for a function z(y) and hence find the profile of the river.
    I have also done measurements of the velocity on different points in the river, so now I wish to find a function v(y,z) that describes the speed at every point in the cross-section. This will give me a way to find the average speed by integrating over the function for the cross-area.

    2. Relevant equations
    Basically, through a number of assumptions and approximations, the problem can be reduced to a set of two differential equations:

    [itex]\eta_{y}\frac{\partial^{2} v(y,z)}{\partial y^2} + \eta_{z}\frac{\partial^{2} v(y,z)}{\partial z^2} = -\rho g Sin(\alpha)[/itex]

    By assuming by turns [itex]\eta_{y} = 0[/itex] and [itex]\eta_{z} = 0[/itex] I can find functions v(y) and v(z) that describe a kind of an average velocity for any given y or given z.
    The problem now is the integration of these two into one function v(y,z).

    Given boundary conditions:
    * v(0,z) = 0
    * v(W,z) = 0 -> W represents max width
    * v(y,0) = 0 -> (y,0) is a point on the bottom
    * [itex]\frac{\partial v}{\partial z}[H] = 0 [/itex] -> change of speed is zero in for the maximum height H.

    3. The attempt at a solution
    I've found solved the equations for the assumptions [itex]\eta_{y} = 0[/itex] and [itex]\eta_{z} = 0[/itex] and using the boundary conditions found:

    [itex]v(y) = \frac{\rho g Sin[\alpha]}{2\eta_{y}} (W-y)y[/itex]

    [itex]v(z) = \frac{\rho g Sin[\alpha]}{2\eta_{z}} (2H-z)z[/itex]

    The question now is whether there is an easier way to find v(y,z) by directly solving it mathematically, possibly by combining solutions for v(y) and v(z)? I can see that a linear combination of the two would work, but I can't see how I can find the general solution of the function.

    I have tried to find a solution using Mathematica, but by filling in the boundary conditions it reduced to the function v(y) that is independent of z. Here is what my code in Mathematica more or less looks like:

    Code (Text):

    DSolve[{a*D[f[y, z], {y, 2}] + b *D[f[y, z], {z, 2}] == -\[Rho]*g*
        Sin[\[Alpha]]}, f, {y, z}]

    {{f -> Function[{y, z}, -((g y^2 \[Rho] Sin[\[Alpha]])/(2 a)) +
         C[1][(Sqrt[-a b] y)/a + z] + C[2][-((Sqrt[-a b] y)/a) + z]]}}

    v[y_, z_] := -((g y^2 \[Rho] Sin[\[Alpha]])/(2 a)) +
       c1 ((Sqrt[-a b] y)/a + z) + c2 (-((Sqrt[-a b] y)/a) + z);
    r1 := v[0, z];
    r2 := v[W, z];
    r3 := v[y, 0];
    r4 := D[v[y, z], z];
    {r1, r2, r3, r4}

    {c1 z + c2 z,
     c2 (-((Sqrt[-a b] W)/a) + z) + c1 ((Sqrt[-a b] W)/a + z) - (
      g W^2 \[Rho] Sin[\[Alpha]])/(2 a), (Sqrt[-a b] c1 y)/a - (Sqrt[-a b] c2 y)/
      a - (g y^2 \[Rho] Sin[\[Alpha]])/(2 a), c1 + c2}

    Clear[c1, c2]
    Solve[{r1 == 0, r2 == 0}, {c1, c2}]
    Solve[{r3 == 0, r4 == 0}, {c1, c2}]
    Solve[{r1 == 0, r3 == 0}, {c1, c2}]
    Solve[{r2 == 0, r3 == 0}, {c1, c2}]

    {{c1 -> (g W \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]),
      c2 -> -((g W \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]))}}

    {{c1 -> (g y \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]),
      c2 -> -((g y \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]))}}

    {{c1 -> (g y \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]),
      c2 -> -((g y \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]))}}

    {{c1 -> -((-b g W^2 \[Rho] Sin[\[Alpha]] + b g W y \[Rho] Sin[\[Alpha]] +
         Sqrt[-a b] g y z \[Rho] Sin[\[Alpha]])/(4 a b z)),
      c2 -> -((-Sqrt[-a b] g W^2 \[Rho] Sin[\[Alpha]] +
         Sqrt[-a b] g W y \[Rho] Sin[\[Alpha]] + a g y z \[Rho] Sin[\[Alpha]])/(
        4 a Sqrt[-a b] z))}}

    Clear[c1, c2]
    c1 := (g W \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]);
    c2 := -((g W \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]));
    Simplify[v[y, z]]

    (g (W - y) y \[Rho] Sin[\[Alpha]])/(2 a)

    c1 := -((-b g W^2 \[Rho] Sin[\[Alpha]] + b g W y \[Rho] Sin[\[Alpha]] +
        Sqrt[-a b] g y z \[Rho] Sin[\[Alpha]])/(4 a b z));
    c2 := -((-Sqrt[-a b] g W^2 \[Rho] Sin[\[Alpha]] +
        Sqrt[-a b] g W y \[Rho] Sin[\[Alpha]] + a g y z \[Rho] Sin[\[Alpha]])/(
       4 a Sqrt[-a b] z));
    Simplify[v[y, z]]

    (g W (W - y) \[Rho] Sin[\[Alpha]])/(2 a)
    In both cases I obtain exactly the solution for v(y).
  2. jcsd
  3. Sep 1, 2011 #2
    Have not gone through it in detail but the first thing I notice is that you seem to be mis-interpreting the line:

    [tex]\left\{\left\{f\to \text{Function}\left[\{y,z\},-\frac{g y^2 \rho \text{Sin}[\alpha ]}{2 a}+C[1]\left[\frac{\sqrt{-a b} y}{a}+z\right]+C[2]\left[-\frac{\sqrt{-a b} y}{a}+z\right]\right]\right\}\right\}[/tex]

    The C[1] and C[2] are function names and not constants c1 and c2 multiplied by those values in brackets. You realize that or no? So the first one is an arbritrary function of a single variable named C1 at the value given in brackets. Or am I insulting your intelligence by saying that? Sorry if so.
  4. Sep 2, 2011 #3
    Ah, thanks, I actually didn't know that! Indeed, now the solution seems to make more sense to me. It's identical to the wave function except that it is non-homogeneous. Now I presume I can find the functions C[1] and C[2] in the same way as the wave function and then add the factor [itex] -\frac{gy^2\rho Sin[\alpha]}{2\eta_{y}} [/itex] to it?

    EDIT: I have simply no idea how to solve for the functions C[1] and C[2] using the boundary conditions. I looked up the solution to the wave equation but the approach there is quite different since it is possible to use the boundary conditions of amplitude and speed at t=0 for x=0. I lack knowledge of differential equations to solve for the functions. Also, if I try to insert the boundary conditions in DSolve in Mathematica, it doesn't give me any solution.
    Last edited: Sep 2, 2011
  5. Sep 2, 2011 #4
    Yeah, me neither. How about solving it numerically? First consider a simplified version:

    [tex]\frac{\partial^2 v}{\partial y^2}+\frac{\partial^2 v}{\partial z^2}=-1/2,\quad 0\leq y\leq W,\quad 0\leq z\leq H[/tex]

    [tex]\begin{array}{cc} v(0,z)=0 & v(W,z)=0 \\
    v(y,0)=0 & \frac{\partial v}{\partial z}\biggr|_{z=H}=0

    Or something simple like that. So I'll set up a grid say for now 10x10, do all the difference equations on it with a spacing of one, then solve for all the internal points so that's what 64 linear equations in 64 variables. Alright, I did this quick so not entirely confident of the code but the back-substitution on the results using a sample point in the grid are close to -1/2 which is the first numeric result and the second result is the derivative with respect to z along the top edge at a select point. It's -0.07. Not zero but pretty close to what we want. So then maybe the code is workable so if you want, you can start adding complexity to the PDE and make the grid spacing smaller. When you're done with that, it's not then hard to integrate the interpolation function as you indicated above.

    Code (Text):

    maxv = 10;

    myary = Array[Subscript[v, #1, #2] & , {maxv, maxv}];

    For[j = 1, j <= maxv, j++, {Subscript[v, 1, j] = 0; Subscript[v, j, 1] = 0;
         Subscript[v, j, maxv] = 0; Subscript[v, maxv, j] = Subscript[v, maxv - 2, j]; }];

    myeqns = Flatten[Table[Subscript[v, n, j + 1] + Subscript[v, n, j - 1] +
           Subscript[v, n - 1, j] + Subscript[v, n + 1, j] - 4*Subscript[v, n, j] == -0.5,
         {n, 2, maxv - 1}, {j, 2, maxv - 1}]];

    myvars = Flatten[Table[Subscript[v, n, j], {n, 2, maxv - 1}, {j, 2, maxv - 1}]];

    mysol = First[Solve[myeqns, myvars]];

    For[n = 2, n <= maxv - 1, n++, For[j = 2, j <= maxv - 1, j++,
         Subscript[v, n, j] = Subscript[v, n, j] /. mysol; ]; ];

    mydatapoints = Flatten[Table[{j, n, myary[[n,j]]}, {n, 1, maxv}, {j, 1, maxv}], 1];

    they = 2.6;
    thez = 3.8;
    thev = Interpolation[mydatapoints];
    vy2[y_, z_] = D[thev[y, z], {y, 2}]
    vz2[y_, z_] = D[thev[y, z], {z, 2}]
    N[vy2[they, thez] + vz2[they, thez]]

    vz1[y_, z_] = D[thev[y, z], z]
    vz1[6.5, 10]

    InterpolatingFunction[][y, z]

    InterpolatingFunction[][y, z]


    InterpolatingFunction[][y, z]


    I can do a plot of the solution via:


    which is below. Note at the top edge with respect to z, it's almost flat which I assume is the first partial set to zero there.

    Attached Files:

    Last edited: Sep 2, 2011
  6. Sep 2, 2011 #5
    Thanks a lot. The graph indeed shows the expected result, with the velocity highest midway between the shores and decreasing in the neighbourhood of the shores and bottom.
    I'm going to complete it by adding extending the formula and making the grid larger. Also, since I have real data for a number of points (around 20 I believe), I could insert those in my grid and then try to find an Interpolation through all the points. I'll show some results if I have some more.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook