diff --git a/input b/input new file mode 100644 index 0000000000..5ed36c65f5 --- /dev/null +++ b/input @@ -0,0 +1,69 @@ +alamo.program = mechanics +alamo.program.mechanics.model = finite.adhesion + +plot_file = output + +# this is not a time integration, so do +# exactly one timestep and then quit +timestep = 0.01 +stop_time = 2.0 + +# amr parameters +amr.plot_int = 1 +amr.max_level = 6 +amr.n_cell = 32 32 32 +amr.blocking_factor = 8 +amr.regrid_int = 1 +amr.grid_eff = 1.0 +amr.cell.all = 1 +amr.thermo.plot_int = 1 + +# geometry +geometry.prob_lo = -8 -8 -8 +geometry.prob_hi = 8 8 8 +geometry.is_periodic = 0 0 0 + +#ic.type = ellipse +#ic.ellipse.a = 1.0 0.75 0.5 # ellipse radii +#ic.ellipse.x0 = 0 0 0 # location of ellipse center +#ic.ellipse.eps = 0.1 # diffuse boundary + +ic.type = expression +ic.expression.constant.eps = 0.3 +ic.expression.constant.R = 4.0 +ic.expression.region0 = " (max( erf((x^2 + y^2 - R)/eps),0.0))^2" +ic.expression.region1 = " (max(-erf((x^2 + y^2 - R)/eps),0.0))^2" +ic.expression.region2 = "1.0 - (max( erf((x^2 + y^2 - R)/eps),0.0))^2 - (max(-erf((x^2 + y^2 - R)/eps),0.0))^2" +#ic.expression.region2 = "0.0" + +# elastic moduli +nmodels = 2 +model1.mu = 2.0 +model1.kappa = 3.0 +model2.mu = 10.0 +model2.kappa = 30.0 +model3.mu = 2.0 +model3.kappa = 3.0 +model1.d = 0.0 +model2.d = 0.0 +model3.d = 0.0 +model1.zeta = 0.0 +model2.zeta = 0.0 +model3.zeta = 1.0 +model1.L = 0.0 +model2.L = 0.0 +model3.L = 0.2 +#model3.mu = 210 +#model3.kappa = 0.3 + +solver.verbose = 3 +solver.nriters = 10 +#solver.max_iter = 20 +solver.fixed_iter = 100 +#elasticop.small = 0.1 + +print_model = 1 + +bc.type = tensiontest +bc.tensiontest.type = uniaxial_stress +bc.tensiontest.disp=(0,1,2:0,8,0) diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H new file mode 100644 index 0000000000..e8c7788642 --- /dev/null +++ b/src/Model/Solid/Finite/Adhesion.H @@ -0,0 +1,242 @@ +// This models a finite-strain, solid with adhesion behavior governed by +// a combination of elastic energy and an adhesion decay term. +// +// The energy and derivatives are: +// +// .. math:: +// :nowrap: +// +// \begin{gather} +// W = d \left[ \frac{\mu}{2J^{2/3}}(F^T F : I - 3) + \frac{\kappa}{2}(J - 1)^2 \right] + \frac{\zeta}{J^n} \\ +// DW = \text{(first Piola-Kirchhoff stress)} \\ +// DDW = \text{(tangent stiffness tensor)} +// \end{gather} +// +// where: +// - :math:`d` The coefficient of degradable part, +// - :math:`\mu` is the shear modulus, +// - :math:`\kappa` is the bulk modulus, +// - :math:`\zeta` is the adhesion decay coefficient or the coefficient of non-degradable part, +// - :math:`n` is the decay exponent, +// - :math:`F` is the deformation gradient, +// - :math:`J = \det(F)` is the volume ratio, +// - :math:`F^T F : I = \text{trace}(F^T F)` is used to capture isochoric deformation. +// +// Notes: +// +// Adhesion energy model W(F): +// Consists of two parts: +// 1. Degradable (elastic shear + volumetric terms) +// 2. Non-degradable (adhesion term decaying with volume expansion) +// This model is rate-independent, reversible, and purely elastic. + +#ifndef MODEL_SOLID_FINITE_ADHESION_H_ +#define MODEL_SOLID_FINITE_ADHESION_H_ + +#include "IO/ParmParse.H" +#include "Model/Solid/Finite/NeoHookean.H" +#include "Model/Solid/Solid.H" +#include "Set/Set.H" +#include // For std::array to represent fourth-order tensors +#include // For std::pow and std::abs +#include // For Eigen::Matrix3d + +namespace Model +{ +namespace Solid +{ +namespace Finite +{ + +class Adhesion : public NeoHookean +{ +public: + static constexpr KinematicVariable kinvar = KinematicVariable::F; + + Adhesion() {} + Adhesion(Solid base) : NeoHookean(base) {} + virtual ~Adhesion() {} + + // Material parameters (declared in order) + Set::Scalar d = NAN; ///< Adhesion strength (degradable part coefficient) + Set::Scalar zeta = NAN; ///< Adhesion decay parameter (non-degradable part) + Set::Scalar n = NAN; ///< Decay exponent + Set::Scalar L = NAN; ///< Mobility + + // Return a zero-parameter adhesion model. + static Adhesion + Zero() + { + Adhesion adhesion; + adhesion.d = 0.0; ///< Adhesion strength (degradable part coefficient) + adhesion.zeta = 0.0; ///< Adhesion decay parameter (non-degradable part) + adhesion.n = 0.0; ///< Decay exponent + adhesion.L = 0.0; + adhesion.mu = 0.0; + adhesion.kappa = 0.0; + return adhesion; + } + + // Return a random-parameter adhesion model. + static Adhesion + Random() + { + Adhesion adhesion; + adhesion.d = Util::Random() ; // Util::Random(); ///< Adhesion strength (degradable part coefficient) + adhesion.zeta = Util::Random(); ///< Adhesion decay parameter (non-degradable part) + adhesion.n = Util::Random(); ///< Decay exponent + adhesion.mu = Util::Random(); + adhesion.kappa = Util::Random();// Util::Random(); + adhesion.L = Util::Random(); + return adhesion; + } + + // Parse function: Reads material parameters from the input file. + // Here, we assume the input file provides: mu, kappa, d, zeta, and n. + static void + Parse(Adhesion &value, IO::ParmParse &pp) + { + pp.queryclass(value); + // Initial value for damage + pp_query_default("d", value.d, 0.0); + // Adhesion decay parameter (non-degradable part) + pp_query_default("zeta", value.zeta, 0.0); + // Decay exponent + pp_query_default("n", value.n, 1.0); + // Mobility + pp_query_default("L", value.L, 0.0); + } + + // Strain energy function: + // W = d * [ 0.5*mu*(trace(F^T F)-3)/J^(2/3) + 0.5*kappa*(J-1)^2 ] + zeta/J^n + Set::Scalar + W(const Set::Matrix &F) const override + { + Set::Scalar w_nh = (1.0 - d) * NeoHookean::W(F); + Set::Scalar J = F.determinant(); + Set::Scalar w_compress = zeta * std::pow(J, -n); + return w_nh + w_compress; + } + + // First Piola-Kirchhoff stress (derivative of W) + Set::Matrix + DW(const Set::Matrix &F) const override + { + Set::Matrix dw_nh = (1.0 - d) * NeoHookean::DW(F); + Set::Scalar J = F.determinant(); + Set::Matrix FinvT = F.inverse().transpose(); + Set::Matrix dw_compress = zeta * (-n * std::pow(J, -n) * FinvT); + return dw_nh + dw_compress; + } + + // Tangent stiffness tensor (second derivative of W) + Set::Matrix4 + DDW(const Set::Matrix &F) const override + { + Set::Matrix4 ddw_nh = + (1.0 - d) * NeoHookean::DDW(F); + + Set::Scalar J = F.determinant(); + Set::Matrix FinvT = F.inverse().transpose(); + + Set::Matrix4 ddw_compress; + + for (int i = 0; i < AMREX_SPACEDIM; i++) + { + for (int j = 0; j < AMREX_SPACEDIM; j++) + { + for (int k = 0; k < AMREX_SPACEDIM; k++) + { + for (int l = 0; l < AMREX_SPACEDIM; l++) + { + // Adhesion-related term (zeta) + + // dFinvT(i,j) / dF(k,l) + // = dFinv(j,i) / dF(k,l) + // = -dFinvT(k,j) * dFinvT(iL) + + Set::Scalar t3 = n * n * std::pow(J, -n) * FinvT(i, j) * FinvT(k, l) + + n * std::pow(J, -n) * FinvT(k, j) * FinvT(i, l); + ddw_compress(i, j, k, l) = zeta * t3; + } + } + } + } + return ddw_nh + ddw_compress; + } + + // Print material parameters + virtual void + Print(std::ostream &out) const override + { + out << "d = " << d << ", mu = " << mu << ", kappa = " << kappa + << ", zeta = " << zeta << ", n = " << n; + } + + virtual void Advance(Set::Scalar dt, Set::Matrix strain, + Set::Matrix P, Set::Scalar /*time*/) override + { + // Util::Message(INFO,dt); + // Util::Message(INFO,L); + // Util::Message(INFO,(P.transpose()*P).trace()); + Set::Scalar dd = std::max(0.0,(P.transpose()*P).trace() * (0.8 - d)); + // Util::Message(INFO,dd); + d += dd * L * dt; + } + + + +#define OP_CLASS Adhesion +#define OP_VARS X(d) X(mu) X(kappa) X(zeta) X(n) X(L) +#include "Model/Solid/InClassOperators.H" +}; + +#include "Model/Solid/ExtClassOperators.H" + +} // namespace Finite +} // namespace Solid +} // namespace Model + +template<> +ALAMO_SINGLE_DEFINITION +int Set::Field::NComp() const +{ + return 4; +} +template<> +ALAMO_SINGLE_DEFINITION +std::string Set::Field::Name(int i) const +{ + if (i == 0) return name + "_mu"; + if (i == 1) return name + "_kapp"; + if (i == 2) return name + "_d"; + if (i == 3) return name + "_L"; + return name; +} +template<> +ALAMO_SINGLE_DEFINITION +void Set::Field::Copy(int a_lev, amrex::MultiFab& a_dst, int a_dstcomp, int a_nghost) const +{ + for (amrex::MFIter mfi(a_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost)); + if (bx.ok()) + { + amrex::Array4 const& src = ((*this)[a_lev])->array(mfi); + amrex::Array4 const& dst = a_dst.array(mfi); + for (int n = 0; n < AMREX_SPACEDIM * AMREX_SPACEDIM; n++) + { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + dst(i, j, k, a_dstcomp + 0) = src(i, j, k).mu; + dst(i, j, k, a_dstcomp + 1) = src(i, j, k).kappa; + dst(i, j, k, a_dstcomp + 2) = src(i, j, k).d; + dst(i, j, k, a_dstcomp + 3) = src(i, j, k).L; + }); + } + } + } +} + + + +#endif // MODEL_SOLID_FINITE_ADHESION_H_ diff --git a/src/mechanics.cc b/src/mechanics.cc index 7638d6f07b..8542b8b262 100644 --- a/src/mechanics.cc +++ b/src/mechanics.cc @@ -12,6 +12,7 @@ #include "Model/Solid/Affine/Isotropic.H" #include "Model/Solid/Linear/Cubic.H" #include "Model/Solid/Affine/Cubic.H" +#include "Model/Solid/Finite/Adhesion.H" #include "Model/Solid/Finite/NeoHookean.H" #include "Model/Solid/Finite/NeoHookeanPredeformed.H" #include "Model/Solid/Finite/PseudoLinear/Cubic.H" @@ -56,6 +57,8 @@ int main (int argc, char* argv[]) pp.select_only>(integrator); else if (model == "finite.neohookeanpre") pp.select_only>(integrator); + else if (model == "finite.adhesion") + pp.select_only>(integrator); else if (model == "finite.pseudolinear.cubic") pp.select_only>(integrator); else if (model == "finite.pseudoaffine.cubic") diff --git a/src/test.cc b/src/test.cc index c5deb8f3df..ca05cf3902 100644 --- a/src/test.cc +++ b/src/test.cc @@ -15,6 +15,7 @@ #include "Model/Solid/Linear/Laplacian.H" #include "Model/Solid/Affine/Isotropic.H" #include "Model/Solid/Affine/Cubic.H" +#include "Model/Solid/Finite/Adhesion.H" #include "Model/Solid/Finite/NeoHookean.H" #include "Model/Solid/Finite/NeoHookeanPredeformed.H" #include "Model/Solid/Finite/PseudoLinear/Cubic.H" @@ -57,7 +58,8 @@ int main (int argc, char* argv[]) MODELTEST(Model::Solid::Finite::PseudoLinear::Cubic); MODELTEST(Model::Solid::Finite::NeoHookeanPredeformed); MODELTEST(Model::Solid::Finite::PseudoAffine::Cubic); - + MODELTEST(Model::Solid::Finite::Adhesion); + Util::Test::Message("Set::Matrix4"); { diff --git a/tests/Solid/input b/tests/Solid/input index fac0d54c40..969d55bdee 100644 --- a/tests/Solid/input +++ b/tests/Solid/input @@ -105,6 +105,25 @@ #@ ignore = model1.E model1.nu #@ coverage=true #@ +#@ [adhesion] +#@ exe=mechanics +#@ dim=3 +#@ check = true +#@ check-file = reference/adhesion.dat +#@ args = timestep=0.01 +#@ args = alamo.program.mechanics.model=finite.adhesion +#@ args = model1.mu=3.0 +#@ args = model1.kappa=6.5 +#@ args = model1.d=1.0 +#@ args = model1.zeta=0.5 +#@ args = model1.n=2.0 +#@ args = bc.tensiontest.disp=(0,1:0,1) +#@ args = solver.nriters=10 +#@ ignore = model1.E model1.nu +#@ benchmark-beaker = 10.73 +#@ benchmark-statler = 8.17 +#@ benchmark-github = 11.35 +#@ alamo.program = mechanics diff --git a/tests/Solid/reference/adhesion.dat b/tests/Solid/reference/adhesion.dat new file mode 100644 index 0000000000..6dc6910767 --- /dev/null +++ b/tests/Solid/reference/adhesion.dat @@ -0,0 +1,101 @@ +time disp_xhi_x disp_xhi_y disp_yhi_x disp_yhi_y trac_xhi_x trac_xhi_y trac_yhi_x trac_yhi_y +0 0 -0.123824 0 0 1.15602 -4.01326e-11 0 0 +0.01 0.01 -0.12456 0 0 1.2205 -2.86443e-10 0 0 +0.02 0.02 -0.125279 0 0 1.28357 -3.93683e-10 0 0 +0.03 0.03 -0.125982 0 0 1.34527 -3.4169e-10 0 0 +0.04 0.04 -0.126672 0 0 1.40568 -3.48891e-10 0 0 +0.05 0.05 -0.127348 0 0 1.46483 -3.40385e-10 0 0 +0.06 0.06 -0.128012 0 0 1.52278 -2.15224e-10 0 0 +0.07 0.07 -0.128666 0 0 1.57958 -2.15809e-10 0 0 +0.08 0.08 -0.129309 0 0 1.63528 -2.85579e-10 0 0 +0.09 0.09 -0.129944 0 0 1.68991 -4.1995e-10 0 0 +0.1 0.1 -0.130571 0 0 1.74353 -1.97542e-12 0 0 +0.11 0.11 -0.13119 0 0 1.79616 -7.2879e-12 0 0 +0.12 0.12 -0.131804 0 0 1.84785 -4.12464e-11 0 0 +0.13 0.13 -0.132411 0 0 1.89863 -1.15129e-13 0 0 +0.14 0.14 -0.133013 0 0 1.94854 -2.45439e-10 0 0 +0.15 0.15 -0.133611 0 0 1.99761 -2.35392e-10 0 0 +0.16 0.16 -0.134205 0 0 2.04587 -4.76978e-11 0 0 +0.17 0.17 -0.134795 0 0 2.09335 1.49388e-10 0 0 +0.18 0.18 -0.135382 0 0 2.14007 -4.91033e-12 0 0 +0.19 0.19 -0.135967 0 0 2.18607 -6.66379e-10 0 0 +0.2 0.2 -0.136549 0 0 2.23137 -8.69189e-11 0 0 +0.21 0.21 -0.137129 0 0 2.27599 4.86789e-10 0 0 +0.22 0.22 -0.137708 0 0 2.31996 -5.14368e-10 0 0 +0.23 0.23 -0.138285 0 0 2.3633 -1.30087e-12 0 0 +0.24 0.24 -0.138861 0 0 2.40602 -1.76863e-10 0 0 +0.25 0.25 -0.139436 0 0 2.44816 -3.11804e-10 0 0 +0.26 0.26 -0.14001 0 0 2.48972 -3.65531e-10 0 0 +0.27 0.27 -0.140583 0 0 2.53073 -2.95396e-10 0 0 +0.28 0.28 -0.141156 0 0 2.5712 3.44428e-11 0 0 +0.29 0.29 -0.141729 0 0 2.61115 -3.22309e-11 0 0 +0.3 0.3 -0.142301 0 0 2.6506 -1.88363e-10 0 0 +0.31 0.31 -0.142873 0 0 2.68956 -2.5663e-10 0 0 +0.32 0.32 -0.143445 0 0 2.72804 -2.13702e-12 0 0 +0.33 0.33 -0.144017 0 0 2.76606 2.24091e-10 0 0 +0.34 0.34 -0.144588 0 0 2.80363 -2.44375e-10 0 0 +0.35 0.35 -0.14516 0 0 2.84077 1.04304e-11 0 0 +0.36 0.36 -0.145731 0 0 2.87748 -2.56168e-10 0 0 +0.37 0.37 -0.146303 0 0 2.91379 -2.05313e-10 0 0 +0.38 0.38 -0.146874 0 0 2.94969 -2.77199e-10 0 0 +0.39 0.39 -0.147446 0 0 2.9852 7.41536e-12 0 0 +0.4 0.4 -0.148017 0 0 3.02033 -8.64127e-12 0 0 +0.41 0.41 -0.148588 0 0 3.0551 -3.00506e-10 0 0 +0.42 0.42 -0.149159 0 0 3.0895 4.60683e-11 0 0 +0.43 0.43 -0.14973 0 0 3.12355 -3.36724e-10 0 0 +0.44 0.44 -0.150301 0 0 3.15726 -5.79441e-11 0 0 +0.45 0.45 -0.150872 0 0 3.19063 -2.86736e-10 0 0 +0.46 0.46 -0.151442 0 0 3.22368 2.43842e-10 0 0 +0.47 0.47 -0.152012 0 0 3.25642 2.86899e-12 0 0 +0.48 0.48 -0.152582 0 0 3.28884 2.16865e-10 0 0 +0.49 0.49 -0.153151 0 0 3.32096 3.02515e-12 0 0 +0.5 0.5 -0.15372 0 0 3.35278 1.91218e-13 0 0 +0.51 0.51 -0.154288 0 0 3.38431 2.53284e-11 0 0 +0.52 0.52 -0.154856 0 0 3.41556 -2.02127e-12 0 0 +0.53 0.53 -0.155423 0 0 3.44654 2.64056e-11 0 0 +0.54 0.54 -0.15599 0 0 3.47724 -2.6267e-10 0 0 +0.55 0.55 -0.156556 0 0 3.50768 1.40424e-10 0 0 +0.56 0.56 -0.157122 0 0 3.53787 -1.07387e-10 0 0 +0.57 0.57 -0.157686 0 0 3.5678 -1.00735e-10 0 0 +0.58 0.58 -0.15825 0 0 3.59748 -1.29752e-10 0 0 +0.59 0.59 -0.158813 0 0 3.62692 -8.31481e-12 0 0 +0.6 0.6 -0.159375 0 0 3.65612 -9.11869e-12 0 0 +0.61 0.61 -0.159936 0 0 3.68509 -9.81089e-12 0 0 +0.62 0.62 -0.160497 0 0 3.71384 -1.25632e-11 0 0 +0.63 0.63 -0.161056 0 0 3.74236 -1.38536e-11 0 0 +0.64 0.64 -0.161614 0 0 3.77066 -1.15191e-11 0 0 +0.65 0.65 -0.162171 0 0 3.79875 -1.08418e-11 0 0 +0.66 0.66 -0.162727 0 0 3.82663 -1.18489e-11 0 0 +0.67 0.67 -0.163282 0 0 3.8543 -8.40617e-12 0 0 +0.68 0.68 -0.163835 0 0 3.88177 -1.031e-11 0 0 +0.69 0.69 -0.164388 0 0 3.90904 8.89774e-12 0 0 +0.7 0.7 -0.164939 0 0 3.93612 -1.94457e-11 0 0 +0.71 0.71 -0.165489 0 0 3.96301 -1.72192e-11 0 0 +0.72 0.72 -0.166037 0 0 3.98971 3.39921e-11 0 0 +0.73 0.73 -0.166584 0 0 4.01623 1.45924e-11 0 0 +0.74 0.74 -0.16713 0 0 4.04257 -2.67873e-12 0 0 +0.75 0.75 -0.167675 0 0 4.06873 -3.08766e-10 0 0 +0.76 0.76 -0.168217 0 0 4.09471 2.23796e-10 0 0 +0.77 0.77 -0.168759 0 0 4.12053 3.09863e-11 0 0 +0.78 0.78 -0.169299 0 0 4.14618 1.1446e-11 0 0 +0.79 0.79 -0.169837 0 0 4.17166 -2.76966e-10 0 0 +0.8 0.8 -0.170374 0 0 4.19698 3.90644e-11 0 0 +0.81 0.81 -0.170909 0 0 4.22215 -2.30649e-13 0 0 +0.82 0.82 -0.171443 0 0 4.24715 1.16604e-10 0 0 +0.83 0.83 -0.171975 0 0 4.27201 9.44755e-11 0 0 +0.84 0.84 -0.172506 0 0 4.29671 -4.09357e-13 0 0 +0.85 0.85 -0.173035 0 0 4.32126 -3.04696e-12 0 0 +0.86 0.86 -0.173562 0 0 4.34567 1.19393e-10 0 0 +0.87 0.87 -0.174088 0 0 4.36993 1.39926e-10 0 0 +0.88 0.88 -0.174611 0 0 4.39406 1.81844e-10 0 0 +0.89 0.89 -0.175134 0 0 4.41804 -6.50169e-11 0 0 +0.9 0.9 -0.175654 0 0 4.44188 3.89623e-11 0 0 +0.91 0.91 -0.176173 0 0 4.4656 -1.05935e-11 0 0 +0.92 0.92 -0.17669 0 0 4.48917 -1.63238e-11 0 0 +0.93 0.93 -0.177205 0 0 4.51262 -2.55951e-12 0 0 +0.94 0.94 -0.177719 0 0 4.53594 7.32922e-12 0 0 +0.95 0.95 -0.17823 0 0 4.55914 -1.55842e-11 0 0 +0.96 0.96 -0.17874 0 0 4.58221 -8.00096e-12 0 0 +0.97 0.97 -0.179248 0 0 4.60515 -1.05624e-11 0 0 +0.98 0.98 -0.179755 0 0 4.62798 1.43395e-10 0 0 +0.99 0.99 -0.180259 0 0 4.65069 -9.58419e-12 0 0