How do I incorporate PDEmodel for uncoupled thermo-elastic analysis?

I've been trying to solve the Navier-Lame equation for thermoelasticity in a reinforced concrete beam by incorporating PDEmodel. It's inconvinient, but I want my Young's modulus to be a function of temperature. My procedure is solving a thermal transient problem and then taking the obtained temperature field into my c and f coefficients. The geometry is created via multicuboid, then transformed into fegeometry and fed to the thermal transient model. Then the same geometry (before conversion to fegeometry) is assigned to the general pde model (createpde(3)). The mesh is seperate between models.

The solution I get is nonsense: displacements for a 2m beam, fixed on both ends and subjected to temperature gradient of 1500K/m is around 1.2m (looks like torsion). I cannot find any mathematical errors. The attached code shows my definition of both c and f coefficients (it is just an overview). Can I use the interpolateTemperature and evaluateTemperatureGradient effectively inside my coefficients? Perhaps there is another problem with my approach or some PDE Toolbox limitations?

Thanks for any potential help on this.

Best regards

KS

%Body of the structural problem, Rt are the thermal transient reults
model = createpde(3);
model.Geometry = g;
generateMesh(model, 'Hmax', 0.01);

ccoeff1 = @(location,state) myfunWithAdditionalArgs1(location,state, Rt, length(tlist));
ccoeff2 = @(location,state) myfunWithAdditionalArgs2(location,state, Rt, length(tlist));
fcoeff1 = @(location,state) myfunWithAdditionalArgs3(location, state, Rt, length(tlist));
fcoeff2 = @(location,state) myfunWithAdditionalArgs4(location,state, Rt, length(tlist));
gcoef = @(location,state) myfunWithAdditionalArgs5(location,state, Rt, length(tlist));

specifyCoefficients(model,'m',0,'d',0,'c', ccoeff1,'a',0,'f',[0;0;0], 'Cell', 1);
specifyCoefficients(model,m=0,d=0,c=ccoeff2,a=0,f=[0;0;0], Cell=[2, 3]);

applyBoundaryCondition(model, "neumann", g=gcoef, Face=[3, 4, 5, 6]);
applyBoundaryCondition(model,"dirichlet",Face=2,u=[0,0,0]);
applyBoundaryCondition(model,"dirichlet",Face=1,u=[0,0,0]);

result = solvepde(model)




%Defining of the coefficients-just the concrete
ccoeff1 = @(location,state) myfunWithAdditionalArgs1(location,state, Rt, length(tlist));
function cmatrix1 =  myfunWithAdditionalArgs1(location, state, T, id)
    n1 = 45;
    nr = numel(location.x);
    cmatrix1 = zeros(n1, nr);
    E = 30e9.*(1-0.0001*interpolateTemperature(T, location.x, location.y, location.z, id));
    nu = 0.2;
    mu = E./(2.*(1+nu));
    lambda = E.*nu./((1-2.*nu).*(1+nu));

    cmatrix1 (1,:) = 2.*mu+lambda;  
    cmatrix1 (3,:) = mu;
    cmatrix1 (6,:) = mu;
    cmatrix1 (8,:) = mu;
    cmatrix1 (10,:) = lambda;
    cmatrix1 (16,:) = mu;
    cmatrix1 (18,:) = 2.*mu+lambda;
    cmatrix1 (21,:) = mu;
    cmatrix1 (24,:) = mu;
    cmatrix1 (28,:) = lambda;
    cmatrix1 (36,:) = mu;
    cmatrix1 (38,:) = lambda;
    cmatrix1 (40,:) = mu;
    cmatrix1 (42,:) = mu;
    cmatrix1 (45,:) = 2.*mu+lambda;
end

fcoeff1 = @(location,state) myfunWithAdditionalArgs3(location, state, Rt, length(tlist));
function fmatrix1 = myfunWithAdditionalArgs3(location, state, T, id)
    n1=3;
    nr = numel(location.x);
    fmatrix1 = zeros(n1, nr);
    alpha = 1.2e-5;
    nu = 0.2;
    E =30e9.*(1-0.0001*interpolateTemperature(T, location.x, location.y, location.z, id));
    [gradTx, gradTy, gradTz] = evaluateTemperatureGradient(T, location.x, location.y, location.z, id);
    fmatrix1(1,:) = E.*alpha.*gradTx./(1-2.*nu);
    fmatrix1(2,:) = E.*alpha.*gradTy./(1-2.*nu);
    fmatrix1(3,:) = E.*alpha.*gradTz./(1-2.*nu);
end

Answer 1

[Anumeha]

Hi Karol,

I understand you are trying to use the "interpolateTemperature" and "evaluateTemperatureGradient" effectively inside your coefficients and getting unexpected output. When using inbuilt functions, there’s often some asuumption already in place which effects the output for a tailored problem. I suggest you go thorugh the following documentation page:

https://www.mathworks.com/help/pde/ug/createpde.html

Also check out the topics within “See Also” section at the end of the page. You can also use this command to understand thoroughly how it is implemented:

edit createpde

Hope this helps!