Line data Source code
1 : #ifndef MODEL_SOLID_FINITE_NEOHOOKEAN_H_
2 : #define MODEL_SOLID_FINITE_NEOHOOKEAN_H_
3 :
4 : #include "IO/ParmParse.H"
5 : #include "Model/Solid/Solid.H"
6 : #include "Unit/Unit.H"
7 :
8 : namespace Model
9 : {
10 : namespace Solid
11 : {
12 : namespace Finite
13 : {
14 : class NeoHookean : public Solid<Set::Sym::Major>
15 : {
16 : public:
17 13460 : NeoHookean() = default;
18 : NeoHookean(Solid<Set::Sym::Major> base) : Solid<Set::Sym::Major>(base) {};
19 : ~NeoHookean() = default;
20 :
21 960 : Set::Scalar W(const Set::Matrix& a_F) const
22 : {
23 : #if AMREX_SPACEDIM==2
24 380 : Eigen::Matrix3d F = Eigen::Matrix3d::Identity();
25 380 : F(0, 0) = a_F(0, 0);
26 380 : F(0, 1) = a_F(0, 1);
27 380 : F(1, 0) = a_F(1, 0);
28 380 : F(1, 1) = a_F(1, 1);
29 : #elif AMREX_SPACEDIM==3
30 580 : Eigen::Matrix3d F = a_F;
31 : #endif
32 :
33 960 : Set::Scalar J = F.determinant();
34 960 : Set::Scalar J23 = std::pow(fabs(J), 2. / 3.);
35 960 : Set::Scalar w = 0.0;
36 960 : w += 0.5 * mu * ((F * F.transpose()).trace() / J23 - 3.);
37 960 : w += 0.5 * kappa * (J - 1.0) * (J - 1.0);
38 960 : return w;
39 : }
40 123488 : Set::Matrix DW(const Set::Matrix& a_F) const
41 : {
42 : #if AMREX_SPACEDIM==2
43 120228 : Eigen::Matrix3d F = Eigen::Matrix3d::Identity();
44 120228 : F(0, 0) = a_F(0, 0);
45 120228 : F(0, 1) = a_F(0, 1);
46 120228 : F(1, 0) = a_F(1, 0);
47 120228 : F(1, 1) = a_F(1, 1);
48 : #elif AMREX_SPACEDIM==3
49 3260 : Eigen::Matrix3d F = a_F;
50 : #endif
51 :
52 123488 : Set::Scalar J = F.determinant();
53 123488 : Set::Scalar J23 = std::pow(fabs(J), 2. / 3.);
54 123488 : Eigen::Matrix3d FinvT = F.inverse().transpose();
55 :
56 123488 : Eigen::Matrix3d dw = Eigen::Matrix3d::Zero();
57 :
58 123488 : dw += mu * (F / J23 - (F * F.transpose()).trace() * FinvT / (3. * J23));
59 123488 : dw += kappa * (J - 1) * J * FinvT;
60 :
61 : #if AMREX_SPACEDIM==2
62 120228 : Set::Matrix r_dw;
63 120228 : r_dw(0, 0) = dw(0, 0);
64 120228 : r_dw(0, 1) = dw(0, 1);
65 120228 : r_dw(1, 0) = dw(1, 0);
66 120228 : r_dw(1, 1) = dw(1, 1);
67 240456 : return r_dw;
68 : #elif AMREX_SPACEDIM==3
69 6520 : return dw;
70 : #endif
71 : }
72 71074 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Major> DDW(const Set::Matrix& a_F) const
73 : {
74 : #if AMREX_SPACEDIM==2
75 71054 : Eigen::Matrix3d F = Eigen::Matrix3d::Identity();
76 71054 : F(0, 0) = a_F(0, 0);
77 71054 : F(0, 1) = a_F(0, 1);
78 71054 : F(1, 0) = a_F(1, 0);
79 71054 : F(1, 1) = a_F(1, 1);
80 : #elif AMREX_SPACEDIM==3
81 20 : Eigen::Matrix3d F = a_F;
82 : #endif
83 :
84 71074 : Set::Matrix4<3, Set::Sym::Major> ddw;
85 71074 : Set::Scalar J = F.determinant();
86 71074 : Set::Scalar J23 = std::pow(fabs(J), 2. / 3.);
87 71074 : Eigen::Matrix3d FinvT = F.inverse().transpose();
88 284296 : for (int i = 0; i < 3; i++)
89 852888 : for (int j = 0; j < 3; j++)
90 2558664 : for (int k = 0; k < 3; k++)
91 7675992 : for (int l = 0; l < 3; l++)
92 : {
93 5756994 : ddw(i, j, k, l) = 0.0;
94 :
95 5756994 : Set::Scalar t1 = 0.0, t2 = 0.0;
96 :
97 5756994 : if (i == k && j == l) t1 += 1.0;
98 5756994 : t1 -= (2. / 3.) * F(i, j) * FinvT(k, l);
99 5756994 : t1 -= (2. / 3.) * FinvT(i, j) * F(k, l);
100 5756994 : t1 += (2. / 9.) * (F * F.transpose()).trace() * FinvT(i, j) * FinvT(k, l);
101 5756994 : t1 += (1. / 3.) * (F * F.transpose()).trace() * FinvT(i, l) * FinvT(k, j);
102 :
103 5756994 : t2 += (2. * J - 1.) * FinvT(i, j) * FinvT(k, l);
104 5756994 : t2 += (1. - J) * FinvT(i, l) * FinvT(k, j);
105 :
106 11513988 : ddw(i, j, k, l) = (mu / J23) * t1 + kappa * J * t2;
107 : }
108 : #if AMREX_SPACEDIM==2
109 71054 : Set::Matrix4<2, Set::Sym::Major> r_ddw;
110 213162 : for (int i = 0; i < 2; i++)
111 426324 : for (int j = 0; j < 2; j++)
112 852648 : for (int k = 0; k < 2; k++)
113 1705296 : for (int l = 0; l < 2; l++)
114 2273728 : r_ddw(i, j, k, l) = ddw(i, j, k, l);
115 142108 : return r_ddw;
116 : #elif AMREX_SPACEDIM==3
117 40 : return ddw;
118 : #endif
119 : }
120 : void Print(std::ostream& out) const
121 : {
122 : out << "mu = " << mu << " kappa = " << kappa;
123 : }
124 :
125 : public:
126 : Set::Scalar mu = NAN, kappa = NAN;
127 : static constexpr KinematicVariable kinvar = KinematicVariable::F;
128 :
129 : public:
130 8747 : static NeoHookean Zero()
131 : {
132 8747 : NeoHookean ret;
133 8747 : ret.mu = 0.0;
134 8747 : ret.kappa = 0.0;
135 8747 : return ret;
136 : }
137 124 : static NeoHookean Random()
138 : {
139 124 : NeoHookean ret;
140 124 : ret.mu = Util::Random();
141 124 : ret.kappa = Util::Random();
142 124 : return ret;
143 : }
144 4 : static void Parse(NeoHookean& value, IO::ParmParse& pp)
145 : {
146 :
147 24 : std::pair<std::string, Set::Scalar> moduli[2];
148 :
149 16 : pp.forbid("lame","Use 'lambda' instead for lame constant");
150 16 : pp.forbid("shear","Use 'mu' instead for shear modulus");
151 12 : pp.forbid("bulk","Use 'K' instead for bulk modulus");
152 :
153 : // Lame constant \(\lambda\),
154 : // shear modulus \(\mu\), Young's modulus \(E\), Poisson's ratio \(\nu\),
155 : // bulk modulus \(K\).
156 : // You can currently specify (mu and kappa), (lambda and mu), or (E and nu).
157 : // std::vector{Unit::Pressure, Unit::Pressure, Unit::Pressure, Unit::Less, Unit::Pressure};
158 16 : pp.query_exactly<2>({"lambda","mu","E","nu","kappa"}, moduli, std::vector<Unit>{Unit::Pressure(), Unit::Pressure(), Unit::Pressure(), Unit::Less(), Unit::Pressure()});
159 :
160 :
161 4 : if (moduli[0].first == "mu" && moduli[1].first == "kappa")
162 : {
163 4 : value.mu = moduli[0].second;
164 4 : value.kappa = moduli[1].second;
165 : }
166 0 : else if (moduli[0].first == "lambda" && moduli[1].first == "mu")
167 : {
168 0 : Set::Scalar lambda = moduli[0].second;
169 0 : value.mu = moduli[1].second;
170 0 : value.kappa = lambda + (2.0 * value.mu) / 3.0;
171 : }
172 0 : else if (pp.contains("E") && pp.contains("nu")) {
173 0 : Set::Scalar E = moduli[0].second;
174 0 : Set::Scalar nu = moduli[1].second;
175 0 : value.kappa = E / (3.0 - 6.0 * nu);
176 0 : value.mu = E / (2.0 + 2.0 * nu);
177 : }
178 : else
179 : {
180 0 : Util::Exception(INFO,"Haven't implemented",moduli[0].first," and ",moduli[1].first," yet (sorry!)");
181 : }
182 16 : }
183 :
184 : #define OP_CLASS NeoHookean
185 : #define OP_VARS X(kappa) X(mu)
186 : #include "Model/Solid/InClassOperators.H"
187 : };
188 : #include "Model/Solid/ExtClassOperators.H"
189 :
190 : }
191 : }
192 : }
193 :
194 : #endif
|
