Rectify the interface traction term

Change the sign in front of the pressure in the interface traction term and use the normal on the solid (as described in the tutorial documentation) instead of the normal on the fluid for the computation of the interface traction contribution to the system matrix. Previously, the pillar was bending in the wrong direction due to the wrong sign.

diff --git a/doc/news/changes/minor/20190515SebastianStark b/doc/news/changes/minor/20190515SebastianStark
index 71a7dd9e9b524fb96518639ac59976e11939535d..cbfbabc6a7d58a704c488ece58df451c1de56d7c 100644 (file)
--- a/doc/news/changes/minor/20190515SebastianStark
+++ b/doc/news/changes/minor/20190515SebastianStark
@@ -1,3 +1,3 @@
-Fixed: Fix a bug in tutorial step-46, which caused the y-displacement (z-displacement in 3d) to be slightly unsymmetric.
+Fixed: Fix a bug in tutorial step-46, which caused the y-displacement (z-displacement in 3d) to be slightly unsymmetric and rectify the interface traction term.
 <br>
 (Sebastian Stark, 2019/05/15)

diff --git a/examples/step-46/doc/intro.dox b/examples/step-46/doc/intro.dox
index 77305ac75a3d8317eefeb4dae7354c1699c8ad4a..a0fa0fdac7e63b9c4119e4362bc5156b35bd5543 100644 (file)
--- a/examples/step-46/doc/intro.dox
+++ b/examples/step-46/doc/intro.dox
@@ -36,7 +36,7 @@ where you may want to read up on the individual equations):
   while we assume free-flow conditions on the remainder of the external
   boundary,
   @f{align*}
-    (2\eta \varepsilon(\mathbf v) + p \mathbf 1) \cdot \mathbf n = 0
+    (2\eta \varepsilon(\mathbf v) - p \mathbf 1) \cdot \mathbf n = 0
      \qquad\qquad
      \text{on}\ \Gamma_{f,2} = \partial\Omega \cap \partial\Omega_f \backslash
      \Gamma_{f,1}.
@@ -69,9 +69,11 @@ where you may want to read up on the individual equations):
   Secondly, the forces (traction) on the solid equal the normal stress from the fluid,
   @f{align*}
     (C \varepsilon(\mathbf u)) \mathbf n =
-    (2 \eta \varepsilon(\mathbf v) + p \mathbf 1) \mathbf n \qquad\qquad
-     \text{on}\ \Gamma_{i} = \partial\Omega_s \cap \partial\Omega_f.
+    (2 \eta \varepsilon(\mathbf v) - p \mathbf 1) \mathbf n \qquad\qquad
+     \text{on}\ \Gamma_{i} = \partial\Omega_s \cap \partial\Omega_f,
   @f}
+  where $\mathbf{n}$ is the normal vector on $\Gamma_{i}$ pointing from
+  the solid to the fluid.
 
 We get a weak formulation of this problem by following our usual rule of
 multiplying from the left by a test function and integrating over the
@@ -86,7 +88,7 @@ H^1(\Omega_s)^d$ such that
        + (\varepsilon(\mathbf b), C \varepsilon(\mathbf u))_{\Omega_s} &
        \\
        - (\mathbf b,
-           (2 \eta \varepsilon(\mathbf v) + p \mathbf 1) \mathbf n)_{\Gamma_i}
+           (2 \eta \varepsilon(\mathbf v) - p \mathbf 1) \mathbf n)_{\Gamma_i}
        &=
        0,
 @f}
@@ -308,7 +310,7 @@ discrete level we recall to be
        + (\varepsilon(\mathbf b_h), C \varepsilon(\mathbf u_h))_{\Omega_s} &
        \\
        - (\mathbf b_h,
-           (2 \eta \varepsilon(\mathbf v_h) + p \mathbf 1) \mathbf n)_{\Gamma_i}
+           (2 \eta \varepsilon(\mathbf v_h) - p \mathbf 1) \mathbf n)_{\Gamma_i}
        &=
        0,
 @f}
@@ -323,8 +325,8 @@ in the @ref vector_valued module (rather than the one from step-8).
 
 The term that is of more interest is the interface term,
 @f[
-       (\mathbf b_h,
-           (2 \eta \varepsilon(\mathbf v_h) + p \mathbf 1) \mathbf n)_{\Gamma_i}.
+       -(\mathbf b_h,
+           (2 \eta \varepsilon(\mathbf v_h) - p \mathbf 1) \mathbf n)_{\Gamma_i}.
 @f]
 Based on our assumption that the interface $\Gamma_i$ coincides with
 cell boundaries, this can in fact be written as a set of face
@@ -334,8 +336,8 @@ notation $\psi_i[\mathbf v],\psi_i[p], \psi_i[\mathbf u]$, then the
 term above yields the following contribution to the global matrix
 entry $i,j$:
 @f[
-       \sum_K (\psi_i[\mathbf u],
-           (2 \eta \varepsilon(\psi_j[\mathbf v]) + \psi_j[p] \mathbf 1)
+       -\sum_K (\psi_i[\mathbf u],
+           (2 \eta \varepsilon(\psi_j[\mathbf v]) - \psi_j[p] \mathbf 1)
           \mathbf n)_{\partial K \cap \Gamma_i}.
 @f]
 Although it isn't immediately obvious, this term presents a slight
@@ -558,7 +560,7 @@ prescribe Dirichlet conditions for the flow at the top so that we get
 inflow at the left and outflow at the right. At the left and right
 boundaries, no boundary conditions are imposed explicitly for the
 flow, yielding the implicit no-stress condition $(2\eta
-\varepsilon(\mathbf v) + p \mathbf 1) \cdot \mathbf n = 0$.
+\varepsilon(\mathbf v) - p \mathbf 1) \cdot \mathbf n = 0$.
 The conditions on the interface between the two domains has already been
 discussed above.
 

diff --git a/examples/step-46/doc/results.dox b/examples/step-46/doc/results.dox
index 6c20e128d163a5c8f103c433e1fbcca3c1da2c49..02ef00fd494ae681350a53d8596c567af8324dca 100644 (file)
--- a/examples/step-46/doc/results.dox
+++ b/examples/step-46/doc/results.dox
@@ -21,29 +21,29 @@ Refinement cycle 1
    Writing output...
 
 Refinement cycle 2
-   Number of active cells: 412
-   Number of degrees of freedom: 3667
+   Number of active cells: 436
+   Number of degrees of freedom: 3723
    Assembling...
    Solving...
    Writing output...
 
 Refinement cycle 3
-   Number of active cells: 1216
-   Number of degrees of freedom: 9999
+   Number of active cells: 1096
+   Number of degrees of freedom: 7737
    Assembling...
    Solving...
    Writing output...
 
 Refinement cycle 4
-   Number of active cells: 2788
-   Number of degrees of freedom: 18537
+   Number of active cells: 2656
+   Number of degrees of freedom: 15177
    Assembling...
    Solving...
    Writing output...
 
 Refinement cycle 5
-   Number of active cells: 6496
-   Number of degrees of freedom: 35985
+   Number of active cells: 5992
+   Number of degrees of freedom: 29061
    Assembling...
    Solving...
    Writing output...

diff --git a/examples/step-46/step-46.cc b/examples/step-46/step-46.cc
index 12fa63c0c44a1182fe699cf4a43b20f44d0f99d3..8f14d0355f967f038696cdb26bda79f7bd453a0d 100644 (file)
--- a/examples/step-46/step-46.cc
+++ b/examples/step-46/step-46.cc
@@ -482,20 +482,20 @@ namespace Step46
     FEFaceValues<dim>    stokes_fe_face_values(stokes_fe,
                                             common_face_quadrature,
                                             update_JxW_values |
-                                              update_normal_vectors |
                                               update_gradients | update_values);
     FEFaceValues<dim>    elasticity_fe_face_values(elasticity_fe,
                                                 common_face_quadrature,
-                                                update_values);
+                                                update_normal_vectors |
+                                                  update_values);
     FESubfaceValues<dim> stokes_fe_subface_values(stokes_fe,
                                                   common_face_quadrature,
                                                   update_JxW_values |
-                                                    update_normal_vectors |
                                                     update_gradients |
                                                     update_values);
     FESubfaceValues<dim> elasticity_fe_subface_values(elasticity_fe,
                                                       common_face_quadrature,
-                                                      update_values);
+                                                      update_normal_vectors |
+                                                        update_values);
 
     // ...to objects that are needed to describe the local contributions to
     // the global linear system...
@@ -796,7 +796,7 @@ namespace Step46
     for (unsigned int q = 0; q < n_face_quadrature_points; ++q)
       {
         const Tensor<1, dim> normal_vector =
-          stokes_fe_face_values.normal_vector(q);
+          elasticity_fe_face_values.normal_vector(q);
 
         for (unsigned int k = 0; k < stokes_fe_face_values.dofs_per_cell; ++k)
           {
@@ -813,7 +813,7 @@ namespace Step46
              ++i)
           for (unsigned int j = 0; j < stokes_fe_face_values.dofs_per_cell; ++j)
             local_interface_matrix(i, j) +=
-              -((2 * viscosity * (stokes_symgrad_phi_u[j] * normal_vector) +
+              -((2 * viscosity * (stokes_symgrad_phi_u[j] * normal_vector) -
                  stokes_phi_p[j] * normal_vector) *
                 elasticity_phi[i] * stokes_fe_face_values.JxW(q));
       }