From fb0b3e8adfa31ca97af41f5f1bb3188db108a8db Mon Sep 17 00:00:00 2001 From: ralberd Date: Wed, 11 Sep 2024 11:55:42 -0600 Subject: [PATCH] fixes bug in the dissipation potential and updates test to better check --- optimism/material/HyperViscoelastic.py | 21 +++++++------ test/material/test_HyperVisco.py | 41 ++++++++++---------------- 2 files changed, 25 insertions(+), 37 deletions(-) diff --git a/optimism/material/HyperViscoelastic.py b/optimism/material/HyperViscoelastic.py index 1dae91c6..dea869f9 100644 --- a/optimism/material/HyperViscoelastic.py +++ b/optimism/material/HyperViscoelastic.py @@ -29,7 +29,7 @@ def compute_state_new(dispGrad, state, dt): return state def compute_material_qoi(dispGrad, state, dt): - return _compute_dissipation(dispGrad, state, dt, props) + return _compute_dissipated_energy(dispGrad, state, dt, props) return MaterialModel(compute_energy_density = energy_density, compute_initial_state = compute_initial_state, @@ -63,9 +63,11 @@ def _energy_density(dispGrad, state, dt, props): Ee = Ee_trial - delta_Ev W_neq = _neq_strain_energy(Ee, props) - Psi = _incremental_dissipated_energy(Ee, dt, props) + + Dv = delta_Ev / dt + Psi = _dissipation_potential(Dv, props) - return W_eq + W_neq + Psi + return W_eq + W_neq + dt * Psi def _eq_strain_energy(dispGrad, props): K, G = props[PROPS_K_eq], props[PROPS_G_eq] @@ -81,22 +83,19 @@ def _neq_strain_energy(elasticStrain, props): G_neq = props[PROPS_G_neq] return G_neq * TensorMath.norm_of_deviator_squared(elasticStrain) -def _incremental_dissipated_energy(elasticStrain, dt, props): +def _dissipation_potential(Dv, props): G_neq = props[PROPS_G_neq] tau = props[PROPS_TAU] eta = G_neq * tau - Me = 2. * G_neq * elasticStrain - M_bar = TensorMath.norm_of_deviator_squared(Me) - - return 0.5 * dt * M_bar / eta + return eta * TensorMath.norm_of_deviator_squared(Dv) -def _compute_dissipation(dispGrad, state, dt, props): +def _compute_dissipated_energy(dispGrad, state, dt, props): Ee_trial = _compute_elastic_logarithmic_strain(dispGrad, state) delta_Ev = _compute_state_increment(Ee_trial, dt, props) - Ee = Ee_trial - delta_Ev - return _incremental_dissipated_energy(Ee, dt, props) + Dv = delta_Ev / dt + return dt * _dissipation_potential(Dv, props) def _compute_state_new(dispGrad, stateOld, dt, props): Ee_trial = _compute_elastic_logarithmic_strain(dispGrad, stateOld) diff --git a/test/material/test_HyperVisco.py b/test/material/test_HyperVisco.py index 19c3a547..090aa318 100644 --- a/test/material/test_HyperVisco.py +++ b/test/material/test_HyperVisco.py @@ -21,7 +21,8 @@ def setUp(self): G_neq_1 = 5.0 tau_1 = 25.0 self.G_neq_1 = G_neq_1 - self.tau_1 = tau_1 # needed for analytic calculation below + self.tau_1 = tau_1 + self.eta = G_neq_1 * tau_1 self.props = { 'equilibrium bulk modulus' : K_eq, 'equilibrium shear modulus' : G_eq, @@ -52,6 +53,7 @@ def test_loading_only(self): state_old = self.compute_initial_state() energies = np.zeros(n_steps) states = np.zeros((n_steps, state_old.shape[0])) + dissipated_energies = np.zeros(n_steps) # numerical solution for n, F in enumerate(Fs): @@ -59,19 +61,15 @@ def test_loading_only(self): energies = energies.at[n].set(self.energy_density(dispGrad, state_old, dt)) state_new = self.compute_state_new(dispGrad, state_old, dt) states = states.at[n, :].set(state_new) + dissipated_energies = dissipated_energies.at[n].set(self.compute_material_qoi(dispGrad, state_old, dt)) state_old = state_new - Fvs = jax.vmap(lambda Fv: Fv.reshape((3, 3)))(states) Fes = jax.vmap(lambda F, Fv: F @ np.linalg.inv(Fv), in_axes=(0, 0))(Fs, Fvs) Evs = jax.vmap(lambda Fv: log_symm(Fv))(Fvs) Ees = jax.vmap(lambda Fe: log_symm(Fe))(Fes) - Dvs = jax.vmap(lambda Ee: (1. / self.tau_1) * deviator(Ee))(Ees) - Mes = jax.vmap(lambda Ee: 2. * self.G_neq_1 * deviator(Ee))(Ees) - dissipations = jax.vmap(lambda D, M: np.tensordot(D, M), in_axes=(0, 0))(Dvs, Mes) - # analytic solution e_v_11 = (2. / 3.) * strain_rate * times - \ (2. / 3.) * strain_rate * self.tau_1 * (1. - np.exp(-times / self.tau_1)) @@ -86,30 +84,19 @@ def test_loading_only(self): [0., 0., e_22]] ), in_axes=(0, 0) )(e_e_11, e_e_22) + Me_analytic = jax.vmap(lambda Ee: 2. * self.G_neq_1 * deviator(Ee))(Ee_analytic) - Dv_analytic = jax.vmap(lambda Me: 1. / (2. * self.G_neq_1 * self.tau_1) * deviator(Me))(Me_analytic) - dissipations_analytic = jax.vmap(lambda Me, Dv: np.tensordot(Me, Dv), in_axes=(0, 0))(Me_analytic, Dv_analytic) + Dv_analytic = jax.vmap(lambda Me: 1. / (2. * self.eta) * deviator(Me))(Me_analytic) + dissipated_energies_analytic = jax.vmap(lambda Dv: dt * self.eta * np.tensordot(deviator(Dv), deviator(Dv)) )(Dv_analytic) # test self.assertArrayNear(Evs[:, 0, 0], e_v_11, 3) self.assertArrayNear(Ees[:, 0, 0], e_e_11, 3) self.assertArrayNear(Ees[:, 1, 1], e_e_22, 3) - self.assertArrayNear(dissipations, dissipations_analytic, 3) + self.assertArrayNear(dissipated_energies, dissipated_energies_analytic, 3) if plotting: plt.figure(1) - plt.plot(times, np.log(Fs[:, 0, 0])) - plt.xlabel('Time (s)') - plt.ylabel('F 11') - plt.savefig('times_Fs.png') - - plt.figure(2) - plt.plot(times, energies) - plt.xlabel('Time (s)') - plt.ylabel('Algorithmic energy') - plt.savefig('times_energies.png') - - plt.figure(3) plt.plot(times, Evs[:, 0, 0], marker='o', linestyle='None', markevery=10) plt.plot(times, e_v_11, linestyle='--') plt.xlabel('Time (s)') @@ -117,7 +104,7 @@ def test_loading_only(self): plt.legend(["Numerical", "Analytic"], loc=0, frameon=True) plt.savefig('times_viscous_stretch.png') - plt.figure(4) + plt.figure(2) plt.plot(times, Ees[:, 0, 0], marker='o', linestyle='None', markevery=10, label='11', color='blue') plt.plot(times, Ees[:, 1, 1], marker='o', linestyle='None', markevery=10, label='22', color='red') plt.plot(times, e_e_11, linestyle='--', color='blue') @@ -127,11 +114,13 @@ def test_loading_only(self): plt.legend(["11 Numerical", "22 Numerical", "11 Analytic", "22 Analytic"], loc=0, frameon=True) plt.savefig('times_elastic_stretch.png') - plt.figure(5) - plt.plot(times, dissipations, marker='o', linestyle='None', markevery=10) - plt.plot(times, dissipations_analytic, linestyle='--') + plt.figure(3) + plt.plot(times, np.cumsum(dissipated_energies), marker='o', linestyle='None', markevery=10) + plt.plot(times, np.cumsum(dissipated_energies_analytic), linestyle='--') + plt.xlabel('Time (s)') + plt.ylabel('Dissipated Energy Density (MPa)') plt.legend(["Numerical", "Analytic"], loc=0, frameon=True) - plt.savefig('times_dissipation.png') + plt.savefig('times_dissipated_energy.png')