From: Sebastian Stark Date: Fri, 17 May 2019 11:03:19 +0000 (+0200) Subject: Rectify the interface traction term X-Git-Tag: v9.2.0-rc1~797^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=68edb3bb1097260a762068ffb3489fb2106bfd23;p=dealii.git 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 71a7dd9e9b..cbfbabc6a7 100644 --- 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.
(Sebastian Stark, 2019/05/15) diff --git a/examples/step-46/doc/intro.dox b/examples/step-46/doc/intro.dox index 77305ac75a..a0fa0fdac7 100644 --- 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 6c20e128d1..02ef00fd49 100644 --- 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 12fa63c0c4..8f14d0355f 100644 --- a/examples/step-46/step-46.cc +++ b/examples/step-46/step-46.cc @@ -482,20 +482,20 @@ namespace Step46 FEFaceValues stokes_fe_face_values(stokes_fe, common_face_quadrature, update_JxW_values | - update_normal_vectors | update_gradients | update_values); FEFaceValues elasticity_fe_face_values(elasticity_fe, common_face_quadrature, - update_values); + update_normal_vectors | + update_values); FESubfaceValues stokes_fe_subface_values(stokes_fe, common_face_quadrature, update_JxW_values | - update_normal_vectors | update_gradients | update_values); FESubfaceValues 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)); }