Since you want to both check mean and variance I'd recommend glmmTMB with its dispformula. Familiarize yourself with what that actually does. You might start here: https://stats.stackexchange.com/a/615211/341520

There are also more advanced options with proper correlation structures. You don't have that many devices though,but check it out here: https://cran.r-project.org/web/packages/glmmTMB/vignettes/covstruct.html

Anyway here's a simulation with huge effects an looking at everything that is not an interaction. You might start from here:

library(glmmTMB)
library(tidyverse)
set.seed("0815")
n_device <- 10
device_random_effect <- rnorm(n_device*2, sd = 1)
dat <- data.frame(articulations = rep(as_factor(c(0, 50, 100, 120)), 
           n_device), 
           angle = as_factor(rep(rep(c(16, 130), each = 4), 
           n_device)),
           device_id = rep(1:(2*n_device), each = 8),
           type = rep(c("X", "Y"), each = 8*n_device))

dat$y <-  5 - as.integer(dat$articulations) + 1*(dat$angle == "16") + 
  device_random_effect[dat$device_id] + 
  rnorm(8*2*n_device, sd = 1 + 2*(dat$angle == "16") + 
  3*(dat$type == "Y"))
# the plot does not show much except the higher variance for y
ggplot(dat, aes( y = y, x = articulations, color = type)) +
  geom_boxplot() +
  geom_jitter() + 
  facet_wrap(vars(angle))
# straight forard model
my_model <- glmmTMB(y ~ articulations + angle + type+ (1|device_id), 
        dispformula = ~ articulations + angle + type,
        data = dat)

summary(my_model)
# advanced and not properly converging 
my_model2 <- glmmTMB(y ~ articulations + angle + type + 
                      us(articulations + angle|device_id), 
                    dispformula = ~ type - 1,
                    data = dat)
summary(my_model2)