On Tue, 22 Jan 2008, David Knezevic wrote:
> I'd like to implement the following boundary condition for the outflow
> of a channel:
> du/dn - n*p = 0,
> where n is the outward normal and p is the pressure. It seems to me that
> this BC would result in an extra boundary term $-\int_{\partial\Omega} n
> p \ds$ being added to the momentum equations?
That seems right.
> In my case the outward normal on the outflow boundary is n = (1,0). I
> tried to implement this BC using the following code in side_constraint:
>
> if (boundary_id == outflow_id)
> Fu(i) -= JxW_side[qp] * p * phi_side[i][qp];
>
> if (boundary_id != outflow_id)
> Fv(i) += JxW_side[qp] * penalty *
> (v - v_value) * phi_side[i][qp];
>
> if (boundary_id != outflow_id)
> Kuu(i,j) += JxW_side[qp] * penalty *
> phi_side[i][qp] * phi_side[j][qp];
> if (boundary_id != outflow_id)
> Kvv(i,j) += JxW_side[qp] * penalty *
> phi_side[i][qp] * phi_side[j][qp];
>
> where p = side_value(p_var, qp), plus there are other BCs not shown
> here for inflow and no-slip. This doesn't appear to work; the Newton
> convergence fails with these BCs. If anyone can point out how to
> properly implement this BC, I'd be most appreciative :)
For the new Fu term, you'll need a corresponding Jacobian term. Try:
if (boundary_id == outflow_id)
Kup(i,j) -= JxW_side[qp] * psi_side[j][qp] * phi_side[i][qp];
---
Roy
-------------------------------------------------------------------------
This SF.net email is sponsored by: Microsoft
Defy all challenges. Microsoft(R) Visual Studio 2008.
http://clk.atdmt.com/MRT/go/vse0120000070mrt/direct/01/
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
