System Modeler
and
Finite Elements

Exploring Integration

System Modeler -> FEM

FEM -> System Modeler

FEM <-> System Modeler

FEM <-> System Modeler

Not considered

Integration

System Modeler -> FEM

Room Heating

FEM -> System Modeler

Spring Constant

Beam Model (proof of concept)

System Modeler -> FEM

Room Heating

Room Heating

model=SystemModel["CollegeThermal.RoomHeating.Exercises.RoomModelRome"];
model["Diagram"]

Manipulate[
SystemModelPlot[SystemModelSimulate["CollegeThermalFEM.RoomHeating.Exercises.RoomModelRome",36000,<|"ParameterValues"->{"heatSupplied"->x1,"heatCapacityRoom"->x2*1000,"roomRefTemp"->x3}|>],ImageSize->Medium,AxesLabel->{"Time[s]","Temp[°C]"}],{{x2,40,"Heat capacity of the room [kJ/K]"},Row[{Slider[#,{0.1,100,1},ContinuousAction->False],Spacer[5],#}]&},{{x3,20,"Room reference temperature [°C]"},Row[{Slider[#,{0,30,1},ContinuousAction->False],Spacer[5],#}]&},{{x1,1200,"Power of heat source [W]"},Row[{Slider[#,{0,5000,1},ContinuousAction->False],Spacer[5],#}]&},ContinuousAction->False,SynchronousUpdating->False,SynchronousInitialization->False,Initialization:>ImportString[,"MO"]]

Room Heating-PDE

HeatTransfer monograph

∇·(-k ∇T(X))=Q

Room Heating-PDE

vars={T[x,y],{x,y}};

width=4;
length=4;
Ωroom=Rectangle[{0,0},{length,width}];
Ωradiator=Rectangle[{1,1/10},{3,3/20}];
Ω=RegionDifference[Ωroom,Ωradiator];
mesh=ToElementMesh[Ω,"RegionHoles"->None];
mesh["Wireframe"]

Room Heating-PDE

(Subscript[k, air])=ThermodynamicData["Air","ThermalConductivity"];
pars=<|"ThermalConductivity"->(Subscript[k, air])|>;

Parameter Q

pars["HeatSource"]=With[{rmf=RegionMember[Ωradiator]},Piecewise[{{Q,rmf[{x, y}]}},0]]

Boundary Conditions

System Modeler variable

sim=SystemModelSimulate[model,40000]

systemModelTime=30000;
srfTemp=sim[{"roomWithWaterRadiator.room.temp"},systemModelTime][[1]]

(Subscript[Γ, wall])=HeatTemperatureCondition[x==0||x==length||y==0||y==width,vars,pars,<|"SurfaceTemperature"->srfTemp|>]

Solve

pde={HeatTransferPDEComponent[vars,pars]==0,(Subscript[Γ, wall])};
pfun=ParametricNDSolveValue[pde,T,{x,y}∈mesh,Q];

Tfun=pfun[200]

Show[
ContourPlot[Tfun[x,y],{x,y}∈Ω,],
ToBoundaryMesh[mesh]["Wireframe"]
]

Average temperature

average[Q_?NumericQ]:=Module[{Tfun},
Tfun=pfun[Q];
NIntegrate[Tfun[x,y],{x,y}∈mesh]/Area[Ω]
]

average[200]

Average temperature (20^o)

root=FindRoot[average[Q]==20,{Q,200}]

TfunFinal=pfun[Q]/.root

Show[
ContourPlot[TfunFinal[x,y],{x,y}∈Ω,],
ToBoundaryMesh[mesh]["Wireframe"]
]

Learnings so far

System Modeler provides info needed for FEM: Boundary conditions and initial conditions

System Modeler may make assumptions: Uniformity of temperature across room, walls, etc

System Modeler provides a multi-domain system model

FEM model provides internal dynamics

FEM -> System Modeler

FEM Beam modeled as Mass-Spring System Model

Spring Constant

length=0.5;
height=0.1;
thickness=height;
Ω=Rectangle[{0,0},{length,height}];
mesh=ToElementMesh[Ω];
mesh["Wireframe"]

pars=<|"ShearModulus"->441000,"Thickness"->thickness,"PoissonRatio"->0.49|>;
vars={{u[x,y],v[x,y]},{x,y}};

(Subscript[Γ, fixed])=SolidFixedCondition[x==0,vars,pars];
(Subscript[Γ, loadParametric])=SolidBoundaryLoadValue[x==length,vars,pars,<|"Force"->{0,-force}|>];

pfun=ParametricNDSolveValue[{SolidMechanicsPDEComponent[vars,pars]==(Subscript[Γ, loadParametric]),(Subscript[Γ, fixed])},vars[[1]],vars[[-1]]∈mesh,{force}]

parametricForces={10,25,50,75,100};
AbsoluteTiming[displacements=pfun/@parametricForces;]

VectorDisplacementPlot[Evaluate[displacements[[-1]]],{x,y}∈MeshRegion[mesh],VectorSizes->Full]

Spring Constant

point={length,height/2};

data=Transpose[{Sqrt[Total[(displacements/.Thread[Rule[{x,y},point]])^2,{2}]],parametricForces}];
ListLinePlot[data,]

dx=D[Interpolation[data][x],x];
springConstant=Mean[Head[dx]["ValuesOnGrid"]]

Spring Constant

comp={"fixed"∈"Modelica.Mechanics.Translational.Components.Fixed","spring"∈"Modelica.Mechanics.Translational.Components.Spring","mass"∈"Modelica.Mechanics.Translational.Components.Mass"};

conn={"fixed.flange"\[UndirectedEdge]"spring.flange_a","spring.flange_b"\[UndirectedEdge]"mass.flange_a"};

linearSpringModel=ConnectSystemModelComponents["MassSpring",comp,conn,<|"InitialValues"->{"mass.m"->5.5,"spring.c"->springConstant,"mass.s"->0.3}|>]

Spring Constant

linearSpringDisplacementPlot=SystemModelPlot[linearSpringModel,{"mass.s"},2]

Spring Constant

parsNeoHook=pars;
parsNeoHook["SolidMechanicsMaterialModel"]="NeoHookean";

pfun=ParametricNDSolveValue[{SolidMechanicsPDEComponent[vars,parsNeoHook]==(Subscript[Γ, loadParametric]),(Subscript[Γ, fixed])},vars[[1]],vars[[-1]]∈mesh,{force}];

parametricForces={10,25,50,75,100};
AbsoluteTiming[displacements=pfun/@parametricForces;]

Spring Constant

data=Transpose[{Sqrt[Total[(displacements/.Thread[Rule[{x,y},point]])^2,{2}]],parametricForces}];
ListLinePlot[data,]

dx=D[Interpolation[data][x],x];

Plot[dx,Evaluate[Flatten[Join[{x},if["Domain"]]]]]

Spring Constant

modelNonlinear=ImportString[,"MO"];
modelNonlinear["Diagram"]

CreateDataSystemModel["SpringNonlinear.nonlinearSpring.forceDisplacement",data]

sim=SystemModelSimulate[modelNonlinear,5,<|"InitialValues"->{"mass.m"->5.5,"mass.s"->0.3}|>];

Spring Constant

nonlinearSpringDisplacementPlot=SystemModelPlot[sim,{"mass.s"},2,PlotStyle->Red];
nonlinearSpringForcePlot=SystemModelPlot[sim,{"mass.flange_a.f"},2,PlotStyle->Red];
Show[linearSpringDisplacementPlot,nonlinearSpringDisplacementPlot]

FEM -> System Modeler

Beam Model (proof of concept)

Jib Crane

craneRigid=ImportString[,"MO"]; 
craneRigid[{"Diagram",}]

simRigid=SystemModelSimulate[craneRigid];
SystemModelPlot[simRigid,]

SystemModelSimulationData[]

Beam Model

Ω=Rectangle[{0,0},{10,1}];
mesh=ToElementMesh[Ω,"MaxCellMeasure"->4,"MeshOrder"->1];

Show[
 mesh["Wireframe"],
 mesh["Wireframe"["MeshElement" -> "PointElements", "MeshElementIDStyle" -> Black]]
 ]

vars={{u[t,x,y],v[t,x,y]},t,{x,y}};
pars=<|"YoungModulus"->210*10^9,"PoissonRatio"->0.3,"MassDensity"->7850,"Thickness"->0.3|>;
(Subscript[Γ, wall])=SolidFixedCondition[x==0,vars,pars];
(Subscript[Γ, force])=SolidBoundaryLoadValue[x==10,vars,pars,<|"Force"->{0,-Sin[2π t]*100000}|>];
op=SolidMechanicsPDEComponent[vars,pars];

Beam Model

ics={u[0,x,y]==0,v[0,x,y]==0,Derivative[1,0,0][u][0,x,y]==0,Derivative[1,0,0][v][0,x,y]==0};

tEnd=5;
displacement=NDSolveValue[{op==(Subscript[Γ, force]),(Subscript[Γ, wall]),ics},{u[t,x,y],v[t,x,y]},{x,y}∈mesh,{t,0,tEnd}];

frames=Table[Show[ToBoundaryMesh[mesh]["Wireframe"],ElementMeshDeformation[{t,mesh},displacement,"ScalingFactor"->200]["Wireframe"[
"ElementMeshDirective"->Directive[EdgeForm[],FaceForm[LightGray]]]],Boxed->True,BoxStyle->Transparent,PlotRange->{All,{-1,2}}],
{t,0,tEnd,tEnd/50}];
ListAnimate[frames]

Beam Model

data=NDSolve`FEM`ProcessPDEEquations[,{u[t,x,y],v[t,x,y]},{x,y}∈mesh,{t,0,tEnd}];

{loadV,stiffnessM,dampingM,massM}=data[[1]]["SystemMatrices"];

MatrixPlot[stiffnessM]

Beam Model

CreateSystemModel["BeamFEM",{n==24,mm.D[v[t],{t,1}]+dd.v[t]+kk.x[t]==f[t],v[t]==D[x[t],t],x[0]==QuantityVariable["zeros",IndependentPhysicalQuantity[""]][n],
mm==massNew,
kk==stiffnessNew,
dd==0.001*stiffnessNew},
  t,,<|"LocalComponentDescriptions"->{"x"->"Field variables, displacements."}|>]

craneFEM=ImportString[,"MO"]; 
craneFEM["Diagram"]
beam2D=SystemModel["ControlledCraneFEM.Beam2D"]

SystemModelPlot[craneFEM,"Internal Displacements",ImageSize->Automatic,PlotRange->All]

Learnings so far

Use FEM to find model parameters for System Modeler

It is possible to get FEM dynamics to System Modeler

System Modeler and NDSolve rock!

Summary:

SystemModeler -> FEM
FEM -> SystemModeler
Community posts: [1] [2]
Documentation

Feel free to contact us:

vedats@wolfram.com
ruebenko@wolfram.com

Thanks

Init

Needs["NDSolve`FEM`"]
$HistoryLength=0;