Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <limits>
7 : #include <pup.h>
8 :
9 : #include "DataStructures/CachedTempBuffer.hpp"
10 : #include "DataStructures/DataBox/Prefixes.hpp"
11 : #include "DataStructures/DataVector.hpp"
12 : #include "DataStructures/TaggedTuple.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "Elliptic/Systems/Elasticity/Tags.hpp"
15 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
16 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
17 : #include "Options/String.hpp"
18 : #include "PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.hpp"
19 : #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
20 : #include "Utilities/ContainerHelpers.hpp"
21 : #include "Utilities/Serialization/CharmPupable.hpp"
22 : #include "Utilities/TMPL.hpp"
23 :
24 0 : namespace Elasticity::Solutions {
25 :
26 : namespace detail {
27 : template <typename DataType>
28 : struct BentBeamVariables {
29 : using Cache = CachedTempBuffer<
30 : Tags::Displacement<2>,
31 : ::Tags::deriv<Tags::Displacement<2>, tmpl::size_t<2>, Frame::Inertial>,
32 : Tags::Strain<2>, Tags::MinusStress<2>, Tags::PotentialEnergyDensity<2>,
33 : ::Tags::FixedSource<Tags::Displacement<2>>>;
34 :
35 : const tnsr::I<DataType, 2>& x;
36 : const double length;
37 : const double height;
38 : const double bending_moment;
39 : const ConstitutiveRelations::IsotropicHomogeneous<2>& constitutive_relation;
40 :
41 : void operator()(gsl::not_null<tnsr::I<DataType, 2>*> displacement,
42 : gsl::not_null<Cache*> cache,
43 : Tags::Displacement<2> /*meta*/) const;
44 : void operator()(gsl::not_null<tnsr::iJ<DataType, 2>*> deriv_displacement,
45 : gsl::not_null<Cache*> cache,
46 : ::Tags::deriv<Tags::Displacement<2>, tmpl::size_t<2>,
47 : Frame::Inertial> /*meta*/) const;
48 : void operator()(gsl::not_null<tnsr::ii<DataType, 2>*> strain,
49 : gsl::not_null<Cache*> cache, Tags::Strain<2> /*meta*/) const;
50 : void operator()(gsl::not_null<tnsr::II<DataType, 2>*> minus_stress,
51 : gsl::not_null<Cache*> cache,
52 : Tags::MinusStress<2> /*meta*/) const;
53 : void operator()(gsl::not_null<Scalar<DataType>*> potential_energy_density,
54 : gsl::not_null<Cache*> cache,
55 : Tags::PotentialEnergyDensity<2> /*meta*/) const;
56 : void operator()(
57 : gsl::not_null<tnsr::I<DataType, 2>*> fixed_source_for_displacement,
58 : gsl::not_null<Cache*> cache,
59 : ::Tags::FixedSource<Tags::Displacement<2>> /*meta*/) const;
60 : };
61 : } // namespace detail
62 :
63 : /*!
64 : * \brief A state of pure bending of an elastic beam in 2D
65 : *
66 : * \details This solution describes a 2D slice through an elastic beam of length
67 : * \f$L\f$ and height \f$H\f$, centered around (0, 0), that is subject to a
68 : * bending moment \f$M=\int T^{xx}y\mathrm{d}y\f$ (see e.g.
69 : * \cite ThorneBlandford2017, Eq. 11.41c for a bending moment in 1D). The beam
70 : * material is characterized by an isotropic and homogeneous constitutive
71 : * relation \f$Y^{ijkl}\f$ in the plane-stress approximation (see
72 : * `Elasticity::ConstitutiveRelations::IsotropicHomogeneous`). In this scenario,
73 : * no body-forces \f$f_\mathrm{ext}^j\f$ act on the material, so the
74 : * \ref Elasticity equations reduce to \f$\nabla_i T^{ij}=0\f$, but the bending
75 : * moment \f$M\f$ generates the stress
76 : *
77 : * \f{align}
78 : * T^{xx} &= \frac{12 M}{H^3} y \\
79 : * T^{xy} &= 0 = T^{yy} \text{.}
80 : * \f}
81 : *
82 : * By fixing the rigid-body motions to
83 : *
84 : * \f[
85 : * \xi^x(0,y)=0 \quad \text{and} \quad \xi^y\left(\pm \frac{L}{2},0\right)=0
86 : * \f]
87 : *
88 : * we find that this stress is produced by the displacement field
89 : *
90 : * \f{align}
91 : * \xi^x&=-\frac{12 M}{EH^3}xy \\
92 : * \xi^y&=\frac{6 M}{EH^3}\left(x^2+\nu y^2-\frac{L^2}{4}\right)
93 : * \f}
94 : *
95 : * in terms of the Young's modulus \f$E\f$ and the Poisson ration \f$\nu\f$ of
96 : * the material. The corresponding strain \f$S_{ij}=\partial_{(i}\xi_{j)}\f$ is
97 : *
98 : * \f{align}
99 : * S_{xx} &= -\frac{12 M}{EH^3} y \\
100 : * S_{yy} &= \frac{12 M}{EH^3} \nu y \\
101 : * S_{xy} &= S_{yx} = 0
102 : * \f}
103 : *
104 : * and the potential energy stored in the entire infinitesimal slice is
105 : *
106 : * \f[
107 : * \int_{-L/2}^{L/2} \int_{-H/2}^{H/2} U dy\,dx = \frac{6M^2}{EH^3}L \text{.}
108 : * \f]
109 : */
110 1 : class BentBeam : public elliptic::analytic_data::AnalyticSolution {
111 : public:
112 0 : using constitutive_relation_type =
113 : Elasticity::ConstitutiveRelations::IsotropicHomogeneous<2>;
114 :
115 0 : struct Length {
116 0 : using type = double;
117 0 : static constexpr Options::String help{"The beam length"};
118 0 : static type lower_bound() { return 0.0; }
119 : };
120 0 : struct Height {
121 0 : using type = double;
122 0 : static constexpr Options::String help{"The beam height"};
123 0 : static type lower_bound() { return 0.0; }
124 : };
125 0 : struct BendingMoment {
126 0 : using type = double;
127 0 : static constexpr Options::String help{
128 : "The bending moment applied to the beam"};
129 0 : static type lower_bound() { return 0.0; }
130 : };
131 0 : struct Material {
132 0 : using type = constitutive_relation_type;
133 0 : static constexpr Options::String help{
134 : "The material properties of the beam"};
135 : };
136 :
137 0 : using options = tmpl::list<Length, Height, BendingMoment, Material>;
138 0 : static constexpr Options::String help{
139 : "A 2D slice through an elastic beam which is subject to a bending "
140 : "moment. The bending moment is applied along the length of the beam, "
141 : "i.e. the x-axis, so that the beam's left and right ends are bent "
142 : "towards the positive y-axis. It is measured in units of force."};
143 :
144 0 : BentBeam() = default;
145 0 : BentBeam(const BentBeam&) = default;
146 0 : BentBeam& operator=(const BentBeam&) = default;
147 0 : BentBeam(BentBeam&&) = default;
148 0 : BentBeam& operator=(BentBeam&&) = default;
149 0 : ~BentBeam() override = default;
150 0 : std::unique_ptr<elliptic::analytic_data::AnalyticSolution> get_clone()
151 : const override {
152 : return std::make_unique<BentBeam>(*this);
153 : }
154 :
155 : /// \cond
156 : explicit BentBeam(CkMigrateMessage* m)
157 : : elliptic::analytic_data::AnalyticSolution(m) {}
158 : using PUP::able::register_constructor;
159 : WRAPPED_PUPable_decl_template(BentBeam); // NOLINT
160 : /// \endcond
161 :
162 0 : BentBeam(double length, double height, double bending_moment,
163 : constitutive_relation_type constitutive_relation)
164 : : length_(length),
165 : height_(height),
166 : bending_moment_(bending_moment),
167 : constitutive_relation_(std::move(constitutive_relation)) {}
168 :
169 0 : double length() const { return length_; }
170 0 : double height() const { return height_; }
171 0 : double bending_moment() const { return bending_moment_; }
172 :
173 0 : const constitutive_relation_type& constitutive_relation() const {
174 : return constitutive_relation_;
175 : }
176 :
177 : /// Return potential energy integrated over the whole beam material
178 1 : double potential_energy() const {
179 : return 6. * length_ * square(bending_moment_) /
180 : (cube(height_) * constitutive_relation_.youngs_modulus());
181 : }
182 :
183 : template <typename DataType, typename... RequestedTags>
184 0 : tuples::TaggedTuple<RequestedTags...> variables(
185 : const tnsr::I<DataType, 2>& x,
186 : tmpl::list<RequestedTags...> /*meta*/) const {
187 : using VarsComputer = detail::BentBeamVariables<DataType>;
188 : typename VarsComputer::Cache cache{get_size(*x.begin())};
189 : const VarsComputer computer{x, length_, height_, bending_moment_,
190 : constitutive_relation_};
191 : return {cache.get_var(computer, RequestedTags{})...};
192 : }
193 :
194 : template <typename... RequestedTags>
195 0 : tuples::TaggedTuple<RequestedTags...> variables(
196 : const tnsr::I<DataVector, 2>& x, const Mesh<2>& /*mesh*/,
197 : const InverseJacobian<DataVector, 2, Frame::ElementLogical,
198 : Frame::Inertial>& /*inv_jacobian*/,
199 : tmpl::list<RequestedTags...> meta) const {
200 : return variables(x, meta);
201 : }
202 :
203 : /// NOLINTNEXTLINE(google-runtime-references)
204 1 : void pup(PUP::er& p) override {
205 : elliptic::analytic_data::AnalyticSolution::pup(p);
206 : p | length_;
207 : p | height_;
208 : p | bending_moment_;
209 : p | constitutive_relation_;
210 : }
211 :
212 : private:
213 0 : double length_{std::numeric_limits<double>::signaling_NaN()};
214 0 : double height_{std::numeric_limits<double>::signaling_NaN()};
215 0 : double bending_moment_{std::numeric_limits<double>::signaling_NaN()};
216 0 : constitutive_relation_type constitutive_relation_{};
217 : };
218 :
219 0 : bool operator==(const BentBeam& lhs, const BentBeam& rhs);
220 0 : bool operator!=(const BentBeam& lhs, const BentBeam& rhs);
221 :
222 : } // namespace Elasticity::Solutions
|
