Skip to content
Snippets Groups Projects
Commit 5b839b29 authored by Vinayak Gholap's avatar Vinayak Gholap
Browse files

Changes for large deformation in weak coupling scheme: to be corrected

parent bc7b431d
Branches
No related tags found
1 merge request!386Major update: Computing average IPVar over time period, rewriting last time...
......@@ -260,6 +260,16 @@ void mlawLinearElecMagTherMech::constitutive(
// here magnetic, electric field and temp dependent
const double ctime = getTime();
const double dh = getTimeStep();
double invdh;
if (dh < 0.0)
Msg::Error("Negative time step size in mlawLinearEMTM::constitutive()");
else if (dh == 0.0)
invdh = 0.0;
else
invdh = 1.0/dh;
// if true then we solve EM problem and thus reevaluate Magnetic potential as well
// else do not consider Magnetic contributions -> TM problem
if (this->_evaluateCurlField)
......@@ -332,15 +342,7 @@ void mlawLinearElecMagTherMech::constitutive(
}
// Add contribution from j_p and j_m too???
const double ctime = getTime();
const double dh = getTimeStep();
double invdh;
if (dh < 0.0)
Msg::Error("Negative time step size in mlawEMTM::constitutive()");
else if (dh == 0.0)
invdh = 0.0;
else
invdh = 1.0/dh;
// phase angle for sin wave function
// phi = M_PI / 2.0 if static analysis in EM
// phi = 0.0 if time domain dynamic analysis
......@@ -523,7 +525,11 @@ void mlawLinearElecMagTherMech::constitutive(
dThermalEMFieldSourcedGradV(m) += (-dedgradV(m,n) * (An(n) - A0(n)) * invdh);
if (_useFluxT)
{
dThermalEMFieldSourcedGradV(m) += (- dedgradV(m,n) * gradV(n));
dThermalEMFieldSourcedA(m) += (-dfluxDdA(m,n) * gradV(n));
}
}
dThermalEMFieldSourcedA(m) += (-fluxD(m) * invdh);
dThermalEMFieldSourcedB(m) += /*(H(m) * invdh)*/ 0.0;
......@@ -564,6 +570,51 @@ void mlawLinearElecMagTherMech::constitutive(
STensorOperation::zero(dmechFieldSourcedgradV);*/
}
/*
if(!this->_evaluateCurlField)
{
STensor3 Ft;
STensorOperation::unity(Ft);
//Ft = Fn;
const double Jac = Ft.determinant();
const double Jacinv = 1.0/Jac;
STensor3 Ftinv;
STensorOperation::zero(Ftinv);
STensorOperation::inverseSTensor3(Ft, Ftinv);
STensor3 I;
STensorOperation::unity(I);
if(stiff)
{
if(!isInductor())
{
STensorOperation::zero(dThermalEMFieldSourcedA);
STensorOperation::zero(dThermalEMFieldSourcedGradV);
for (unsigned int m = 0; m < 3; ++m)
{
for (unsigned int n = 0; n < 3; ++n)
{
for (unsigned int k = 0; k < 3; ++k)
{
for (unsigned int l = 0; l < 3; ++l)
{
dThermalEMFieldSourcedA(m) += (-Jacinv * Ft(n, m) * dfluxDdA(n, k) *
(Ftinv(k, l) * (An(l) - A0(l)) * invdh));
if (this->_useFluxT)
dThermalEMFieldSourcedA(m) += (-Jacinv * Ft(n,m) * dfluxDdA(n,k) * Ftinv(k,l) *
gradV(l));
}
}
dThermalEMFieldSourcedA(m) += (-Jacinv * Ft(n,m) * fluxD(n) * invdh);
if (_useFluxT)
dThermalEMFieldSourcedGradV(m) += (-Jacinv * Ft(n,m) * fluxD(n));
}
}
}
}
}*/
}
double mlawLinearElecMagTherMech::deformationEnergy(STensor3 defo,const SVector3 &gradV,const SVector3 &gradT,
......
......@@ -389,13 +389,19 @@ if(evaluateCurlField)
STensorOperation::zero(dqdV);
STensorOperation::zero(dqdF);
STensor3 Ft;
STensorOperation::unity(Ft);
//if (!evaluateCurlField)
//Ft = Fn;
for(int i=0;i<3;i++)
{
fluxT(i)=0.;
for(int j=0;j<3;j++)
//for(int k = 0; k < 3; ++k)
{
fluxT(i)+=(_k(i,j)+_l0(i,j)*_seebeck*_seebeck*T)*gradT(j)+_l0(i,j)*_seebeck*T*gradV(j);
fluxjy(i)+=(_k(i,j)+_l0(i,j)*_seebeck*_seebeck*T+_l0(i,j)*_seebeck*V)*gradT(j)+(_l0(i,j)*V+_l0(i,j)*_seebeck*T)*gradV(j);
fluxT(i)+=(_k(i,j)+_l0(i,j)*_seebeck*_seebeck*T)*gradT(j)+_l0(i,j)*_seebeck*T/**Ft(k,j)*/*gradV(j);
  • Author Developer

    Note: Need to pull back of EM fields when computing in TM problem? a-v computed in EM and then employed in TM problem on reference config. So: a = F^{-T} A, e = F^{-T} E. Also EM material properties? Need to check.

  • Please register or sign in to reply
fluxjy(i)+=(_k(i,j)+_l0(i,j)*_seebeck*_seebeck*T+_l0(i,j)*_seebeck*V)*gradT(j)+(_l0(i,j)*V+_l0(i,j)*_seebeck*T)/**Ft(k,j)*/*gradV(j);
}
}
......@@ -403,12 +409,13 @@ if(evaluateCurlField)
{
dqdT(i)=0.;
for(int j=0;j<3;j++)
//for(int k = 0; k < 3; ++k)
{
dqdT(i)+=dkdT(i,j)*gradT(j)+2*_seebeck*dseebeckdT*_l0(i,j)*T*gradT(j)+_seebeck*_seebeck*dldT(i,j)*T*gradT(j)+_seebeck*_seebeck*_l0(i,j)*gradT(j)
+dseebeckdT*_l0(i,j)*T*gradV(j)+_seebeck*dldT(i,j)*T*gradV(j)+_seebeck*_l0(i,j)*gradV(j);
+dseebeckdT*_l0(i,j)*T*/*Ft(k,j)**/gradV(j)+_seebeck*dldT(i,j)*T/**Ft(k,j)*/*gradV(j)+_seebeck*_l0(i,j)/**Ft(k,j)*/*gradV(j);
djydT(i)+=_seebeck*_seebeck*_l0(i,j)*gradT(j)+_seebeck*_l0(i,j)*gradV(j);
djydV(i)+=_seebeck*_l0(i,j)*gradT(j)+_l0(i,j)*gradV(j);
djydT(i)+=_seebeck*_seebeck*_l0(i,j)*gradT(j)+_seebeck*_l0(i,j)/**Ft(k,j)*/*gradV(j);
djydV(i)+=_seebeck*_l0(i,j)*gradT(j)+_l0(i,j)*/*Ft(k,j)**/gradV(j);
}
}
STensorOperation::zero(dqdgradT);
......
......@@ -3799,8 +3799,15 @@ void nonLinearMechSolver::fillIPDataOnPhysical(const double curtime, const int n
{
// get IPVariable value to store
const double val = ips[j]->getState(IPStateBase::current)->get(IPField::Output(ipval));
const double Jac = ips[j]->getState(IPStateBase::current)->get(IPField::JACOBIAN);
// pull-back operation for the IPVariable
// W_em = Jac * w_em
const double RefStateVal = Jac * val;
// insert GP->num, IPField value
gpNum_IPVariableValue.emplace(j, val);
gpNum_IPVariableValue.emplace(j, RefStateVal);
}
// insert ele_num, GP->num (0,1,2,3), IPVariable value
eleNum_gpNum_IPVariableValue.emplace(elnum, gpNum_IPVariableValue);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment