elasticmatrixproblem.hh 9.37 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/**
 * \file
 * \brief Definition of a problem, for the linear elasticity problem:
 * Problem definition for the deformation of an elastic solid.
 */
#ifndef DUMUX_ELASTICMATRIXPROBLEM_HH
#define DUMUX_ELASTICMATRIXPROBLEM_HH

#include <dune/grid/io/file/dgfparser/dgfyasp.hh>

#include <dumux/geomechanics/elastic/elasticmodel.hh>
#include <dumux/implicit/common/implicitporousmediaproblem.hh>

#include "elasticspatialparams.hh"

namespace Dumux
{

template <class TypeTag>
class ElasticMatrixProblem;

namespace Properties
{
NEW_TYPE_TAG(ElasticMatrixProblem, INHERITS_FROM(BoxElastic,ElSpatialParams));

// Set the grid type
45
SET_TYPE_PROP(ElasticMatrixProblem, Grid, Dune::YaspGrid<3>);
46
47

// Set the problem property
48
SET_TYPE_PROP(ElasticMatrixProblem, Problem, Dumux::ElasticMatrixProblem<TypeTag>); 
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

}

/*!
 * \ingroup ElasticBoxProblems
 * \ingroup ImplicitTestProblems
 *
 * \brief Problem definition for the deformation of an elastic matrix.
 *
 * The problem defined here leads to the following linear distribution of the
 * solid displacement: u(x,y,z) = 1/E * (x,0,-nu*z) which for the given grid
 * The numerical results can be verified analytically.
 *
 * The 3D domain given in linearelastic.dgf spans from (0,0,0) to (10,1,10).
 *
 * Dirichlet boundary conditions (u=0.0) are applied for the displacement in y-direction
 * in the whole domain, for the displacement in x-direction on the left boundary (x < eps)
 * and for the displacement in z-direction for the lower left edge (x<eps && z<eps).
 * On the remaining boundaries Neumann boundary conditions are applied.
 * The normal stress applied on each boundary results from solving the momentum balance
 * analytically for the solid displacement function u(x,y,z) = 1/E * (x,0,-nu*z).
 * This leads to the normal stresses: \f$ \boldsymbol{\sigma_{xx}} = 2\,\frac{\mu}{E} + \frac{\lambda}{E} ( 1-\nu)\f$,
 *                                       \f$ \boldsymbol{\sigma_{yy}} = \frac{\lambda}{E} ( 1-\nu)\f$,
 *                                       \f$ \boldsymbol{\sigma_{zz}} = -2\,\frac{\mu \nu}{E} + \frac{\lambda}{E}\f$.
 * The shear stresses are set to zero.
 *
75
 * This problem uses the \ref ElasticModel model.
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
 *
 * To run the simulation execute the following line in shell:
 * <tt>./test_elastic -parameterFile ./test_elastic.input</tt>
 */

template <class TypeTag>
class ElasticMatrixProblem: public ImplicitPorousMediaProblem<TypeTag>
{
    typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;

    // copy some indices for convenience
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
    enum {
        // Grid and world dimension
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld,
    };
    enum {
        // indices of the primary variables
            uxIdx = Indices::uxIdx,
            uyIdx = Indices::uyIdx,
            uzIdx = Indices::uzIdx,
    };
    enum {
        // indices of the equations
            momentumXEqIdx = Indices::momentumXEqIdx,
            momentumYEqIdx = Indices::momentumYEqIdx,
            momentumZEqIdx = Indices::momentumZEqIdx,
    };

    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<dim>::Entity Vertex;
    typedef typename GridView::Intersection Intersection;

    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;

public:
    ElasticMatrixProblem(TimeManager &timeManager, const GridView &gridView)
        : ParentType(timeManager, gridView)
    {}

    /*!
     * \name Problem parameters
     */
    // \{

    /*!
     * \brief The problem name.
     *
     * This is used as a prefix for files generated by the simulation.
     */
    const char *name() const
    {    return "elasticmatrix";}
    // \}

    /*!
     * \name Boundary conditions
     */
    // \{

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     *
     * \param values The boundary types for the conservation equations
     * \param vertex The vertex for which the boundary type is set
     */
149
150
    void boundaryTypesAtPos(BoundaryTypes &values,
                            const GlobalPosition &globalPos) const
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
    {
        values.setAllNeumann();
        values.setDirichlet(uyIdx, momentumYEqIdx);

        if(globalPos[0] < eps_)
        {
            values.setDirichlet(uxIdx, momentumXEqIdx);
            if(globalPos[2] < eps_)
                values.setDirichlet(uzIdx, momentumZEqIdx);
        }
    }

    /*!
     * \brief Evaluate the boundary conditions for a Dirichlet
     *        boundary segment.
     *
     * \param values The Dirichlet values for the primary variables
     * \param vertex The vertex for which the boundary type is set
     *
     * For this method, the \a values parameter stores primary variables.
     */
    void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
    {
        values = 0.0;
    }

    /*!
     * \brief Evaluate the boundary conditions for a Neumann
     *        boundary segment.
     *
     * For this method, the \a values parameter stores the mass flux
     * in normal direction of each phase. Negative values mean influx.
     */
    void neumann(PrimaryVariables &values,
            const Element &element,
186
187
            const FVElementGeometry &fvGeometry,
            const Intersection &intersection,
188
189
190
191
192
193
            int scvIdx,
            int boundaryFaceIdx) const
    {
        values = 0.0;
        // get Lame parameters

194
195
196
197
        Scalar lambda = this->spatialParams().lameParams(element, fvGeometry, scvIdx)[0];
        Scalar mu = this->spatialParams().lameParams(element, fvGeometry, scvIdx)[1];
        Scalar E = this->spatialParams().E(element, fvGeometry, scvIdx);
        Scalar nu = this->spatialParams().nu(element, fvGeometry, scvIdx);
198
199
200
201
202
203
204
205
206
207
208

        // calculate values of sigma in normal direction
        Dune::FieldMatrix<Scalar, dim, dim> sigma(0);
        sigma[0][0] = 2.0*mu + lambda*(1 - nu);
        sigma[1][1] = lambda*(1 - nu);
        sigma[2][2] = -2.0*mu*nu + lambda*(1 - nu);

        sigma *= -1.0/E;

        // determine normal vector of current face
        Dune::FieldVector<Scalar, dim-1> localDimM1(0);
209
        Dune::FieldVector<Scalar,dim> normal = intersection.unitOuterNormal(localDimM1);
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265

        // use stress in normal direction as boundary condition
        sigma.mv(normal, values);
    }
    // \}

    /*!
     * \name Volume terms
     */
    // \{

    /*!
     * \brief Evaluate the source term for all phases within a given
     *        sub-control-volume.
     *
     * For this method, the \a priVars parameter stores the rate momentum
     * is generated or annihilate per volume
     * unit. Positive values mean that momentum is created, negative ones
     * mean that it vanishes.
     */
    void sourceAtPos(PrimaryVariables &priVars,
                     const GlobalPosition &globalPos) const
    {
        priVars = Scalar(0.0);
    }

    /*!
     * \brief Evaluate the initial value for a control volume.
     *
     * \param values The initial values for the primary variables
     * \param globalPos The position for which the initial condition should be evaluated
     *
     * For this method, the \a values parameter stores primary
     * variables.
     */
    void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
    {
        initial_(values, globalPos);
    }

    // \}


private:
    // the internal method for the initial condition
    void initial_(PrimaryVariables &priVars,
                  const GlobalPosition &globalPos) const
    {
        priVars = 0.0; // initial condition for the solid displacement
    }

    static constexpr Scalar eps_ = 3e-6;
};
} //end namespace

#endif