Re-introduction of viscosity term

Author: Federica Ninno

scripts/ChamberSphere.yaml

@@ -24,7 +24,8 @@ derivatives:
 constants:
   - n
   - volume0
-  - gamma_W1_over_n # k = gamma_W1/n  (identifiable effective stiffness; gamma_W1 = k*n)
+  - gamma_W1
+  - gamma_eta
   - gamma_sigma_max # gamma * sigma_max
   - act
   - act_plus
@@ -57,7 +58,7 @@ helper_functions: |
 
 residuals:
   - stretch(volume) * stress - Pout * CG(volume)
-  - -stress + tau + 4 * (1 - CG(volume) ** (-3)) * gamma_W1_over_n * n + prestress
+  - -stress + tau + 4 * gamma_W1 * (1 - CG(volume) ** (-3)) + prestress + gamma_eta * dCG(volume, dvolume_dt) * (1 + 2 * CG(volume) ** (-6))
   - dtau_dt + act * tau - gamma_sigma_max * act_plus
   - Qin - Qout - dvolume_dt
   - Pin - Pout

src/model/ChamberSphere.cpp

@@ -12,6 +12,10 @@ void ChamberSphere::setup_dofs(DOFHandler& dofhandler) {
 
 void ChamberSphere::update_constant(SparseSystem& system,
                                     std::vector<double>& parameters) {
+  const double gamma_eta = parameters[global_param_ids[ParamId::gamma_eta]];
+  const double volume0 = parameters[global_param_ids[ParamId::volume0]];
+  const double n = parameters[global_param_ids[ParamId::n]];
+  system.E.coeffRef(global_eqn_ids[1], global_var_ids[6]) = 2*gamma_eta/(n*volume0);
   system.E.coeffRef(global_eqn_ids[2], global_var_ids[5]) = 1;
   system.E.coeffRef(global_eqn_ids[3], global_var_ids[6]) = -1;
   system.F.coeffRef(global_eqn_ids[0], global_var_ids[2]) = -1;
@@ -36,19 +40,22 @@ void ChamberSphere::update_solution(
     const Eigen::Matrix<double, Eigen::Dynamic, 1>& y,
     const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy) {
   const double prestress = parameters[global_param_ids[ParamId::prestress]];
-  const double gamma_W1_over_n = parameters[global_param_ids[ParamId::gamma_W1_over_n]];
+  const double gamma_W1 = parameters[global_param_ids[ParamId::gamma_W1]];
+  const double gamma_eta = parameters[global_param_ids[ParamId::gamma_eta]];
   const double volume0 = parameters[global_param_ids[ParamId::volume0]];
   const double gamma_sigma_max = parameters[global_param_ids[ParamId::gamma_sigma_max]];
   const double n = parameters[global_param_ids[ParamId::n]];
   const double volume = y[global_var_ids[6]];
+  const double dvolume_dt = dy[global_var_ids[6]];
   const double stress = y[global_var_ids[4]];
   const double Pout = y[global_var_ids[2]];
   system.C.coeffRef(global_eqn_ids[0]) = -Pout*pow((volume + volume0)/volume0, (2.0/3.0)/n) + Pout + stress*pow((volume + volume0)/volume0, (1.0/3.0)/n) - stress;
-  system.C.coeffRef(global_eqn_ids[1]) = 4*gamma_W1_over_n*n - 4*gamma_W1_over_n*n*pow(volume/volume0 + 1, -2/n) + prestress;
+  system.C.coeffRef(global_eqn_ids[1]) = -2*dvolume_dt*gamma_eta/(n*volume0) + (2.0/3.0)*dvolume_dt*gamma_eta*pow(volume/volume0 + 1, -(5.0/3.0)/n)*pow(volume/volume0 + 1, -1 + (7.0/3.0)/n)/(n*volume0) + (4.0/3.0)*dvolume_dt*gamma_eta*pow(volume/volume0 + 1, -(17.0/3.0)/n)*pow(volume/volume0 + 1, -1 + (7.0/3.0)/n)/(n*volume0) + 4*gamma_W1 - 4*gamma_W1*pow(volume/volume0 + 1, -2/n) + prestress;
   system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[2]) = 1 - pow((volume + volume0)/volume0, (2.0/3.0)/n);
   system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[4]) = pow((volume + volume0)/volume0, (1.0/3.0)/n) - 1;
   system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[6]) = (1.0/3.0)*pow((volume + volume0)/volume0, (1.0/3.0)/n)*(-2*Pout*pow((volume + volume0)/volume0, (1.0/3.0)/n) + stress)/(n*(volume + volume0));
-  system.dC_dy.coeffRef(global_eqn_ids[1], global_var_ids[6]) = 8*gamma_W1_over_n*pow((volume + volume0)/volume0, -2/n)/(volume + volume0);
+  system.dC_dy.coeffRef(global_eqn_ids[1], global_var_ids[6]) = (2.0/9.0)*pow((volume + volume0)/volume0, -(17.0/3.0)/n)*(dvolume_dt*gamma_eta*pow((volume + volume0)/volume0, (19.0/3.0 - n)/n)*(2 - 3*n) - 2*dvolume_dt*gamma_eta*pow((volume + volume0)/volume0, (1.0/3.0)*(7 - 3*n)/n)*(3*n + 10) + 36*gamma_W1*n*volume0*pow((volume + volume0)/volume0, (11.0/3.0)/n))/(pow(n, 2)*volume0*(volume + volume0));
+  system.dC_dydot.coeffRef(global_eqn_ids[1], global_var_ids[6]) = (2.0/3.0)*gamma_eta*pow((volume + volume0)/volume0, -(22.0/3.0)/n)*(-3*pow((volume + volume0)/volume0, (22.0/3.0)/n) + 2*pow((volume + volume0)/volume0, (4 - n)/n) + pow((volume + volume0)/volume0, (8 - n)/n))/(n*volume0);
 
   // active stress
   system.C.coeffRef(global_eqn_ids[2]) = -act_plus*gamma_sigma_max;

src/model/ChamberSphere.h

@@ -139,14 +139,15 @@ class ChamberSphere : public Block {
   enum ParamId {
     n = 0,
     volume0 = 1,
-    gamma_W1_over_n = 2,
-    gamma_sigma_max = 3,
-    prestress = 4,
-    alpha_max = 5,
-    alpha_min = 6,
-    tsys = 7,
-    tdias = 8,
-    steepness = 9
+    gamma_W1 = 2,
+    gamma_eta = 3,
+    gamma_sigma_max = 4,
+    prestress = 5,
+    alpha_max = 6,
+    alpha_min = 7,
+    tsys = 8,
+    tdias = 9,
+    steepness = 10
   };
 
   /**
@@ -159,7 +160,8 @@ class ChamberSphere : public Block {
       : Block(id, model, BlockType::chamber_sphere, BlockClass::vessel,
               {{"n", InputParameter()},
                {"volume0", InputParameter()},
-               {"gamma_W1_over_n", InputParameter()},
+               {"gamma_W1", InputParameter()},
+               {"gamma_eta", InputParameter()},
                {"gamma_sigma_max", InputParameter()},
                {"prestress", InputParameter()},
                {"alpha_max", InputParameter()},