From c5d4aceb97c8ae540e1e4a641f38e5537bdb36ef Mon Sep 17 00:00:00 2001
From: "vinugholap@gmail.com" <Vinugitlab$123>
Date: Tue, 9 Mar 2021 18:57:10 +0100
Subject: [PATCH] Updates in Terms for force and stiffness with new coupling
terms from emFieldSource
---
dG3D/src/dG3DTerms.cpp | 52 ++++++++++++++++++++++++++++++++++++++++++
1 file changed, 52 insertions(+)
diff --git a/dG3D/src/dG3DTerms.cpp b/dG3D/src/dG3DTerms.cpp
index d46e436d5..a0f1d2f80 100644
--- a/dG3D/src/dG3DTerms.cpp
+++ b/dG3D/src/dG3DTerms.cpp
@@ -770,6 +770,11 @@ void dG3DForceBulk::get(MElement *ele,int npts,IntPt *GP,fullVector<double> &vFo
if (getConstitutiveExtraDofDiffusionAccountMecaSource()){
w += ipv->getConstRefToMechanicalSource()(extraDOFField);
}
+ double w_AV = 0.0;
+ if (getConstitutiveCurlAccountFieldSource() && extraDOFField == 0)
+ {
+ w_AV += ipv->getConstRefToEMFieldSource()(extraDOFField);
+ }
int indexField=3+getNumNonLocalVariable()+extraDOFField;
int nbFFRow = nlspace->getNumShapeFunctions(ele,indexField);
int nbFFTotalLastRow = nlspace->getShapeFunctionsIndex(ele,indexField);
@@ -782,6 +787,7 @@ void dG3DForceBulk::get(MElement *ele,int npts,IntPt *GP,fullVector<double> &vFo
if(extraDOFField==0)
{ // corresponds to temperature, but when considering electro-mechanics, ???
vFor(j+nbFFTotalLastRow) += getConstitutiveExtraDofDiffusionEqRatio()*ratio*(w*Vals[j+nbFFTotalLastRow]);// related to cp
+ vFor(j+nbFFTotalLastRow) += getConstitutiveExtraDofDiffusionEqRatio()*ratio*(w_AV*Vals[j+nbFFTotalLastRow]); // related to eddy current + hysteresis loss
}
}
}
@@ -1510,6 +1516,21 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
dwdV += ipv->getConstRefTodMechanicalSourcedField()(extraDOFField,extraDOFField2);
}
}
+ SVector3 dw_AVdGradT(0.0);
+ SVector3 dw_AVdGradV(0.0);
+ double dw_AVdT = 0.0;
+ if (getConstitutiveCurlAccountFieldSource())
+ {
+ if (extraDOFField == 0 && extraDOFField2 == 0)
+ {
+ dw_AVdT += ipv->getConstRefTodEMFieldSourcedField()(extraDOFField,extraDOFField2);
+ dw_AVdGradT += ipv->getConstRefTodEMFieldSourcedGradField()[extraDOFField][extraDOFField2];
+ }
+ if (extraDOFField == 0 && extraDOFField2 == 1)
+ {
+ dw_AVdGradV += ipv->getConstRefTodEMFieldSourcedGradField()[extraDOFField][extraDOFField2];
+ }
+ }
for(int j=0;j<nbFFRow; j++)
{
for(int k=0;k<nbFFCol;k++)
@@ -1528,11 +1549,24 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
{
// related to cp
mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol)+=getConstitutiveExtraDofDiffusionEqRatio()*ratio*dwdt*Vals[j+nbFFTotalLastRow]*Vals[k+nbFFTotalLastCol];
+ mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol) -= getConstitutiveExtraDofDiffusionEqRatio()*ratio*dw_AVdT*Vals[j+nbFFTotalLastRow]*Vals[k+nbFFTotalLastCol];
+
+ for (unsigned int m = 0; m < 3; ++m)
+ {
+ mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol) -= getConstitutiveExtraDofDiffusionEqRatio()*ratio*Vals[j+nbFFTotalLastRow]*
+ dw_AVdGradT.operator()(m)*Grads[k+nbFFTotalLastCol][m];
+ }
}
if (extraDOFField == 0 && extraDOFField2 == 1)
{
mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol)+=getConstitutiveExtraDofDiffusionEqRatio()*ratio*dwdV
*Vals[j+nbFFTotalLastRow]*Vals[k+nbFFTotalLastCol];
+
+ for (unsigned int m = 0; m < 3; ++m)
+ {
+ mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol) -= getConstitutiveExtraDofDiffusionEqRatio()*ratio*Vals[j+nbFFTotalLastRow]*
+ dw_AVdGradV.operator()(m)*Grads[k+nbFFTotalLastCol][m];
+ }
}
}
}
@@ -1683,6 +1717,8 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
const STensor3 & dFluxdVectorCurl = ipv->getConstRefTodFluxdVectorCurl()[extraDOFField][curlField];
SVector3 dwTdVectorPotential(0.0);
+ SVector3 dw_AVdVectorPotential(0.0);
+ SVector3 dw_AVdVectorCurl(0.0);
if (extraDOFField == 0)
{
@@ -1695,6 +1731,12 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
{
dwTdVectorPotential += ipv->getConstRefTodMechanicalSourcedVectorPotential()[extraDOFField][curlField];
}
+
+ if (getConstitutiveCurlAccountFieldSource())
+ {
+ dw_AVdVectorPotential += ipv->getConstRefTodEMFieldSourcedVectorPotential()[extraDOFField][curlField];
+ dw_AVdVectorCurl += ipv->getConstRefTodEMFieldSourcedVectorCurl()[extraDOFField][curlField];
+ }
}
for (unsigned int j = 0; j < nbFFRow; ++j)
@@ -1723,6 +1765,16 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
(Vals[j+nbFFTotalLastRow] *
dwTdVectorPotential.operator()(l) *
curlVals[k+numEdgeShapeFuncTotalLastCol][l]);
+ mStiff(j+nbFFTotalLastRow,
+ k+numEdgeShapeFuncTotalLastCol+numStandardDof) -= getConstitutiveExtraDofDiffusionEqRatio() * ratio *
+ (Vals[j+nbFFTotalLastRow] *
+ dw_AVdVectorPotential.operator()(l) *
+ curlVals[k+numEdgeShapeFuncTotalLastCol][l]);
+ mStiff(j+nbFFTotalLastRow,
+ k+numEdgeShapeFuncTotalLastCol+numStandardDof) -= getConstitutiveExtraDofDiffusionEqRatio() * ratio *
+ (Vals[j+nbFFTotalLastRow] *
+ dw_AVdVectorCurl.operator()(l) *
+ CurlcurlVals[k+numEdgeShapeFuncTotalLastCol][l]);
}
}
}
--
GitLab