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