right so in the Stokes subsystem of my equations, in simplified form I am
trying to solve:
div tau - grad p = rhs1
div v = rhs2
where my tau is 2 epsilon as in step-22
The boundary conditions I now want to implement to solve the 'real' problem is:
on the sides (boundary 0): zero tangential stress and no normal flux, so we
have v dot n = 0 and n.(-pI + tau)t = 0.
Here, you would use VectorTools::compute_no_normal_flux_constraints(), which
takes care of the v.n=0 condition. The tangential stress would lead to an
additional term in the right hand side of the weak form, but because it's
zero, there is nothing additional you need to do.
at the top (boundary 1): zero tangential stress and prescribed normal
component of normal stress, so that n.(-pI+tau)n= top stress value,
and similarly, on the bottom (boundary 2): zero tangential stress and
prescribed normal component of normal stress, so that n.(-pI+tau)n= bottom
stress value.
n and t here are the normal and tangential vectors relative to the surface at
the respective boundaries, respectively.
The test code works with Dirichlet conditions around the boundary it seems,
but i need to implement the problem as above.
Since the last message, I realised that I can just expand n.(-pI+tau) into
normal and tangential components, which means that for the above real problem,
I can just implement Neumann conditions with top stress value*normal vector to
the surface, as the tangential stress are zero and therefore in the
expansions, just simply equal zero. ie. I have n.(-pI+tau) =
stress_value*normal vector. This is for the top and bottom boundaries.
Yes, exactly -- normal and tangential stress are just two components of one
vector, and they can be treated together exactly as you suggest.
I am now having a small issue with this, and tried to look up more information
about fe_face_values, but I am a little stuck.
I am doing (eg. for boundary id 1, and i have the same for 2):
for (unsigned int face_number=0;
face_number<GeometryInfo<dim>::faces_per_cell; ++face_number)
if (cell->face(face_number)->at_boundary()
&&
(cell->face(face_number)->boundary_id() == 1))
{
fe_face_values.reinit (cell, face_number);
for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
{
for (unsigned int i=0; i<dofs_per_cell; ++i)
local_rhs(i) +=
data::topstress*(fe_face_values.normal_vector(q_point) *
fe_face_values.shape_value(i,q_point) *
fe_face_values.JxW(q_point));
}
}
I can't seem to find info about the fe_face_values.normal_vector, but it seems
I am struggling with the types (tensor/double/vector) in the calculation - I
am getting the error no match for óperator+= '(operand types are 'double'' and
'dealii::Tensor<1, 2, double>')
Here, local_rhs(i) is just a number (double) whereas the compiler tells you
that the right hand side is a Tensor<1,2>. I suspect that's because you
declare 'topstress' as a Tensor<2,2>? If so, the expression on the right would be
Tensor<2,2>
* Tensor<1,2> // the normal vector
* double // shape_value
* double // JxW
which obviously is a Tensor<1,2>. I think you forget that the top stress
should really be just the normal component of the stress (a "force") because
that's all you can prescribe. In other words, you either need to write this as
normal_vector * topstress * normal_vector * ...
if you want 'topstress' to be a Tensor<2,2>, or as
topforce * normal_vector * ...
if you correctly only describe the *force* on the top boundary.
Best
W.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: bange...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
